Population viability analysis of the endangered Dupont’s Lark Chersophilus duponti in Spain

Steppe lands in Europe are critically affected by habitat loss and fragmentation, and hold over 50% of IUCN Red List bird species in Europe. Dupont’s Lark is a threatened steppe-specialist passerine whose European geographic range is restricted to Spain, with less than 2000 pairs and an annual population decline of − 3.9%. Its strongly fragmented habitat leads to a metapopulation structure in the Iberian Peninsula that includes 24 populations and 100 subpopulations. We present an updated Population Viability Analysis based on the latest scientific knowledge regarding distribution, population trends, breeding biology and connectivity. Our results predict metapopulation extinction in 2–3 decades, through a centripetal contraction process from the periphery to the core. The probability of extinction in 20 years was 84.2%, which supports its relisting to Endangered in Spain following IUCN criteria. We carried out a sensitivity analysis showing that some parameters, especially productivity and survival of adults and juveniles, help to increase metapopulation viability. Simulation of management scenarios showed that habitat restoration in a subset of key subpopulations had a positive effect on the overall metapopulation persistence. Translocations of a limited number of individuals from source to recipient locations may help to rescue the most endangered subpopulations without reducing the global time to extinction of the metapopulation. In addition, we identified the most critical areas for action, where local populations of the species are prone to extinction. This work suggests that the viability of the Dupont’s Lark metapopulation could be improved and its risk of extinction reduced if urgent and localized conservation measures are applied. In the short-term, habitat loss and fragmentation due to ploughing, reforestation and infrastructures implementation in Dupont’s Lark habitat must be avoided. Habitat restoration and translocations could help to avoid imminent extinction of critical subpopulations. Restoration of extensive grazing is recommended as the most effective way to achieve the long-term conservation of Dupont’s Lark in Spain.

