Low genetic diversity and shallow population structure in the endangered vulture, Gyps coprotheres

Globally, vulture species are experiencing major population declines. The southern African Cape vulture (Gyps coprotheres) has undergone severe population collapse which has led to a listing of Endangered by the IUCN. Here, a comprehensive genetic survey of G. coprotheres is conducted using microsatellite markers. Analyses revealed an overall reduction in heterozygosity compared to other vulture species that occur in South Africa (Gypaetus barbatus, Necrosyrtes monachus, and Gyps africanus). Bayesian clustering analysis and principal coordinate analysis identified shallow, subtle population structuring across South Africa. This provides some support for regional natal philopatry in this species. Despite recent reductions in population size, a genetic bottleneck was not detected by the genetic data. The G. coprotheres, however, did show a significant deficiency of overall heterozygosity. This, coupled with the elevated levels of inbreeding and reduced effective population size, suggests that G. coprotheres is genetically depauperate. Given that genetic variation is considered a prerequisite for adaptation and population health, the low genetic diversity within G. coprotheres populations is of concern and has implications for the future management and conservation of this species.

Based on collection locality samples were grouped into three geographic regions (Western, Middle and Northern regions see Fig. 1). Genetic diversity estimates varied among the three geographic regions ( Table 2). The number of alleles were the lowest in the Western region (45 alleles) and highest in the Middle region (134 alleles). Although the mean number of alleles differed among the three geographic regions, with the Middle region having the highest mean number of alleles (Ā = 10.308; SE = 1.082), when allelic richness was estimated using the rarefaction index (accounting for the differences in sample sizes among the regions) little difference in allelic richness was observed. The number of private alleles ranged from 0.41 alleles (Western) to 0.79 alleles (Northern). Similar levels of observed heterozygosity were observed across the three regions (observed heterozygosity <0. 39). No signs of excess heterozygosity were observed in any of the three regions. The Middle region and the Northern region showed signs of heterozygote deficiency (F > 0. 15).
In addition to the regional grouping, 266 samples collected from six breeding colonies were analysed separately. Genetic diversity also varied among the six breeding colonies ( Table 2). The number of alleles ranged from 45 alleles in Potberg to 58 alleles in Msikaba. Allelic richness estimates showed little difference across the six colonies. Similar levels of heterozygosity were observed in all six colonies (observed heterozygosity <0. 39). No signs of excess heterozygosity were observed in any of the six colonies. The inbreeding coefficient ranged from −0.045 in Collywobbles to 0.062 in Potberg.
South African representatives of G. barbatus (n = 54), N. monachus (n = 54), and G. africanus (n = 68) were also included and genotyped using the same microsatellite loci, to allow for direct comparison of genetic diversity values. The mean observed heterozygosity in G. coprotheres (Ho = 0.38) was much lower than that observed in G. barbatus (Ho = 0.50), N. monachus (Ho = 0.71), and G. africanus (Ho = 0.65). Gyps coprotheres and the much smaller, isolated South African G. barbatus population showed elevated levels of inbreeding (mean F = 0.17). Inbreeding was not detected in the other two vulture species N. monachus F = −0.07 and G. africanus F = 0.07 (Supplementary Table 2).  Table 2. Genetic diversity estimates for the 605 Gyps coprotheres grouped by geographic region and for the 266 Gyps coprotheres individuals grouped by breeding colony. Number of individuals (N), total number of alleles (A), mean number of alleles (Ā), allelic richness (A R ), private alleles (A P ), observed heterozygosity (Ho), unbiased expected heterozygosity (uHe), and inbreeding coefficient (F). Allelic richness for each region was based on the minimum number of gene copies. Standard errors are shown in parentheses.
www.nature.com/scientificreports www.nature.com/scientificreports/ population structure. Two optimal genetic clusters were recovered from the Bayesian clustering analysis conducted on the 605 G. coprotheses, including all 24 sampling localities (K = 2, delta K = 42.16; Supplementary Table 3). Figure 2 shows the structure barplots for all 605 G. coprotheses samples grouped by geographic region for K = 2 and K = 3. No distinct geographic structure was observed in the data across both K-values. However, significant genetic differentiation was observed in all region pairs shown in the pairwise F ST values reported in Supplementary Table 4.
The second Bayesian clustering analysis, conducted using only the 266 G. coprotheses collected at the six breeding colonies, recovered three genetic clusters (K = 3, delta K = 11.05; Fig. 2 Table 3). All three genetic clusters were again present across the six breeding colonies. In this case, the two colonies (Skeerpoort and Kransberg) are distinct from other colonies. Surprisingly, the geographically isolated Potberg colony in Western Cape, South Africa is not distinguishable from colonies in Eastern Cape and Kwazulu-Natal provinces. However, significant pairwise F ST values (p-value < 0.003) were recovered between the Potberg colony and all colony pairs and the Skeerpoort colony and all colony pairs as shown by pairwise F ST estimates reported in Supplementary  Table 5. Pairwise F ST analysis showed no clear genetic partitioning when individuals were assigned to localities (Supplementary Table 6). Mantel tests showed no correlation between pairwise geographic distances and pairwise genetic distances (individuals grouped by locality: R = 0.011, p-value = 0.317; individuals grouped by colony: R = 0.052, p-value = 0.105).

