Genetic signatures of population bottlenecks, relatedness, and inbreeding highlight recent and novel conservation concerns in the Egyptian vulture

Study species and study area

The Egyptian vulture is a small (∼2 kg) obligate scavenger living in a variety of open, arid, montane, and rugged landscapes in Southern Europe, Africa, and Asia (Cramp & Simmons, 1980; Del Hoyo, Elliott & Sargatal, 1994). Breeding adults nest on cliffs, defend territories against conspecifics, and behave as solitary or social foragers depending on the abundance, predictability, and concentration of food resources. Carcasses of wild animals (e.g., ungulates, lagomorphs, birds, and reptiles) and livestock represent the main food sources during the breeding season (Donázar, 1993; Blanco, Junza & Barrón, 2017; Blanco et al., 2019). Egyptian vultures are monogamous, highly philopatric, and territorial during the breeding season (Donázar, 1993; Sanz-Aguilar et al., 2017). This species shows a slow pace of life, with long generation time (14 years according to BirdLife International, 2019), average age of first breeding in the 7th year of age, low fecundity and reproductive rate (clutch size: 2 eggs; 0–2 fledglings per breeding attempt), and long lifespan (maximum age recorded in the wild: 24 years) (Sanz-Aguilar et al., 2017). This species is categorized globally as Endangered due to severe long-term declines through much of the distribution range (BirdLife International, 2019). The bulk of the European population is concentrated in Spain (about 1,500 pairs), where the breeding population has showed a strong decline in the last decades (BirdLife International, 2019).

The study was conducted in the province of Segovia (Castilla y Léon, Central Spain). This area holds a small and declining population of Egyptian vultures mostly nesting in cliffs located in two major systems of large gorges in two protected areas (Natural Parks of Hoces del Río Duratón and Hoces del Río Riaza) and their surroundings, while isolated pairs nest in small and scattered cliffs. The breeding population has suffered a sharp decline in the study area since the 1990s (from 48 pairs in 1993 to 27 pairs in 2019; Del Moral, 2009; Sanz-Aguilar et al., 2017). This population has been genetically characterized as forming part of the central-northern cluster from the distribution range in Iberia, including breeding sub-nuclei in Navarra and Aragón, north-eastern Spain (Agudo et al., 2011a). The Egyptian vulture pairs nesting in the area between Segovia province (the study area) and the northern distribution range in Iberia are also expected to belong to the same genetic cluster (Agudo et al., 2011a), as indicated by natal dispersal involving emigrants and immigrants recruited as breeders throughout this area (Serrano et al., 2021). Overall, for 19 Egyptian vultures recruited in or from the study area between 5 and 10 years after being marked as nestlings since 2003, we recorded a median dispersal distance of 41.2 km (range: 5.5–201). Most of the individuals recruited as breeders in the study area correspond to those also born in the study area (13 of 16). The remaining individuals were recruited in neighbouring provinces (Burgos and Soria respectively), while another individual was recruited in a more distant province (Zaragoza, Ebro Valley, north-eastern Spain). In addition, only two immigrant individuals born in the province of Zaragoza were recruited in the study area, while another individual born in Zaragoza was recruited in a neighbouring province (Guadalajara) to the study area. The study was carried out in accordance with the permission of the regional government of Castilla y Léon, Dirección General del Medio Natural, Servicio de Espacios Naturales.

Fieldwork and sampling procedures

An intensive monitoring program, including searching for territories, ringing nestlings, and identifying banded breeding individuals at territories was conducted at the Segovia province and its surroundings during the last four decades (Martinez & Blanco, 2002; Sanz-Aguilar et al., 2017). The nests were accessed by climbing when nestlings were feathered (at ages of about 50 days) but without risk of leaving the nests. The clutch size in this species is of two eggs, but in a proportion of nests only one nestling survives, mostly the first-hatched nestling (Donázar et al., 2020). Therefore, sampled nestlings can correspond to the single nestling in the nest (brood size = 1) or two siblings in the same nest (brood size = 2). Nestlings (n = 180) were measured, weighed, and banded with metal and plastic rings with alphanumeric codes allowing individual identification at long distance with terrestrial telescopes. A sample of blood as source of DNA was taken from the brachial wing vein and preserved in absolute ethanol.

