Genetic diversity of the Griffon vulture population in Serbia and its importance for conservation efforts in the Balkans

The Griffon vulture was once a widespread species across the region of Southeast Europe, but it is now endangered and in some parts is completely extinct. In the Balkan Peninsula the largest Griffon vulture inland population inhabits the territory of Serbia. We present, for the first time, the genetic data of this valuable population that could be a source for future reintroduction programs planned in South-eastern Europe. To characterize the genetic structure of this population we used microsatellite markers from ten loci. Blood samples were collected from 57 chicks directly in the nests during the ongoing monitoring program. We performed a comparative analysis of the obtained data with the existing data from three native populations from French Pyrenees, Croatia, and Israel. We have assessed the genetic differentiation between different native populations and determined the existence of two genetic clusters that differentiate the populations from the Balkan and Iberian Peninsulas. Furthermore, we analysed whether the recent bottleneck events influenced the genetic structure of the populations studied, and we found that all native populations experienced a recent bottleneck event, and that the population of Israel was the least affected. Nevertheless, the parameters of genetic diversity suggest that all analysed populations have retained a similar level of genetic diversity and that the Griffon vulture population from Serbia exhibits the highest value for private alleles. The results of this study suggest that the Griffon vulture populations of the Balkan Peninsula are genetically differentiated from the populations of the Iberian Peninsula, which is an important information for future reintroduction strategies.