www.nature.com/scientificreports/ Subpopulations' extinction occurred radially, showing a progressive range contraction from the most peripheral subpopulations to the metapopulation core (Supplementary File S1). During the first decade the model predicted a strong decline in population size with a 100% persistence probability of the metapopulation, with a steep decline in persistence probability after year 12 (Supplementary Figure S1). Critical areas. Nineteen subpopulations became extinct in the first two years (T mean < 2, Supplementary   Table S2). These subpopulations had low initial population sizes (mean N 0 = 1.29; min. = 0; max. = 4 individuals) and a mean linear distance to the metapopulation core of 451.72 km (min. = 42.4, max. = 449.1 km). Those subpopulations overcoming a T mean of 10 years (n = 29, Supplementary Table S1) had higher initial population sizes (mean N 0 = 156.90; min. = 1; max. = 1,125 individuals) and were located closer to the metapopulation centroid (mean = 79.47, min. = 10.2, max. = 154.1 km). The Pearson's test showed a significant correlation between T mean and N 0 (r = 0.46, p < 0.001), and T mean and the distance to the metapopulation centroid (r = -0.44, p < 0.001).
According to population decline rates, two thirds of the subpopulations (69%) showed negative values of population growth (r). Among the subpopulations with a stronger decline were those with higher initial population size and located in the core of the distribution (see Supplementary Table S2 for a summary of the 10 subpopulations with the strongest annual decline, and Supplementary Table S1 for the complete list). Sensitivity analysis. All the parameters evaluated markedly affected the metapopulation viability. A range in productivity from 1.2 to 1.9 offspring/clutch produced a variation in probability of extinction in 20 years from 1 to 0, respectively, with a stronger effect for the intermediate values of 1.4-1.7 (Fig. 1A). In the base model with a value of productivity of 1.5 offspring/brood we obtained a probability of extinction in 20 years of 84.2%. An increase to 1.6 offspring/brood showed a reduction of that probability to 60.6%, while a reduction to 1.4 offspring/brood would increase the probability of extinction to 96.8%. The effects on r, T mean , and T max were similar to P 0 (20) (see Supplementary Table S3 and Supplementary Figs. S2-S4).
In relation to the percentage of breeding females, the metapopulation showed a high sensitivity. Our base value was of 100% breeding females due to the species sex ratio, obtaining 84.2% probability of extinction in  www.nature.com/scientificreports/ 20 years. A reduction to 90% breeding females predicted an increase in that probability to 95%, and reaching the threshold of 75% breeding females would drive the species to extinction in 20 years (Fig. 1B). The metapopulation showed a similar sensitivity to reductions in survival of the different sexes and age classes. We used the values of 48% adult male survival, and 31% for females and juveniles, obtaining 84.2% probability of extinction. Reaching the lowest thresholds of 30% (adult males), 20% (adult females) and 25% (juveniles) predicted a 100% probability of metapopulation extinction in 20 years (Fig. 1C-E). The effect of increasing survival rates differed among classes. Increasing male survival had a positive effect reducing extinction probability, though always above 0% (P 0 (20) = 45.4% with a 100% survival of adult males; Fig. 1C, Supplementary Table S3), while increases in female and juvenile survival (respectively, thresholds of 60% and 40%) avoided metapopulation extinction in 20 years (Fig. 1D-E).
Finally, according to dispersal survival, we used a value of 50% in the base model to obtain a 84.2% probability of extinction in 20 years. Without dispersal survival, the model predicted a probability of extinction of 99.4%, while a survival of 100% would result in a P 0 (20)  Habitat restoration. The simulation of habitat management had a positive effect on the metapopulation overall values of T mean , r, P 0 (20) and T max (Supplementary Tables S4, S5). All scenarios of habitat management increased the viability of the metapopulation compared to the base model. Mean time to extinction increased with the number of subpopulations and years of management ( Fig. 2A). Two-way ANOVA showed a stronger effect of increasing the number of subpopulations (F(3, 4391) = 225.328, p < 0.001) than the time period (F(2,1) = 0.077, p = 0.926), and the Tukey's post hoc test showed significant differences (p < 0.001) among all levels of the factor (0, 5, 10 and 15 subpopulations managed). The intensity of management increased the metapopulation mean time to extinction (Fig. 2B). One-way ANOVA showed a significant effect (F(5, 18,416) = 617.2, p < 0.001) and the post hoc test revealed differences between all levels: no management, 1, 5, 10, 15 and 20% increase in productivity and K.
Similarly, habitat restoration simulations increased metapopulation growth rate (Supplementary Figure S5), reduced the probability of extinction in 20 years (Supplementary Figure S6) and increased the maximum time to extinction (Supplementary Figure S7).
Translocation. The translocation of individuals produced a slight decrease in the overall metapopulation mean time to extinction (Fig. 3A). Two-way ANOVA showed a significant effect of harvesting (F(6, 138) = 4.372, p < 0.001). The post-hoc Tukey test revealed a significant difference with the 'before translocation' scenario appeared when harvesting 10 males and 10 females (diff. = − 0.414, p = 0.001). The rest of the harvesting levels showed no significant differences with the 'before translocation' scenario.
Harvesting also showed a significant effect in donor subpopulations (F(5, 7118) = 199.6, p < 0.001), which in general reduced their mean time to extinction (Fig. 3B). Significant differences with the 'before translo-   Factors affecting population viability. The GLM showed a significant influence of isolation and climate on metapopulation persistence (Table 1). In the first case, the distance to the metapopulation centroid had a negative effect on the mean time to extinction (Table 1), meaning shorter extinction times as the subpopulation is further from the core, supporting the aforementioned range contraction process. According to climate, the model revealed shorter times to extinction with increases in annual precipitation trends ( Table 1), suggesting that subpopulations with increases in aridity along the period 2000-2018 could persist over more years.

