From local monitoring to a broad-scale viability assessment : a case study for the Bonelli ’ s Eagle in western Europe

Population viability analysis (PVA) has become a basic tool of current conservation practice. However, if not accounted for properly, the uncertainties inherent to PVA predictions can decrease the reliability of this type of analysis. In the present study, we performed a PVA of the whole western European population (France, Portugal, and Spain) of the endangered Bonelli’s Eagle (Aquila fasciata), in which we thoroughly explored the consequences of uncertainty in population processes and parameters on PVA predictions. First, we estimated key vital rates (survival, fertility, recruitment, and dispersal rates) using monitoring, ringing, and bibliographic data from the period 1990–2009 from 12 populations found throughout the studied geographic range. Second, we evaluated the uncertainty about model structure (i.e., the assumed processes that govern individual fates and population dynamics) by comparing the observed growth rates of the studied populations with model predictions for the same period. Third, using the model structures suggested in the previous step, we assessed the viability of both the local populations and the overall population. Finally, we analyzed the effects of model and parameter uncertainty on PVA predictions. Our results strongly support the idea that all local populations in western Europe belong to a single, spatially structured population operating as a source– sink system, whereby the populations in the south of the Iberian Peninsula act as sources and, thanks to dispersal, sustain all other local populations, which would otherwise decline. Predictions regarding population dynamics varied considerably, and models assuming more constrained dispersal predicted more pessimistic population trends than models assuming greater dispersal. Model predictions accounting for parameter uncertainty revealed a marked increase in the risk of population declines over the next 50 years. Sensitivity analyses indicated that adult and pre-adult survival are the chief vital rates regulating these populations, and thus, the conservation efforts aimed at improving these survival rates should be strengthened in order to guarantee the long-term viability of the European populations of this endangered species. Overall, the study provides a framework for the implementation of multi-site PVAs and highlights the importance of dispersal processes in shaping the population dynamics of long-lived birds distributed across heterogeneous landscapes.