Results
Parameters of genetic diversity. Parameters of genetic diversity per population are presented in Table 1, whereas the parameters of genetic diversity per locus per population are presented in Tables S1, S2. The obtained results demonstrated lack of polymorphism for locus BV13 in the native Griffon vulture populations inhabiting Serbia and Croatia, as well as in the introduced individuals in Navacelles population (Table S1). All other loci were polymorphic in all analysed populations, and the most polymorphic locus overall was BV12, with the highest number of alleles and greatest allelic range (Table S1). The average number of alleles per locus varied from 2.0 to 13.4 (Table S2) while the average number of alleles per population varied from 4.8 to 6.7 (Table 1  and Table S2). Among native populations the highest number of alleles as well as the highest mean number of alleles based on minimal sample size was detected in the Griffon vulture population inhabiting Israel (Table 1  and Table S1). Considering the mean number of private alleles based on minimal sample size the highest number among native populations was observed in Serbia, although the actual number of detected private alleles Scientific Reports | (2020) 10:20394 | https://doi.org/10.1038/s41598-020-77342-1 www.nature.com/scientificreports/ was highest in Israel and Causses ( Table 1). The average effective number of alleles was similar in all native populations and it varied from the lowest detected in the Pyrenees (3.36) to the highest detected in Israel (3.93) ( Table 1). The effective number of alleles per locus per population is presented in Table S1. Values for the observed (H O ) and expected (H E ) heterozygosity were found to be very close in all populations, ranging from 0.53 to 0.60 for H O and from 0.55 to 0.60 for H E . No deviation from Hardy-Weinberg equilibrium was detected in any of the analysed populations (Table S1). The Wilcoxon signed-rank test before the Bonferroni correction demonstrated that all native populations, except the one in Israel, exhibit significant heterozygosity excesses. After the correction the population from Serbia was also without significant heterozygosity excesses ( Table 1). The values for G-W index were high for all analysed populations and, among the native populations, the one from Israel exhibited the highest value. The mean number of pairwise differences among native populations was the highest in the Pyrenees ( Table 1). The average number of pairwise differences between and within populations is visualized in Fig. S1 together with Nei's distances. The lowest value for the RMP was detected in the population of the Pyrenees ( Table 1). The population with the lowest number of loci with linkage disequilibrium was the population from Israel and out of all native populations the one in the Pyrenees had the highest number of loci with linkage disequilibrium (Table S3). Estimated effective population sizes for all analysed populations are presented in Table 1 and Table S4.
Differentiation of analysed Griffon vulture populations. The AMOVA performed across all loci revealed that 2.24% of the genetic variance can be contributed to the variation among the populations while the rest of genetic variance can be explained with the variation within populations ( Table 2). The AMOVA performed on different groups of populations showed the existence of generally low values of genetic variance among the groups of populations and the highest observed value, 3.85%, was observed between the group of populations from the Balkan Peninsula (Serbia and Croatia) and the group of populations derived from the Iberian Peninsula (the Pyrenees with reintroduced populations from Baronnies, Causses, Verdon, Navacelles and Diois) ( Table S5). The lowest observed value of genetic variance between the groups of analysed populations was observed between the group of populations derived from the Iberian Peninsula (Table S5). The pairwise population F ST values between the Griffon vulture population of Serbia and other analysed populations are generally low but statistically significant for all pairs of native populations (Table 3). Pairwise population F ST is visualized in Fig. S2 and visualization by MDS plot shows the positioning of populations in two dimensions (Fig. 1). Griffon vulture populations of Serbia and Croatia occupy outlying position compared to other populations and together they form a distinct cluster. The Griffon vulture population from Israel is www.nature.com/scientificreports/ positioned close to the population of the Pyrenees, but between Pyrenees and populations from the Balkan Peninsula. In addition, all introduced populations cluster together with the native population of the Pyrenees from which they originate. Analysis performed by the DAPC method ( Fig. 2) showed that populations cluster in a pattern similar to the one presented in Fig. 1. The birds from Serbia and Croatia form two different clusters of their own, while the cluster that consists of the birds from Israel occupies an intermediary position between all the clusters. The birds from the Pyrenees and the birds from introduced populations are grouped in cluster overlaps, indicating the similarity of these populations. The populations of Baronies, Verdon and Diois contain the birds of Spanish origin, and these birds were used as the sample population representative for Spain in Le Gouar 25 . Additionally, in the same paper, the authors decided to pool the data for French Pyrenees (Ossau) and Spain, since there were no detectable genetic differences between these two populations. These results indicate that the population of the Pyrenees, together with the reintroduced populations, can be used to represent the population of the Iberian Peninsula.
STRU CTU RE Harvester showed that K = 2 is the most likely scenario ( Fig. 3a and Fig. S3). The distribution of two genetic clusters shows that in the Balkans (Serbia and Croatia) the second cluster is dominant, which is opposite from the Pyrenees where the first cluster is dominant, while the population that inhabits the Middle East (Israel) has almost equal contributions from both genetic clusters (Fig. 3b,c). STRU CTU RE analysis was performed on native and introduced populations as well, and the results are shown in Fig. S4.   Table 3.