To evaluate temporal trends in genetic relatedness, bottlenecks, and inbreeding in the population, we selected nestlings sampled in 2006 (n = 19 from 15 territories) and 2018 (n = 14 from 10 territories). The sampled nestlings correspond to about 40–50% of the active territories each year in the population, as a variable proportion of breeding pairs failed in the reproduction each year (Martinez & Blanco, 2002; Sanz-Aguilar et al., 2017). Each territory can include several nesting sites, but only a particular nesting site is used each year in each territory. Therefore, we sampled a relatively well-represented proportion of the territories available for nestling sampling each study year. Overall, our study includes the sampling of 42 Egyptian vulture individuals (see Appendix S1).

In the breeding season of 2017, we located a novel breeding pair composed of two individuals born in the same nest but in different years. Because breeders can use the same nest and territory year to year (Sanz-Aguilar et al., 2017), we assessed whether these individuals were genetically related or if they simply born in the same territory used by different breeders in different years. To this aim, we evaluated the relatedness among the nestlings (n = 10) born in 2004 (n = 1), 2005 (n = 1), 2006 (n = 1), 2007 (n = 2), 2010 (n = 2), 2018 (n = 2) and 2020 (n = 1) in the two nesting sites (distant <200 m) included in this single territory (a single nesting site used each year). This allowed us to evaluate the validity of the relatedness estimates in reflecting full-sib relationships. To determine whether the same breeders used the same territory across the years, we collected feathers of the breeders attending this particular territory during the nest access to band the nestlings in 2018. We wanted to determine whether these breeders were the parents of the nestlings raised in the same territory the previous years, and then to establish long-term nest re-use and changes in the composition of the pair mate.

DNA extraction, molecular sexing and microsatellite genotyping

The DNA was isolated by using the Quick-DNA Miniprep Plus Kit (Zymo Research) according to manufacturer’s protocols. Molecular sexing was performed using specific primers for partial amplification of the CHD1 gene (Fridolfsson & Ellegren, 1999).

Seventeen microsatellite markers previously characterized for Egyptian vulture were analysed: BV9, BV13, BV14, BV20, NP39, NP51, NP93, NP140, NP141, NP155, NP163, NP166, NP229, NP249, NP257, NP259, and NP296 (Kretzmann et al., 2003; Agudo et al., 2008). Microsatellite forward primers were labelled with 6-FAM, NED, PET, and VIC fluorescent dyes to determine the genotypes using capillary electrophoresis. PCR amplifications were performed in a total volume of 10 µl containing 5 µl of 2x MyTaq HS Mix (Bioline), 2 µM of each primer, and 5 ng DNA. PCR thermal conditions were as follows: initial denaturation at 95 °C for 5 min, followed by 40 cycles of 95 °C for 30 s, 58 °C for 1 min, 72 °C for 30 s, and a final extension at 60 °C for 10 min. Amplified fragments were electrophoresed on an ABI PRISM 3130xl Genetic Analyzer (Applied Biosystems) using the GeneScan 500 ROX size standard. Allele sizes were determined using Peak Scanner v.3.0 (Thermo Fisher Cloud) and are presented as supplementary data (Appendix S1).

Evaluation of microsatellite loci quality and genetic diversity

The screening of genotyping errors, large allele dropout, null alleles, and stuttering were performed with MICRO-CHECKER v.2.2.3 (Van Oosterhout et al., 2004) using Bonferroni (Dunn-Sidak) adjusted confidence intervals (95%) obtained after 10,000 Monte Carlo simulations. The significance of linkage disequilibrium (LD) and deviations from Hardy-Weinberg equilibrium were assessed using GENEPOP v.4.4.2 (Rousset, 2008) applying the Markov chain method (10,000 dememorisation steps, 1000 batches and 5000 iterations/batch), and a sequential Bonferroni correction. The number of alleles per locus (N_A), observed and expected heterozygosities (H_O and H_E) were calculated with GENALEX v.6.5 (Peakall & Smouse, 2012). Allelic richness (A_R) values were obtained using FSTAT v.2.9.3.2 (Goudet, 1995).

Analysis of genetic differentiation