Discussion
This work offers updated information on the viability of the Spanish metapopulation of the Dupont's Lark Spanish metapopulation and highlights the dramatic conservation status of the species in Europe. In summary, our work suggests that, if there is no change in present-day driving forces, the species will reach extinction in a few decades, through a process of centripetal contraction from the periphery to the core subpopulations of the . Effects in mean time to extinction (T mean and 95% confidence interval) in the metapopulation, donor and recipient subpopulations after the translocation program. Donor subpopulations were those with more than 100 males (n = 7). Recipient subpopulations were those with the most unfavorable situation in the PVA results (shortest mean time to extinction), two scenarios were simulated: seven recipient subpopulations close to the metapopulation core and seven recipients distant from the core. For each scenario, six different harvest alternatives are offered: 1 + 1, 2 + 2, 3 + 3, 5 + 5, 10 + 10 and 10 + 6 males/females, this last one adjusted to the species sex ratio. In all of them, movements were carried out during 3 consecutive years. Harvested individuals were introduced randomly in recipient subpopulations (Supplementary Tables S6, S7).  (Fig. 4, Supplementary File S1). However, the species is sensitive to management measures, especially to increases in productivity and survival of adults and juveniles, which offers opportunities for its conservation. Finally, subpopulations at high risk of extinction may be rescued via translocation of wild individuals from source subpopulations. Our model shows that small and more isolated subpopulations have a shorter mean time to extinction. Subpopulations facing an imminent extinction (T mean < 2 years) according to the model predictions had extremely low population sizes and longer distances to the nearest subpopulation and metapopulation centroid (Supplementary  Table S2). These patches are more vulnerable to stochastic demographic and environmental processes, as well as genetic drift 40,64 . The core areas (Layna, Altos de Barahona, Altiplano de Teruel, amongst others; Supplementary  Table S2) showed longer times to extinction but also the most pronounced population decline (higher negative   Table 1. Summary of the GLM coefficients to evaluate the effect of geographic and climatic variables on the mean time to extinction (T mean ) of the metapopulation. Variables initially included were: distance to the metapopulation centroid (Dist. C), distance to the nearest subpopulation (Dist. N), minimum January (T jan ) and maximum August (T aug ) temperatures (mean values of the period 2000-2018), annual precipitation (mean 2000-2018) and the trends of the three climatic variables (as the slope of a linear regression of the period 2000-2018). Only the significant variables remaining after sequential model simplification are shown (see details in text).  www.nature.com/scientificreports/ r values). An important fraction of the population inhabits small patches with a high risk of extinction, where local reproduction may be insufficient to compensate for mortality and whose persistence would rely on the immigration of individuals from nearby source patches 65 .
The available data reveals that recolonizations are an extremely rare event in Dupont's Lark, and local extinctions seem to be irreversible. Our monitoring in the Iberian Range has revealed the gradual extinction of small patches with only one recolonization in the post-2000 period 26,36 . At a metapopulation scale, among the 23 subpopulations extinct since the year 2000 only one recolonization has been recorded 66 . Extinct or unoccupied areas may play an important role as stepping stones, which were included in the dispersion matrix of the PVA model, and could play a relevant role in the maintenance of the overall metapopulation, as they increase the connectivity among patches 49 . Our results showed that large and connected patches are crucial in the conservation of the Dupont's Lark, which is in line with the conceptual framework of metapopulation theory. We may have expected higher recolonization and persistence probability 17 in larger and connected patches 40,64 however, our results do not support a Levins model of colonization-extinction balance since recolonizations are extremely rare events, while extinctions are permanent. The process seems to be a radial contraction of the distribution range towards the metapopulation core, probably associated with an increase in landscape fragmentation as the distance to the center increases in a source-sink pattern revealing a centripetal contraction process from the periphery to the core subpopulations of the Iberian Range and Ebro Valley 49 ).
We suggest that the range contraction process, acting first on smaller and peripheral patches, is occurring both at local and metapopulation scales, which offers little chance for a patch to be reoccupied once it has gone extinct. This is not incompatible with recent findings of new settlements of the species (explained by the use of new monitoring methods 67 ), which we attribute to previous insufficient sampling effort rather than true recolonizations. Our results at the metapopulation level are similar to those obtained in a recent analysis of Dupont's Lark occupancy at a smaller scale 55 , which detected a hierarchical pattern in habitat use, with metapopulation-scale factors (connectivity and patch size) explaining most of the variance versus landscape and microhabitat-scale variables. While genetics suggests scarce gene flow 39 among Spanish subpopulations, individual movements could be occurring to some extent according to both the data available about juvenile dispersion 68 and connectivity analysis 49 . Several efforts to analyze gene flow and structure in Dupont's Lark have been made in the recent years 39,40 . However, the rapid extinction process of the species requires more genetic research to quantify effects of isolation and trends toward homozygosity, especially in small and peripheral populations.
Management and conservation measures should be guided to enhance three crucial aspects: interpatch dispersion, minimum patch size and habitat quality. We highlight the convenience of a two-level management approach, concordant with both types of critical areas (Supplementary Table S2). First, the core subpopulations showed strong declines, with high negative r values. In these areas, measures directed at maintaining or increasing habitat quality and carrying capacity should benefit productivity and slow down this decline. Second, the small and most isolated subpopulations must be protected against imminent extinction, to urgently reduce the contraction process of the metapopulation. Here, management should aim at guaranteeing a minimum patch size and interpatch connectivity. In this sense, connectivity analysis 49 could help identifying key subpopulations with a potential effect as stepping stones and connectors for their neighbors.
As a general consideration, actions that improve productivity and adult and young survival should be implemented, at least, in the critical subpopulations. Our results showed that increasing the number of managed subpopulations offered better results than lengthening the time period of management (Supplementary Table S4, Fig. 2A). Predator control has been implemented as a useful and cost-effective tool to increase breeding success of ground nesters 69,70 (but see also [69][70][71][72][73] ) and could be evaluated for Dupont's Lark. Food supplementation enhances population viability in endangered species 74,75 and could be introduced as a short-term action in critical subpopulations. Particularly, sown sheep dung could simulate grazing in terms of coprophagous availability 76 . Habitat management showed a high potential for improving the viability of the metapopulation, as the action over a small subset of subpopulations had an effect on the whole metapopulation persistence (Fig. 2). Habitat restoration has been successfully carried out in Dupont's Lark peripheral subpopulations successfully 50 and may be implemented in the short term in critical areas, including shrub clearing and tree removal. In the medium term, the reintroduction of extensive grazing is the most recommended measure, as it preserves habitat quality and food availability in shrub steppes.
Translocations offer a good opportunity for management and should be explored as a useful tool for shorttime rescue of small subpopulations at imminent extinction risk, but must be applied with caution, as the simulations showed that increasing the number of yearly harvested individuals up to 20 (10 males and 10 females) could promote a significant reduction in the metapopulation time to extinction. Moving a low number of individuals can balance the conservation of critical subpopulations with metapopulation persistence. Special attention must be paid to donor subpopulations, as the effects in the reduction of their time to extinction were significant in almost all scenarios of harvesting. The effects were slight up to 10 individuals (5 males, 5 females) harvested yearly (Fig. 3B), but with strong improvements in recipient subpopulations. Thus, this measure should be considered as a short-term way to reinforce subpopulations facing imminent extinction or rescue recently extinct ones. More research is needed regarding translocations, as they could play an important role in mitigating inbreeding and increasing allele diversity in such a fragmented population, considering that more than 50% of the Spanish subpopulations are below 5 individuals 50 and gene flow is lacking 39 .
As a medium-to-long term objective, management actions should be considered only as part of a broader framework of self-viability. A reliable conservation strategy to ensure the persistence of Dupont's Lark in Europe should include the restoration of extensive grazing in natural steppes 56 , as it shapes plant structure and provides food through its positive effect on arthropod availability 55,56 . Habitat loss and fragmentation, both due to abandonment or land use changes (ploughing, reforestation or infrastructure installation), must be urgently stopped to slow down species decline. The results of modeling geographic and climate effects on the metapopulation persistence also highlighted the radial contraction process towards the core regions. Our PVA results forecast a future Dupont's Lark distribution restricted to the Iberian Range and the Ebro Valley in about one decade if the process is not reverted (Fig. 4). The higher resistance to extinction of those subpopulations with a more pronounced trend toward aridity ( Table 1) could suggest that this steppe species would benefit from a future climate change scenario, but this must be considered with caution as it may be associated with a geographic effect 77 . While steppes and other sub-desert systems could take advantage of a future increase in aridity (as some records support 78 ), steppe bird species are usually strict specialists, as is the case of Dupont's Lark, making them extremely susceptible to habitat changes, whatever the direction. Thus, more research about effects of future climate change scenario is needed in this sense.
Finally, according to the IUCN Red List Categories and Criteria (IUCN 2012), our results point to the necessity of updating the species status to Endangered in Spain, following the criterion E (quantitative analysis indicating the probability of extinction in the wild to be ≥ 20% in 20 years or 5 generations, whichever is longer). Also, management measures aimed at improving the conservation status of Dupont's Lark may also benefit other declining steppe birds, such as Little bustard 9 , Stone curlew 10 , Northern wheatear 11 or Ortolan bunting 12 , among others, which have experienced a general negative trend in recent years in Spain and Europe [13][14][15] , as well as the steppe habitat itself.