Discussion
In the present paper we present for the first time genetic data of the Griffon vulture population from Serbia that inhabits parts of the Balkan Peninsula with a continental climate. We have also performed a comparative analysis with the data obtained on the Griffon vulture populations from Croatia, Israel (Mediterranean climate), and the Pyrenees in France, published in Le Gouar 25 . The endangered Griffon vulture species of the Balkan Peninsula is the last inland population adapted to the continental climate. It could thus be an important source for the reintroduction programs in Southeast Europe. In the present paper we have assessed the genetic variability using the highly variable microsatellite loci described and developed for conservation and reintroduction programs 35,36 . Presented data demonstrate that Griffon vulture population from Serbia exhibits similar values of genetic diversity as other native populations. Similar values of the parameters of genetic diversity were observed for the Griffon vulture populations of Cyprus and Spain 26 although a slightly different set of microsatellite loci were used in that study. Compared to other endangered vulture species, the Griffon vulture population of Serbia exhibits higher values of genetic diversity than the ones observed in the populations of the Indian vulture (G. indicus) from Pakistan 26 , Cape vulture (G. coprotheres) from South Africa 33 , Cinereous vulture (A. monachus) from Spain, Georgia, Armenia and Kazakhstan 32 , and Whitetailed eagles (Haliaeetus albicilla) from Norway, Estonia, and Germany 37 . Interestingly the population of Griffon vulture from Serbia exhibits similar values for parameters of genetic diversity when compared with the populations of the Bearded vulture (Gypateus barbatus) from the Pyrenees 38 , White-rumped vulture (G. bengalensis) from Pakistan 26 , and the Golden eagle (Aquila chrysaetos canadensis) from North America 39 , while it exhibits lower values in comparison to the White-backed vulture (G. africanus) from Namibia 26 . Although the population of Griffon vulture from Serbia experienced a serious bottleneck during the last decade of the twentieth century, it retained its rich genetic diversity which can be observed even when compared to the biggest and historically most stable population of Europe, the population from the Iberian Peninsula i.e. Pyrenees 40 . This finding corresponds with the observation that the similar long-lived species like White-tailed eagles retained relatively high levels of genetic diversity even after they had experienced a bottleneck event 37 . Some parameters, like the number of private alleles and the mean effective number of alleles, even suggest that the Balkan Peninsula, i.e. populations from Croatia and Serbia, exhibit higher genetic diversity than the ones observed in the Pyrenees. The highest number of private alleles based on a sample of 12 individuals was detected in the Griffon vulture population of Serbia, suggesting that this population could act as an important source of genetic variation 34 .
Considering the bottleneck in the Griffon vulture population of Serbia, the observed G-W index is the lowest observed among all the populations analysed, but is still relatively high and close to the observed values for other analysed native and introduced populations. According to the G-W index, the populations with the values closest to value 1 are considered to be the most stable populations, which is, according to our analysis, the population of Israel. Wilcoxon test, which evaluates the probability of recent bottleneck events based on heterozygosity excess, even suggests that the highest probability for this event is detected for the population of the Pyrenees, and not the population of Serbia. In addition, even the linkage disequilibrium between the pairs of loci suggests a higher probability of the recent bottleneck event for the population of the Pyrenees than the population of Serbia. It has been shown that the reduction in population size leads to a non-random association of alleles for unlinked loci [41][42][43] and, further, that the detection of linkage disequilibrium for unlinked loci suggests a recent reduction in effective population size 44 . The effective population sizes estimated for the analysed populations showed that Israel has the biggest population, followed by Serbia and the Pyrenees.
In addition, the discrepancy between the expected and observed bottleneck effects for the population of Serbia can be explained by possible immigration of the birds from nearby related colonies/populations. It is possible that after war conflicts in the region in the early '90 s the four known Griffon vulture colonies in Eastern Herzegovina (Bosnia and Herzegovina) didn't truly disappear, instead, they might have resettled in Serbia 45,46 forming a new colony in 1995 in the Mileševka Gorge 45,46 . These birds might have helped to maintain the high genetic diversity in the population of Serbia, thus reversing the negative effects of the previous rapid population decline.
Population structure analysis revealed the existence of two genetic clusters, and the distribution of these clusters can clearly distinguish two groups that exhibit different proportions, with the addition of one admixed population. Populations from Serbia and Croatia form one group specific to the Balkan Peninsula, while another group would be specific to the Pyrenees, and all the populations derived from the Iberian Peninsula ( Fig. 3 and Fig. S4). The third group, i.e. the admixed population, with almost equal proportions of both genetic clusters, is the Israeli population from the Middle East. In addition, all introduced populations are remarkably similar, almost identical, to the population from the Pyrenees. This is expected, since all introduced populations are formed from the birds that are caught and brought from the Spanish and French sides of the Pyrenees 25 . This confirms that, although all introduced populations exhibit high genetic variations, they form one continuous and uniform population with their source population (Fig. S4), which is a pattern also observed in Le Gouar 25 .
Additionally, the F ST values, as well as MDS, clearly demonstrate that the populations from the Balkan Peninsula represent one genetically differentiated group and should be viewed as a distinct population (Fig. 1). Even though Le Gouar 25 presented the existence of genetic differences between the Balkan population from Croatia and those originating from the Pyrenees and Spain, the Griffon vulture population from Spain was used as a source population for the reintroduction program in Bulgaria 47 . Le Gouar 25 explained the observed genetic differentiation as the result of the recent isolation, and not as a measure of significant long term differentiation among two distinct populations. It should be noted that the Griffon vultures which inhabit the Balkan Peninsula exhibit some specific morphological differences compared to their counterparts from the Iberian Peninsula. These could be adaptations to different climate conditions 48,49 . On average, the Griffon vultures from the Balkan Peninsula exhibit a greater body mass and have later hatching time, 1 month after the Mediterranean populations, most likely as an adaptation to a colder climate and delays of herds pasturing due to longer periods of cold 50 .
Another interesting result for the Griffon vulture population from Israel suggests that the Middle East could be recognized as a hybridization/admixture zone 26,51 , or, more likely, as the region from which European populations have originated 6 . It has been reported that young birds from the Balkan Peninsula migrate to the Middle East, where they stay between 3 and 4 years until they are sexually mature 52 , then return to the region where they have hatched, to form nests 53,54 . According to our results, it is possible that some of the birds choose to stay in the Middle East where they interbreed with the birds from the region. Another possibility is that the observed equal contribution of both genetic clusters in the Griffon vulture population from Israel suggests that it was the source population from which this species colonized the Mediterranean. It is suggested that during the Last Glacial Maximum (LGM) in Europe the European Griffon vulture populations retreated to refugia in North Africa and the Arabian Peninsula 55 . After the end of LGM the recolonization of Europe occurred in two directions: across Gibraltar into the Iberian Peninsula and across Bosporus into the Balkan Peninsula. This hypothesis is supported by the results which show that the population of Israel is the only population without a recent bottleneck. It is plausible that during the initial colonization of Europe from the Middle East the populations of the Iberian and the Balkan Peninsulas went through a founder effect and successive bottleneck, which led to the prevalence of two different genetic clusters, thus differentiating them.
Even though the two distinct groups can be recognized, AMOVA results suggest that the highest genetic variation is observed within populations and less between them. Since the within-population variance based on microsatellites is generally high, it can lead to an underestimation of variance among the populations. However, STRU CTU RE analysis supports the finding that these populations highly overlap in their genetic composition, since they share two genetic clusters.
Griffon vulture species is considered to be socially monogamous and once formed couples stay together 56,57 , although rare cases of extra-pair copulation within the same colony have been observed 58 but never confirmed by genetic analyses 59 . Furthermore, it has a wide dispersion and migrates throughout the Mediterranean. In their migrations, the birds from Serbia can be found almost exclusively in the eastern Mediterranean towards the Arabian Peninsula, but when the nesting period approaches, these birds return to the region in which they have hatched 14 . This natal philopatric behaviour is important because it can act as a mechanism for increasing genetic differences between the populations that nest in distinct regions, especially in combination with somewhat different climatic conditions. In conservation biology, this effect of philopatry is well recognized as it can be one of the causes of phylogeographic structuring due to geographic isolation caused by this kind of behaviour. This long-term isolation can be used as a major criterion to identify population, i.e. evolutionary significant unit, that deserves separate management or priority for conservation due to its high distinctiveness 60 . The effect of natal philopatry on genetic differentiation of Cinereous vulture populations from the Iberian and Balkan peninsulas was used to explain the observed genetic differences, resulting in the fact that they are now recognized as evolutionary significant units worthy of protetction 32 . New data presented in our paper show that the Griffon vulture populations from the Balkan and Iberian Peninsulas are genetically differentiated as well, and that they may have been adapted to different climatic conditions. This raises the question whether the introduction of foreign birds to the Balkan Peninsula is justified. The morphological differences, as well as the observed genetic differences 25 , should be a rationale to pay closer attention to the local populations as the source for reintroduction programs. The fact that a small percentage of introduced birds from the Iberian Peninsula to Bulgaria survived suggests that the continental Balkan population deserves different conservation strategy and that the reintroduction should not be performed with foreign birds but the ones from the already existing and viable native populations of the Balkan Peninsula. The observed differences between the populations of the Balkan and Iberian Peninsulas could be the result of a non-adaptive processes, like a bottleneck and/or a founder effect, as well as the result of the natal philopatry. In the present paper, we used selectively neutral markers which limit the conclusions about adaptive differences, enabled from the loci under selection or the morphometric data. Nevertheless, the experience of reintroduction with the birds from the Iberian Peninsula into the Bulgaria region may suggest that these differences could be the result of local adaptations as well.

