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).