Methods
The Species and study area. The Dupont´s Lark is a relatively well monitored species. Broad scientific information has been produced in the last two decades regarding many aspects of its ecology and biology. Our study area includes the totality of the known European distribution range, which is restricted to continental Spain, currently comprising 24 populations and 100 subpopulations 49 . Criteria for the definition of these entities are based on the species' movements 49 and on the map of georeferenced locations and associated distance buffers. As a result, subpopulations are defined by clusters of observations separated by less than 5 km (considered resident movements) while populations are groups of subpopulations separated by a maximum distance of 20 km (considered the probable juvenile dispersive distance), see details in García-Antón et al. 49 .
The core of the metapopulation comprises two large geographic regions: the Iberian Range and the Ebro Valley (Fig. 5), which are merged into one, single population following the definition criteria. We consider as a starting point the metapopulation structure recently published in the National Conservation Strategy of the species 50 . Some of the 100 subpopulations considered have extremely low population sizes and are expected to reach extinction during the following years.
Risk of extinction has been previously assessed at a local scale in the Ebro Valley 79 and in the Iberian Range 80 , and for the whole Spanish population 81 . However, new data have been published since then, and monitoring programs have increased capture-recapture data and updated the basics on distribution 48 , connectivity and dispersal 49,68 , breeding biology 82 and population trends 36,83,84 . This allowed us to update critical parameters such as longevity and mortality, productivity, initial population size and carrying capacity, and re-evaluate the population viability of Dupont's Lark at a national scale (Table 2). We used data from different monitoring programs carried out during the last 20 years, including the II National Census (2004-2007 85 ) and from monitoring and research projects performed during the last 10 years, which included seasonal censuses, territory mapping, capture-recapture programs, radiotracking and record of behavioral data. We gathered 16,676 georeferenced locations of Dupont's Lark from 2000 to 2018, including our own data and that from public administrations, other research institutions and private ornithologists (see also 48 ).