The possible genetic differentiation between offspring from 2006 and 2018 was assessed using the F_ST estimator of Weir & Cockerham (1984), as implemented in GENALEX v.6.5. Allele frequency differentiation between samples of 2006 and 2018 was tested using Fisher’s exact test in GENEPOP v.4.4.2. We carried out a multi-locus spatial autocorrelation analysis (Smouse & Peakall, 1999) for both sampling years pooled and independently by using GENALEX. The spatial autocorrelation coefficient of genetic distance (r) was estimated for a maximum of 9 distance classes (a total distance of 90 km). Significance tests for each distance class were performed using 9,999 random permutations and 9,999 bootstrap replicates. In addition, we performed a Bayesian clustering analysis using the program STRUCTURE v.2.3.4 (Pritchard, Stephens & Donnelly, 2000). The analysis consisted of 10 independent runs for K values ranging from 1 to 3, each with 1,000,000 MCMC iterations after a burn-in of 100,000 interactions, under a model of admixture and correlated allele frequencies. The outputs were processed in the CLUMPAK server (Kopelman et al., 2015).

Genetic relatedness and inbreeding measures

Various studies have been reported that genome-wide single nucleotide polymorphisms (SNPs) may outperform microsatellite markers in kingship assignment and relatedness estimations (e.g., Hauser et al., 2011; Weinman, Solomon & Rubenstein, 2015; Thrasher et al., 2018). However, the microsatellites could be also a reliable choice if (i) there is a microsatellite panel available for the species under study, (ii) the markers have suitable polymorphic levels, (iii) they have a good resolution power and (iv) DNA samples have low-quality to generate reliable genome wide SNPs (Flanagan & Jones, 2019). We selected a well-known microsatellite panel previously developed and applied in population genetic studies of Egyptian vultures for our relatedness approaches, taking into account all these factors, our target species and the biological samples.

Relatedness (r) among individuals was calculated using different estimators to evaluate the robustness and reliability of our inferences, since the performance of the different relatedness algorithms are dependent of the dataset used (Anderson & Weir, 2007; Wang, 2014). The COANCESTRY application (Wang, 2011) was used to test the performance of five relatedness estimators reported by Wang (2002), Lynch & Li (Lynch, 1988; Li, Weeks & Chakravarti, 1993), Lynch & Ritland (1999), Ritland (1996) and Queller & Goodnight (1989). The pairwise relatedness values (dyads) were calculated and the best estimators were evaluated using the reference values obtained for the know full-sib relationships (full-sibs 1V5, 1V9 and G15749; full-sibs 290, 291, 32U, 32T, 9MA and 9LX; full-sibs 9MC and 9MF; see Appendix S1). The mean pairwise relatedness estimated from these individuals were compared with the theoretically expected value for full-siblings (0.5). Relatedness were also estimated for the group of all Egyptian vultures included in the study (a value less than 0.5 is expected). GENALEX v.6.5 (Smouse & Peakall, 2012) was applied to generate the tridiagonal matrix of pairwise estimates of relatedness using the Queller and Goodnight’s estimator (QG-r), which minimizes a downward bias for small sample sizes (Queller & Goodnight, 1989). The ML-RELATE software (Kalinowski, Wagner & Taper, 2006) was used to estimate the genetic relationships between all individual vultures. This program calculates the maximum likelihood estimates of relatedness (ML-r) and relationship categories between pairs of individuals from genotypic data. The likelihood values of the relationships available from the software (U, Unrelated; HS, Half sib; FS, Full sib; PO, Parent-Offspring) were calculated to determine the relationship with highest likelihood for each pair of individuals. When individuals born in the same year were categorized as PO, kinship was assumed to be FS. The inbreeding coefficients (F_IS) were estimated using GENEPOP v.4.4.2 (Rousset, 2008) to assess the level of inbreeding in the groups of individuals analysed. The program CERNICALÍN v.1 (Aparicio, Ortego & Cordero, 2006) was used to calculate the homozygosity by loci (HL), estimating the indices of inbreeding considering the proportion of loci within the genotypes of an individual that are homozygous. Statistical differences between years in the parameters of genetic diversity, relatedness, and inbreeding were assessed by Mann–Whitney U-test.

The spatial pattern of relatedness was tested by two-tailed Mantel test with 9,999 permutations considering the geographic distances (in km) between nests (Appendix S2) and the QG-r and ML-r values of each nestling from 2006 and 2018 (a single nestling per nest selected at random was considered in the case of broods with two siblings to avoid pseudoreplication).

Demographic events and effective population sizes

The detection of potential bottleneck events that lead to severe population size reductions in the recent (contemporary) past over the last decades was performed using the program BOTTLENECK version 1.2.02 (Piry, Luikart & Cornuet, 1999). The stepwise mutation model (SMM), infinite allele model (IAM), and two-phase mutation (TPM) model and mode-shift test were used. This analysis was complemented with the M-ratio model implemented in the software M_P_VAL (Garza & Williamson, 2001), which allows the detection of bottlenecks that occurred in the recent and distant (historical) past over centuries or millennia (Piry, Luikart & Cornuet, 1999; Williamson-Natesan, 2005). All parameters used in the analysis were previously described (see Blanco et al., 2021).