Materials and methods
Population location and sampling. In the present study we have analysed new samples from the Griffon vulture population that inhabits the region of South-western Serbia, the gorge of the river Uvac. Uvac belongs to the protected area of the "Special nature reserve Uvac" (43°24′N 19°56′E 25 . In the comparative analyses, we used the raw data (amplicon lengths) provided by Prof. Le Gouar (personal communication). In order to be able to compare genetic data generated in two different laboratories, six samples from the Griffon vulture population of Croatia used in the study of Le Gouar 25 were reanalysed. The feathers from the same birds used in Le Gouar 25 were chosen for the DNA extraction in order to be processed in the same way as the new samples from Serbia. For all analyses performed in the present work, we used data for individual birds which had at least 5 genotyped loci, thus creating the differences in the sample sizes for the same populations used in our analysis and Le Gouar 25 . In a previous study, sample sizes for Baronies, Verdon, and Diois were smaller than the sizes used in our analysis, because the authors excluded the birds of known Spanish origin from the group of these samples, which they used to create a representative group for Spain. Considering that no genetic structuring between the introduced groups and between introduced groups and the native population in Ossau was observed in Le Gouar 25 we included these birds in the samples for Baronies, Verdon, and Diois. Furthermore, due to the fact that in the same work the low level of genetic differentiation between the samples from Ossau and Spain was observed, we dubbed the Ossau as the population of the French Pyrenees derived from the Iberian Peninsula.
DNA extraction and amplification of microsatellite loci. The DNA extraction was performed using the GeneJET Genomic DNA Purification Kit (Thermo Fischer Scientific Cat.No. K0721). The DNA concentration and quality were checked both by a spectrophotometer (NanoPhotometer, IMPLEN, Germany) and agarose gel electrophoresis. We used ten microsatellite loci. Five were taken from Mira 36 and the other five were chosen from Gautschi 35 (Table 4). Microsatellite loci were amplified in PCR reactions in which forward primers were labelled with a fluorescent dye (Table 4). PCR was performed in five reactions that differed in annealing temperature (Table 4) using the following program: one cycle of initial denaturation at 94 °C for 5 min, after which there were 30 cycles of 35 s at 94 °C, 35 s at the annealing temperature (Table 4) and 35 s at 72 °C. The step of final elongation was performed at 72 °C for one hour. Loci were amplified in a 2720 Thermal Cycler (Applied Biosys- Table 4. List of loci used in genetic analyses, primers used for their amplification, fluorescent dyes used for tagging the forward primers, annealing temperature (Tm) used in each reaction for the amplification of specific microsatellite loci with the combination of loci amplified together in multiplex reactions I, II, III, and individual reactions.
Fragment analysis. For fragment analysis, the first and second multiplex reactions were multipooled by the GF9C1 reaction, while the third multiplex reaction was multipooled by the GF3F3 reaction. All reactions were mixed in equal volumes and plated as one reaction in a volume of 0.5 µL. Each amplification mix contained five different loci. GeneScan 600 LIZ Size Standard was used to score alleles (Applied Biosystems, UK). Fragment analysis was performed on the 3130 Genetic Analyzer (Applied Biosystems, UK). Data were analysed using Gene Mapper Software (Life Technologies, USA).
Population genetics. In order to assess the genetic diversity of analysed populations, we determined the following parameters for each locus and population: number of alleles, number of alleles based on a minimal sample size (obtained by rarefaction), number of private alleles based on a minimal sample size (obtained by rarefaction), observed (H O ) and expected (H E ) heterozygosity, random match probability and mean number of pairwise differences. Random match probability represents the parameter that expresses the probability that two randomly picked individuals from a population have a matching genotype, and is calculated as the sum of square frequencies 62 . The mean number of pairwise differences represents the measure of differences between all pairs of haplotypes in the sample. Calculations were performed using Arlequin ver. 3.5.2.2 63 and HP-Rare 1.1 64 . The effective number of alleles per locus was calculated using the following formula: A e = 1/ p 2 i , where pi is the frequency of an ith allele. Arlequin software was used for assessing genetic differentiation among populations through Analysis of molecular variances (AMOVA) and estimating pairwise population and overall F ST values. Statistical significance of all performed tests was assessed with 10,000 permutations. The matrix of pairwise population F ST values was visualized both by R functions and two-dimensional scaling (non-metric MDS) using PAST 3.25 65 . Hardy-Weinberg equilibrium was tested using Arlequin software with 1,000,000 number of steps in MC and 100,000 dememorization steps.
In order to detect whether a significant number of targeted loci featured heterozygosity excess, we used the Wilcoxon signed-rank test implemented in the software Bottleneck v.1.2.02 66 this analysis was used, since excessed heterozygosity would imply a recent bottleneck event. The two-phase mutation model (TPM) 67 was used and parameterized to 25% of stepwise mutation. The effect of the recent bottleneck was also evaluated based on the Garza-Williamson index (G-W) which represents the ratio between the number of alleles and allelic range and was calculated using Arlequin software. Linkage disequilibrium between the pairs of loci was estimated using the likelihood ratio test in Arlequin software with 10,000 numbers of steps in MC and 10,000 dememorization steps. Effective population sizes for the analysed populations were calculated with NeEstimator 2.1 68 using the linkage disequilibrium method and monogamy mating model. In order to correct the probabilities when multiple tests were done simultaneously, we performed a sequential Bonferroni test for all performed analyses 69 .
Population structure. Clustering of the analysed populations and the observed distances among samples were presented using discriminant analysis of principal components (DAPC) 70 using the R package adegenet 71 . This method consists of performing linear discriminant analysis (LDA) on the principal components analyses (PCA) transformed matrix. The number of retained PCs was estimated using randomly repeated cross-validation (100 iterations) which consisted of performing DAPC on 90% of randomly sampled training set observations (stratified sampling was used so that the training set consisted of 90% of the observations from each population) after retaining 10-83 PCs and using the obtained models to predict the groups (populations) in the remaining 10% of samples (test set). Average prediction success per group was used as the metric to choose the optimal model. In order to estimate the number of genetic clusters represented in the native populations of the Griffon vulture we used the STRU CTU RE v 2.3.4 software [72][73][74][75] . We also performed the same analysis including the reintroduced populations. The analysis was performed using the admixture model with a burn length of 10,000 and a Markov Chain Monte Carlo (MCMC) of 100,000 randomizations. The range of the possible number of clusters (K) was from 1 to 10, with a series of 10 runs for each K. STRU CTU RE harvester (Web v0.6.94, https ://taylo r0.biolo gy.ucla.edu/struc tureH arves ter/ 76 ) was used to analyse the results obtained by STRU CTU RE. Based on the results generated by STRU CTU RE this software creates a plot of the mean likelihood value per K value and calculates the highest value of the second-order rate of change (Delta K) according to the method of Evanno 77 in order to detect the number of K groups that best fits the data set. The model choice criterion, LnP(D), implemented in the STRU CTU RE which detects the true K as an estimate of the posterior probability of the data for a given K was evaluated as well. The most likely scenario was chosen and used to graphically plot both the analysed individuals and populations.