Population viability analysis.
We used the stochastic simulation software VORTEX version 10.2.5.0 86 and run individual-based models with 500 iterations to account for demographic, environmental and genetic stochasticity. All values in VORTEX were kept at their default setting except those indicated below. To evaluate long-term metapopulation persistence we projected the model to 50 years, although we evaluated the species status after the first 20 years to fit IUCN criteria. Individual-based models rely on independent outcomes of the fates of individuals, so they include stochasticity in demography (random variation in births and deaths), environment (variation in disease, predation, food availability, weather, natural catastrophes) and genetics (decrease of fitness due to inbreeding depression and loss of diversity due to random genetic drift 87 ). All PVA models were carried out at the subpopulation level (n = 100) using the above mentioned current Spanish metapopulation structure 49 (Fig. 5). We defined the species extinction as the absence of at least one sex.
We included inbreeding depression to introduce evolutionary processes in the model. Here, we use the default value of 6.29 lethal equivalents suggested by Lacy and Pollak 86 , which represents the combined effect of inbreeding on fecundity and first year survival reported by O'Grady et al. 88 . The correlation in the environmental variability among populations was fixed at an intermediate value of 50% (see also 81 ). The value of the correlation between breeding and survival was also maintained at the default value of 50% 86 . Similar to Suárez and Carriles 81 , we included one regional catastrophe with a frequency of 5% of the years in which 5% of females would not breed and survival would suffer a reduction of 5%.
We built a base model considering the most plausible value of each population parameter regarding the available current information (see below and Table 2) and estimated the probability of metapopulation extinction after 20 years. Then, we ran alternative models to assess the impact on the species of critical parameters and potentially Scientific Reports | (2021) 11:19947 | https://doi.org/10.1038/s41598-021-99125-y www.nature.com/scientificreports/ handled through management plans and habitat restoration: productivity; breeding females, survival of adults and juveniles, and survival during dispersion (see below). For each of these sensitivity models we varied the specified parameter and kept the rest constant, to avoid interactions and obtain comparative results, and all the simulations were run with 500 iterations. We used the variation range that provided a response from 0 to 100% in the probability of metapopulation extinction in 20 years. For each model we estimated the mean year growth rate (r), mean population size each 10-year period (N t ), probability of extinction each 10-year period (P 0 (t): equivalent to the proportion of the 500 iterations in which the population is extinct or remains extant), median time to extinction (T med ), mean time to extinction (T mean ) and maximum time to extinction (T max , equivalent to the year in which all 500 iterations result in metapopulation extinction).