The contemporary effective population size (N_e) was estimated using the program LDNE v.1.31 (Waples & Do, 2008). The random mating model was used. Alleles with frequencies lower than 0.02 were excluded and the confidence intervals were obtained using jack-knife statistics.

Screening of infectious agents

The inbred offspring (nestling with feather malformations, ref: 9MC from 2018; and nestling apparently healthy, ref: 9MF from 2020) were analysed for the presence of three potential infectious agents that can cause plumage alterations (beak and feather disease virus, BFDV; avian polyomavirus, APV; and Salmonella). Specific molecular markers previously reported were used in the PCR screening of BFDV (Amery-Gale et al., 2017) and APV (Johne & Müller, 1998). PCR was performed in a reaction mixture of 10 µl containing 5 µl of 2x MyTaq HS Mix (Bioline), 2.5 µM of each primer and ∼20 ng of template DNA. The amplification protocol was composed of the following steps: 95 °C for 5 min followed by 40 cycles of 95 °C for 30 s, annealing at 58 °C for 1 min, 72 °C for 30 s, and a final extension at 60 °C for 10 min. To determine the presence of Salmonella, nestlings were sampled for cloacal mucosa using sterile microbiological swabs with Amies transport medium. The samples were analysed by conventional microbiological culture (ISO 6579:2002/Amd. 1:2007) and by using the real-time PCR IQ-CHECK™ Salmonella II kit (Bio-Rad) and Easy I protocol, according to the manufacturer’s instructions (see details of both methods in Blanco & Díaz de Tuesta, 2018; Blanco & Díaz de Tuesta, 2021).

Ethics statement

Our study followed ethical guidelines proposed for the Spanish Royal Decree 1205/2005 on the protection of animals used in experiments and scientific research. The study was carried out in accordance with permits from the Spanish Bird Ringing Centre (Permit Number: 530115), and the regional government of Castilla y Léon (Expte: EP/CyL/298/2016).

Genetic data quality, genetic diversity and differentiation

The probability of null alleles was negligible for all loci except for NP140 (dataset of all individuals and offspring 2018). Therefore, the subsequent analyses were carried out with and without this locus. The results were similar when NP140 was included and excluded from the data set. Therefore, we present the results including this locus. Only the loci NP140 (offspring 2018) showed a significant deviation from Hardy-Weinberg equilibrium (HWE) after Bonferroni correction. When all individuals were analysed in the same dataset no significant deviation from HWE was observed (Table S1). Exact tests for genotypic LD confirmed the absence of physical linkage at most loci. Only two significant linkages were found; however, these significances were lost when all siblings were removed from the dataset. In this context, all microsatellite loci were included in the multilocus analysis.

All microsatellite loci were polymorphic, with a total of 75 different alleles amplified. The number of alleles per locus ranged from two (BV14, BV20 and NP257) to eight (NP39), with an average of four alleles (Table S1). In the offspring of 2018, only one allele for the locus NP51was detected. Mean values of A_R, H_O, and H_E did not differ between nestlings from 2006 and 2018 (Table 1).

Genetic differentiation was not significant between offspring from 2006 and 2018 (F_ST = 0.02, p = 0.38). No significant differences were detected when Fisher’s exact tests of genetic differentiation were calculated from allele frequency data (Table S2). We found no evidence of genetic spatial autocorrelation when all individuals were analysed together and when the offspring of 2006 and 2018 were analysed as independent groups (Fig. S1). The results of Bayesian clustering analysis also show no evidence of genetic structure and differentiation between years (Fig. S2).

Relatedness

The evaluation of the five genetic relatedness estimators showed a good performance for most of the estimators (Table S3; Appendix S3). The estimators that allowed a better identification of the true full-sibs using the relatedness values were Wang (0.59 ± 0.13), Lynch & Li (0.60 ± 0.12) and Queller & Goodnight (0.60 ± 0.12) (Table S3). The mean relatedness values when all individuals were analysed were also very similar for all estimators, ranging from −0.02 ± 0.24 (Queller & Goodnight) to −0.01 ± 0.26 (Wang and Lynch & Li) (Table S3). In addition, most of the pairwise values obtained for the five estimators were highly correlated (average correlation value of 0.87 ± 0.07) (see Table S4 for pairwise correlation values). Therefore, for simplification purposes, we showed only the relatedness values for Queller & Goodnight estimates (QG-r) in the following analysis.

