3 RESULTS

3.1 Sequencing and genotyping quality

We succeeded to sequence a 390 bp fragment of the mtDNA CR from 82 out of the 86 golden eagles. Four individuals were excluded due to low sequence quality. With the addition of five Finnish individuals from Kylmänen et al. (2023), we created a mtDNA dataset of 87 golden eagles. All 86 individuals were genotyped with 12 polymorphic loci (Table S6). Including the five Finnish individuals (Kylmänen et al. , 2023), the total microsatellite dataset consisted of 91 individuals. Mean genotyping success was 87.6%, being the lowest for loci Aa35 (51.6%) and Aa11 (67.0%). However, since most samples were from old museum specimens that usually have lower genotyping success because of the degraded DNA (Tsai et al. , 2020), and since these loci had high PIC values (0.74 and 0.59, respectively), we included them into the analyses. The number of alleles varied between 2 and 13 per locus and was on average 7.8. The error rate was 0.036 per allele and 0.078 per locus. The mean ADO rate was 0.056 and the mean FA rate was 0.022. MicroChecker suggested the presence of null alleles in loci Aa04, Aa11, Aa35, Aa36, Aa39, and NVHfr142, and stuttering in loci Aa35, Aa36, and Aa39. However, these results were not consistent when the four geographical groups were analyzed separately, indicating that the excess of homozygotes was not a genotyping artefact and was rather attributed to undetected population structure (e.g., Wahlund effect; Garnier-Géré & Chikhi, 2013). The mean null allele frequency over all loci was 0.076. No consistent significant deviations from HWE and linkage equilibrium were observed when groups of different datasets were tested separately. Therefore, all loci were kept for the downstream analyses.

3.2 Spatial genetic variation

3.2.1 Genetic diversity

Genetic diversity of Eurasian golden eagles is summarized in Table 1. From the dataset of 434 individuals, that included the GenBank sequences, we identified 40 haplotypes. Overall, nucleotide diversity was higher than theta (π = 0.0164, θ = 0.0079), and this pattern remained for Central and Eastern Europe and Central Asia and Caucasus. These two groups also had the highest nucleotide diversity. Far East exhibited the highest haplotype diversity and the highest number of haplotypes (Figure S2A), while its nucleotide diversity was relatively low. The lowest mitochondrial diversity was in Northern Europe, with notably lower nucleotide diversity compared to theta (π = 0.0034, θ = 0.0087). Between the lineages, each mitochondrial genetic diversity parameter was lower in the Mediterranean group than in the Holarctic group, including the number of haplotypes after the rarefaction-extrapolation analysis (Figure S2B).
When we estimated mitochondrial diversity using only our data (N = 87), the sample sizes among the groups became more balanced, yet smaller (Table S7). With these, we identified 23 haplotypes (GenBank accession numbers: OR635080–OR635102; Table S8). In the entire Eurasian population, nucleotide diversity was lower than theta (π = 0.0067, θ = 0.0112), and the pattern remained in all groups except for Central Asia and Caucasus. The latter group also had the highest haplotype and nucleotide diversities (h = 0.752, π = 0.0141), while the lowest nucleotide diversity was in Northern Europe (π = 0.0024). Interestingly, Far East, now represented only by individuals from continental Eurasia (i.e., excluding the Japanese population), exhibited the lowest haplotype diversity (h = 0.571). Individuals carrying the Mediterranean lineage had notably higher haplotype and nucleotide diversities compared to the Holarctic lineage, but these groups had unbalanced sample sizes (N = 8 and 79 for Mediterranean and Holarctic, respectively).
From the microsatellite data, the HO in Eurasia was 0.504, and HE 0.633. HO was lower than HE in all analyzed groups (Table 1). Among the four geographical groups, Far East had the highest HO (0.597), and the second highest HE (0.651; the highest HE = 0.659 was in Central Asia and Caucasus). Inbreeding coefficient (FIS) was positive in all groups and varied between 0.098 (Far East) and 0.242 (Central and Eastern Europe). The FIS significantly deviated from zero in all groups except in Far East (Northern Europe: V = 49, p = 0.032; Central and Eastern Europe: V = 68, p = 0.021; Central Asia and Caucasus: V = 67, p = 0.027; Far East: V = 62, p = 0.077). When comparing genetic diversity among the groups, significant difference was found only in HO between Central and Eastern Europe and Far East (W = 36.5, p = 0.04; Table S9). Allelic richness was similar across the four groups (on average 4.23). The number of private alleles per locus (i.e., private allelic richness, PAR) in Central Asia and Caucasus was nearly twice as high as in the other groups, while Northern Europe had the lowest PAR (Table 1, Figure S3A). The two mitochondrial lineages showed similar levels of heterozygosity and inbreeding coefficients (Table 1). As these groups had a large difference in sample sizes, the results must be interpreted with caution. Nevertheless, no statistical difference was observed in the estimates of A, HO, HE, and FIS between the mitochondrial lineages (Table S9). On the other hand, the Holarctic group had lower allelic and private allelic richness compared to the Mediterranean group (AR = 3.81 and 4.39, PAR = 1.08 and 1.34, respectively, corrected for sample size; Figure S3B). FIS significantly deviated from zero (V = 74, p = 0.003) in the Holarctic group, but not in the Mediterranean group (V = 52, p = 0.100).