Demographic parameters.
Regarding dispersal, Dupont's Lark adults appear to be resident 85,87-92 while juveniles are assumed to be the dispersive fraction of the metapopulation 49,68 . Our base model assumed that only first year individuals of both sexes disperse and included 50% of dispersal survival. To avoid the use of a fixed arbitrary dispersal rate, we used the probability of connection between paired subpopulations considering juvenile movements of 20 km 49 and a proportion of 10% of dispersers 79 .
In relation to longevity, our capture-recapture work during 2000-2018 in the Iberian Range showed a maximum lifespan of five years, despite some rare events of longer longevity (one 8-year-old bird recaptured 85 ; 15 years used in Laiolo et al. 79 ). Despite recent evidence of extra-pair copulation in this 82 and other alaudids 93 , monogamy is considered the reproductive system in the Dupont's Lark [77][78][79][80][81] . The age of reproduction was considered from year 1 to 5, concordant with the lifespan. Number of clutches per year was three, with a maximum of 5 eggs per clutch 82 . Sex ratio at birth was 50%.
As adult sex ratio is strongly male biased (0.79 following Vögeli et al. 51 ; 0.61 following Suárez et al. 52 ), the base model considered that all females were available to breed 85 . We introduced a random 10% SD in percent of www.nature.com/scientificreports/ breeding females to account for environmental variation. In this way, the value varied in each of the 500 iterations of the model within that range. We considered 80% of breeding males, due to the 20% of floating males registered by radiotracking 81 . We used distribution of broods as in previous works 79,81 , which estimated this parameter using data available on the Skylark 94 , in which 6, 29 and 65% of females lay 1, 2 and 3 clutches, respectively. Productivity of 1.5 ± 0.5 offspring/brood was based on the most recent study available 82 .
Regarding mortality rates we averaged the two data available, mean male mortality used in Laiolo et al. 79 and in Vögeli et al. 90 (54% and 50% respectively; mean = 52%). We estimated adult female mortality at 69% using the same method as Suárez and Carriles 81 . We considered that juveniles of both sexes suffer higher mortality rates than adults, and accordingly used the highest of both values available (69%). To account for stochasticity, as well as for percentage of breeding females, we introduced a variability of 10% SD in mortality of adults and juveniles.
Regarding initial population size, we used the most recent census data for each of the 100 subpopulations. We inferred the current number of males in 2020 applying the annual change rate of − 3.9% 36 . We then applied Table 2. Values of the population parameters used to build the base model in VORTEX. The base model was built considering the most plausible values given the available current information. Alternative scenarios and iterations were built to carry out the sensitivity analysis and simulations of translocation and habitat restoration programs (details in the text).  85 to obtain the number of females and total population size per subpopulation in 2020 (see initial population size per subpopulation in Supplementary Table S1). The carrying capacity (K) of each subpopulation was estimated based on the mean density of 1 individual/10 ha obtained by radiotracking 43 . To calculate the surface of adequate habitat within each subpopulation we first intersected the 16,676 observations from 2000 to 2018 with all available CORINE Land Cover layers in this period (2000, 2006, 2012 and 2018), and assigned to each point the most updated land category according to the observation date. We selected as the preferred categories of the species those accounting for 95% of the observations (Supplementary Table S8). We then extracted those categories from the last CORINE layer (2018) to obtain the current distribution of the species' habitat and clipped them with the 100 subpopulations. This method implies a risk of surface overestimation, as it can include CORINE categories with little probability of occupation by the species (low % of total observations, see Supplementary Table S8), but with high surface within the subpopulation. To control for this, we applied a correction weighting the obtained surface by its habitat quality (mean probability of presence of the 1 × 1 km cells included in it 48 ), reducing carrying capacity estimates in subpopulations with a large occupation of such land categories with low probability of presence. Finally, we discarded the surface with a slope over 15%, which is strongly avoided by the Dupont's Lark [39][40][41][42][43][44][45][46][47] and those habitat patches smaller than 20 ha, matching the minimum area requirements of the species 91 . For a similar procedure, see García-Antón et al. 48 . We included a 10% SD in K to account for environmental stochasticity. Finally, we considered the current habitat loss described recently for the Spanish metapopulation, which is considered one of the main threats for the species, and introduced in the base model an annual change rate of -3.9% in K 36 .
Sensitivity analysis. We varied the following parameters to evaluate their effect on the metapopulation probability of extinction. For each of them, we used the range of values that allowed an observation of a 0-100% response in the probability of metapopulation extinction in 20 years: productivity (variation range of 1.2-1. Critical areas. We used the base model to identify subpopulations of a higher concern and in need of special attention, following two criteria: imminent extinction and fast population decline. In the first case, we considered critical those subpopulations reaching extinction in the first 2 years (T mean < 2 across the 500 iterations). Regarding population growth, we selected the subpopulations with a fastest annual decline (higher negative mean r values). We used a Pearson's correlation test to evaluate the association between the mean time to extinction of each subpopulation and its initial population size and distance to the metapopulation centroid.
Habitat restoration. We tested for the effects that a habitat management program applied in a subset of key subpopulations would have on metapopulation persistence. Assuming that the program would increase habitat suitability (more available foraging areas and/or safer nesting sites), we increased carrying capacity and productivity to simulate successful habitat restoration. We first generated nine different scenarios to evaluate the effect of spatial (number of subpopulations managed) and temporal (number of years of application) scale of the program on the overall metapopulation parameters. The selected subpopulations were those with a higher initial population size, and the nine scenarios tested were: (i) the 5 largest subpopulations during 3 consecutive years; (ii) 5 subpopulations/5 years; (iii) 5 subpopulations/10 years; (iv) 10 subpopulations/3 years; (v) 10 subpopulations/5 years; (vi) 10 subpopulations/10 years; (vii) 15 subpopulations/3 years; (viii) 15 subpopulations/5 years; and (ix) 15 subpopulations/10 years. We used two-way ANOVA test to evaluate whether the number of subpopulations and years of management had a significant effect on the metapopulation persistence, using mean time to extinction as the response variable and number of subpopulations and years as factors. Tukey's Honest Significant Difference test was employed to find differences among groups.
Based on the scenario that we considered most reliably applicable (10 subpopulations managed for 3 years), we tested the effect of increasing the intensity of the restoration on the selected subset of subpopulations. We assessed five different alternatives: 1, 5, 10, 15 and 20% increases in both productivity and carrying capacity. We evaluated the effects with one-way ANOVA and a Tukey's post-hoc test.
Translocation. We simulated the effects of a translocation program on the metapopulation viability, as well as on donor and recipient subpopulations. We selected donor subpopulations as those with more than 100 males (n = 7, Supplementary Table S6), all of them in the core area of the metapopulation. As recipients, we selected those seven subpopulations (same number as donors) with the most unfavorable situation based on the PVA results (shortest mean time to extinction). We simulated two different scenarios: recipients close to and distant from the metapopulation core (Supplementary Table S6).
We used a realistic time period of translocating individuals during three consecutive years. For each scenario, we simulated different harvesting alternatives varying the number of individuals moved from donor to recipient subpopulations yearly: 1 + 1, 2 + 2, 3 + 3, 5 + 5, 10 + 10 and 10 + 6 males/females, this last one adjusted to the species' sex ratio. Movements were directed from donor to recipient subpopulations and individuals were randomized to increase allele diversity.
We carried out an ANOVA test to evaluate differences in the mean time to extinction of the metapopulation before and after the translocation, for the two different scenarios. We used harvesting alternatives (6 classes + no translocation) and close/distant scenario as factors. The response variable was the mean time to extinction Scientific Reports | (2021) 11:19947 | https://doi.org/10.1038/s41598-021-99125-y www.nature.com/scientificreports/ obtained in the simulation. We tested differences in donor and recipient subpopulations using harvesting alternatives and subpopulation ID as factors, and T mean as a response variable.
Factors affecting population viability. We assessed the effect of geographic and climatic variables on the mean time to extinction (T mean ) of the 100 subpopulations using a GLM. As geographic predictors, we included the linear distances to the nearest subpopulation (border to border) and to the metapopulation centroid (mean X and Y values of the centroids of all subpopulations). As climate predictors we included several variables obtained from the WorldClim database 95 with 1 km resolution from the time series 2000-2018, concordant with our data set: minimum January temperature (to account for the coldest temperature of the year), maximum August temperature (as the hottest peak) and annual precipitation (cumulative rainfall of the 12 months). We obtained the mean value of the raster cells that intersected with the habitat patches within each subpopulation. For each of the three climate variables we calculated the overall average value and the 18-year trend, estimated as the slope of the linear regression for the whole time period. As the response variable, we used the value of mean time to extinction of the 100 subpopulations under the base model obtained in VORTEX. We used a Gaussian link, and model fit was assessed via examination of residuals. Starting from the saturated model and using the drop1 function in R Chambers 19 we removed those variables that showed the highest p-value, subsequently running a new reduced model and comparing it to the previous one using a Log-Likelihood ratio Chi-squared test for the significance of the model. This procedure was carried out until all remaining variables showed a significant effect, and a significant difference between consecutive models tested using ANOVA was found. To avoid circularity, model parameters previously used to build the PVA, and therefore not independent of mean time to extinction, were not included in the GLM (such as initial population size, habitat surface and habitat quality).

Data availability
All data generated or analyzed during this study are included in this published article (and its Supplementary Information files).