Measures of QG-r (Appendix S4) and ML-r (Appendix S5) were highly correlated (Mantel test: r = 0.835; p < 0.0001). The strategy used for categorizing the genetic relationships allowed the accurate assignment of all direct relationships in individuals with known pedigree, and only two deviations in the kinship assignation were observed for inferred relationships. Overall, we correctly assigned 93% (26 of the 28 pair-wise relationships) obtained for the group of individuals with known pedigree. Therefore, it is expected that relationship categories estimated for individuals with unknown pedigree reflect the actual relatedness with accuracy. For nestlings from 2006, 33 half-sib (19.8%) and three full-sib (1.8%) relationships were estimated, excluding the four full-sibs from broods of two nestlings reared in the same nest and year (Appendix S6). For nestlings from 2018, 13 half-sibs (14.9%) and four full-sibs (4.6%) relationships were estimated, excluding the four full-sibs from the same nest (Fig. 1A; Appendix S6). The geographic location of the territories sampled in 2018 is shown in Fig. 1B. The proportion of nestlings in each kinship category did not differ between years (Fisher exact test, p = 0.19), being about 20% those estimated being closely related individuals (half and full-sibs) when pooled years; for this analysis, we considered a single nestling selected at random from broods of two nestlings to avoid pseudoreplication. Mean values of relatedness and inbreeding were similar between nestlings from 2006 and 2018 (Table 1). The HL values for offspring of 2006 ranged from 0.19 to 0.52 (mean = 0.33 ± 0.09), while these values were slightly higher (mean = 0.39 ± 0.10) in offspring from 2018 (range = 0.26–0.59), indicating that on average 35–50% of an individuals’ loci are homozygous (Table 1; Table S5).

All territories showed at least one closely related individual (full or half sib) among the nestlings sampled in 2006 (100%, n = 15 territories), and most territories in 2018 (90%, n = 10). The number of territories with closely related nestlings ranged from one to seven in 2006, and from zero to four in 2018 (Fig. 2; see also Figs. 1A and 1B for the specific and geographical representations of these relationships in 2018). Relatedness between nestlings from different nests (Fig. 1A for 2018) showed no apparent spatial pattern (Fig.‘1B), as inter-nest distance was not significantly correlated with QG-r and ML-r values of each nestling from 2006 (r =  − 0.06, p = 0.62, n = 105; r =  − 0.08, p = 0.43, n = 105, respectively) and 2018 (r = 0.13, p = 0.35, n = 45; r = 0.08, p = 0.57, n = 45, respectively).

Demographic history and effective population sizes

The analysis with BOTTLENECK suggests population contractions in the recent past by analysing the samples from 2006 and 2018 populations (Table 2). Evidence of recent bottlenecks was found through SMM and the mode-shift test in 2006 and through all models in 2018 (Table 2). The M-ratio test does not support the occurrence of bottleneck events in a more distant (historical) past (values of M >Mc) (Table S6). Similar results were attained when nestlings from different territories and years were pooled to increase sample size (Table 2; Table S6).

Estimates of effective population size (N_e) for each sampled year were 100 (CI = 24-inf) for 2006 and 48 (CI = 13-inf) for 2018. The value of N_e considering all sampled individuals (2006 + 2018) was 53 (CI = 26–314).

Long-term genetic relationships and inbreeding effects in a familiar lineage

The retrospective analysis of relatedness in the sampled familiar linage showed a replacement of the mated pair present in 2004–2006 by another pair in 2007–2018 in this particular territory. This was confirmed by several facts: (i) nestlings from 2004–2006 were full sibs, and therefore shared the same progenitors; (ii) nestlings from 2007, 2010, and 2018 were also full sibs (Fig. 3) but unrelated to nestlings from 2004–2006 except for a single half-sib relationships, and therefore were offspring of a different breeding pair; (iii) nestlings from 2007–2018 showed parent–offspring relationships with their putative parents sampled for feathers in 2018 (Fig. 3); (iv) nestlings from 2004–2006 were unrelated with the female breeder but showed half-sib relationships with the male breeder in 2018 (Appendix S6).