3.2.2 Population structure

A median-joining haplotype network of 663 golden eagles from across the Northern Hemisphere revealed 56 haplotypes clustered into two mitochondrial lineages: Mediterranean and Holarctic (Figure 1). The Mediterranean lineage had one central haplotype M1, while the Holarctic lineage had two central haplotypes differing by one SNP: H1 in Eurasia and CR1 mainly in North America. We discovered eight new haplotypes from Russia and Central Asia (RUS1, RUS2, RUS3, KAZ1, UZB1, IRN1, KYR1, and KYR2). Haplotype CR4, previously reported only in North America, was also found in our dataset from Kamchatka. Moreover, we identified several haplotypes from the North American cluster (H6–H9) in golden eagles sampled from Far East and Central Asia and Caucasus. Thus, of the 56 haplotypes thus far detected in golden eagles, geographically three were truly Holarctic (CR1, CR4, N12), 16 Nearctic, and 37 Palearctic. Golden eagles with Mediterranean haplotypes were found exclusively in southern Eurasia, spanning between Spain and eastern Kazakhstan, except for one individual from Northern Finland (Figure 2). On the other hand, Holarctic haplotypes were found across Eurasia (Figure 2).
From the STRUCTURE results, the most likely number of genetic clusters was two when geographical locations were used as LOCPRIOR (Table S10, Figure S4). Central Asia and Caucasus was genetically distinct from the other groups (Figure 3A). On the contrary, DAPC with the pre-defined geographical groups did not reveal nuclear genetic structure (Figure 3B). Nevertheless, Northern Europe and Central Asia and Caucasus were somewhat separated along the Discriminant function 1. Also, PCA showed some level of differentiation of Central Asia and Caucasus (Figure 3C). The uniqueness of Central Asia and Caucasus and Northern Europe were also consistent with log genotype probability (LGP) plots of pairwise comparisons of the four geographical groups (Figure S5). Individuals with Mediterranean and Holarctic haplotypes (N = 87) were subdivided into two nuclear clusters by STRUCTURE when the lineage was used as LOCPRIOR (Figure 3D, Table S10, Figure S6). Nuclear differentiation of the Mediterranean and the Holarctic groups was also supported by DAPC and GenePlot results (Figure 3E, F). No population differentiation was found using the de novo STRUCTURE (Table S10, Figure S7). In thede novo DAPC, K = 2 was selected as the most likely number of clusters, but the assignment was done with only 10 PCs (Figure S8A). Besides, no spatial pattern was observed for the suggested grouping (Figure S8B). No IBD pattern was detected in Eurasia with the Mantel test (Figure S9) nor was there any spatial autocorrelation (Figure S10).
AMOVA analyses with the GenBank sequences estimated 53.8% of mitochondrial genetic variation to be among the five geographical groups (p < 0.001). Furthermore, the pairwise ɸSTvalues were significant and high, pointing to existing genetic differentiation among these groups (Table 2). However, when the same analyses were done using only our data, only 22.8% of mitochondrial variation was among the groups (p < 0.001), and the pairwise ɸST ranged from -0.012 to 0.329, being significant only for Central Asia and Caucasus (Table S11). On the other hand, from the microsatellite data, AMOVA assigned only 1.1% of nuclear variation to among the geographical groups (p = 0.07), and the pairwise microsatellite FST values were low and varied from -0.003 to 0.038 (Table 2). Only Northern Europe was significantly, yet weakly, differentiated (Table 2). The mitochondrial lineages were highly differentiated according to pairwise ɸST estimated with the GenBank sequences (ɸST = 0.934, p < 0.001) and with our data only (ɸST = 0.925, p < 0.001). However, the nuclear FST between the Mediterranean and Holarctic lineages was low and not significant (FST = 0.020, p = 0.120).