INTRODUCTION
Population viability analysis (PVA), the use of quantitative methods to predict the likely future status of a population, has become a basic, widely used tool in current conservation practice (Gilpin and Soule´1986, Boyce 1992, Brook et al. 2000, Morris and Doak 2002, Lindenmayer et al. 2003). One of the major contributions of this type of mathematical method is its ability to provide technical and scientific information for guiding decision makers working with populations of conservation target species. However, the uncertainties inherent to PVA predictions may decrease the reliability of predictions if not accounted for properly and thus may limit their applied value , McGowan et al. 2011.
Two crucial elements in the formulation of a useful PVA are, firstly, an accurate estimation of the vital rates on which the PVA is based and, secondly, a model that mimics the real world well enough to carry out its stated purpose (Reed et al. 2002). In particular, vital rates should ideally be estimated using large sample sizes derived from long-term data sets, thereby ensuring representative information for the studied population, as well as a reduction in sampling variances around parameter estimates and the reliable capture of the effect of environmental variation Kikkawa 1998, Crone et al. 2011). However, the demographic data required for a detailed simulation of PVA models is often scarce in the case of rare or endangered species, and therefore, accounting for uncertainties on the estimation of these parameters is nearly always necessary in order to increase the reliability of PVA predictions (Doak et al. 2005).
Uncertainty about model structure, the assumed processes that govern individual fates and population dynamics, may be much harder to evaluate, but can also be a great source of uncertainty in model predictions (Burgman et al. 2005, Valle et al. 2009. A reliable method for testing the validity of the assumptions implicit in the model structure is to compare past population trends and model predictions over the same period of time (Brook et al. 2000, Lindenmayer et al. 2001, Schtickzelle and Baguette 2004. To do so, it is necessary to have information on viability measures based on past population trends and estimated independently of the vital rates used to make model predictions. Unfortunately, information on past population trends is rarely available for endangered species, and therefore, most published PVAs do not fulfill this criterion. When data is available, however, model validity evaluation was conducted using several viability metrics such as quasi-extinction risk, patch occupancy, and population size (e.g., Brook et al. 2000, Coulson et al. 2001; nevertheless, few studies have made use of population growth rates, one of the most basic measures of mean population performance (Bierzychudek 1999, Crone et al. 2011). More than other viability metrics, this measure may be a better indicator of possible future problems than an extinction probability if the short-term risk of extinction is low (Morris and Doak 2002). Additionally, this measure may be reliably estimated with incomplete or short-term population size data. In order to validate model structure, observed trends are usually compared to the predictions of just a single model (e.g., Schtickzelle and Baguette 2004). Models can be also ''calibrated'' by comparing the predictions of the initially assumed model and subsequently modifying model assumptions in order to improve the match between predicted and observed population trends (e.g., Mc-Carthy et al. 2000). Despite the practical usefulness of this procedure, it is preferable to select the most plausible of a set of a priori-defined models with different assumptions, which should be based on general theory and the characteristics of the study species (Hilborn and Mangel 1997).
Model structure may be inaccurate if knowledge of the species' life-history traits or of the ecological drivers behind the variation in its vital rates is incomplete, or if the spatial structure of the population is either not taken into account or is inaccurately portrayed in the model (Lindenmayer et al. 2001, Schtickzelle et al. 2005). However, these inaccuracies may be of little relevance if they do not compromise PVA predictions. For example, the effect of ecological drivers may be obviated if past variation in vital rates correctly accounts for their future variation. Similarly, in a collection of populations distributed over a geographical area, very low or very high movement rates between local populations are expected to have little effect on PVA predictions if all the local populations are treated as a single population (Morris and Doak 2002). However, low to medium movement rates may explicitly need to take into account the spatial structure of the populations if they are to properly predict their overall demographic fate. On theoretical grounds, the metapopulation paradigm has played a significant role in the development of general theory in population ecology (Levins 1969, Hanski 1999. Essentially, this paradigm assumes that species occur in habitable patches, some of which may be occupied by local populations connected by dispersal. In heterogeneous habitats, the population growth rate of these local populations differs and thus may generate source-sink dynamics, characterized by the existence of sink populations (in which deaths outnumber births) whose persistence depends on dispersal from source populations (in which births outnumber deaths) (Pulliam 1988). Although empirical data have provided strong support for the main assumptions behind the metapopulation paradigm, metapopulation viability analyses still only represent a minority of the PVAs that are carried out (but see McCarthy et al. 2000, Lindenmayer et al. 2003, Schtickzelle et al. 2005) and rarely have they accounted for metapopulations ranging over broad geographical areas (Gonzalez-Suarez et al. 2006, Bonnot et al. 2011). This gap is due in part to the difficulties inherent in modeling complex spatially structured populations, but also in part to the difficulty of obtaining detailed demographic information from multiple local populations arrayed across broad spatial scales.
In the present study, we performed a viability assessment of a spatially structured population of an endangered avian predator, in which we thoroughly evaluated model assumptions with empirical data on past population trends and then assessed how uncertainties affect the model's predictions. The study focused on the western European population of Bonelli's Eagle (Aquila fasciata; see Plate 1), which corresponds to a fairly discrete demographic unit within a much larger distribution range extending from southeast Asia through the Middle East to the western Mediterranean (del Hoyo et al. 1992). The European Bonelli's Eagle population is estimated at 920-1100 pairs, of which ;80% are found in the Iberian Peninsula (see Fig. 1 for ranges and population codes; BirdLife International 2004). This species has undergone a dramatic decline in numbers and range in recent decades and is now listed as an endangered species in Europe (79/409/EEC; Rocamora 1994, BirdLife International 2004, Ontiveros et al. 2004. The main factor behind this decline is thought to be a demographic imbalance related to decreased fertility and, above all, increased adult and pre-adult mortality (Real and Mañosa 1997, Carrete et al. 2002, Soutullo et al. 2008. Other factors such as the loss of suitable habitat caused by changes in land use (Balbontı´n 2005), the decline in key prey populations (Real 1991, Moleo´n et al. 2009, and interspecific competition with the Golden Eagle Aquila chrysaetos (Ferna´ndez and Insausti 1990, Lo´pez-Lo´pez et al. 2004, Carrete et al. 2006 have also had negative effects on Bonelli's Eagle populations. Interestingly, however, trends at local population scale within western Europe have differed markedly over the last two decades: whereas most populations in central and northern Iberia have declined, those in southern Iberia have remained stable or even increased (Real and Mañosa 1997, del Moral 2006, Beja and Palma 2008, Soutullo et al. 2008 Previous demographic assessments of this species have afforded researchers valuable information on the dynamics of these populations and have suggested that the survival of adults and pre-adults is the chief vital rate regulating this eagle's populations (Real and Mañosa 1997, Carrete et al. 2005, Soutullo et al. 2008. Due to the scarcity of information on dispersal patterns, these previous studies assumed either that local populations were closed or that the probability of dispersal between populations was independent of distance. These simplifications could have seriously affected results given that empirical data indicates that, although Bonelli's Eagles tend to breed within or close to their population of origin, medium-and long-distance natal dispersal movements are not uncommon , suggesting that distance-dependent dispersal patterns may play an important role in overall population dynamics. Additionally, the multi-site model used by Soutullo et al. (2008) does not account for the role played by the non-studied local populations in the overall trend of the whole population. Finally, uncertainty in model structure or parameter estimates was not accounted for in these studies and thus their effects on viability predictions remain unknown.
The aims of the present study were to estimate the key vital rates of 12 local populations located throughout the distribution and demographic range of Bonelli's Eagle in western Europe (France, Portugal, and Spain) FIG. 1. Distribution area and population codes of Bonelli's Eagle (Aquila fasciata) in western Europe (shown in the areas of any gray darker than the light-gray background). The distribution area was drawn by buffering a 7-km area around the exact location of nests (FRA, CAT, SPT, and ARR-NPT) and around the center of the occupied 10 3 10 UTM quadrants given by Martı´and Moral (2003) anddel Moral (2006). In panel (A), the darker gray areas shows the local populations that were monitored, and in panel (B), the darker gray areas shows those considered in the viability assessment. Shades of gray highlight the different local populations. and use this information to develop a multi-site PVA to evaluate the likely future status of the whole western European population of this eagle. We used PVA predictions regarding past local population trends to test the model's assumptions. We evaluated viability by predicting the future demographic fate of populations and their probabilities of declining as far as certain population thresholds. We conducted a thorough analysis of the effects of both model and parameter uncertainties on the model predictions. Finally, we performed sensitivity and elasticity analyses to identify the key vital rates on which conservation efforts should be strengthened to improve the conservation status of this endangered species.

Definition of local populations
Bonelli's Eagle has a generally continuous distribution within western Europe (longitude 9830 0 W to 5855 0 E; latitude 36800 0 to 44824 0 N), although its range becomes more patchy at its western and northern edges (Fig. 1). Therefore, local populations cannot be defined straightforwardly on a basis of spatial discontinuities in the distribution of individuals. Nevertheless, available information highlights the existence of marked differences across the species' range in key demographic parameters such as survival, fertility, and local demographic trends (Real and Mañosa 1997, del Moral 2006, Carrascal and Seoane 2009. Therefore, we used as a general criterion for defining local populations the existence of groups of individuals sharing the same geographical area that have similar demographic characteristics. We delineated local populations by combining existing demographic and environmental data (Real and Mañosa 1997, Carrascal and Seoane 2009, Galicia et al. 2010 with data from our own field experience with the species. Despite not corresponding to highly isolated patches, as is generally used in metapopulation terminology, the local populations we identified did allow for spatial structure based on demographically relevant units (also see Carrete et al. 2008, Bonnot et al. 2011. To address the questions posed in our study, the following local populations were considered: (1) 12 monitored populations located across the western European range of Bonelli's Eagle used to estimate vital rates (Fig. 1A), (2) 18 populations used to evaluate uncertainty in the model structure and the natal dispersal kernel, and (3) 10 populations used for the viability assessment. In all our simulations we considered the whole population within the studied range. To do so, we initially split the whole population into 10 local populations based on their geographical proximity and distribution gaps, as well as on their demographic and habitat characteristics (Fig. 1B, Table 1). In a second step, we split the whole population into 18 local populations: the 12 monitored local populations (Fig.  1A) and the 6 non-monitored populations, which were defined by redrawing the boundaries of the 10 initial populations and the monitored populations ( Fig.   1A, B). This was done to keep the 12 monitored populations separate to be able to compare the observed and predicted trends in the model uncertainty evaluation. In our models, to estimate dispersal probabilities between populations, we considered the location of individual territories so that the shape of the boundaries of the populations did not affect our ability to properly estimate this probability (see Appendix B).

Life-history traits and life cycle of Bonelli's Eagle
Bonelli's Eagle is a long-lived territorial bird with delayed maturity and low breeding rates. Like other territorial raptors, Bonelli's Eagles pass through a transient nomadic phase (i.e., a dispersal period) after the post-fledging dependence period and before territorial recruitment. During the dispersal period, birds perform medium-to long-distance movements and often settle for extended periods in dispersal areas (geometric mean is 101 km, range 1-1020 km; , where they apparently suffer high human-induced mortality rates , Moleo´n et al. 2007, Cadahı´a et al. 2010. Territorial recruitment mostly occurs between three and four years of age . During the territorial phase, individuals exhibit strong pair-bonding behavior and establish home ranges of ;50 km 2 , with strong fidelity to the breeding area throughout the year and from one year to the next , Herna´ndez-Matı´as et al. 2011a. Most eagles in western Europe nest on cliffs, although in southern Portugal, most pairs nest in trees. Egg-laying (one or two, rarely three, eggs) usually takes place in January-March, and mean productivity rates range between 0.6 and 1.4 fledglings per pair and year (del Moral 2006).
The life cycle considered in our models was based on a post-breeding census of females in five age classes (Fig.  2). We chose to structure the population into age classes because age appears to be strongly correlated with the main vital rates in Bonelli's Eagle; this allowed us to model appropriately the differing contribution to population growth of individuals from different classes. Specifically, this eagle progressively increases its survival and fertility rates, and propensity for territorial behavior until the age of four, at which point birds acquire their adult plumage and these vital rates stabilize (Martı´nez et al. 2008, 2011b. The main vital rates considered were survival (S1, S23, S4, and SA, corresponding, respectively, to first-year, second-and third-year, fourth-year, and adult survival rates), fertility per female per year (F2, F3, FA, corresponding, respectively, to second-year, third-year, and adult fertility), and recruitment rates (R2, R3, R4, corresponding, respectively, to second-year, third-year, and fourth-year recruitment rates), which reflect the increasing propensity of eagles to become territorial as they age (see next section). We assumed an offspring sex ratio of 1:1 Mañosa 1997, Soutullo et al. 2008).

Vital-rate estimations
The estimation of vital rates was mainly based on: (1) monitoring data from territorial pairs in all studied populations (population size, adult survival, and fertilities), (2) data on individually known eagles from longterm ringing programs in certain populations (dispersal patterns from birth to territorial recruitment), and (3) previously published information on this species (preadult survival and recruitment rates).
Monitoring (1990-2009; Table 2) consisted of repeated visits during the breeding season (January-July) to a representative sample of territories (from 65% to 100% of all known territories) in 12 of the western European  Fig. 1 for the location of populations. Vital rates are coded with numbers to illustrate which values are shared by populations. The same code corresponds to the same estimated parameter in columns with the same heading: S1, S2, S3, S4, and SA (first-, second-, third-, fourth-year, and adult survivals, respectively) and F2, F3, and FA (second-, third-year, and adult fertility, respectively). Environmental characteristics are based on the classification proposed by Galicia et al. (2010), which describes the main climatic and lithological features of the Iberian Peninsula; these are as follows: A, mountain ranges moderated by proximity to oceans; B, rugged limestone country in the northern part of the central Spanish plateau and around the Ebro depression; C, Baetic mountains moderated by proximity to oceans and plains around Gibraltar; D, rugged siliceous landscapes and depressions in the central Spanish plateau and siliceous mountains in the eastern half of the peninsula; E, the Ebro depression, southern half of the central Spanish plateau and calcareous Baetic mountains with continental environments; F, mid-course of the river Guadiana, basin of the river Guadalquivir, and littoral and sub-littoral areas located south of the river Ebro. The degree of humanization is divided into three categories based on the density of people in the regions where eagle populations are located: the values of densities range from 16.0 to 43.6 inhabitants/km 2 (low), from 56.5 to 92.4 inhabitants/km 2 (medium), and from 129.9 to 272.1 inhabitants/km 2 (high) (data sources for human population density: http://www.ine.pt, http://www.ine.es, and http://www.insee.fr/). The solid lines represent the transitions between age classes, while the dashed lines represent the connection between age classes due to reproduction. SR and F represent the sex ratio and fertilities, respectively, for any given age class (F2, F3, and FA being fertility of two-and three-year olds and adults, respectively). S represents the survival rates of individuals during their first (S1), second and third (S23), fourth (S4), and fifth and subsequent (adult) years of life (SA). R represents the propensity to become territorial in immature birds (R2), subadults (R3), and first-year adults (R4) (denoted as a i in the text).
populations (see Fig. 1A). This allowed us to obtain the occupation rates of territories, identity of individuals (if marked), the plumage-age and sex of territorial birds, and the number of fledged chicks in each territory. We classified territorial birds according to their plumageages as juvenile (first year), immature (second year), subadult (third year), or adult (fourth year or older; Parellada 1984). Populations from North Africa were not monitored and were not considered in our models as available data suggest that migration across the Straits of Gibraltar is very uncommon and the few birds observed during migration periods are mostly territorial or dispersing birds that do not actually cross the Straits (Programa Migres 2009; G. M. Arroyo, personal communication).
Annual fertility was estimated for each local population as the mean productivity of monitored territories occupied by territorial pairs (i.e., mean number of fledglings per pair and year). We estimated fertility separately for immature, subadult, and adult females. We estimated the mean fertility during the study period to be the arithmetic mean of annual estimates. In cases where the sample size for non-adult fertilities was small and provided unreliable estimates, we used estimates from other populations with reliable data for this parameter (Tables 1 and 2). For adult fertility, we obtained estimates of the environmental variance to be taken into account in the models by applying White's method based on the code proposed by Morris and Doak (2002), which allows the observation error due to sampling variation to be discounted (White 2000).
Adult survival was estimated (except in France) from turnover rates of territorial birds based on plumage ages, but was corrected by the proportion of nonobservable replacements (NOR; i.e., the replacement of an adult-plumaged territorial eagle by another adultplumaged eagle; see details in Herna´ndez-Matı´as et al. 2011a and Appendix B). We estimated the adult survival for each population as the geometric mean of annual survival probabilities. Environmental variance on survival was estimated by White's method (White 2000). In France, we used the adult survival estimates given by Herna´ndez-Matı´as et al. (2011b) that were estimated using capture-recapture multistate models on a longterm data set of individually ringed birds (1990-2008; n ¼ 423 birds).
Survival during the first three years of life (pre-adult survival) was also estimated for the French population using Herna´ndez-Matı´as et al. (2011b), and we assumed this to be the same for all populations except those in southern Spain and southern Portugal. This assumption was based on data showing that dispersing birds from several western European populations share the same dispersal areas (Cheylan et al. 1996, Cadahı´a et al. 2010, Moleo´n et al. 2011) and thus are exposed to similar environmental conditions. In southern populations, however, published and unpublished data suggests that birds tend to use nearby dispersing areas intensively and that their survival rates are considerably higher than in other populations (Balbontı´n et al. 2003, Balbontı´n 2004, Pais 2005, Consejerı´a de Medio Ambiente 2008, Balbontı´n and Ferrer 2009).
Recruitment rates were estimated as the proportion of territorial birds (a i ) and were calculated on the basis of the age-specific probability that a non-territorial bird will attempt to recruit (a i ), as estimated by Herna´ndez-Matı´as et al. (2010) for the French population using the formulas linking a i to a i in Pradel and Lebreton (1999). It was necessary to use these two rates (a i and a i ) for our models in order to correctly account for the number of potential breeders in the event of saturated populations. We assumed that a i applied to the territorial propensity of birds of a given age from any given population. We Notes: Pre-adult survival rates were S1 ¼ 0.687, S2 and S3 ¼ 0.720 for CAD, GRA, and SWP; and S1 ¼ 0.480, S2 and S3 ¼ 0.574 for the other local populations. ''Var'' expresses temporal variance after discounting sampling variance (see Methods: Vital-rate estimations). Population codes are as in Fig. 1.
See details in Appendix A. à Estimated from Catalonia. defined territorial propensity as the probability of recruiting in a population below its carrying capacity. In our models, propensity to recruit may not correspond to actual recruitment in populations that have already reached their carrying capacities (see next section). We assumed that the observed a i values in France were suitable for our purposes since recruitment is apparently not constrained in this population by territorial availability (i.e., there is relatively high adult mortality and a high proportion of non-adult territorial birds).
Natal dispersal patterns were analyzed from a sample of 52 individuals from the French and Catalan populations. These individuals correspond to all records of birds equipped with darvic bands when fledglings that are known to have been recruited into a territory during the period 1986-2010; thus, for these birds, we knew the distance between their natal and their recruitment territories. Using maximum likelihood procedures we estimated the dispersal kernel by fitting three different distributions: negative exponential, gamma, and truncated power law (see Appendix B).
Correlations among and between different vital rates can have a substantial effect on predicted population viability (Morris and Doak 2002). To analyze this possibility, we used Pearson's correlations at three levels: spatial and temporal autocorrelation, and temporal correlations between vital rates. Correlations were only computed using vital-rate estimates that had been independently estimated. These correlations were always low (see Appendix C) and so we did not account for any of these correlations in our models.

Definition of the model structure
The formulation of all our models followed the life cycle described in Fig. 2: after each breeding cycle, surviving females pass onto the next age class. According to recruitment rates, a fraction of the individuals in age classes 2, 3, and !4 attempt to recruit as territorial birds (i.e., they actually recruit if there are vacant territories available). Density dependence was incorporated into the model by setting an absolute maximum for the number of breeding pairs that each population can sustain (i.e., the carrying capacity). We decided against including declines in fertility or survival with increasing population size, since the highest survival and productivities were found in populations close to their carrying capacity (also see Carrete et al. 2006, Martı´nez et al. 2008. The carrying capacity for each population was established as the historical maximum number of pairs and, in the case of increasing populations, as 10% above the current population size. We chose this conservative threshold in our viability assessment given that growing populations act as sources in our system and available information suggests that these populations will not be able to grow at high rates in the future (see Discussion). In density-dependent models and in populations with fewer vacant territories than eagles seeking to recruit (based on recruitment rates), older eagles were dominant over younger ones. This meant that in our models, when a population reaches its carrying capacity, the age of territorial recruitment may be delayed just as has been seen to occur in natural populations of long-lived territorial birds when the number of breeding sites is limited. In models considering that local populations are connected by dispersal, individuals were allowed to disperse at the time they attempted to recruit. The probabilities of natal dispersal between populations were calculated considering the location and numbers of territories in each local population, as well as the natal dispersal kernels we estimated from observed dispersal events (see Appendix B). In the models we used these probabilities to randomly assign a dispersing individual from any given population of origin to a given population of destination. Individuals that had dispersed from one population to another (either birds attempting to recruit or already-recruited birds) were subject to the vital rates of the destination population. Once a bird becomes territorial, no further dispersal was allowed for since empirical data suggest that such behavior is extremely rare in Bonelli's Eagle (Herna´ndez-Matı´as et al. 2011a). All models accounted for demographic stochasticity in all vital rates and environmental stochasticity in adult survival and fertility, which was estimated from their observed temporal variance. Model scripts were developed in MATLAB code (see Supplement), while scripts provided by Morris and Doak (2002) were used to simulate survivals (beta distribution) and productivities (stretched beta distribution).
Evaluating uncertainty in the model structure and the natal dispersal kernel We compared the population trends observed in the studied populations over the period 1990-2005 with those predicted by six different models (Table 3, Fig. 3) accounting for several sources of uncertainty: (1) whether local populations were closed or connected by dispersal (i.e., open); (2) if closed, whether or not local populations were regulated by density dependence and, if open, (3) the shape of the natal dispersal kernel in multisite models (i.e., the exponential or power law); and (4) the behavioral ''rules'' of dispersing birds in multi-site models. The latter possibility was considered because we still lack knowledge regarding how Bonelli's Eagles decide where to recruit, and so we explicitly modeled two possible scenarios that we will refer to hereafter as the ''constrained'' and ''unconstrained'' dispersal scenarios. Under constrained dispersal, once an individual ''decides'' to settle in a given local population (on the basis of recruitment rates and dispersal kernels), it remains there, independently of whether or not a site is available. Under unconstrained dispersal, eligible breeders that cannot settle because there is no site available remain in that population for one year and if they survive, they then have priority over eagles of the same age attempting to breed for the first time when a territory becomes vacant; however, if there are still no sites, they are allowed to disperse (based again on the dispersal kernel). This unconstrained dispersal scenario implicitly assumes density-dependent dispersal as it allows individuals that do not find a territory to settle in nearby areas in which territories are available.
Observed population growth rates were estimated from census data from the studied populations at the beginning and the end of the study period. These Notes: Under constrained dispersal, once an individual ''decides'' to settle in a given local population (on the basis of recruitment rates and dispersal kernels), it is deemed to remain there, independently of whether a site is available or not. Under unconstrained dispersal, individuals attempting to breed that cannot settle because there is no site available are allowed to disperse again (see Methods: Evaluating uncertainty in the model structure and the natal dispersal kernel for details). We also used models assuming natal dispersal kernels with either exponential or power law distributions. censuses consisted of counts of territorial pairs using the methods described in Real and Mañosa (1997) and in del Moral (2006), and were independent of vital-rate estimations. We estimated model growth rates from 10 000 replicates. We calculated the expected number of breeding pairs as the mean and the corresponding 95% interpercentile range between 2.5 and 97.5 percentiles (hereafter referred to as the percentile 95% confidence interval, 95% CI) of all the replicates at the end of the period (McGowan et al. 2011). The population growth rate, or lambda, was estimated as the geometric mean of the population growth rates each year for each replicate.
Our approach to model choice relies on how likely the observed lambdas are given each model's set of predictions. For each model and for each population with an observed lambda, we binned the simulated lambda values for that population using fairly fine bins (0.01 bins). Then, we obtained the probability of the simulated lambdas being within each bin as a means of estimating the probability of seeing the observed lambda for that population, P(k i ). Finally, we obtained the loglikelihood of the observed lambdas for each model as the sum of ln (P(k i )). Models were chosen by comparing the negative log likelihoods for each considered model (Hilborn and Mangel 1997;D. Doak, personal communication).
At this point in our analyses, we split the whole population into 18 local populations (see Definition of local populations) to be able to compare the observed trends with the trends predicted by the models in the monitored populations. For those populations lacking monitoring data, we initially assumed that survival (either pre-adult or adult) was the same as in nearby monitored populations with similar environmental characteristics (Table 1). Given the low temporal variance observed in the studied populations, productivities could be based on those provided by del Moral (2006) for 2005.

Evaluation of the assumptions regarding vital rates with little available information
Based on the results from the previous section, we used the most parsimonious model in Table 3 to test whether or not our initial assumptions regarding vital rates for which little information is available were plausible. We thus considered models that assumed four additional scenarios: (1) pre-adult survival in the southern populations (CAD, GRA, SWP) is the same as that in the remaining populations, (2) adult survival in non-monitored populations is the same as the average survival in monitored populations, (3) the carrying capacity of all populations is the same as the 2005 census, and (4) a combination of scenarios (2) and (3). Again, models were run 10 000 times and their predicted lambda values were compared to those observed in the local populations based on the negative log likelihoods for each model.

Assessing population viability and the effect of the sources of uncertainty
We assessed the viability of both the whole western European population and the studied local populations. To do so, we ran models predicting the expected number of breeding pairs over the next 50 years (10 000 runs for each model) and estimated the predicted lambdas (and their corresponding percentile 95% confidence intervals), as well as the probability of population decline to thresholds of 75%, 50%, and 20% of the current population sizes. Based on our results, we only accounted for model uncertainty for natal dispersal behavior (constrained or unconstrained) when making predictions regarding the fate of these populations.
We also evaluated the effect of parameter uncertainty on model predictions. On the basis of the sample variance and the corresponding 95% confidence intervals of the vital rates (see Appendix B), we generated a set of 100 values using the beta and the stretched beta distributions (see Morris and Doak 2002) for survival and fertility, respectively. Next, we performed simulations in which each trajectory (over a period of 50 years) took into account randomly chosen values of the 100 values previously generated for each vital rate (McGowan et al. 2011). Random values of adult survival and fertility were different for each local population. Nevertheless, as in our previous models, pre-adult survival values were established as being the same for two groups of local populations, AND and SPT, on the one hand, and the remaining populations, on the other. The effects of environmental stochasticity were kept at the best estimated values (Table 4).

Identifying conservation targets
We used four simulation-based methods to identify the vital rates that had the strongest effects on population dynamics. First, we simulated for each vital rate 1000 replicates in which the values of the rates in every local population were increased by 5%, and then compared to either the absolute (sensitivity) or the relative (elasticity) increment in lambda in relation to the increment in the vital rate (Morris and Doak 2002). Second, we generated a set of simulations to evaluate the effect of the relative change in the vital rates on the absolute change in the growth rate (see also Soutullo et al. 2008). This is a measure that combines the advantages of using the same units to compare the vital rates of different magnitudes (like elasticity) with the ability to visualize the absolute effects on the growth rate of these changes. The increments in the values of the vital rates considered were À15, À10, À5, À4, À3, À2, À1, 1, 2, 3, 4, 5, 10, and 15%. Third, we performed a set of simulations to investigate the sensitivity values across the entire range of values that each vital rate may take (values ranging from 0.05 to 0.95 at 0.05 intervals). We did this separately for S1, S23, SA, and FA. In the case of FA, the number of fledged females was considered in the Results (FA 3 0.5). In this analysis, the vital rate under investigation was established as the same for all populations. Although this may seem at first unrealistic, we believe that it is a useful way of identifying critical values operating at the whole population scale that may be a threat throughout the whole population of this eagle. Finally, we carried out a simulation similar to the previous one, but with the values reached by the vital rate under analysis bounded to their actual observed range. In this case, we compared the values of lambda across the entire observed range of the vital rates represented as a percentage, where 0% corresponded to the lower limit and 100% to the upper limit of the estimated vital rates. Here, our aim was to explicitly address how the observed range of each vital rate might affect the population trend. Table 2 provides a summary of the main vital rates estimated for 12 local populations. Annual adult survival ranged from 0.868 to 0.940 (mean ¼ 0.904, SE ¼ 0.006, n ¼ 12 populations from 5851 bird-years of observations), with temporal variance ranging from 4.43 3 10 À4 to 3.18 3 10 À3 (mean ¼ 1.67 3 10 À3 , SE ¼ 0.28 3 10 À3 , n ¼ 12 populations from 181 population-years). Based on previously published information, pre-adult survival was estimated at 0.480 (SE ¼ 0.064) for fledglings (to one-year old), 0.574 (SE ¼ 0.055) for one-and two-year olds, and 0.821 (SE ¼ 0.105) for three-year olds in the French population (S1, S23, and S4, respectively; Herna´ndez-Matı´as et al. 2011b), and 0.69 (SE ¼ 0.114) for fledglings and 0.72 (SE ¼ 0.137) for one-and two-year olds in Andalusia (S1, and S23, respectively; based on Balbontı´n et al. [2003] and Balbontı´n [2004]).
Recruitment rates (a i ) were estimated at 0.161 for two-year olds, 0.680 for three-year olds, 0.934 for fouryear olds, and 1 for five-year olds and older.
In general, we found no correlation between adult survival and productivity, no temporal autocorrelation in these vital rates, and no geographical correlation of the vital rates between local populations. In the few cases in which we did detect significant correlations, significance was lost after applying Bonferroni's correction (see Appendix C).
Evaluating uncertainty in the model structure and the natal dispersal kernel The predictions of the models considering closed populations and no density dependence revealed marked differences between the predicted and the observed lambdas during the period 1990-2005 (Fig. 3). Populations located close to the Mediterranean coast (from France to southeast Spain) declined at a lower rate than predicted by these models, thereby suggesting that these populations received significant numbers of individuals from other populations. The population in northwest Spain in ARR exhibited a similar pattern. Most interior populations in central Spain behaved in a broadly similar fashion to the model predictions. Likewise, the increasing population in SWP showed size variations that were similar to the model predictions. Contrast- Notes: The vital rates presented are the assumed carrying capacity (K ), adult survival (SA), second-year (F2), third-year (F3), and adult (FA) fertility. Environmental variances are also shown for SA and FA as Var(SA) and Var(FA), respectively. Pre-adult survival rates were S1 ¼ 0.687, S2 and S3 ¼ 0.720 for AND and SPT, and S1 ¼ 0.480, S2 and S3 ¼ 0.574 for the rest. Also shown are the predicted lambdas for the next 50 years under constrained (kC) and unconstrained (kU) dispersal, with 95% confidence intervals in parentheses. See Fig. 1 for population codes. ingly, the models revealed that the populations in southern Spain in CAD and GRA increased at a lower rate than predicted by the models. In these southern areas, the models assuming negative density dependence are better fitted to the observed lambdas.
Our results strongly support the idea that the local populations studied constitute parts of a single spatially structured population, and shows source-sink dynamics (Table 3). Models assuming unconstrained dispersal were less well supported than those assuming constrained dispersal. The results also provided better support for power-law than for negative exponential dispersal kernels. The mean lambdas predicted by the best supported model in Table 3 explained 79% of the observed variance, supporting the idea that our multisite model is sufficiently realistic to be able to make reliable predictions regarding the fate of the study populations.

Evaluation of the assumptions regarding vital rates with little available information
The model assuming that adult survival in nonmonitored populations is the same as the average survival in monitored populations (scenario 2) performed slightly better than the model assuming that adult survival in non-monitored populations is the same as in similar nearby populations (scenario 0; Table 5). By contrast, the models considering that the carrying capacity of all local populations is the same as the 2005 census (scenarios 3 and 4), and the models assuming that pre-adult survival in the southern populations is equal to that of the rest of local populations (scenario 1) performed less well.

Assessing population viability and the effect of the sources of uncertainty
Predictions regarding future population dynamics varied considerably in terms of the assumptions regarding natal dispersal behavior in the eagles. Models assuming constrained dispersal gave more pessimistic population trends than models assuming unconstrained dispersal, both at whole population (lambda of 0.997 [95% CI ¼ 0.996-0.998] vs. 1.001 [95% CI ¼ 1.000-1.002]) and at local population levels (Table 4,Figs. 4 and 5). Models with constrained dispersal predicted an mean reduction of 14.1% in the number of territorial pairs for the whole population, while models with unconstrained dispersal predicted an mean increase of 6.8% over the next 50 years. The main differences between these two models appear in local populations with population growth rates below 1 if they were closed, which act as sink populations since their viability depends on the number of immigrants arriving from the southern local populations. The most isolated populations in France (FRA) and north and northwest Iberia (ANR-BUR and ARR-NPT, respectively) are predicted to decrease markedly under the scenario of constrained dispersal, but remain stable or even increase under an assumption of unconstrained dispersal. ''Sink'' populations located closer to ''source'' populations (e.g., ARA-ECLM or MUR) remain almost stable under the constrained dispersal scenario, but are expected to increase rapidly to their carrying capacity under unconstrained dispersal. Eastern populations (PVAL and CAT) are intermediate in pattern, decreasing under constrained dispersal, but increasing to near carrying capacity under unconstrained dispersal. By contrast, the populations in the south and southwest Iberian Peninsula have a positive demographic balance and behave as source populations, and the numbers of territorial pairs in these areas remain almost constant since they are close to or have already reached their carrying capacity ( Table 4).
The cumulative probability of population reduction was highly dependent on both the local demographic balance between births and deaths and the distance from source populations (Fig. 6). The populations at greatest risk are ARR-NPT, FRA, and ANR-BUR, all of which under constrained dispersal have high probabilities of reduction to 75% (0.99, 0.93, and 0.73, respectively), 50% (0.77, 0.60, and 0.40), and 20% (0.07, 0.03, and 0.05), respectively, of their current size over the next 25 years. In PVAL and CAT there is also a high probability of a reduction of 25% (0.89 and 0.63, respectively) or, to a lesser extent, of 50% (0.22 and 0.07, respectively) in the next 25 years. On the other hand, in the other local populations there is a very low risk of population reduction, either as a result of a positive demographic balance (AND and SPT) or due to their proximity to source populations (EXT-WCLM, MUR, and ARA-ECLM). Southern populations are CAD, GRA, Rest of AND, SWP, and Rest of SPT. à The most plausible model in Table 3. As expected, parameter uncertainty provided predictions that were slightly more pessimistic and revealed a marked increase in the risk of overall population decline given the wide percentile 95% confidence intervals of the estimated growth rates (model assuming constrained dispersal: 0.995 [95% CI ¼ 0.958-1.003]; Fig. 7). Indeed, the probability of overall population decline to 75%, 50%, and 20% of its current size was estimated at 0.20, 0.09, and 0.04, respectively. Similarly, the wide ampli-tude of the predicted trajectories increased the predicted risk of population decline in all local populations, except in ANR-BUR, ARR-NPT, PVAL, and FRA (Fig. 6).

Identifying conservation targets
Sensitivity analyses show that adult survival had the strongest effect on the overall population growth rate (0.085), followed by pre-adult survival (0.055), adult fertility (0.010), and carrying capacity (,0.001). Elas- FIG. 4. Predicted Bonelli's Eagle population trends for local populations (codes) over the next 50 years. Plots show the mean (solid line) and percentile 95% confidence intervals (dashed lines) of 10 000 population trajectories. We assumed in the models two types of possible natal dispersal behavior: constrained or unconstrained (see Methods: Evaluating uncertainty in the model structure and the natal dispersal kernel for details). See Fig. 1 for population code locations. ticity analyses showed a similar pattern in the case of adult (0.078) and pre-adult survival (0.034), although contrastingly gave higher values for carrying capacity (0.020) than for adult fertility (0.011; Fig. 8A).
Analyses of the relationships between the relative change in vital rates (from À15% to 15%) and the absolute change in the whole population lambda reveal that these populations would crash if there was a reduction in adult survival of 5% (Fig. 8B). Additionally, the curve representing the variation in lambda with adult survival flattens out noticeably for survival increments .5%. This can be explained by the fact that, on the one hand, we constrained the values of adult survival up to a maximum of 0.98 and that, on the other, the value to which we fixed the carrying capacity of the system does not allow the population to grow any further. By contrast, pre-adult survival and adult fertility showed a fairly constant elasticity (i.e., slope in Fig. 8B) over the range of relative changes in the vital rates studied. While the slopes of the carrying capacity were similar to adult fertility for relative decreases in its values, for increases in its values it had steeper slopes.
Our third analyses, in which we studied the sensitivities of each vital rate that was set as constant for all the local populations, showed marked declines in growth rate for values in adult survival under 0.95 and for second-and third-year survival under 0.65 (Fig. 8C). Additionally, the whole population crashed for adult survival under 0.90, for second-and third-year survival under 0.50, for first-year survival under 0.35, and for adult fertility under 0.60 fledglings per year (i.e., 0.3 females per year). On the other hand, the absolute change in lambda across the entire observed range of the vital rates illustrates much more subtle differences between vital rates (Fig. 8D). By increasing vital rates to their maximum observed value, the overall population growth rate reached values of 1.0008, 1.0004, 1.0001, and 0.9994 for adult survival, second-and third year survival, adult fertility, and first-year survival, respectively. On the other hand, by decreasing vital rates to their minimum observed value the growth rate reached values of 0.9927, 0.9952, 0.9918, and 0.9954 for the same vital rates, respectively.

DISCUSSION
Our study focusing on the western European population of Bonelli's Eagle illustrates the presence of a large-scale spatially structured population with sourcesink dynamics, in which the fate of local populations located hundreds of kilometers apart are linked through natal dispersal, and in which the survival of adults and pre-adults are the most important vital rates governing the viability of both the overall and local populations. Achieving these results was possible thanks to the implementation of a complex model, which required large amounts of demographic data collected across a large geographical area over nearly two decades; this fact highlights the importance of large-scale and long-term monitoring when attempting to assess the viability of populations of long-lived species with complex spatial structures.
Despite the extensive fieldwork carried out in this study, uncertainty regarding model structure and vital rates still exists, mainly in relation to aspects such as dispersal behavior, migration rates, and pre-adult survival that are all difficult to estimate. As in other PVAs, we thus had to make a number of assumptions, which could seriously affect model predictions (Burgman et al. 2005). We tested thoroughly the validity of our assumptions using a novel approach (but see also Brook et al. 2000, Lindenmayer et al. 2001, Schtickzelle and Baguette 2004, which consisted of (1) predefining a set of candidate models with different realistic assumptions, (2) evaluating the likelihood of observed population growth rates given model predictions, and (3) selecting the most plausible models for the PVA. Furthermore, we evaluated the effects of uncertainty in model structure and vital rates on population viability, which revealed both the importance of detailed information on dispersal processes as a means of providing more reliable estimates of the viability of spatially structured populations and the fact that accounting for parameter uncertainty markedly increases the risk of population declines predicted by the models. Our study thus highlighted the importance of thoroughly testing model assumptions and evaluating the effect of uncertainties on model predictions, thereby providing an analytic framework that can be widely used for generating reliable conservation prescriptions for a range of species and ecological contexts.

Population dynamics of Bonelli's Eagle in western Europe
Our results strongly support the idea that local Bonelli's Eagle populations in western Europe belong to a single, spatially structured population that exhibits source-sink dynamics (Pulliam 1988). Although the local populations we studied were not totally isolated patches (in the typical metapopulation sense), they did allow for spatial structure based on demographically relevant units, which we defined on the basis of their geographical proximity and environmental characteristics (see Table 1). Essentially, the populations from southern Iberia act as sources (births outnumber deaths) and, thanks to dispersal, sustain all other local populations that would otherwise decline. This mechanism occurred in all our multisite models, including the model assuming a common pre-adult survival rate for all populations (scenario 1 in Table 5), and was mainly explained by differences in adult survival between populations. The dependence on immigrating eagles appears particularly strong in the most isolated populations in north, northwest, and eastern Spain, and in France. Nevertheless, the conservation problems in these populations vary from one region to another.
Over the last three decades, the northern population in BUR has suffered the most severe decline reported for any raptor in Europe (.90% from the 1970s to the present day) and stands today on the brink of extinction (Real and Mañosa 1997). However, survival is relatively high in this population and its problems can be traced back to habitat loss due to an abandonment of traditional land practices and a crash in the numbers of its main prey species, which have led to very low productivity and high rates of territorial abandonment (Ferna´ndez et al. 1998). To a lesser extent, the population from ARR has also undergone a severe decline, likewise related to low productivity; neverthe-less, in this latter case, the fall in numbers can also be attributed to low adult survival.
By contrast, the French and eastern Iberian populations (FRA, CAT, PVAL, and MUR) have fairly good levels of productivity, but with the lowest survival rates (Carrete et al. 2005, del Moral 2006, Herna´ndez-Matı´as et al. 2011a, which agrees with the view that, due to both electrocution and direct persecution (the chief causes of mortality in this species), the highest mortality risk is to be found in highly humanized areas ). The populations from central Iberia (ARA and GUA) are more self-sustainable because they have better survival values than eastern populations, although productivity is slightly lower given the poorer local environmental conditions (Real 2004, Carrascal andSeoane 2009).
In contrast to all the others, the three southern populations (CAD, GRA, and SWP) have a notable positive balance between births and deaths, mainly because they exhibit the highest adult and pre-adult survival values. The local population in CAD has remained almost constant over the last two decades, suggesting that this population has reached its carrying capacity. GRA has increased its size, but at a much lower rate than predicted by the model that does not account for density dependence; additionally, in this latter population, most recently occupied territories are found in suboptimal areas for the species (i.e., in areas that are highly humanized, have a low availability of cliffs and are relatively cool), which suggests that the population is approaching its carrying capacity . Finally, SWP behaves rather differently since eagle numbers have increased steadily over the past two decades. This is the only population throughout the species' range in Europe to exhibit such a positive trend, which has occurred at the same time as an increase in its breeding range. Interestingly, this trend can be attributed both to high survival values and to the tendency of eagles from this population to nest in trees, which allows them to occupy new territories even in areas without cliffs (Palma et al. 2006, Beja andPalma 2008).

Effect of uncertainty on Bonelli's Eagle viability assessment
There was a marked degree of variation in viability assessments associated with uncertainty in dispersal behavior (i.e., constrained vs. unconstrained dispersal). In general, our predictions were optimistic for most populations from central and southern Iberia, which fits fairly well with the observed stable-positive trend in these populations in recent years (del Moral 2006, Moleo´n and. By contrast, populations from France and eastern and northwestern Iberia are predicted to decline under the constrained dispersal scenario, but to remain stable or even increase under unconstrained dispersal, which highlights the strong dependence of these populations on immigration and their ''true'' poor conservation status, mostly associated with high levels of human-induced mortality in territorial eagles . Indeed, recent censuses suggest that these populations have remained stable despite the fact that, in some, the main vital rates such as productivity and, especially, adult survival have progressively worsened. The most likely explanation for this previous decline in these populations is that, in the past, the number of immigrants was not sufficient to support the size of the population, as our models suggest. Additionally, it is also possible that environmental conditions in the dispersal areas (where non-territorial eagles are found) have improved in recent years as a result of the concerted efforts to correct dangerous power lines in the main dispersal areas in central (Guil et al. 2011) and southern Spain (Gil-Sa´nchez et al. 2005, Moleo´n et al. 2007) and the recent increase in rabbit numbers (the preferred prey species of this eagle) after their dramatic decline in the 1990s (Delibes-Mateos et al. 2007, Moleo´n et al. 2009, Resano et al. 2011. Cumulative probabilities of Bonelli's Eagle population decline to 75%, 50%, and 20% (black solid, black dashed, and gray solid lines, respectively) of current population sizes over the next 50 years under the scenario of constrained dispersal (black lines). The cumulative probability of 50% decline under parameter uncertainty in the constrained scenario is shown by gray dashed lines. See Fig. 1 for local population codes.
Despite these optimistic signs, when taking into account uncertainty in parameter estimates, our models predict a relatively high probability of population decline. Interestingly, although the overall population viability decreased markedly, this effect was more complex at the level of each local population. In both source populations and in sink populations located close to source populations, the risk of population decline was always greater under the uncertainty scenario. By contrast, more isolated sink populations showed a steeper increase in the short-term cumulative probability of population decline, but an earlier flattening out of the curve describing this probability. This means that in some cases the probability of decline on a 50-year horizon in models taking into account parameter uncertainty was similar (PVAL) or lower than (ANR-BUR, ARR-NPT, and FRA) in models that did not take it into account. This pattern is due to the fact that the confidence intervals of adult and, especially, pre-adult survival were broad, thereby generating more optimistic futures for a greater proportion of the simulated trajectories. Indeed, for this eagle, we still lack robust pre-adult survival estimations for most populations and have little information on the yearly variation or the spatial correlation of this important vital rate. Thus, it is worth highlighting the need to maintain current monitoring schemes, to extend them to include unmonitored populations, and to implement large-scale ringing schemes (similar to that used in France).
Besides the sources of uncertainty considered in our models, a number of other factors may also potentially affect our model predictions. We assumed that populations were structured in age classes rather than according to any other individual characteristics. Age appears to be strongly correlated with the main vital rates in Bonelli's Eagle , 2011b and thus makes possible the modeling of the different contributions of individuals in different classes to population growth. In our models, the territorial stage determined both the dispersal probabilities of individuals and their fertility rates given that only territorial eagles can breed. The territorial stage may also have an effect on the survival probability, although available data suggest that age is a better determinant of survival than the territorial stage in Bonelli's Eagle (Herna´ndez-Matı´as et al. 2011b). We also assumed that the carrying capacity of the populations was the same as the historical maximum number of pairs. Nevertheless, it is possible that environmental conditions have changed due to, for example, a decrease in prey abundance or habitat loss as a result of land abandonment or the construction of infrastructures, thereby reducing the amount of habitat available (Carrete et al. 2005). In this sense, the results suggest that our assumptions regarding the carrying capacity for the period 1990-2005 were correct. Our models did not take into account the presence of competitor species such as the Golden Eagle. Although demographic analyses suggest that competition with this larger eagle does not threaten all Bonelli's Eagle populations (Carrete et al. 2005), this output varies depending on the habitat availability for each species (Carrete et al. 2005) and on the size of the Golden Eagle population, which has been shown to decrease productivity in Bonelli's Eagle (Gil-Sa´nchez et al. 2004, Carrete et al. 2006. In fact, our elasticity analyses showed that the carrying capacity FIG. 7. Predicted population trends for the western European Bonelli's Eagle population over the next 50 years under the scenario of constrained dispersal. The mean (solid lines) and percentile 95% confidence intervals (dashed lines) of 10 000 population trajectories under parameter uncertainty (black lines) and no parameter uncertainty (gray lines) are represented.
(above all, in the case of source populations) is an important demographic parameter in determining the fate of the whole population; thus, uncertainty regarding this parameter may increase considerably the uncertainty of the model predictions.

Target vital rates in Bonelli's Eagle conservation
Adult survival had the highest values in sensitivity and elasticity (to be expected given the life strategy of this species), which thus suggests that conservation efforts should concentrate on improving this vital rate Mañosa 1997, Carrete et al. 2005). Bonelli's Eagle territories show notable environmental heterogeneity both within and between populations, which can lead to high variability in the causes and levels of mortality between territories (e.g., . Additionally, the territories with the greatest risks in the most used areas of home ranges have the highest mortality rates (Rollan et al. 2010). Thus, in terms of conservation, the best strategy would be to identify territories with the highest mortality rates, analyze the causes, and then implement measures aimed at mitigating this mortality (e.g., Tinto´et al. 2010).
Our analyses disagree with a previous study by Soutullo et al. (2008), who suggested that non-adult survival was the chief vital rate regulating this eagle's populations. This discrepancy is probably due to the fact that these authors used high natal dispersal levels within local populations based on a small sample size of dispersing birds, which does not agree with the dispersal kernels used in our study. Additionally, Soutullo et al. (2008) did not consider the role of all local populations FIG. 8. Sensitivity and elasticity analyses. (A) Sensitivity and elasticity of adult (SA) and pre-adult (S123) survival, adult fertility (FA), and carrying capacity (K ). (B) Absolute change in population growth rates (lambda) in relation to relative changes in vital rates. Lines were constructed by combining the values estimated as relative changes of minus and plus 1, 2, 3, 4, 5, 10, and 15%: SA, solid line; S123, dash-dotted line; K, dashed line; and FA, dotted line. Survival rates were not allowed to reach values .0.98. (C) Sensitivity analysis assuming that the overall metapopulation reaches the same value for the perturbed vital rate, while other vital rates vary according to their local populations. Lines show: SA, solid; S23, dashed; S1, dash-dot; and FA 3 0.5, dotted. (D) Absolute change in lambda in relation to relative changes in vital rates. Each vital-rate range is represented in a 0-100% scale (with 0% being the minimum limit of estimated values for each vital rate, and 100% being the maximum limit of estimated values); as in panel (C), we assumed that the overall metapopulation reaches the same value for the perturbed vital rate. Lines show: SA, solid (range ¼ 0.868-0.940); S23, dashed (range ¼ 0.574-0.720); S1, dash-dot (range ¼ 0.480-0.663); and FA, dotted (range ¼ 0.591-1.422).
in western Europe, which probably also accentuated the importance of survival during the dispersal period in determining the overall fate of the population. Nevertheless, in agreement with these authors, we have shown that the ratio of adult survival to non-adult survival sensitivities is considerably lower than estimated in previous studies in which no spatial structure was considered, i.e., in which local populations were assumed to be isolated (see Real and Mañosa 1997), and therefore, our results reveal that an increase in juvenile and immature survival will also have a very positive effect on the dynamics of Bonelli's Eagle populations. In this sense, we estimated under the constrained dispersal scenario that the whole population would remain stable under the scenario of a 3% increase in adult survival or a 10% increase in pre-adult survival.
Our results also emphasize the fact that natal dispersal plays a key role in guaranteeing the persistence of sink populations. Given that all individuals pass through a nomadic phase (i.e., a dispersal period) in which they use specific dispersal areas intensively (e.g., , Cadahı´a et al. 2010, Moleo´n et al. 2011, the implementing of suitable conservation measures in these areas aimed at mitigating mortality (above all, avoiding electrocution by power lines and illegal persecution) would appear to be a key issue in securing suitable survival rates amongst dispersers and in guaranteeing the long-term persistence of these populations. Ideally, conservation measures should be able to cope with the fact that the location of these dispersal areas may change over time.
Our final analysis revealed that differences between vital rates in their impact on the overall fate of the population are smaller when the elasticity analyses are restricted to their observed range of variation (also see Wisdom et al. 2000). Again, increasing adult survival to the maximum observed values had the greatest impact on the overall growth rate of the population; nevertheless, in contrast to conventional sensitivity and elasticity analyses, decreasing productivity to the lowest observed values had the most negative impact on the population. Therefore, it appears necessary for Bonelli's Eagle conservation strategies to work towards mitigating low productivity levels, which in our populations are mainly caused by habitat loss, decreases of main prey species, and disturbances in the breeding areas caused by either human outdoor activities or the construction of new infrastructures.
Finally, and given that the carrying capacity may play an important role in the fate of a population, it would be desirable to implement conservation measures aimed at preserving or restoring the suitability of territories as a PLATE 1. Female Bonelli's Eagle ringed as a chick in Catalonia (CAT) in spring 2008 and recruited as territorial in 2011 in southern France (FRA), where this photo was taken. This bird and its mate disappeared from the territory in spring 2012 and have not been observed since. The fact that five territorial eagles have been replaced in the past three years at this territory suggests that human-induced mortality, possibly persecution, lies behind the (likely) death of this eagle. Photo credit: Alain Marmasse. means of maintaining or even improving the carrying capacity of populations. In practice, this means allowing southern populations to produce a suitable number of dispersing eagles to sustain other populations and, in northern populations, encouraging unused territories to be reoccupied while avoiding the abandonment of currently used territories.

General implications for population viability analysis
Our study highlights the challenges to be confronted when undertaking viability analysis of long-lived species occurring in spatially structured populations over broad geographical areas, and thus provides a framework for multi-site PVA application that may be useful under a wide set of circumstances. In particular, our study revealed the importance of testing model assumptions and evaluating the consequences of uncertainty, which are critical when generating useful conservation guidelines for endangered species.
The approach adopted herein for testing different model assumptions using past population trends appears to be reliable for evaluating the plausibility of considered scenarios (Brook et al. 2000, Lindenmayer et al. 2001, Arnold et al. 2006, and therefore, we believe that this step should always be performed if data on past population trends are available. We used the observed and predicted population growth rates instead of other measures of viability, which may be particularly useful for species with a low short-term risk of extinction but vulnerable current population sizes, as is the case of many endangered long-lived species. To do so, growth rate estimates for populations and the vital rate data used in the models must be estimated independently. When knowledge of the species' biology permits, we also advocate predefining a set of feasible models (considering different assumptions), and then selecting the most plausible model using likelihood-based methods, thereby providing a better theoretical basis and understanding of the assumed processes governing population dynamics. Even so, uncertainty in the model assumptions may still be high and these sources of uncertainty must be accounted for in order to increase the soundness of the PVA predictions.
In our analyses, uncertainty in dispersal components had a strong effect on model predictions. Despite being a crucial element in most population models (MacArthur and Wilson 1967, Levins 1969, Pulliam 1988, Hanski 1999, there are still significant deficiencies in our understanding of behavioral and spatial variation in dispersal (MacDonald andJohnson 2001, Clobert et al. 2009). This gap in our knowledge means that only rarely in population viability assessments has much attention been paid to dispersal (but see McCarthy et al. 2000, Lindenmayer et al. 2003, Schtickzelle et al. 2005, Bonnot et al. 2011). In the case of long-lived species, in which dispersal may occur across large geographical ranges and over relatively long periods of time, the quantification of dispersal is particularly challenging and, conse-quently, this factor has been misrepresented in empirical studies in the literature, illustrating the role of dispersal in shaping the spatial structure of populations (Oro et al. 2004), especially over broad spatial scales. Here, we show that, along with survival, dispersal may have a profound effect on the population dynamics of longlived species such as Bonelli's Eagle. Unlike in metapopulations of short-lived or highly mobile long-lived species, where local extinctions are prone to happen in habitable patches over observable timescales, our sink populations remain stable (or decrease slowly) and their predicted extinction probabilities are relatively low in the mid-term (at least when the uncertainty parameter is not taken into account); in other words, the ''rescue effect'' works continuously and local extinctions are rarely seen to occur.
Our results also emphasize the fact that the overall dynamics of our population are powerfully determined by both the shape of the dispersal kernel and the behavior of dispersing individuals. The dispersal kernel that we fitted with observed natal dispersal distances was better matched to a negative exponential function, although the evaluation of past variation in local population size strongly supports the idea that dispersal followed a power-law function; this is understandable if we bear in mind the fact that long-distance dispersers may be misrepresented in our sample because the monitoring effort was lower in more distant populations. Furthermore, these results highlight the crucial importance of long-distant dispersers in rescuing sink populations and exemplify the fact that assumptions made on the spatial variation of dispersal may substantially change the predictions deriving from population models. In addition, the behavioral rules governing dispersal decisions may also have a strong effect on population dynamics (Clobert et al. 2009). If floaters (i.e., potential breeders) change their dispersal behavior during recruitment based on a perception of territory availability, they are not expected to remain in the same patch if they have little chance of recruiting, as occurs in saturated populations. The aim of our unconstrained dispersal was to recreate this scenario, whose consequences are that the bulk of floaters, that would otherwise remain in the source population under constrained dispersal, now move into the surrounding sink populations in order to recruit. Thus, they also increase the number of floaters in these populations and thus even allow very distant sink populations to receive a significant number of immigrants. In other words, in our system, density-dependent dispersal increases the persistence probabilities of sink populations and ultimately improves the overall viability of the whole population.
Finally, we have illustrated that accounting for parameter uncertainty may also cause a marked increase in the probability of population decline, even when PVAs rely on fairly detailed demographic data , McGowan et al. 2011. Assessing the effects of uncertainty sources on the viability of target populations enhances the detection of those parameters for which available data is insufficient and on which future research and/or monitoring should be concentrated, and also provides more reliable information on the risk associated under different monitoring and management scenarios ). In conclusion, we emphasize that it will always be necessary to take into account the main sources of uncertainty when looking to increase the reliability of population viability analysis, particularly when the status of a species of conservation concern is being evaluated. In the case of endangered species showing a fragmented distribution, it will be of great importance to account for uncertainties in dispersal processes if we are to guide management and conservation decision-making appropriately.