Mean relatedness values for the familiar group from 2007 to 2018 was high (QG-r = 0.62 ± 0.11; ML-r = 0.53 ± 0.12), as is expected for dyads of closely related individuals. Two of these full siblings mated between them (QG-r = 0.59; ML-r = 0.43) in 2017: a female (ref. 291) born in 2007 and a male (ref. 32U) born in 2010 (Fig. 3), and nested in a territory distant 59 km from their natal nest (Fig. 1B). We found two half-sib relationships in this familiar lineage (Fig. 3): (i) between the male breeder (ref. M) and his granddaughter (ref. 9MC) from the inbred pair (refs. 291-32U); (ii) between a nestling from 2017 (ref. 290) and the son (ref. 9MF) of the inbred pair (refs. 291-32U). A single half-sib relationship was recorded by molecular methods between two actually full-sibs from different cohorts in this territory, which suggest slight deviations in the kinship assignation due to allelic frequency. Relatedness of all offspring with their actual parents and grandfathers were confirmed by pedigree and molecular methods. Note, however, that the ML-RELATE software may not differentiate true full-sibs and parent–offspring relationships in familiar cases with many closely related individuals, but these relationships were ascertained in individuals with known pedigree in this familiar lineage.

The inbred pair failed to reproduce in 2017. In 2018, the same pair reared a female offspring (ref. 9MC) that was confirmed to be daughter of their putative parents by microsatellite allelic patterns (Appendix S6) and relatedness values (QG-r = 0.72 ± 0.18; ML-r = 0.64 ± 0.19). This nestling showed alterations in feather structure and appearance across all the plumage (Fig. 4A). These acute malformations mainly consisted in greasy plumage, constricted calamus in all flight feathers with barbs appearing only in the distal part, and pinched appearance in a proportion of feathers (Fig. 4B). As a consequence of these severe and generalized malformations, this individual was unable to fly, and was retired by competent authorities and admitted in a wildlife recovery centre. In 2019, the same inbred pair failed in the reproduction apparently due to nestling predation. In 2020, the same pair produced an apparently healthy offspring (ref. 9MF) that was confirmed as son of their putative parents by microsatellite allelic patterns (Appendix S6) and relatedness values (QG-r = 0.72 ± 0.07; ML-r = 0.60 ± 0.07). This nestling was confirmed as full-sib of the nestling with feather malformations in 2018 (QG-r = 0.73; ML-r = 0.64), and as full-sibs with all other nestlings from this family lineage in 2007–2018 (Fig. 3). Values of HL were high for both inbreed individuals (Ref. 9MC = 0.48; Ref. 9MF = 0.54) as expected for offspring from parental full-sibs (Table S5).

The screenings of infectious agents performed in these inbreed offspring showed a BFDV positive result in the nestling with malformations (Ref. 9MC) and negative results for APV and Salmonella. All results were negative in the nestling apparently healthy (Ref. 9MF).

Genetic diversity, relatedness and inbreeding

A growing body of evidence has related reduction in genetic diversity with declining population size of wild species (Frankham, 1996; Willoughby et al., 2015). Our data show that the genetic diversity in a small open population of Egyptian vulture remained relatively stable despite of its decline in census and effective population size in the last decades. A relatively high proportion of nestlings from different and distant territories showed high relatedness in each study year. In addition, we found an increasing impact of recent demographic bottlenecks and the first evidence of mating between full siblings coinciding with lethal malformations in offspring.

We found similar levels of genetic variability and non-significant differences between two breeding seasons twelve years apart. This lack of inter-year differences is not surprising considering that the time period between both study years (12 years) is shorter that the generation time in the Egyptian vulture (14 years; BirdLife International, 2019). Some studies have suggested that F_IS values tend to be lower (more outbred) when estimates are calculated using only samples from juveniles, and higher (more inbred) in samples from adults (Basset, Ballouxf & Perrin, 2001; Parreira & Chikhi, 2015). This is a potential consequence of a heterozygote excess increased by reproductive events, which are later counteracted by dispersal (Basset, Ballouxf & Perrin, 2001). Thus, we cannot exclude a potential underestimation bias in our estimates of inbreeding coefficients, since they were obtained from nestlings. Overall, the values of genetic diversity estimates were similar to those reported for several declining vulture species (Gautschi et al., 2003; Arshad et al., 2009; Çakmak et al., 2019), higher than those observed in other declining species (Kleinhans & Willows-Munro, 2019), and even higher than in vulture species showing larger population size and stable or increasing demographic trends (Le Gouar et al., 2008; Arshad et al., 2009). Therefore, these results should be interpreted with caution because different markers were used to estimate genetic diversity indices within and among species (see Väli et al., 2008). In addition, low genetic diversity is not always the result of a reduction in effective population size, as the former can be high in species with very small, isolated, and declining populations (Kekkonen, Wikström & Brommer, 2012; Woolaver et al., 2013). A low genetic diversity may be due to demographic bottlenecks and population rebounds, as well as due to founder events with a reduced number of individuals, and absence of random mating (Thévenon et al., 2003; Kekkonen, Wikström & Brommer, 2012). Current levels of neutral genetic diversity, relatedness, and inbreeding are similar to those reported for this species in the study area and other regions in continental Spain using the same suite of microsatellite loci (Agudo et al., 2011a). However, the values reported previously for the study area are not strictly comparable with those presented here because the former partially include full-grown individuals of unknown geographical origin (natal nest), age, and cohort. These individuals can include immigrants from other genetic clusters and geographical areas due to natal dispersal and non-breeding nomadic movements (IBERIS, 2007; Sanz-Aguilar et al., 2017).