3.3 Temporal genetic variation and demographic history

Genetic diversity of the temporal groups is presented in Table 1. The Bottleneck group exhibited slightly higher haplotype and nucleotide diversities and a slightly lower theta (h = 0.754, π = 0.0141, θ = 0.0121) than the Post-bottleneck group (h = 0.800, π = 0.0174 θ = 0.0101). On the contrary, rarefaction-extrapolation analyses suggested a substantially higher number of haplotypes in the Bottleneck period compared to Post-bottleneck (Figure 4A). When we subdivided the temporal groups into Mediterranean and Holarctic lineages, we found that mitochondrial genetic diversity was higher in the Bottleneck subgroups (Table 1). With only our data (N = 77), nucleotide diversity and theta were higher and haplotype diversity was lower in the Bottleneck group (h = 0.587, π = 0.0066, θ = 0.0103) than the Post-bottleneck group (h = 0.859, π = 0.0042, θ = 0.0074; Table S7).
Observed and expected heterozygosities were on a similar level in both groups (HO = 0.510 and 0.518, HE = 0.628 and 0.610 for the Bottleneck and Post-bottleneck groups, respectively). Moreover, there were no significant differences in nuclear genetic diversity (A, HO, HE, and FIS) between the temporal groups (Table S9), yet the inbreeding coefficient was lower in the Post-bottleneck group (FIS = 0.146 for Post-bottleneck and 0.195 for Bottleneck). In both temporal groups, the inbreeding coefficients were positive and significantly deviated from zero (Bottleneck: V = 65, p = 0.005; and Post-bottleneck: V = 57, p = 0.032). Allelic and private allelic richness were higher in the Bottleneck group (AR = 4.86 and 4.63, PAR = 1.08 and 0.95; Figure S3C).
A temporal haplotype network of 393 golden eagles from across Eurasia (Figure 4B) showed that the Post-bottleneck group was missing 14 haplotypes in comparison to the Bottleneck group. Nine of these haplotypes were from the Holarctic, and five from the Mediterranean lineage. Meanwhile, ten haplotypes not found in the Bottleneck group were discovered in the Post-bottleneck group, of which seven were Holarctic and three Mediterranean. Central haplotypes in both Mediterranean and Holarctic lineages remained the same throughout time.
From the demography analyses, signs of a previous population expansion were observed in the Bottleneck Mediterranean subgroup as indicated by the significantly negative Tajima’s D and Fu’s Fs, a small significant R2 (Table 1) and the shape of the mismatch distributions (Figure S11). The BSP analyses revealed a decline in the effective female population size (Nef) in Eurasia around 1975 (Figure 5). The BSP results for the Mediterranean lineage also showed a slight declining trend during that time, while in the Holarctic lineage, the Nef remained constant over the last four centuries (Figure 5). Wilcoxon’s sign rank test results were not indicative of a bottleneck in any of the analyzed groups, and only the Mediterranean group had signs of population bottleneck based on the mode-shift test (Table S12).
Among the five geographical groups, the bimodal and ragged shapes of mismatch distributions were detected for Central and Eastern Europe, Western Europe, Central Asia and Caucasus, and Far East, suggesting admixed and stable populations. Only Northern Europe showed a unimodal distribution, suggesting past population expansion (Figure S12), further supported by significantly negative Fu’s Fs and R2, also found in Far East (Table 1).
While both DAPC and GenePlot showed some level of nuclear genetic differentiation between the temporal groups (Figure S13), the STRUCTURE results suggested K = 1 as the most likely number of clusters with temporal groups as LOCPRIOR (Figure S14). When we used only our data (N = 77), neither nuclear pairwise FST nor the mitochondrial ɸST between the temporal groups were significant (FST = 0.011, p = 0.158; ɸST= 0.011, p = 0.265). With the addition of the GenBank sequences, the ɸST between Bottleneck and Post-bottleneck was small yet significant (ɸST= 0.075, p = 0.0001). When the temporal groups were further subdivided according to the mitochondrial lineage, the ɸST values were high and significant only between the lineages, and not between the temporal periods (Table S13).