, Supplementary
Principle coordinate analysis (PCoA) was performed on G. coprotheses individuals and G. coprotheses samples grouped by populations (where population is either locality, region or colony groupings). Figure 3 shows the PCoA for all 605 G. coprotheres sampled. When all 605 G. coprotheres individuals were analysed ( Fig. 3 graph a) no correlation between individual genetic distance and locality was observed. Similarly, no correlation was observed when samples were grouped by the 24 sampling localities ( Fig. 3 graph b). In contrast, when individuals were grouped by the three geographic regions (Fig. 3 graph c), all three regions were genetically distinguishable.
In PCoAs for the 266 G. coprotheres individuals collected at the six breeding colonies (Fig. 4) no correlation between individual genetic distances and colonies were observed ( Fig. 4 graph a). The PCoA showing individuals grouped by colony ( Fig. 4 graph b), support the STRUCTURE results by grouping together the three colonies from the Eastern Cape and Kwazulu-Natal provinces (Collywobbles, Msikaba and Umzimkulu). Interestingly, the PCoA was able to separate the Skeerpoort and Kransberg colonies.
AMOVA results reported that the global F ST values for G. coprotheses individuals grouped by sampling locality, G. coprotheses assigned to geographic regions and only G. coprotheses collected at the six breeding colonies deviated significantly from zero (F ST > 0.01; p-value = 0.00). The majority of genetic variation, however, occurred within individuals (64-77%; Supplementary Table 7). population connectivity. Migration rates below 0.10 were used to indicate demographic independence 30 . The migration rates (m) for G. coprotheses grouped by geographic region showed that the highest gene flow occurred from the Middle region to the Northern region (m = 0.31; confidence interval (CI) = 0.30 to 0.32; Table 3). This is unsurprising as these birds are highly mobile and travel long distances to forage 14 and the majority of these samples were collected from feeding sites. Interestingly, the rates of migration between the Western and Northern regions are low (m < 0.04), indicating that these two regions are demographically independent www.nature.com/scientificreports www.nature.com/scientificreports/ from each other. The results suggest that the Middle region may act as a source population for both the Western and Northern regions.
The migration rates between G. coprotheses collected at the six breeding colonies showed that the highest migration rates were observed out of the Collywobbles colony to all other colonies (m > 0.20; Table 4). Potberg, Msikaba, Umzimkulu, Skeerpoort and Kransberg colonies are all demographically independent from each other (m < 0.10).
population bottleneck. Table 5 shows the results for the bottleneck analysis for all 605 G. coprotheses grouped by geographic region. The heterozygote excess method (Hx) for both mutation models showed no signs of recent bottleneck (p-value > 0.003). No significant deviation from the normal L-shaped distribution (Mode-shift) was observed in any of the regions. A similar result was observed when G. coprotheses from the six www.nature.com/scientificreports www.nature.com/scientificreports/ breeding colonies were analyzed ( Table 5). The heterozygote excess method for both models showed no sign of a recent bottleneck (p-value > 0.05); neither did the Mode-shift test. Both the Wilcoxon test for heterozygote excess and the Mode-shift provide strong evidence that G. coprotheres populations, when considering both individuals grouped by region and only those individuals collected at the six breeding colonies, have not undergone a recent bottleneck.
Interestingly, two regions (Middle and Northern) showed significant heterozygosity deficiency (p-value < 0.003) for the Wilcoxon test under both mutation models. Additionally, the sign tests for the Middle region and Northern region indicated significant heterozygosity deficiency (p-value < 0.003) under both mutation models. Both regions significantly deviate from the expected ratio (1:1 for heterozygous deficiency to heterozygous excess) for non-bottlenecked, mutation equilibrium populations. When G. coprotheres was analyzed as a single population, a similar result was observed with significant heterozygosity deficiency (p-value < 0.003) for the Wilcoxon test under both mutation models. G. coprotheres population deviated from the expected ratio (1:1 for heterozygous deficiency to heterozygous excess) for non-bottlenecked, mutation equilibrium populations. None of the colonies deviated from the expected ratio for any of the tests. However, when all individuals collected from the six breeding colonies were observed as a single colony unit, significant deviation from the expected Wilcoxon test for heterozygous deficiency and sign test under SMM was observed. The heterozygous deficiency indicates that when the breeding colonies are considered as a single colony unit, they do not behave in mutation-drift equilibrium.
Effective population size (Ne) was estimated assuming both monogamous and random mating for G. coprotheses. The monogamous mating model estimated Ne = 409 individuals (confidence interval (CI): 318.5; 537.2) the ratio of effective population size to census population size (Ne/N) 31 was 0.044. When the random mating model was selected Ne = 208 individuals (CI: 161.4; 274.6) Ne/N = 0.022. This is much lower than that reported number of G. coprotheses in recent census data (approximate 9400 individuals) 19 . www.nature.com/scientificreports www.nature.com/scientificreports/ Discussion G. coprotheres has a relatively narrow distribution in comparison to other African vultures, which have continental scale distributions. As such this species provides a unique opportunity for us to intensively sample across the distribution of the species and thoroughly assess the genetic diversity of the species. This study represents the first comprehensive genetic assessment of the endangered G. coprotheres and provides important genetic baseline data for the conservation and management of the species.
Low levels of observed heterozygosity (mean Ho = 0.38) in G. coprotheres were recorded in this study. Genetic diversity estimates of G. coprotheres were also lower than that observed in the three other ecologically similar vulture species (G. barbatus, N. monachus and G. africanus). Although the sample sizes of the other three species were much smaller than that used for G. coprotheres the same genetic markers were amplified across all vulture species examined, allowing for direct comparison of levels of genetic diversity. This study found higher allelic diversity (number of alleles) in G. coprotheres probably as a result of differences in sample size, but, alarmingly N. monachus and G. africanus showed higher mean observed heterozygosity estimates. N. monachus revealed significantly higher levels of heterozygosity and low levels of inbreeding (F = −0.07). The genetic diversity values for G. africanus and G. barbatus estimated in this study are similar to the values reported in the literature for the Namibian G. africanus population 7 and European G. barbatus populations 32 . This provides support for the results recovered in the present study not being biased by sample size or marker choice, but rather reflect the true levels of genetic variation present in populations of these African vultures. The lower than expected genetic diversity in G. coprotheres is of concern, given that both N. monachus and G. africanus are currently at higher extinction risk than G. coprotheres. N. monachus and G. africanus are currently listed as Critically Endangered by the IUCN 10 . The genetic diversity observed in G. coprotheres was lower than that observed in other threatened raptors such as the Spanish imperial eagle (Aquila adalberti) 3 ; and the Eurasian black vulture (Aegypius monachus) 33 .
The Bayesian assignment tests indicated very shallow genetic structure across the range of the species. However, the pairwise F ST values suggested genetic partitioning at a regional scale. Vultures from the two northern colonies (Skeerpoort and Kransberg) were genetically distinct from the other colonies. Given the geographic isolation and small size of the Potberg colony, we expected this colony to be the most effected by genetic drift. Although no marked allele frequency differences were observed in the STRUCTURE analysis, this colony did show the highest private allelic richness which is indicative of restricted gene flow. This is supported by the pairwise F ST values in Supplementary Table 5 where Potberg is significantly dissimilar to the other five breeding colonies.
The shallow genetic structuring seen in this species was not unexpected; given that birds often show shallower genetic structuring than that of other vertebrate species [34][35][36][37] . The lack of genetic structuring in avian populations is attributed to high dispersal capabilities [38][39][40] . A similar pattern of low genetic structure was observed in Egyptian vulture (Neophron percnopterus) 41 and griffon vulture (G. fulvus) 42 . In contrast, other wide-ranging vulture and raptor species such as the bearded vulture G. barbatus 43 , the Eurasian black vulture (Aegypius monachus) 33 ; and the white-tailed eagle (Haliaeetus albicilla) 44 have shown high levels of population differentiation.
The degree of natal philopatry (the likelihood that individuals breed at or near their place of origin) can influence the extent of genetic structuring in animal populations. Data from radio tracking and ringing studies have suggested that G. coprotheres exhibit natal philopatry 27,28 . Analysis of individuals collected from the six breeding colonies across South Africa did not support strong natal philopatry in this species, rather the genetic data suggests regional philopatry. The migration and pairwise F ST results showed that the Potberg colony was demographically independent from all other colony pairs. The Skeerpoort and Kransberg colonies were also distinct, while the Collywobbles, Msikaba and Umzimkulu colonies seems to be operating as a single large regional unit.
G. coprotheres are listed as endangered by IUCN due to declines in overall population numbers 19 . Some recent surveys have however, reported increasing population number in some areas 45 . Populations that have experienced a recent bottleneck event are usually characterized by an excess of heterozygotes, as opposed to a population at mutation-drift equilibrium 46,47 . Over time however, inbreeding and genetic drift following a severe bottleneck event would lead to populations containing more homozygotes (heterozygote deficient). In the data generated in the present study, no signs of heterozygosity excess were observed in South African G. coprotheres populations. This could be because G. coprotheres may be currently experiencing a bottleneck and not enough generations have passed for the genetic signal of this event to be detected. In contrast, population genetic theory also suggests that G. coprotheres populations may have reduced genetic variation due to older bottleneck events 2,48 . This latter hypothesis was supported by results of Wilcoxon signed rank tests which found significant heterozygote deficiency in vultures from all three geographic regions and four colonies (Potberg, Collywobbles, Msikaba and Kransberg).
Identifying and prioritizing G. coprotheres populations, which have higher genetic variability, is an important step forward in conserving these birds. Genetically diverse populations could be valuable source populations for future translocation programmes. Six of the 24 localities showed elevated heterozygosity (F < 0; (Commando Drift and Elliot in Eastern Cape; Natal Midlands in KwaZulu-Natal; Soetdoring in Freestate; Blouberg in Limpopo; and Mala Mala in Mpumalanga), whereas all other localities showed reduced heterozygosity. The isolated Potberg colony in the Western Cape had the lowest genetic variation of all the breeding colonies. This is not unexpected given the isolated nature of this colony. Careful management and monitoring of this colony is needed to avoid the detrimental effects of inbreeding and genetic drift. Given the shallow genetic structuring observed across the range of the species and the high levels of gene flow among breeding colonies, the present study suggests that the entire South African G. coprotheres population be managed as a single management unit.
The results from this study showed low levels of genetic diversity and variability in South African G. coprotheres populations. This, together with the rapid decline of this species 9,21 , highlights the need for conservation strategies to include the maintenance of genetic diversity as a population management tool. Reductions in genetic diversity can have serious implications on the evolutionary potential of a species 49 . The present study highlights www.nature.com/scientificreports www.nature.com/scientificreports/ the utility of microsatellite markers for the assessment and monitoring of genetic diversity in an African vulture species. The data produced in this study can be used as a baseline for future genetic monitoring and species recovery programmes.

sampling.
A total of 605 G. coprotheres from 24 localities, across the South African distribution of the species, were sampled for this study (Supplementary Table 8). Based on sample localities, samples were grouped into three geographic regions: Western (n = 18), Middle (n = 462) and Northern (n = 125). In addition to the regional grouping, the 266 samples collected from six breeding colonies (Fig. 1) were analysed separately. South African representatives of G. barbatus (n = 54), N. monachus (n = 54), and G. africanus (n = 68) were also included and genotyped using the same microsatellite loci, to allow for direct comparison of genetic diversity values.
Samples consisted of feather samples, collected opportunistically from feeding sites, sites of electrocutions, poisoning events and below nests at six main breeding colonies 18 . Blood samples were also collected when vultures were captured and fitted with global positioning system/global system for mobile transmitters 14 . Bloods were stored on Whatman FTA ® Elute cards (Sanford, USA). Archival museum samples (dried skin snips) were sourced from local South African museums (Supplementary Table 8 Table 9). The 13 loci were amplified in six multiplex reactions (Supplementary Table 9) using the KAPA2G TM Fast Multiplex PCR kit (KAPA Biosystems). The 10 μl reactions consisted of ~2-30 ng template DNA, 5 μl KAPA2G Fast Multiplex mix, 0.2 μM of each primer, 0.1 μl of 1 mg/ml bovine serum albumin (BSA) and purified water. In reactions performed using DNA extracted from feather and archival samples, template DNA was increased to ~20-200 ng to improve amplification success. The cycling parameters for all loci followed the standard KAPA2G TM Fast Multiplex PCR kit protocol except for multiplex 1, where the annealing temperature was reduced to 58 °C. All amplified products were sent to the Central Analytical Facility at Stellenbosch University, South Africa for fragment analysis. The software GeneMarker v2.4.0 (Soft Genetics) was used for genotype scoring. To ensure genotyping consistency, all archival samples were re-amplified, and each locus was genotyped multiple times (up to five times). In addition, 20% of all feather, muscle and blood samples were re-amplified multiple times (up to five times) to verify the reliability of the data.

Data analysis. Assessing genetic variation.
To ensure that duplicated genotypes were not included when feathers were collected, identity analysis was performed in Cervus v3.0.7 52 . Cervus was also used to estimate polymorphic information content (PIC) for each locus. Null allele frequencies for each locus were estimated in FreeNA 53 using the expectation maximization (EM) algorithm 54 . Because null alleles can bias population structure analysis 55 , uncorrected global F ST were compared to F ST values corrected using the excluding null alleles (ENA) method 53 using a paired t-test. Linkage disequilibrium was tested using Genepop v4.2 56 . Deviations from Hardy-Weinberg equilibrium were estimated in GenAlEx v6.502 57 . The effective population size (Ne) was estimated in NeEstimator v2.01 58 implementing both random and monogamous mating.
Genetic diversity estimates (number of alleles, observed heterozygosity, unbiased expected heterozygosity and inbreeding coefficient) were calculated in GenAlEx. Tests for deviation from Hardy-Weinberg equilibrium were performed in Genepop. Allelic richness was estimated in FSTAT v2.9.3.2 59 . The observed number of alleles in a population is dependent on sample size and given the opportunistic nature of sampling used in this study sample size varied across localities. Consequently, allelic richness was calculated using a rarefaction index estimated in FSTAT and the number of private alleles within populations were estimated using the rarefaction method implemented in HP-RARE v1.0 software 60 . All multiple comparisons were adjusted using the Bonferroni correction. population structure. Bayesian assignment tests were performed in STRUCTURE v2.3.4. 61 . For each analysis ten independent runs were performed consisting of 100000 Markov chain Monte Carlo (MCMC) replicates with a burn-in of 10000 and the proposed number of genetic clusters (K) ranging from one to ten. The admixture ancestry model with correlated allele frequencies was selected for all runs. Sampling locality information was incorporated using the LOCPRIOR model. STRUCTURE Harvester 62 was used to estimate the most probable number of genetic clusters 63 . Membership probabilities (Q-values) for each individual and for each genetic cluster were estimated using ClumpAK 64 . Bar plots from STRUCTURE runs were created in Pophelper 65 .
The correct identification of the number of genetic clusters can be challenging if there is unbalanced sampling or complex phylogeographic structure i.e. when there are unequal genetic distances among subpopulations or unbalanced sample sizes 66 . For this reason, Principle Coordinate Analysis (PCoA) was also performed in GenAlEx using pairwise codominant genotypic genetic distance. PCoA does not rely on a model but rather multidimensional scaling to visualize similarity in a dataset.
Analysis of molecular variance (AMOVA) was performed in GenAlEx. Three AMOVA's were conducted, one grouping individuals by sampling locality, a second grouping individuals by geographic region. A third AMOVA www.nature.com/scientificreports www.nature.com/scientificreports/ was conducted only on individuals collected from the six breeding colonies. Pairwise F ST values were estimated in FSTAT. Mantel tests 67 were performed in GenAlEx, to test for correlation between geographic distance and genetic distance. population connectivity. Migration rates (gene flow) between the three geographic regions and the six breeding colonies were estimated in BayesAss v1.3 68 . Two analyses were performed consisting of 10000000 iterations, with a burn-in of 1000000 iterations and a sampling frequency of 90000 iterations. Multiple runs were performed with varying run lengths to ensure MCMC chain convergence. The delta values for each parameter were adjusted to achieve a 20-60% acceptance rate 68 . In analysis to estimate gene flow among the three geographic regions, the final delta values were delta allele frequency = 0.30, delta migration rate = 0.10 and delta inbreeding coefficient = 0.50. In analysis to estimate gene flow among the six breeding colonies, the final delta values were delta allele frequency = 0.50, delta migration rate = 0.10 and delta inbreeding coefficient = 0.70. For both analyses, migration rates below 0.10 were used to indicate demographic independent populations 30 . population bottleneck. BOTTLENECK v1.2.02 69 was used to test for evidence of heterozygosity excess (Hx) 5,47 in G. coprotheres populations. Two analyses were performed, one grouping individuals by geographic region, and a second only including individuals collected at the six breeding colonies. Each analysis ran two mutation models. The conservative stepwise mutation model (SMM) and the two-phase model (TPM) with 90% stepwise mutation 70,71 and a variance of 12 69 to include the observed range of multistep mutations in natural populations 72 . Two statistical tests (Wilcoxon sign-rank test and the sign test) as well as a qualitative test (Mode-shift test), were performed. The Wilcoxon sign-rank test is effective in detecting recent declines in effective population size (Ne) and assumes that populations that have recently undergone a decline in effective population size will have a higher level of heterozygosity compared to a population at mutation-drift equilibrium 5 . The sign test was used to identify the number of microsatellite loci that have either heterozygosity excess (Hx) or heterozygosity deficiency (Hd). Changes in heterozygosity (excess or deficit) can occur after a recent change in the effective population size or if heterozygotes have a selective advantage or disadvantage 69 . In populations experiencing a bottleneck event, alleles are usually lost faster than heterozygosity, thus producing a heterozygote excess 47 . The Mode-shift is a qualitative indicator that can differentiate between bottlenecked and stable (non-bottlenecked) populations. This test can detect the genetic changes caused by a population decline within a few dozen generations 73 . Therefore, only recent bottleneck events will be identified using this method.

Data Availability
The data analysed in this study will be made available as Supplementary information upon publication.