Our results show a relatively high proportion (about 20%) of dyads comprising half- and full-sibs relationships among nestlings born in different territories in each study year. Indeed, a proportion of nestlings from distant territories in each year were genetic full-siblings. However, we found no evidence of a relationship between pairwise geographic inter-territory distances and pairwise genetic relatedness between nestlings from different territories in each cohort. This suggests a low contribution of extra-pair paternity to these intra- and inter-nest kinship relationships. Although extra-pair copulation has been reported at low frequency in the Egyptian vulture (Donázar, Ceballos & Tella, 1994), the contribution to extra-pair paternity remains generally unexplored in vultures (Le Gouar et al., 2011). Our results showed exclusive full-sib relationships between nestlings from the same nest and cohort (broods of two siblings). However, we found two half-sib relationships in the familiar lineage with known genealogy. This can be indicative of slight and conservative deviations in the kinship assignation due to allelic frequency, which deserves further assessment including more parents and offspring with known pedigree. The assessment of genetic relatedness using genome-wide SNP loci could improve the estimates (Flanagan & Jones, 2019), although our microsatellite panel with a probability of an accurate kinship assignment higher than 90% seems reliable enough for the aims of this study. Proportions of dyads falling in the kinship categories were higher than those recorded in other vulture species, although these relatedness estimates were calculated for wild and captive full-grown individuals (Gautschi et al., 2003; Ishtiaq et al., 2015). Given that these studies analysed relatedness of individuals of different ages, they can include actual parent–offspring relationships and full and half siblings from different cohorts. Instead, we found closely related individuals among nestlings from the same cohorts but distant territories, which clearly indicate that their parents were also closely related within and between pair mates from different territories. At the scale of the whole population, we found that most territories showed closely related individuals among the nestlings sampled in each study year. Indeed, nestlings of each particular territory were closely related to those from a variable number of other territories in each year. For instance, nestlings from three territories were closely related to those from other seven territories each in 2006, when the sampling was more extensive, while nestlings from one territory were closely related to those from up four territories in 2018. This implies that each one of these nestlings was closely related to those from about half of the territories sampled each year. Given that we sampled a representative proportion of the successful pairs each year, this further indicates a high relatedness in the population. Overall, this highlights the importance of knowing natal origin, age, and cohort in the assessment of genetic diversity and relatedness of wildlife populations.

Demographic history

The analysis of bottlenecks detected a clear signature of genetic contractions during the last decades, coupled with the recent decline according to census and effective population size in the study area and other regions included in the same genetic cluster (Carrete et al., 2007; Agudo et al., 2011a; Sanz-Aguilar et al., 2017). The significance observed in the heterozygosity excess contrasts with the consistent and non-significant results of the M-ratio test in both study years, which indicate that historical bottlenecks for the last centuries or millennia are not evident (Williamson-Natesan, 2005; Peery et al., 2012). However, it is important to highlight that statistical power of bottlenecks tests is highly related to the number of samples and markers used, which in many cases affect the detection of historical population reductions (Peery et al., 2012; Hoban, Gaggiotti & Bertorelle, 2013). Furthermore, the different sensitivity of heterozygosity-excess and M-ratio tests to possible violations of mutation model assumptions may affect the inferences about the timing of population declines (Peery et al., 2012). Therefore, we cannot completely rule out the occurrence of historical bottlenecks (distant past) in this population. The clear signature of recent bottlenecks disagrees with the relatively high values of genetic diversity estimates. This apparent contradiction may arise because of long generation times owing to delayed breeding in this long-lived species (Sanz-Aguilar et al., 2017) can still reflect the retained genetic diversity in the parents of the nestlings analysed. This time lag associated to a slow pace of life could act as an intrinsic buffer to stave off genetic erosion in long-lived species (Hailer et al., 2006). This implies that long-term studies are necessary to adequately detect temporal changes in genetic diversity in these species. Given that Egyptian vultures can reproduce during many consecutive years across their long lifespans (Sanz-Aguilar et al., 2017), some of these breeders can be old enough to still harbour relatively high genetic diversity typical of not bottlenecked populations. In support of this possibility, no evidence of recent bottlenecks was found in the study area sampled prior to our study (TPM model, p = 0.153, from Agudo et al., 2011a). This suggests that the impact of recent bottlenecks has increased in the last two decades to the point that it can be observed in its current progression in parallel to the reduction of the effective and census population size, but still without a clear reflection in genetic diversity estimates (see also Çakmak et al., 2019). However, the interpretation of these estimates should be made with caution when the proportion of closely related individuals in the population is high (Gautschi et al., 2003). In addition, the values of relatedness previously reported for the study area were obtained by pooling data from local nestlings with full-grown individuals of unknown natal origin, age, and cohort captured in foraging areas (IBERIS, 2007; Agudo et al., 2011a). Therefore, the lack of recent bottlenecks previously to our sampling could arise through an increase in genetic diversity by including individuals from other geographic regions and genetic clusters in the relatively continuous distribution range across central-north Iberia (Agudo et al., 2011a). This rationale could also apply to our data because, although the nestlings were sampled in a relatively small area of the distribution range of this genetic cluster, a small proportion of breeders could have been born outside the study area. In this case, the evidence of recent bottlenecks is supported by our conservative data, potentially including a proportion of offspring from immigrant individuals, which should be confirmed with a more extensive assessment over a larger geographical area.

Implications of inbreeding and emerging pathogens in conservation

The increase of close kin among mated breeders in small populations can lead to inbreeding depression. The consequences of a reduction of genetic diversity on offspring health and fitness can be severe for the affected populations (Charlesworth & Willis, 2009). This is a critical factor in the conservation of endangered species when it occurs rapidly enough to preclude the purging of deleterious alleles (Brook et al., 2002; Fox, Scheibly & Reed, 2008). Consequently, host susceptibility to disease agents is expected to increase due to inbreeding, although their effects on health can surface depending on whether the population is exposed to particularly dangerous and contagious pathogens (Smith, Sax & Lafferty, 2006; Gupta, Robin & Dharmarajan, 2020). The evidence of recent bottlenecks and the reduction in census and effective population size coincided in time with the first evidence in this species of close inbreeding, confirmed by pedigree of individually-marked nestlings monitored subsequently, and by genetic analysis. The inbred pair showed low reproductive success (only one apparently healthy nestling raised from four reproductive attempts) and the first evidence in this species of lethal feather malformations in offspring, which suggests a possible cause-and-effect relationship. The affected nestling showed alterations in the overall plumage that precluded fledging and these malformations remained the following two years in captivity. This kind of feather dystrophy resembles those characterizing the pinching off syndrome reported associated to inbreeding in several raptor species in absence of any identified aetiological agent (Müller et al., 2007). Similar malformations have been recorded in parrot and allies (Psittaciformes) affected by the beak and feather disease virus (BFDV), and recently this avian circovirus has been recorded affecting non-psittacine species, especially in small and isolated populations and species with low genetic diversity (Raidal & Peters, 2018). We found that the nestling with generalized feather malformations was positive for this virus, which genetic sequence is currently being characterized to determine its phylogenetic position informing on its potential origin (Harkins et al., 2014; Morinha et al., 2020). In contrast, the apparently healthy sibling (born in 2020) of the affected offspring (born in 2018) from the same inbred pair was negative to BFDV. This suggests a potential interaction between inbreeding and malformation depending on infection by BFDV, although other non-surveyed pathogens or pollutants could be involved in disease causation, which requires further research. In support of an interaction between inbreeding and malformation due to BFDV, we found other BFDV-positive nestlings in the population (n = 12), but no other with malformations has been recorded in previous studies (Blanco, Junza & Barrón, 2017; Blanco et al., 2020; Sanz-Aguilar et al., 2017).