Holarctic phylogeography of golden eagles (Aquila chrysaetos) and evaluation of alternative North American management approaches

Holarctic phylogeography of golden eagles (Aquila chrysaetos) and evaluation of alternative North... Abstract The golden eagle (Aquila chrysaetos) is a long-lived bird of prey with a Holarctic distribution. This species has survived severe anthropogenic stressors that have reduced or eliminated populations in some parts of its range. Despite the ecological and cultural importance of the golden eagle, few attempts have been made to determine the partitioning of genetic variation over large areas of its range. This study generated DNA sequence data from the mitochondrial DNA control region from 115 North American golden eagles and combined these data with existing control region sequences from over 300 Nearctic, Palearctic and Mediterranean golden eagles to provide a clearer holistic picture of the Holarctic phylogeographical patterns of genetic variation in this species and the genetic variation and demographic history of golden eagles in North America. The results support that there are two genetic lineages of golden eagles, one representing the Mediterranean and the other occurring throughout the Holarctic. The Holarctic lineage reveals little to no modern-day gene flow between Nearctic and Palearctic golden eagles. Furthermore, the current distribution of haplotypes in the Nearctic shows a recent population expansion with moderate levels of gene flow. INTRODUCTION Climate perturbations throughout history have had effects on the evolutionary history and phylogeography of many species. The climate fluctuations of the Quaternary (~1.8 Mya to present) have been well studied and shown to have impacted the current distributions and speciation events of many taxa in the Northern Hemisphere (Hewitt, 2000, 2004). Raptors throughout the Holarctic display varying patterns of genetic differentiation with much of this variation being attributed to how they were impacted during the Pleistocene. When considering the modern-day phylogeography of raptors, an east to west differentiation of haplotypes is commonly observed. Cooper’s hawk (Accipiter cooperi; Sonsthagen et al., 2012b) and sharp-shinned hawk (A. striatus; Hull & Girman, 2005) in North America, as well as white-tailed eagle (Haliaeetus albicilla; Hailer et al., 2006), bearded vulture (Gypaetus barbatus; Godoy et al., 2004) and Eurasian black vulture (Aegypius monachus; Poulakakis et al., 2008) in the Palearctic exhibit this trend. This contrasts with northern goshawk (A. gentilis) in North America in which an east–west differentiation was not observed. Instead, the south-western populations differed significantly from the other sampled locations (Bayard de Volo, 2008). When the range of a species stretches across the Holarctic, typically two patterns of genetic structuring between the Palearctic and Nearctic exist. In the first pattern, American and Eurasian lineages live in sympatry with a distinct American haplogroup resulting from recent gene flow across Beringia or the Bering Strait (Nebel et al., 2015). Alternatively, due to a lack of gene flow across Beringia or the Bering Strait, speciation occurs and distinct species are found in both North America and Eurasia (Nebel et al., 2015). Golden eagle (Aquila chrysaetos) is the only eagle species found throughout the Holarctic. Golden eagles are able to survive in diverse environments by adapting to the available prey, landscape and their choice of nesting sites. For example, golden eagles have adapted to nesting on the ground in tree-less steppe and desert habitat as opposed to cliffs and trees in other parts of the world, and their diet includes diverse prey such as ptarmigan, hares, tortoise, snakes and lizards (Watson, Brockie & Watson, 2010). Furthermore, golden eagles are highly mobile at different stages of their lives or during different times of the year. Juvenile eagles are nomadic for approximately 5 years until they become sexually mature, at which time it has been shown, through banding studies, that North American individuals travel back to within approximately 46–175 km of their natal location to establish a breeding territory (Millsap et al., 2014). North American golden eagles that breed in northern latitudes are obligatory migrants, as compared to their southern counterparts which do not migrate, and will move to more southern latitudes in autumn and return to their nesting areas in spring. Although they are thought to be highly philopatric to their natal site (Millsap et al., 2014), these movements may allow for gene flow throughout the range of the species through behaviours such as extra-pair copulations or dispersal events. Little is known about the phylogeographical patterns of genetic variation in this species. Nebel et al. (2015) provided insight into the phylogeographical structure of Palearctic golden eagles through examination of 402 bp of the mitochondrial DNA (mtDNA) control region, but the North American subspecies was poorly represented due to a lack of available samples. Nebel et al. (2015) identified 26 haplotypes, only one of which was restricted to North America, and they speculated that there were two historical refugia. Sonsthagen et al. (2012a) and Craig et al. (2016) also utilized this region of the mtDNA as well as microsatellite loci to analyse golden eagles in the far western portion of their North American range. Sonsthagen et al. (2012a) identified five novel haplotypes and high gene flow when they analysed golden eagles in California and the Channel Islands, while Craig et al. (2016) identified the five hapolotypes previously detected by Sonsthagen et al. (2012a) and three novel haplotypes when analysing samples from Alaska, Idaho and California. Other genetic studies on golden eagles have utilized microsatellites, allozymes and mtDNA, but these studies focused on very limited portions of the species range (Masuda et al., 1998; Suchentrunk, Haller & Ratti, 1999; Wink et al., 2004; Bielikova et al., 2010; Bourke et al., 2010). Throughout their range, golden eagles have faced many anthropogenic stressors including illegal shooting, windfarms, lead poisoning, habitat loss, organic pollutants and game keepers driving the population into bottlenecks or causing local extinctions (Fielding, Whitfield & McLeod,, 2006; Bourke et al., 2010; Stauber et al., 2010; Watson et al., 2010; Pagel et al., 2013). For example, breeding populations of golden eagles have been extirpated east of the Mississippi River in the United States (Morneau et al., 2015), in Ireland (Bourke et al., 2010), and in the Alpine foothills and lowlands of Germany, Austria, the Czech Republic and Poland (Nebel et al., 2015). These bottlenecks and the population fragmentation caused by regional and local extirpations can allow for inbreeding depression and loss of genetic variation, adding to the need of determining the current partitioning of genetic variation for management purposes. In the United States there are two proposed models for developing Eagle Management Units (U.S. Fish and Wildlife Service, 2016). The first management proposal includes the migratory flyways found in North America (Atlantic, Mississippi, Central and Pacific) with the Mississippi and Atlantic flyways being combined. The alternative proposal is based on the currently recognized Bird Conservation Regions (BCRs) with support coming from dispersal from the nesting site (U.S. Fish and Wildlife Service, 2016). Recent studies utilizing nuclear markers (Craig et al., 2016; Doyle et al., 2016; Van Den Bussche et al., 2017) reveal similar groupings of golden eagles across studies, with some support leaning towards the BCR model. Generally, these studies support the genetic uniqueness of golden eagles in northern California, southern Idaho and southern Oregon, and another group consisting of golden eagles from western Canada and Alaska. They differ from the BCR model by grouping most of the central to eastern portion of the western golden eagle range into a single unit (Doyle et al., 2016; Van Den Bussche et al., 2017). Unfortunately, due to a lack of sufficient sampling across the North American range of golden eagles it is unknown if this pattern is an artefact of sampling or if these individuals should be grouped into a single management unit. This study utilized a fragment of the mtDNA control region previously used by Nebel et al. (2015), Sonsthagen et al. (2012a) and Craig et al. (2016) to generate DNA sequence data from 115 North American golden eagles representing Alaska, Alberta, Arizona, California, Colorado, Idaho, Montana, Nebraska, New Mexico, Oklahoma, Oregon, Wyoming, Yukon and British Colombia. Our goals were to fill in geographical gaps across the range of North American golden eagles to provide a clearer picture of the Holarctic phylogeographical patterns of genetic variation of golden eagles. Additionally, we were interested in evaluating whether the mitochondrial sequence data provided insight into the most appropriate geographical structure for managing golden eagles by evaluating levels of genetic differentiation when the samples were partitioned by Migratory Flyways, BCRs or Genetic Regions. Finally, we were interested in elucidating the demographic history of golden eagles in North America to answer questions regarding the impacts of population bottlenecks and subsequent expansions on North America golden eagles. This additional information will allow for a more complete understanding of the historical movements and current patterns of mitochondrial genetic variation for golden eagles throughout the Holarctic. MATERIALS AND METHODS DNA isolation and sequence development Working with permited U.S. Fish and Wildlife biologists and wildlife rehabilitators, we obtained blood samples from 124 North American golden eagles. When nests contained more than one individual, only one sample was included in the study to eliminate over-representation of a haplotype. Individuals that were not collected as a hatchling had natal location determined by GPS transmitters or isotope data. Only two samples (one individual from Colorado and a second individual from Wyoming) may not have a true natal origin represented as they were collected during December and natal location data were unavailable. Blood samples (~ 0.5 mL) were collected, stored in lysis buffer (Longmire, Maltbie & Baker, 1997) and subsequently shipped to our laboratory at Oklahoma State University where DNA was extracted following the protocol of Longmire et al. (1997). DNA quality was assessed by running an aliquot on a 1% agarose gel and quantified using a NanoDrop 3300 spectrophotometer (Thermo Scientific). A 442-bp region of domain I and II of the hypervariable control region-1 (D-loop) was amplified using polymerase chain reaction (PCR) and primers GOEA_CR1L and GOEA_CR595H developed by Sonsthagen et al. (2012b). PCRs contained 6 µL of GoTaq Flexi Buffer (Promega), 25 mM MgCl2, 1 µM of each primer, 10 mg/mL bovine serum albumin, 4 mM dNTPs, 1 unit of GoTaq (Promega) and double distilled (dd)H2O to make a 30-µL reaction. The thermoprofile was 95 °C for 5 min, 35 cycles of 95 °C for 1 min, 40 °C for 1 min, 72 °C for 1 min and a final extension at 72 °C for 7 min. PCR products were visualized by running on a 1% agarose gel and subsequently cleaned using the Promega Wizard Kit. Cleaned products were sequenced using BigDye v3.1 chain terminators (Thermo Fisher), 5× sequencing buffer (Edge Bio), GOEA_CR1L primer and ddH2O. The sequencing thermoprofile consisted of 35 cycles of 96 °C for 30 s, 50 °C for 35 s and 60 °C for 4 min. Sequencing products were cleaned using either Performa DTR V3 96 well short plates or Performa DTR Gel Filtration cartridges (Edge Bio). The cleaned sequencing product (~20 µL) was added to 1.0 µL of HiDi (Life Technologies) followed by electrophoresis using an ABI 3130 Genetic Analyzer (Applied Biosystems). Base calling was performed using Sequencing Analysis 5.2 (Applied Biosystems). Sequences were visually inspected and those sequences that were not long enough or had ambiguous bases were resequenced using the alternative primer (GOEA_CR595H) and de novo aligned in Geneious 7.1.9 to acquire the full 442-bp region. In addition to the sequences we generated, we also downloaded from GenBank the Old World golden eagle sequences generated by Nebel et al. (2015; accession numbers KR259251–KR259276), the southern California sequences generated by Sonsthagen et al. (2012a; accession numbers JQ246417–JQ246421), and the sequences from California, Idaho and Oregon generated by Craig et al. (2016; accession numbers KX687705–KX687711). All sequences were aligned using the ClustalX multiple sequence alignment in Geneious 7.1.9. The resulting alignment was imported into MacClade version 4.06 (Maddison and Maddison, 2005) for visual inspection and to trim our generated sequences as well as those of Sonsthagen et al. (2012a) and Craig et al. (2016) to 402 bp so that we could directly compare these data with those of Nebel et al. (2015). Using the Redundant Taxa option in MacClade we grouped sequences into haplotypes and assigned names to the haplotypes. Further analyses were conducted utilizing the data of Craig et al. (2016) and Nebel et al. (2015) because the Sonsthagen et al. (2012a) haplotype frequencies and locations were not available. Furthermore, approximately half of the golden eagles represented in Sonsthagen et al. (2012a) were from the Channel Islands, a location that only acquired golden eagles in the mid-1990s when bald eagles (Haliaeetus leucocephalus) were extirpated from the island, meaning the natal locations of the eagles that moved into the area are unknown. Data analysis Relationships among control region haplotypes of golden eagles were evaluated through the construction of a minimum spanning network at the 95% confidence level using TCS v2.1 (Clement, Posada & Crandall, 2000). As described below, analyses were conducted on: (1) the entire data set including the samples from our dataset, Craig et al. (2016) and Nebel et al. (2015); (2) data partitioned by geography into Nearctic, Palearctic and Mediterranean; and (3) with the Nearctic samples further broken down by the proposed Migratory Flyways (U.S. Fish and Wildlife Service, 2016), BCRs (U.S. Fish and Wildlife Service, 2016) or Genetic Regions suggested by nuclear markers in previous studies (Craig et al., 2016; Doyle et al., 2016; Van Den Bussche et al., 2017). To assess genetic diversity for each of the three data partitions, we calculated haplotype and nucleotide diversity in Arlequin 3.1 (Excoffier & Lischer 2010). We also calculated the allelic richness using Contrib 1.4 (Petit, El Mousadik & Pons, 1998) in which the population sample size was normalized using a rarefaction correction. The rarefaction value was set to the smallest population size. We performed neutrality tests, Tajima’s D and Fu’s Fs, in Arlequin 3.1 to confirm the selective neutrality for the region. Tajima’s D and Fu’s Fs were further used to determine demographic histories. We used an analysis of molecular variance (AMOVA) and pairwise ΦSTusing the Kimuira two-parameter model (K80) of evolution in Arlequin 3.1 to assess the partitioning of genetic variation within and among the Nearctic, Palearctic and Mediterranean sample localities. For the North American samples, we utilized a hierarchical AMOVA to test the three proposed conservation models (Migratory Flyways, BCRs, Genetic Regions). The hierarchical approach grouped samples by State or Canadian Province and then the States or Provinces were grouped according the Migratory Flyways, BCRs or Genetic Regions. RESULTS We were able to generate DNA sequence data from 115 North American golden eagles (S1) which combined with the samples from Nebel et al. (2015) and Craig et al. (2016) provided us with a total of 417 samples. The 402-bp region contained 27 polymorphic sites and resulted in 44 haplotypes. The 115 Nearctic samples that we generated DNA sequences from provided 19 haplotypes (Table 1). Of these 19 haplotypes, 11 were novel to this study, four previously described by Sonsthagen et al. (2012a) and Craig et al. (2016), three previously described by Craig et al. (2016), and one previously described by Sonsthagen et al. (2012a), Nebel et al. (2015) and Craig et al. (2016). Although Nebel et al. (2015) classified their haplotypes as Holarctic and Mediterranean, based on our increased sampling of North American golden eagles, only haplotype H5 is Holarctic. To provide greater clarity for understanding golden eagle phylogeographical patterns, we used the following haplotypic designations. The single Holarctic haplotype (H5 of Nebel et al., 2015) is represented as ‘H1’, and all other Holarctic haplotypes of Nebel et al. (2015) are Palearctic and identified in this study by the letter ‘P’. We retained the Mediterranean haplotypic designation of Nebel et al. (2015) and assigned it with the prefix ‘M’. Finally, all haplotypes that are found only in North America are designated by the prefix ‘N’. Representative sequences of 11 new haplotypes detected in this study and restricted to North America have been deposited in GenBank (accession numbers MG491978–MG491995). Table 1. Nearctic sample haplotype number with the total number of samples (N), total number of unique haplotypes, and the total number of each haplotype by state; states are classified according to Migratory Flyways [Pacific (P) or Central (C)], Bird Conservation Regions (CR1, 2, 3, 4 or 5) and Genetic Regions (GR1, 2, 3 or 4)   N  Total haplotypes  N1  N2  N3  N4  N5  N6  N7  N8  N9  N10  N11  N12  N13  N14  N15  N16  N17  N18  H1  Alaska (P, CR1, GR1)  24  8  5  2        1  1  1  1    1                12  Alberta (P, CR1, GR1)  2  2                      1                1  British Columbia (P, CR1, GR1)  4  3  2  1                                  1  Yukon (P, CR1, GR1)  2  2  1                                    1  California (P, CR2, GR2)  33  6    9                2  2    1    3        16  Oregon (P, CR2, GR2)  30  9    12      1      2    1  1  2    1      1    9  Idaho (P, CR2, GR2)  26  6    5            3            1  2  1      14  Wyoming (C, CR3, GR3)  16  3  7  1                                  8  Montana (C, CR3, GR3)  7  5  2  1  1                1                2  Nebraska (C, CR4, GR4)  5  4  1  1    1                              2  Oklahoma (C, CR4, GR4)  1  1  1                                      New Mexico (C, CR5, GR4)  4  3    1                1                  2  Colorado (C, CR5, GR4)  14  6  1  2  1                1              1  8  Arizona (C, CR5, GR4)  2  2    1                                  1  Total      20  36  2  1  1  1  1  6  1  4  7  2  1  2  5  1  1  1  77    N  Total haplotypes  N1  N2  N3  N4  N5  N6  N7  N8  N9  N10  N11  N12  N13  N14  N15  N16  N17  N18  H1  Alaska (P, CR1, GR1)  24  8  5  2        1  1  1  1    1                12  Alberta (P, CR1, GR1)  2  2                      1                1  British Columbia (P, CR1, GR1)  4  3  2  1                                  1  Yukon (P, CR1, GR1)  2  2  1                                    1  California (P, CR2, GR2)  33  6    9                2  2    1    3        16  Oregon (P, CR2, GR2)  30  9    12      1      2    1  1  2    1      1    9  Idaho (P, CR2, GR2)  26  6    5            3            1  2  1      14  Wyoming (C, CR3, GR3)  16  3  7  1                                  8  Montana (C, CR3, GR3)  7  5  2  1  1                1                2  Nebraska (C, CR4, GR4)  5  4  1  1    1                              2  Oklahoma (C, CR4, GR4)  1  1  1                                      New Mexico (C, CR5, GR4)  4  3    1                1                  2  Colorado (C, CR5, GR4)  14  6  1  2  1                1              1  8  Arizona (C, CR5, GR4)  2  2    1                                  1  Total      20  36  2  1  1  1  1  6  1  4  7  2  1  2  5  1  1  1  77  View Large Phylogeographical structure Using the computer program TCS we produced a haplotype network in which the 95% confidence limit to connection was eight steps (Fig. 1). Within this network the Mediterranean haplotypes appear to be a distinct lineage with two connections to the North American–Holarctic–Palearctic lineage. Although we detected 18 haplotypes restricted to North America, there is no resolution between the Nearctic and Palearctic haplotypes. Figure 1. View largeDownload slide Network diagram of all 44 control region haplotypes as determined by a minimum spanning network at the 95% confidence level using TCS v2.1. Colours represent lineages as determined in this study with yellow representing the single Holarctic haplotype, black representing the Mediterranean lineage, blue representing the Nearctic lineage, and orange representing the Palearctic lineage. The size of the shape is proportional to the number of individuals represented by that haplotype. Figure 1. View largeDownload slide Network diagram of all 44 control region haplotypes as determined by a minimum spanning network at the 95% confidence level using TCS v2.1. Colours represent lineages as determined in this study with yellow representing the single Holarctic haplotype, black representing the Mediterranean lineage, blue representing the Nearctic lineage, and orange representing the Palearctic lineage. The size of the shape is proportional to the number of individuals represented by that haplotype. Genetic diversity Haplotype and nucleotide diversities for the three lineages were: Palearctic (h = 0.8149, π = 0.01412), Nearctic (h = 0.7354, π = 0.0028) and Mediterranean (h = 0.5835, π = 0.0022) (Table 2). Allelic richness with the rarefaction correction for the Nearctic was 8.4, Palearctic was 7.0 and Mediterranean was 6.7 (Table 2). The Nearctic lineage was statistically significant for Tajima’s D, Fu’s F, raggedness and sum of squared deviations (SSD), while the Palearctic and Mediterranean lineage was not statistically significant for any neutrality measure or the mismatch analysis. Table 2. Statistical analyses of American golden eagles with samples grouped according to geographical lineages, Migratory Flyways, Bird Conservation Regions and Genetic Regions   N  h  π  Allelic richness  Tajima’s D  Fu’s Fs  SSD  Raggedness  Lineages  North American  170  0.81  0.003  8.4  −1.48  13.62  0.01  0.12  Palearctic  217  0.74  0.014  7  1.71  −1.41  0.08  0.7  Mediterranean  27  0.58  0.002  6.7  −0.9  −3.1  0.0001  0.54  Flyways  Central  99  0.73  0.004    −0.75  −3.65  0.202  0.66  Pacific  80  0.72  0.003    −0.95  −7.22  0.14  0.47  Conservation Regions  AK-AB-BC-YT  24  0.72  0.003  2.58  −1.33  −3.42  0.02  0.17  CA-ID-OR  89  0.75  0.002  2.37  −0.89  −5.67  0  0.09  MT-WY  40  0.71  0.002  1.94  −1.3  −2.77  0.02  0.15  NE-OK  6  0.87  0.003  3  −0.68  −0.99  0.01  0.13  AZ-CO-NM  20  0.68  0.002  2.29  −1.33  3.61  0.01  0.09  Genetic Regions  AK-AB-BC-YT  24  0.72  0.003  6.13  −1.33  3.42  0.02  0.17  CA-ID-OR  89  0.75  0.002  5.7  −0.89  5.67  0.01  0.09  MT-WY  40  0.71  0.002  4  −1.3  2.77  0.02  0.15  NE-OK-AZ-CO-NM  26  0.72  0.003  6.42  −1.22  3.88  0.004  0.09    N  h  π  Allelic richness  Tajima’s D  Fu’s Fs  SSD  Raggedness  Lineages  North American  170  0.81  0.003  8.4  −1.48  13.62  0.01  0.12  Palearctic  217  0.74  0.014  7  1.71  −1.41  0.08  0.7  Mediterranean  27  0.58  0.002  6.7  −0.9  −3.1  0.0001  0.54  Flyways  Central  99  0.73  0.004    −0.75  −3.65  0.202  0.66  Pacific  80  0.72  0.003    −0.95  −7.22  0.14  0.47  Conservation Regions  AK-AB-BC-YT  24  0.72  0.003  2.58  −1.33  −3.42  0.02  0.17  CA-ID-OR  89  0.75  0.002  2.37  −0.89  −5.67  0  0.09  MT-WY  40  0.71  0.002  1.94  −1.3  −2.77  0.02  0.15  NE-OK  6  0.87  0.003  3  −0.68  −0.99  0.01  0.13  AZ-CO-NM  20  0.68  0.002  2.29  −1.33  3.61  0.01  0.09  Genetic Regions  AK-AB-BC-YT  24  0.72  0.003  6.13  −1.33  3.42  0.02  0.17  CA-ID-OR  89  0.75  0.002  5.7  −0.89  5.67  0.01  0.09  MT-WY  40  0.71  0.002  4  −1.3  2.77  0.02  0.15  NE-OK-AZ-CO-NM  26  0.72  0.003  6.42  −1.22  3.88  0.004  0.09  Statistics include number of samples (N) haplotype diversity (h), nucleotide diversity (π), allelic richness, Tajima’s D and Fu’s F, SSD and raggedness. Results that are bold indicate significance at the 0.05 level. State and Canadian Province abbreviations are as follows: Alberta (AB), Alaska (AK), Arizona (AZ), British Columbia (BC), California (CA), Colorado (CO), Idaho (ID), Montana (MT), Nebraska (NE), New Mexico (NM), Oklahoma (OK), Oregon (OR), Wyoming (WY), Yukon (YT). View Large Within North America, haplotype diversity did not vary greatly when samples were partitioned by Migratory Flyways, BCRs or Genetic Regions with most groupings possessing haplotype diversities of approximately 0.74 (Table 2). Nucleotide diversity ranged from 0.0019 to 0.0039 depending on how the samples were partitioned (Table 2). Tajima’s D was not statistically significant for any of the three North American groupings while Fu’s F was statistically significant for all populations across groupings (Table 2). Genetic structure When the data were partitioned into Nearctic, Palearctic and Mediterranean, the AMOVA revealed a high level of genetic structuring with 47.62% of the genetic variation partitioned among these three units. Furthermore, the pairwise ΦST comparison revealed statistically significant genetic differentiation between all lineages, with the Palearctic/Nearctic differentiation being 0.517, the Nearctic/Mediterranean differentiation being 0.9139 and the Palearctic/Mediterranean differentiation being 0.3. Within North America, we evaluated the partitioning of genetic variation when individuals were partitioned in accordance with Migratory Flyways, BCRs and Genetic Regions. Figure 3 illustrates the distribution of haplotypes when partitioned by these alternative management scenarios. When samples were grouped by Migratory Flyways, there was statistically significant differentiation among populations within Migratory Flyways (Fsc = 0.0336) and within populations (Fst = 0.0624), but there was no statistically significant differentiation between the two Migratory Flyways (Table 3). When the hierarchical model was removed and samples were grouped only according to Migratory Flyways, the pairwise ΦST was statistically significant (Table 4). The hierarchical AMOVA for the BCRs did not detect statistically significant population structure among populations, among regions or within populations (Table 3), but the pairwise ΦST revealed significant genetic differentiation between golden eagles from Alaska, Alberta, British Columbia and Yukon when compared to individuals from California, Idaho and Oregon (Table 4). Additionally, individuals from California, Idaho and Oregon were significantly differentiated from individuals from Montana and Wyoming (ΦST = 0.12). Finally, low but statistically significant levels of genetic differentiation (ΦST = 0.044) was also detected between individual from Montana and Wyoming when compared to the grouping of individuals from Arizona, Colorado and New Mexico (Table 4). Combining the Nebraska and Oklahoma samples with samples from Arizona, Colorado and New Mexico to reflect the Genetic Region grouping, the hierarchical AMOVA did not reveal significant population differentiation at any level tested (Table 3), but similar to other analyses the pairwise ΦST detected statistically significant pairwise genetic differentiation between the grouping of individuals from Alaska, Alberta, British Columbia and Yukon relative to individuals from California, Idaho, and Oregon (Table 4). In congruence with the analysis based on BCRs, this analysis also detected statistically significant genetic differentiation between individuals from California, Idaho and Oregon when compared to individuals from Montana and Wyoming (Tables 3 and 4). Table 3. Results of hierarchical AMOVA when golden eagle samples are partitioned in accordance with Migratory Flyways, Bird Conservation Regions and Genetic Regions. F-statistics show the molecular variance partitioned among groups (Fct), among populations within groups (Fst), and within populations (Fsc)   Migratory Flyways  Bird Conservation Regions  Genetic Regions    d.f.  Sum of squares  Variance components  Percentage of variation  d.f.  Sum of squares  Variance components  Percentage of variation  d.f.  Sum of squares  Variance components  Percentage of variation  Among groups  1  2.01  0.17  2.99  3  3.55  0.03  4.32  3  3.63  0.023  3.99  Among populations  12  8.95  0.02  3.26  7  3.48  −0.01  −0.88  9  4.70  −0.004  −0.71  Within population  156  84.20  0.54  93.75  137  76.13  0.56  96.55  141  79.33  0.56  96.72  Fsc  0.0336        −0.00917        −0.00737        Fst  0.0624        0.0345        0.03278        Fct  0.0299        0.04325        0.0399          Migratory Flyways  Bird Conservation Regions  Genetic Regions    d.f.  Sum of squares  Variance components  Percentage of variation  d.f.  Sum of squares  Variance components  Percentage of variation  d.f.  Sum of squares  Variance components  Percentage of variation  Among groups  1  2.01  0.17  2.99  3  3.55  0.03  4.32  3  3.63  0.023  3.99  Among populations  12  8.95  0.02  3.26  7  3.48  −0.01  −0.88  9  4.70  −0.004  −0.71  Within population  156  84.20  0.54  93.75  137  76.13  0.56  96.55  141  79.33  0.56  96.72  Fsc  0.0336        −0.00917        −0.00737        Fst  0.0624        0.0345        0.03278        Fct  0.0299        0.04325        0.0399        Numbers identified by bold font represent statistically significant values (0.05). View Large Figure 3. View largeDownload slide Map of North American illustrating haplotype distributions when samples are partitioned according to Migratory Flyways, Bird Conservation Regions and Genetic Regions. Pie charts represent the haplotype frequency for the given Migratory Flyway, Bird Conservation Region and Genetic Region. The colours in each map represent the different proposed management units for each model. Figure 3. View largeDownload slide Map of North American illustrating haplotype distributions when samples are partitioned according to Migratory Flyways, Bird Conservation Regions and Genetic Regions. Pie charts represent the haplotype frequency for the given Migratory Flyway, Bird Conservation Region and Genetic Region. The colours in each map represent the different proposed management units for each model. Population fluctuations A unimodal mismatch distribution was detected for the Mediterranean, Palearctic and Nearctic lineages, indicating a more recent demographic expansion, although the only lineage that had a significant P-value for raggedness was the Nearctic (Fig. 2, Table 2). Within the Nearctic groupings, raggedness and SSD were both statistically significant when the samples were grouped into Migratory Flyways. When grouped into BCRs, raggedness was significant for the grouping of individuals from Alaska, Alberta, British Columbia and Yukon as well as the grouping of individuals from California, Idaho and Oregon. When samples were partitioned in accordance with the Genetic Regions, raggedness was only significant for the grouping of our most northern individuals from Alaska, Alberta, British Columbia and Yukon. SSD was not significant for either BCRs or Genetic Regions, adding support to the recent expansion hypothesis. Figure 2. View largeDownload slide Mismatch distribution for each of the lineages: Mediterranean (M), Palearctic (P) and Nearctic (N). The vertical bars represent the observed frequencies while the line represents the expected frequencies under a model of sudden expansion. Figure 2. View largeDownload slide Mismatch distribution for each of the lineages: Mediterranean (M), Palearctic (P) and Nearctic (N). The vertical bars represent the observed frequencies while the line represents the expected frequencies under a model of sudden expansion. DISCUSSION When the range of a species occupies the Holarctic, one of two phylogeographical patterns typically emerge. North American and Eurasian lineages occur in sympatry with a distinct North American haplogroup resulting from recent gene flow across Beringia or the Bering Strait, with examples such as the raven (Corvus corax) (Omland, Baker & Peters, 2006) and the peregrine falcon (Falco peregrinus) (Bell et al., 2014). Alternatively, due to a lack of gene flow across Beringia or the Bering Strait, speciation occurs and distinct species are found in both North America and Eurasia with examples such as white-tailed eagle and bald eagle (Hailer et al., 2006) or three-toed woodpecker (Picoides tridactylus) and American three-toed woodpecker (Picoides dorsalis) (Zink et al., 2000). Our results support the conclusions of Nebel et al. (2015) that there are two distinct mtDNA lineages: a Mediterranean lineage and a Holarctic lineage (Fig. 1). As also determined by Nebel et al. (2015), the Mediterranean lineage represents a monophyletic group, with individuals characterized by these haplotypes being found only in Mediterranean locations. The Holarctic lineage represents haplotypes found throughout the Palearctic and the Nearctic, with only one haplotype (H1) being truly Holarctic. In North America, H1 is the most numerous and most widespread haplotype, making up 45% of individuals sequenced, as opposed to the Palearctic where it represents few (<1%) individuals with a strong location bias in Europe. The widespread range of H1 together with its central location in the haplotype network suggests that this is an ancestral haplotype. The pattern of having a single shared haplotype between hemispheres was also found when using the mitochondrial control region in ravens (Omland et al., 2006). The other haplotypes defined as Holarctic by Nebel et al. (2015) are only found on a single continent, but our data do not recover reciprocal monophyly between the Palearctic and Nearctic haplotypes (Figs 1, 2). Despite not detecting reciprocal monophyly between the Nearctic and Palearctic haplotypes, pairwise genetic differentiation indicates little to no current gene flow between the Nearctic and Palearctic linages (ΦST = 0.517). Additionally, pairwise comparison of ΦST between the Nearctic and Mediterranean indicates almost complete lack of gene flow (ΦST = 0.9139). Among these three large geographical regions, there does appear to be some gene flow occurring between the Palearctic and Mediterranean lineages with a pairwise ΦST of 0.3. Our data also suggest that the Nearctic population of golden eagles has undergone a recent population expansion. When comparing the Nearctic to the Palearctic and Mediterranean lineages, both Tajima’s D and Fu’s F were statistically significant. Together with the unimodal mismatch analysis, which had significant values for raggedness, this supports a recent population expansion. The first mitochondrial analysis of North American golden eagles examined individuals representing the California Channel Islands and the California mainland (Sonsthagen et al., 2012a). Among the 42 individuals for which they were able to obtain mitochondrial sequence data, Sonsthagen et al. (2012a) detected five haplotypes. More recently, Craig et al. (2016) examined the same mitochondrial region of golden eagles from Alaska, California and Idaho. Among the 49 individuals they examined, they detected eight haplotypes, three unique to birds in their study and the five haplotypes previously detected in southern California golden eagles by Sonsthagen et al. (2012a). Additionally, Craig et al. (2016) reported the possibility of two region-specific haplotypes, one that appeared to be restricted to California and a second that appeared to be restricted to Idaho. To these previous mitochondrial studies of North American golden eagles, we add data from 115 additional individuals. Overall, we detected 19 haplotypes with 11 being newly detected variants. Our measures of haplotype diversity are higher for the Nearctic samples than detected by Nebel et al. (2015) for Old World samples, but our measures of nucleotide and haplotype diversity are similar to the results obtained by Craig et al. (2016). Regarding the conclusion of Craig et al. (2016), our data do not support the possibility of their haplotypes being region-specific because N14, which Craig et al. (2016) suggested was restricted to California, was also detected in Oregon in our analyses. Similarly, haplotype N10, which Craig et al. (2016) suggested was restricted to Idaho, was detected in samples from Oregon and New Mexico. Of the 11 new haplotypes that we identified, eight were region-specific. Haplotypes N5 and N17 were only found in Oregon, N6, N7 and N9 were only found in Alaska, N13 was only found in California, N16 was only found in Idaho and N18 was only found in Colorado (Table 1). Because these sequences were only single occurrences, we suspect that with more thorough sampling, they may be found to be in other areas. The only Holarctic haplotype (H1) was found in every sampling location in the Nearctic except Oklahoma, where only one individual was sampled (Table 1). Due to the wide geographical distribution of this haplotype and that locations bordering Oklahoma had this haplotype, we suspect that if sampling efforts were increased in future studies, this haplotype would be found in Oklahoma as well. Due to the uncertainty in how best to manage North American golden eagles (U.S. Fish and Wildlife, 2016), we partitioned our samples so that the proposed Migratory Flyways and BCRs could be tested as well as the results indicated in studies utilizing nuclear DNA (Genetic Regions). Although the hierarchical AMOVAs failed to reveal statistically significant levels of overall genetic differentiation for any of the potential groupings, the analysis of pairwise ΦST did reveal several patterns of significant genetic differentiation (Table 4). Most notably, for both the BCRs and the Genetic Regions, statistically significant pairwise genetic differentiation was detected between the grouping of northern golden eagles (Alaska and Canada) and golden eagles from California, Idaho and Oregon. Moreover, we detected significant genetic differentiation between the California, Idaho and Oregon birds and those from Montana and Wyoming. Both of these groupings are found in different Conservation Regions and Flyways. Table 4. Pairwise ΦST values for Migratory Flyways, Bird Conservation Regions and Genetic Regions   Migratory Flyways`    Bird Conservation Regions    Genetic Regions    West  East    AK-AB-BC-YT  CA-ID-OR  MO-WY  NE-OK  AZ-CO-NM    AK-AB-BC-YT  CA-ID-OR  MO-WY  NE-OK-AZ- CO-NM  Central  0    AK-AB-BC-YT  0          AK-AB-BC-YT  0        Pacific  0.018  0  CA-ID-OR  0.090  0        CA-ID-OR  0.090  0            MT-WY  −0.026  0.117  0      MT-WY  −0.0262  0.117  0          NE-OK  −0.062  0.068  −0.0532  0    NE-OK-AZ- CO-NM  −0.0002  0.015  0.0178  0        AZ-CO-NM  0.020  −0.003  0.044  0.0072  0              Migratory Flyways`    Bird Conservation Regions    Genetic Regions    West  East    AK-AB-BC-YT  CA-ID-OR  MO-WY  NE-OK  AZ-CO-NM    AK-AB-BC-YT  CA-ID-OR  MO-WY  NE-OK-AZ- CO-NM  Central  0    AK-AB-BC-YT  0          AK-AB-BC-YT  0        Pacific  0.018  0  CA-ID-OR  0.090  0        CA-ID-OR  0.090  0            MT-WY  −0.026  0.117  0      MT-WY  −0.0262  0.117  0          NE-OK  −0.062  0.068  −0.0532  0    NE-OK-AZ- CO-NM  −0.0002  0.015  0.0178  0        AZ-CO-NM  0.020  −0.003  0.044  0.0072  0            Numbers in bold are statistically significant at 0.05. View Large Based on the samples that we analysed, the primary difference between the BCRs and the Genetic Regions is that for the BCRs our samples representing Arizona, Colorado and New Mexico are considered to belong to a group genetically differentiated from our samples from Nebraska and Oklahoma, whereas for the Genetic Regions approach the samples from these five states are combined into a single unit. The additional partitioning of samples in the BCR detected significant genetic differentiation between samples from the Northern Plains (Montana and Wyoming) and the southern Rocky Mountain regions (Arizona, Colorado and New Mexico). These findings based on BCRs and Genetic Regions are similar to results reported by Doyle et al. (2016) in which the Alaska and California eagles were significantly different from the other western eagle population. This was also reported by Van Den Bussche et al. (2017), in which the Alaska and British Columbia eagles created a significant unit and the Idaho, California and Oregon eagles created a significant unit. Finally, when partitioning the Nearctic samples into Migratory Flyways, BCRs or Genetic Regions, the flyway subgroupings as well as the most northern golden eagles (Alaska and Canada) and the California, Idaho and Oregon golden eagles appear to have experienced a recent population expansion (Table 2). CONCLUSION Our study provides another example of intermediate divergence and the importance of studying the continuum of divergence for understanding geographical speciation, species limits and conservation priorities (Omland et al., 2006). This study also provides the first insight into the phylogeographical distribution and partitioning of mitochondrial genetic variation for North American golden eagles. To aid in determining more accurate historical events and to aid in determining proper management units, a more thorough sampling scheme should be carried out utilizing nuclear loci. Currently, single nucleotide polymorphisms have been used to begin to elucidate the genetic structure of North American golden eagles, but with limited numbers of samples from throughout the range (Doyle et al., 2016; Van Den Bussche et al., 2017). Future studies that couple mitochondrial and nuclear loci can provide greater insight into the demography of golden eagles by revealing patterns of genetic variation at both maternally and biparentally inherited loci. SUPPORTING INFORMATION Additional Supporting Information may be found in the online version of this article at the publisher’s web-site Table S1. Locality, age, provider and haplotype information for the 115 samples sequenced for the Nearctic portion of the study. ACKNOWLEDGMENTS We thank Bryan Bedrosian, Rob Domenech, Garth Herring, Gail Kratz, Brian Millsap, Gary Roemer, Victor Roubidoux and Brian Smith for contributing golden eagle blood samples from their own research and rehabilitation efforts. Without their kindness and contribution, we would not have had the geographical coverage represented in this study. We further thank the Iowa Tribe of Oklahoma for funding this study in collaboration with the Grey Snow Eagle House and the Ford Foundation for funding Megan Judkins’ research assistantship while the work for this project was performed. We extend our appreciation to the other members of Megan Judkins’ dissertation committee, Drs Meredith Hamilton, Michi Tobler, Andrew Doust and Jim Lish, whose review and feedback of an earlier version of the manuscript was critical. Finally, we thank the four anonymous reviewers of a previous version of the manuscript for their helpful comments. REFERENCES Bayard de Volo S. 2008. Genetic studies of Northern Goshawks (Accipiter Ggentilis): Genetic tagging and individual identification from feathers, and determining phylogeography, gene flow and population history for goshawks in North America . Unpublished D. Phil. Dissertation, Colorado State University. Bell DA, Griffiths CS, Caballero IC, Hartley RR, Lawson RH. 2014. Genetic evidence for global dispersal in the peregrine falcon (Falco peregrinus) and affinity with the Taita falcon (Falco fasciinucha). Journal of Raptor Research  48: 44– 53. Google Scholar CrossRef Search ADS   Bielikova M, Ficek A, Valkova D, Turna J. 2010. Multiplex PCR amplification of 13 microsatellite loci for Aquila chrysaetos in forensic applications. Biologia  65: 1081– 1088. Google Scholar CrossRef Search ADS   Bourke BP, Frantz AC, Lavers CP, Davison A, Dawson DA, Burke TA. 2010. Genetic signatures of population change in the British golden eagle (Aquila chrysaetos). Conservation Genetics  11: 1837– 1846. Google Scholar CrossRef Search ADS   Clement M, Posada D, Crandall KA. 2000. TCS: A computer program to estimate gene genealogies. Molecular Ecology  9: 1657– 1659. Google Scholar CrossRef Search ADS PubMed  Craig EH, Adams JR, Waits LP, Fuller MR, Whittington DM. 2016. Nuclear and mitochondrial DNA analyses of golden eagles (Aquila chrysaetos canadensis) from three areas in western North America; initial results and conservation implications. PLoS ONE  11: 1– 15. Google Scholar CrossRef Search ADS   Doyle JM, Katzner TE, Roemer GW, CainIII JW, Millsap BA, McIntyre CL, Sonsthagen SA, Fernandez NB, Wheeler M, Bulut Z, Bloom PH, DeWoody A. 2016. Genetic structure and viability selection in the golden eagle (Aquila chrysaetos): a vagile raptor with a Holarctic distribution. Conservation Genetics  17: 1– 16. Google Scholar CrossRef Search ADS   Excoffier L, Lischer HEL. 2010. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources  10: 564– 567. Google Scholar CrossRef Search ADS PubMed  Fielding AH, Whitfield DP, McLeod DRA. 2006. Spatial association as an indicator of the potential for future interactions between wind energy developments and golden eagles Aquila chrysaetos in Scotland. Biological Conservation  131: 359– 369. Google Scholar CrossRef Search ADS   Godoy JA, Negro JJ, Hiraldo F, Donazar JA. 2004. Phylogeography, genetic structure and diversity in the endangered bearded vulture (Gypaetus barbatus, L.) as revealed by mitochondrial DNA. Molecular Ecology  13: 371– 390. Google Scholar CrossRef Search ADS PubMed  Hailer F, Helander B, Folkestad AO, Ganusevich SA, Garstad S, Hauff P, Koren C, Nygard T, Volke V, Vila C, Ellegren H. 2006. Bottlenecked but long-lived: high genetic diversity retained in white-tailed eagles upon recovery from population decline. Biology Letters  2: 316– 319. Google Scholar CrossRef Search ADS PubMed  Hewitt G. 2000. The genetic legacy of the Quaternary ice ages. Nature  405: 907– 913. Google Scholar CrossRef Search ADS PubMed  Hewitt GM. 2004. Genetic consequences of climatic oscillations in the Quaternary. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences  359: 183– 195. Google Scholar CrossRef Search ADS PubMed  Hull JM, Girman DJ. 2005. Effects of Holocene climate change on the historical demography of migrating sharp-shinned hawks (Accipiter striatus velox) in North America. Molecular Ecology  14: 159– 170. Google Scholar CrossRef Search ADS PubMed  Longmire JL, Maltbie M, Baker RJ. 1997. Use of “lysis buffer” in DNA isolation and its implication for museum collections. Occasional Papers, The Museum, Texas Tech University  163: 1– 3. Maddison DR, Maddison WP. 2005. MacClade 4: Analysis of phylogeny and character evolution . Sunderland, MA: Sinauer Associates, Inc., Publishers. Masuda R, Noro M, Kurose N, Nishida-Umehara C, Takechi H, Yamazaki T, Kosuge M, Yoshida MC. 1998. Genetic characteristics of endangered Japanese golden eagles (Aquila chrysaetos japonica) based on mitochondrial DNA D-loop sequences and karyotypes. Zoo Biology  17: 111– 121. Google Scholar CrossRef Search ADS   Millsap BA, Harmata AR, Stahlecker DW, Mikesic DG. 2014. Natal dispersal distance of bald and golden eagles originating in the coterminous United States as inferred from band encounters. Journal of Raptor Research  48: 13– 23. Google Scholar CrossRef Search ADS   Morneau F, Tremblay JA, Todd C, Chubbs TE, Maisonneuve C, Lemaitre J, Katzner T. 2015. Known breeding distribution and abundance of golden eagles in Eastern North America. Northeastern Naturalist  22: 236– 247. Google Scholar CrossRef Search ADS   Nebel C, Gamauf A, Haring E, Segelbacher G, Villers A, Zachos FE. 2015. Mitochondrial DNA analysis reveals Holarctic homogeneity and a distinct Mediterranean lineage in the golden eagle (Aquila chrysaetos). Biological Journal of the Linnean Society  116: 328– 340. Google Scholar CrossRef Search ADS   Omland KE, Baker JM, Peters JL. 2006. Genetic signatures of intermediate divergence: Population history of Old and New World Holarctic ravens (Corvus corax). Molecular Ecology  15: 795– 808. Google Scholar CrossRef Search ADS PubMed  Pagel JE, Kritz KJ, Millsap BA, Robert K. 2013. Bald eagle and golden eagle mortalities at wind energy facilities in the contiguous United States. Journal of Raptor Research  47: 311– 315. Google Scholar CrossRef Search ADS   Petit RJ, El Mousadik A, Pons O. 1998. Identifying populations for conservation on the basis of genetic markers. Conservation Biology  12: 844– 855. Google Scholar CrossRef Search ADS   Poulakakis N, Antoniou A, Mantziou G, Parmakelis A, Skartsi T, Vasilakis D, Elorriaga J, De La Puente J, Gavashelishvili A, Ghasabyan M, Katzner T, McGrady M, Batbayar N, Fuller M, Natsagdorj T. 2008. Population structure, diversity, and phylogeography in the near-threatened Eurasian black vultures Aegypius monachus (Falconiformes; Accipitridae) in Europe: Insights from microsatellite and mitochondrial DNA variation. Biological Journal of the Linnean Society  95: 859– 872. Google Scholar CrossRef Search ADS   Sonsthagen SA, Coonan TJ, Latta BC, Sage GK, Talbot S. 2012a. Genetic diversity of a newly established population of golden eagles on the Channel Islands, California. Biological Conservation  146: 116– 122. Google Scholar CrossRef Search ADS   Sonsthagen SA, Rosenfield RN, Bielefeldt J, Murphy RK, Steward AC, Stout WE, Driscoll TG, Bozek MA, Sloss BL, Talbot SL. 2012b. Genetic and morphological divergence among Cooper’s Hawk (Accipiter cooperii) populations breeding in North-Central and Western North America. The Auk  129: 427– 437. Google Scholar CrossRef Search ADS   Stauber E, Finch N, Talcott PA, Gay JM. 2010. Lead poisoning of bald (Haliaeetus leucocephalus) and golden (Aquila chrysaetos) eagles in the U.S. inland Pacific Northwest region – an 18-year retrospective study: 1991–2008. Journal of Avian Medicine and Surgery  24: 279– 287. Google Scholar CrossRef Search ADS PubMed  Suchentrunk F, Haller H, Ratti P. 1999. Gene pool variability of a golden eagle (Aquila chrysaetos) population from the Swiss Alps. Biological Conservation  90: 151– 155. Google Scholar CrossRef Search ADS   U.S. Fish and Wildlife Service. 2016. Bald and Golden Eagles: Population demographics and estimation of sustainable take in the United States, 2016 update . Washington: Division of Migratory Bird Management. Van Den Bussche RA, Judkins ME, Montague MJ, Warren WC. 2017. A resource of genome-wide single nucleotide polymorphisms (SNPs) for the conservation and management of golden eagles (Aquila chrysaetos). Journal of Raptor Research  51: 368– 377. Google Scholar CrossRef Search ADS   Watson J, Brockie K, Watson D. 2010. The golden eagle, 2nd edn . New Haven: Yale University Press. Wink M, Clouet M, Goar LJ, Barrau C. 2004. Sequence variation in the cytochrome b gene of subspecies of Golden Eagles Aquila chrysaetos. Alauda  72: 153– 157. Zink RM, Barrowclough GF, Atwood JL, Blackwell-Rago RC. 2000. Genetics, taxonomy, and conservation of the threatened California Gnatcatcher. Conservation Biology  14: 1394– 1405. Google Scholar CrossRef Search ADS   © 2017 The Linnean Society of London, Biological Journal of the Linnean Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Biological Journal of the Linnean Society Oxford University Press

Holarctic phylogeography of golden eagles (Aquila chrysaetos) and evaluation of alternative North American management approaches

Loading next page...
 
/lp/ou_press/holarctic-phylogeography-of-golden-eagles-aquila-chrysaetos-and-4EDnbTg0wV
Publisher
The Linnean Society of London
Copyright
© 2017 The Linnean Society of London, Biological Journal of the Linnean Society
ISSN
0024-4066
eISSN
1095-8312
D.O.I.
10.1093/biolinnean/blx138
Publisher site
See Article on Publisher Site

Abstract

Abstract The golden eagle (Aquila chrysaetos) is a long-lived bird of prey with a Holarctic distribution. This species has survived severe anthropogenic stressors that have reduced or eliminated populations in some parts of its range. Despite the ecological and cultural importance of the golden eagle, few attempts have been made to determine the partitioning of genetic variation over large areas of its range. This study generated DNA sequence data from the mitochondrial DNA control region from 115 North American golden eagles and combined these data with existing control region sequences from over 300 Nearctic, Palearctic and Mediterranean golden eagles to provide a clearer holistic picture of the Holarctic phylogeographical patterns of genetic variation in this species and the genetic variation and demographic history of golden eagles in North America. The results support that there are two genetic lineages of golden eagles, one representing the Mediterranean and the other occurring throughout the Holarctic. The Holarctic lineage reveals little to no modern-day gene flow between Nearctic and Palearctic golden eagles. Furthermore, the current distribution of haplotypes in the Nearctic shows a recent population expansion with moderate levels of gene flow. INTRODUCTION Climate perturbations throughout history have had effects on the evolutionary history and phylogeography of many species. The climate fluctuations of the Quaternary (~1.8 Mya to present) have been well studied and shown to have impacted the current distributions and speciation events of many taxa in the Northern Hemisphere (Hewitt, 2000, 2004). Raptors throughout the Holarctic display varying patterns of genetic differentiation with much of this variation being attributed to how they were impacted during the Pleistocene. When considering the modern-day phylogeography of raptors, an east to west differentiation of haplotypes is commonly observed. Cooper’s hawk (Accipiter cooperi; Sonsthagen et al., 2012b) and sharp-shinned hawk (A. striatus; Hull & Girman, 2005) in North America, as well as white-tailed eagle (Haliaeetus albicilla; Hailer et al., 2006), bearded vulture (Gypaetus barbatus; Godoy et al., 2004) and Eurasian black vulture (Aegypius monachus; Poulakakis et al., 2008) in the Palearctic exhibit this trend. This contrasts with northern goshawk (A. gentilis) in North America in which an east–west differentiation was not observed. Instead, the south-western populations differed significantly from the other sampled locations (Bayard de Volo, 2008). When the range of a species stretches across the Holarctic, typically two patterns of genetic structuring between the Palearctic and Nearctic exist. In the first pattern, American and Eurasian lineages live in sympatry with a distinct American haplogroup resulting from recent gene flow across Beringia or the Bering Strait (Nebel et al., 2015). Alternatively, due to a lack of gene flow across Beringia or the Bering Strait, speciation occurs and distinct species are found in both North America and Eurasia (Nebel et al., 2015). Golden eagle (Aquila chrysaetos) is the only eagle species found throughout the Holarctic. Golden eagles are able to survive in diverse environments by adapting to the available prey, landscape and their choice of nesting sites. For example, golden eagles have adapted to nesting on the ground in tree-less steppe and desert habitat as opposed to cliffs and trees in other parts of the world, and their diet includes diverse prey such as ptarmigan, hares, tortoise, snakes and lizards (Watson, Brockie & Watson, 2010). Furthermore, golden eagles are highly mobile at different stages of their lives or during different times of the year. Juvenile eagles are nomadic for approximately 5 years until they become sexually mature, at which time it has been shown, through banding studies, that North American individuals travel back to within approximately 46–175 km of their natal location to establish a breeding territory (Millsap et al., 2014). North American golden eagles that breed in northern latitudes are obligatory migrants, as compared to their southern counterparts which do not migrate, and will move to more southern latitudes in autumn and return to their nesting areas in spring. Although they are thought to be highly philopatric to their natal site (Millsap et al., 2014), these movements may allow for gene flow throughout the range of the species through behaviours such as extra-pair copulations or dispersal events. Little is known about the phylogeographical patterns of genetic variation in this species. Nebel et al. (2015) provided insight into the phylogeographical structure of Palearctic golden eagles through examination of 402 bp of the mitochondrial DNA (mtDNA) control region, but the North American subspecies was poorly represented due to a lack of available samples. Nebel et al. (2015) identified 26 haplotypes, only one of which was restricted to North America, and they speculated that there were two historical refugia. Sonsthagen et al. (2012a) and Craig et al. (2016) also utilized this region of the mtDNA as well as microsatellite loci to analyse golden eagles in the far western portion of their North American range. Sonsthagen et al. (2012a) identified five novel haplotypes and high gene flow when they analysed golden eagles in California and the Channel Islands, while Craig et al. (2016) identified the five hapolotypes previously detected by Sonsthagen et al. (2012a) and three novel haplotypes when analysing samples from Alaska, Idaho and California. Other genetic studies on golden eagles have utilized microsatellites, allozymes and mtDNA, but these studies focused on very limited portions of the species range (Masuda et al., 1998; Suchentrunk, Haller & Ratti, 1999; Wink et al., 2004; Bielikova et al., 2010; Bourke et al., 2010). Throughout their range, golden eagles have faced many anthropogenic stressors including illegal shooting, windfarms, lead poisoning, habitat loss, organic pollutants and game keepers driving the population into bottlenecks or causing local extinctions (Fielding, Whitfield & McLeod,, 2006; Bourke et al., 2010; Stauber et al., 2010; Watson et al., 2010; Pagel et al., 2013). For example, breeding populations of golden eagles have been extirpated east of the Mississippi River in the United States (Morneau et al., 2015), in Ireland (Bourke et al., 2010), and in the Alpine foothills and lowlands of Germany, Austria, the Czech Republic and Poland (Nebel et al., 2015). These bottlenecks and the population fragmentation caused by regional and local extirpations can allow for inbreeding depression and loss of genetic variation, adding to the need of determining the current partitioning of genetic variation for management purposes. In the United States there are two proposed models for developing Eagle Management Units (U.S. Fish and Wildlife Service, 2016). The first management proposal includes the migratory flyways found in North America (Atlantic, Mississippi, Central and Pacific) with the Mississippi and Atlantic flyways being combined. The alternative proposal is based on the currently recognized Bird Conservation Regions (BCRs) with support coming from dispersal from the nesting site (U.S. Fish and Wildlife Service, 2016). Recent studies utilizing nuclear markers (Craig et al., 2016; Doyle et al., 2016; Van Den Bussche et al., 2017) reveal similar groupings of golden eagles across studies, with some support leaning towards the BCR model. Generally, these studies support the genetic uniqueness of golden eagles in northern California, southern Idaho and southern Oregon, and another group consisting of golden eagles from western Canada and Alaska. They differ from the BCR model by grouping most of the central to eastern portion of the western golden eagle range into a single unit (Doyle et al., 2016; Van Den Bussche et al., 2017). Unfortunately, due to a lack of sufficient sampling across the North American range of golden eagles it is unknown if this pattern is an artefact of sampling or if these individuals should be grouped into a single management unit. This study utilized a fragment of the mtDNA control region previously used by Nebel et al. (2015), Sonsthagen et al. (2012a) and Craig et al. (2016) to generate DNA sequence data from 115 North American golden eagles representing Alaska, Alberta, Arizona, California, Colorado, Idaho, Montana, Nebraska, New Mexico, Oklahoma, Oregon, Wyoming, Yukon and British Colombia. Our goals were to fill in geographical gaps across the range of North American golden eagles to provide a clearer picture of the Holarctic phylogeographical patterns of genetic variation of golden eagles. Additionally, we were interested in evaluating whether the mitochondrial sequence data provided insight into the most appropriate geographical structure for managing golden eagles by evaluating levels of genetic differentiation when the samples were partitioned by Migratory Flyways, BCRs or Genetic Regions. Finally, we were interested in elucidating the demographic history of golden eagles in North America to answer questions regarding the impacts of population bottlenecks and subsequent expansions on North America golden eagles. This additional information will allow for a more complete understanding of the historical movements and current patterns of mitochondrial genetic variation for golden eagles throughout the Holarctic. MATERIALS AND METHODS DNA isolation and sequence development Working with permited U.S. Fish and Wildlife biologists and wildlife rehabilitators, we obtained blood samples from 124 North American golden eagles. When nests contained more than one individual, only one sample was included in the study to eliminate over-representation of a haplotype. Individuals that were not collected as a hatchling had natal location determined by GPS transmitters or isotope data. Only two samples (one individual from Colorado and a second individual from Wyoming) may not have a true natal origin represented as they were collected during December and natal location data were unavailable. Blood samples (~ 0.5 mL) were collected, stored in lysis buffer (Longmire, Maltbie & Baker, 1997) and subsequently shipped to our laboratory at Oklahoma State University where DNA was extracted following the protocol of Longmire et al. (1997). DNA quality was assessed by running an aliquot on a 1% agarose gel and quantified using a NanoDrop 3300 spectrophotometer (Thermo Scientific). A 442-bp region of domain I and II of the hypervariable control region-1 (D-loop) was amplified using polymerase chain reaction (PCR) and primers GOEA_CR1L and GOEA_CR595H developed by Sonsthagen et al. (2012b). PCRs contained 6 µL of GoTaq Flexi Buffer (Promega), 25 mM MgCl2, 1 µM of each primer, 10 mg/mL bovine serum albumin, 4 mM dNTPs, 1 unit of GoTaq (Promega) and double distilled (dd)H2O to make a 30-µL reaction. The thermoprofile was 95 °C for 5 min, 35 cycles of 95 °C for 1 min, 40 °C for 1 min, 72 °C for 1 min and a final extension at 72 °C for 7 min. PCR products were visualized by running on a 1% agarose gel and subsequently cleaned using the Promega Wizard Kit. Cleaned products were sequenced using BigDye v3.1 chain terminators (Thermo Fisher), 5× sequencing buffer (Edge Bio), GOEA_CR1L primer and ddH2O. The sequencing thermoprofile consisted of 35 cycles of 96 °C for 30 s, 50 °C for 35 s and 60 °C for 4 min. Sequencing products were cleaned using either Performa DTR V3 96 well short plates or Performa DTR Gel Filtration cartridges (Edge Bio). The cleaned sequencing product (~20 µL) was added to 1.0 µL of HiDi (Life Technologies) followed by electrophoresis using an ABI 3130 Genetic Analyzer (Applied Biosystems). Base calling was performed using Sequencing Analysis 5.2 (Applied Biosystems). Sequences were visually inspected and those sequences that were not long enough or had ambiguous bases were resequenced using the alternative primer (GOEA_CR595H) and de novo aligned in Geneious 7.1.9 to acquire the full 442-bp region. In addition to the sequences we generated, we also downloaded from GenBank the Old World golden eagle sequences generated by Nebel et al. (2015; accession numbers KR259251–KR259276), the southern California sequences generated by Sonsthagen et al. (2012a; accession numbers JQ246417–JQ246421), and the sequences from California, Idaho and Oregon generated by Craig et al. (2016; accession numbers KX687705–KX687711). All sequences were aligned using the ClustalX multiple sequence alignment in Geneious 7.1.9. The resulting alignment was imported into MacClade version 4.06 (Maddison and Maddison, 2005) for visual inspection and to trim our generated sequences as well as those of Sonsthagen et al. (2012a) and Craig et al. (2016) to 402 bp so that we could directly compare these data with those of Nebel et al. (2015). Using the Redundant Taxa option in MacClade we grouped sequences into haplotypes and assigned names to the haplotypes. Further analyses were conducted utilizing the data of Craig et al. (2016) and Nebel et al. (2015) because the Sonsthagen et al. (2012a) haplotype frequencies and locations were not available. Furthermore, approximately half of the golden eagles represented in Sonsthagen et al. (2012a) were from the Channel Islands, a location that only acquired golden eagles in the mid-1990s when bald eagles (Haliaeetus leucocephalus) were extirpated from the island, meaning the natal locations of the eagles that moved into the area are unknown. Data analysis Relationships among control region haplotypes of golden eagles were evaluated through the construction of a minimum spanning network at the 95% confidence level using TCS v2.1 (Clement, Posada & Crandall, 2000). As described below, analyses were conducted on: (1) the entire data set including the samples from our dataset, Craig et al. (2016) and Nebel et al. (2015); (2) data partitioned by geography into Nearctic, Palearctic and Mediterranean; and (3) with the Nearctic samples further broken down by the proposed Migratory Flyways (U.S. Fish and Wildlife Service, 2016), BCRs (U.S. Fish and Wildlife Service, 2016) or Genetic Regions suggested by nuclear markers in previous studies (Craig et al., 2016; Doyle et al., 2016; Van Den Bussche et al., 2017). To assess genetic diversity for each of the three data partitions, we calculated haplotype and nucleotide diversity in Arlequin 3.1 (Excoffier & Lischer 2010). We also calculated the allelic richness using Contrib 1.4 (Petit, El Mousadik & Pons, 1998) in which the population sample size was normalized using a rarefaction correction. The rarefaction value was set to the smallest population size. We performed neutrality tests, Tajima’s D and Fu’s Fs, in Arlequin 3.1 to confirm the selective neutrality for the region. Tajima’s D and Fu’s Fs were further used to determine demographic histories. We used an analysis of molecular variance (AMOVA) and pairwise ΦSTusing the Kimuira two-parameter model (K80) of evolution in Arlequin 3.1 to assess the partitioning of genetic variation within and among the Nearctic, Palearctic and Mediterranean sample localities. For the North American samples, we utilized a hierarchical AMOVA to test the three proposed conservation models (Migratory Flyways, BCRs, Genetic Regions). The hierarchical approach grouped samples by State or Canadian Province and then the States or Provinces were grouped according the Migratory Flyways, BCRs or Genetic Regions. RESULTS We were able to generate DNA sequence data from 115 North American golden eagles (S1) which combined with the samples from Nebel et al. (2015) and Craig et al. (2016) provided us with a total of 417 samples. The 402-bp region contained 27 polymorphic sites and resulted in 44 haplotypes. The 115 Nearctic samples that we generated DNA sequences from provided 19 haplotypes (Table 1). Of these 19 haplotypes, 11 were novel to this study, four previously described by Sonsthagen et al. (2012a) and Craig et al. (2016), three previously described by Craig et al. (2016), and one previously described by Sonsthagen et al. (2012a), Nebel et al. (2015) and Craig et al. (2016). Although Nebel et al. (2015) classified their haplotypes as Holarctic and Mediterranean, based on our increased sampling of North American golden eagles, only haplotype H5 is Holarctic. To provide greater clarity for understanding golden eagle phylogeographical patterns, we used the following haplotypic designations. The single Holarctic haplotype (H5 of Nebel et al., 2015) is represented as ‘H1’, and all other Holarctic haplotypes of Nebel et al. (2015) are Palearctic and identified in this study by the letter ‘P’. We retained the Mediterranean haplotypic designation of Nebel et al. (2015) and assigned it with the prefix ‘M’. Finally, all haplotypes that are found only in North America are designated by the prefix ‘N’. Representative sequences of 11 new haplotypes detected in this study and restricted to North America have been deposited in GenBank (accession numbers MG491978–MG491995). Table 1. Nearctic sample haplotype number with the total number of samples (N), total number of unique haplotypes, and the total number of each haplotype by state; states are classified according to Migratory Flyways [Pacific (P) or Central (C)], Bird Conservation Regions (CR1, 2, 3, 4 or 5) and Genetic Regions (GR1, 2, 3 or 4)   N  Total haplotypes  N1  N2  N3  N4  N5  N6  N7  N8  N9  N10  N11  N12  N13  N14  N15  N16  N17  N18  H1  Alaska (P, CR1, GR1)  24  8  5  2        1  1  1  1    1                12  Alberta (P, CR1, GR1)  2  2                      1                1  British Columbia (P, CR1, GR1)  4  3  2  1                                  1  Yukon (P, CR1, GR1)  2  2  1                                    1  California (P, CR2, GR2)  33  6    9                2  2    1    3        16  Oregon (P, CR2, GR2)  30  9    12      1      2    1  1  2    1      1    9  Idaho (P, CR2, GR2)  26  6    5            3            1  2  1      14  Wyoming (C, CR3, GR3)  16  3  7  1                                  8  Montana (C, CR3, GR3)  7  5  2  1  1                1                2  Nebraska (C, CR4, GR4)  5  4  1  1    1                              2  Oklahoma (C, CR4, GR4)  1  1  1                                      New Mexico (C, CR5, GR4)  4  3    1                1                  2  Colorado (C, CR5, GR4)  14  6  1  2  1                1              1  8  Arizona (C, CR5, GR4)  2  2    1                                  1  Total      20  36  2  1  1  1  1  6  1  4  7  2  1  2  5  1  1  1  77    N  Total haplotypes  N1  N2  N3  N4  N5  N6  N7  N8  N9  N10  N11  N12  N13  N14  N15  N16  N17  N18  H1  Alaska (P, CR1, GR1)  24  8  5  2        1  1  1  1    1                12  Alberta (P, CR1, GR1)  2  2                      1                1  British Columbia (P, CR1, GR1)  4  3  2  1                                  1  Yukon (P, CR1, GR1)  2  2  1                                    1  California (P, CR2, GR2)  33  6    9                2  2    1    3        16  Oregon (P, CR2, GR2)  30  9    12      1      2    1  1  2    1      1    9  Idaho (P, CR2, GR2)  26  6    5            3            1  2  1      14  Wyoming (C, CR3, GR3)  16  3  7  1                                  8  Montana (C, CR3, GR3)  7  5  2  1  1                1                2  Nebraska (C, CR4, GR4)  5  4  1  1    1                              2  Oklahoma (C, CR4, GR4)  1  1  1                                      New Mexico (C, CR5, GR4)  4  3    1                1                  2  Colorado (C, CR5, GR4)  14  6  1  2  1                1              1  8  Arizona (C, CR5, GR4)  2  2    1                                  1  Total      20  36  2  1  1  1  1  6  1  4  7  2  1  2  5  1  1  1  77  View Large Phylogeographical structure Using the computer program TCS we produced a haplotype network in which the 95% confidence limit to connection was eight steps (Fig. 1). Within this network the Mediterranean haplotypes appear to be a distinct lineage with two connections to the North American–Holarctic–Palearctic lineage. Although we detected 18 haplotypes restricted to North America, there is no resolution between the Nearctic and Palearctic haplotypes. Figure 1. View largeDownload slide Network diagram of all 44 control region haplotypes as determined by a minimum spanning network at the 95% confidence level using TCS v2.1. Colours represent lineages as determined in this study with yellow representing the single Holarctic haplotype, black representing the Mediterranean lineage, blue representing the Nearctic lineage, and orange representing the Palearctic lineage. The size of the shape is proportional to the number of individuals represented by that haplotype. Figure 1. View largeDownload slide Network diagram of all 44 control region haplotypes as determined by a minimum spanning network at the 95% confidence level using TCS v2.1. Colours represent lineages as determined in this study with yellow representing the single Holarctic haplotype, black representing the Mediterranean lineage, blue representing the Nearctic lineage, and orange representing the Palearctic lineage. The size of the shape is proportional to the number of individuals represented by that haplotype. Genetic diversity Haplotype and nucleotide diversities for the three lineages were: Palearctic (h = 0.8149, π = 0.01412), Nearctic (h = 0.7354, π = 0.0028) and Mediterranean (h = 0.5835, π = 0.0022) (Table 2). Allelic richness with the rarefaction correction for the Nearctic was 8.4, Palearctic was 7.0 and Mediterranean was 6.7 (Table 2). The Nearctic lineage was statistically significant for Tajima’s D, Fu’s F, raggedness and sum of squared deviations (SSD), while the Palearctic and Mediterranean lineage was not statistically significant for any neutrality measure or the mismatch analysis. Table 2. Statistical analyses of American golden eagles with samples grouped according to geographical lineages, Migratory Flyways, Bird Conservation Regions and Genetic Regions   N  h  π  Allelic richness  Tajima’s D  Fu’s Fs  SSD  Raggedness  Lineages  North American  170  0.81  0.003  8.4  −1.48  13.62  0.01  0.12  Palearctic  217  0.74  0.014  7  1.71  −1.41  0.08  0.7  Mediterranean  27  0.58  0.002  6.7  −0.9  −3.1  0.0001  0.54  Flyways  Central  99  0.73  0.004    −0.75  −3.65  0.202  0.66  Pacific  80  0.72  0.003    −0.95  −7.22  0.14  0.47  Conservation Regions  AK-AB-BC-YT  24  0.72  0.003  2.58  −1.33  −3.42  0.02  0.17  CA-ID-OR  89  0.75  0.002  2.37  −0.89  −5.67  0  0.09  MT-WY  40  0.71  0.002  1.94  −1.3  −2.77  0.02  0.15  NE-OK  6  0.87  0.003  3  −0.68  −0.99  0.01  0.13  AZ-CO-NM  20  0.68  0.002  2.29  −1.33  3.61  0.01  0.09  Genetic Regions  AK-AB-BC-YT  24  0.72  0.003  6.13  −1.33  3.42  0.02  0.17  CA-ID-OR  89  0.75  0.002  5.7  −0.89  5.67  0.01  0.09  MT-WY  40  0.71  0.002  4  −1.3  2.77  0.02  0.15  NE-OK-AZ-CO-NM  26  0.72  0.003  6.42  −1.22  3.88  0.004  0.09    N  h  π  Allelic richness  Tajima’s D  Fu’s Fs  SSD  Raggedness  Lineages  North American  170  0.81  0.003  8.4  −1.48  13.62  0.01  0.12  Palearctic  217  0.74  0.014  7  1.71  −1.41  0.08  0.7  Mediterranean  27  0.58  0.002  6.7  −0.9  −3.1  0.0001  0.54  Flyways  Central  99  0.73  0.004    −0.75  −3.65  0.202  0.66  Pacific  80  0.72  0.003    −0.95  −7.22  0.14  0.47  Conservation Regions  AK-AB-BC-YT  24  0.72  0.003  2.58  −1.33  −3.42  0.02  0.17  CA-ID-OR  89  0.75  0.002  2.37  −0.89  −5.67  0  0.09  MT-WY  40  0.71  0.002  1.94  −1.3  −2.77  0.02  0.15  NE-OK  6  0.87  0.003  3  −0.68  −0.99  0.01  0.13  AZ-CO-NM  20  0.68  0.002  2.29  −1.33  3.61  0.01  0.09  Genetic Regions  AK-AB-BC-YT  24  0.72  0.003  6.13  −1.33  3.42  0.02  0.17  CA-ID-OR  89  0.75  0.002  5.7  −0.89  5.67  0.01  0.09  MT-WY  40  0.71  0.002  4  −1.3  2.77  0.02  0.15  NE-OK-AZ-CO-NM  26  0.72  0.003  6.42  −1.22  3.88  0.004  0.09  Statistics include number of samples (N) haplotype diversity (h), nucleotide diversity (π), allelic richness, Tajima’s D and Fu’s F, SSD and raggedness. Results that are bold indicate significance at the 0.05 level. State and Canadian Province abbreviations are as follows: Alberta (AB), Alaska (AK), Arizona (AZ), British Columbia (BC), California (CA), Colorado (CO), Idaho (ID), Montana (MT), Nebraska (NE), New Mexico (NM), Oklahoma (OK), Oregon (OR), Wyoming (WY), Yukon (YT). View Large Within North America, haplotype diversity did not vary greatly when samples were partitioned by Migratory Flyways, BCRs or Genetic Regions with most groupings possessing haplotype diversities of approximately 0.74 (Table 2). Nucleotide diversity ranged from 0.0019 to 0.0039 depending on how the samples were partitioned (Table 2). Tajima’s D was not statistically significant for any of the three North American groupings while Fu’s F was statistically significant for all populations across groupings (Table 2). Genetic structure When the data were partitioned into Nearctic, Palearctic and Mediterranean, the AMOVA revealed a high level of genetic structuring with 47.62% of the genetic variation partitioned among these three units. Furthermore, the pairwise ΦST comparison revealed statistically significant genetic differentiation between all lineages, with the Palearctic/Nearctic differentiation being 0.517, the Nearctic/Mediterranean differentiation being 0.9139 and the Palearctic/Mediterranean differentiation being 0.3. Within North America, we evaluated the partitioning of genetic variation when individuals were partitioned in accordance with Migratory Flyways, BCRs and Genetic Regions. Figure 3 illustrates the distribution of haplotypes when partitioned by these alternative management scenarios. When samples were grouped by Migratory Flyways, there was statistically significant differentiation among populations within Migratory Flyways (Fsc = 0.0336) and within populations (Fst = 0.0624), but there was no statistically significant differentiation between the two Migratory Flyways (Table 3). When the hierarchical model was removed and samples were grouped only according to Migratory Flyways, the pairwise ΦST was statistically significant (Table 4). The hierarchical AMOVA for the BCRs did not detect statistically significant population structure among populations, among regions or within populations (Table 3), but the pairwise ΦST revealed significant genetic differentiation between golden eagles from Alaska, Alberta, British Columbia and Yukon when compared to individuals from California, Idaho and Oregon (Table 4). Additionally, individuals from California, Idaho and Oregon were significantly differentiated from individuals from Montana and Wyoming (ΦST = 0.12). Finally, low but statistically significant levels of genetic differentiation (ΦST = 0.044) was also detected between individual from Montana and Wyoming when compared to the grouping of individuals from Arizona, Colorado and New Mexico (Table 4). Combining the Nebraska and Oklahoma samples with samples from Arizona, Colorado and New Mexico to reflect the Genetic Region grouping, the hierarchical AMOVA did not reveal significant population differentiation at any level tested (Table 3), but similar to other analyses the pairwise ΦST detected statistically significant pairwise genetic differentiation between the grouping of individuals from Alaska, Alberta, British Columbia and Yukon relative to individuals from California, Idaho, and Oregon (Table 4). In congruence with the analysis based on BCRs, this analysis also detected statistically significant genetic differentiation between individuals from California, Idaho and Oregon when compared to individuals from Montana and Wyoming (Tables 3 and 4). Table 3. Results of hierarchical AMOVA when golden eagle samples are partitioned in accordance with Migratory Flyways, Bird Conservation Regions and Genetic Regions. F-statistics show the molecular variance partitioned among groups (Fct), among populations within groups (Fst), and within populations (Fsc)   Migratory Flyways  Bird Conservation Regions  Genetic Regions    d.f.  Sum of squares  Variance components  Percentage of variation  d.f.  Sum of squares  Variance components  Percentage of variation  d.f.  Sum of squares  Variance components  Percentage of variation  Among groups  1  2.01  0.17  2.99  3  3.55  0.03  4.32  3  3.63  0.023  3.99  Among populations  12  8.95  0.02  3.26  7  3.48  −0.01  −0.88  9  4.70  −0.004  −0.71  Within population  156  84.20  0.54  93.75  137  76.13  0.56  96.55  141  79.33  0.56  96.72  Fsc  0.0336        −0.00917        −0.00737        Fst  0.0624        0.0345        0.03278        Fct  0.0299        0.04325        0.0399          Migratory Flyways  Bird Conservation Regions  Genetic Regions    d.f.  Sum of squares  Variance components  Percentage of variation  d.f.  Sum of squares  Variance components  Percentage of variation  d.f.  Sum of squares  Variance components  Percentage of variation  Among groups  1  2.01  0.17  2.99  3  3.55  0.03  4.32  3  3.63  0.023  3.99  Among populations  12  8.95  0.02  3.26  7  3.48  −0.01  −0.88  9  4.70  −0.004  −0.71  Within population  156  84.20  0.54  93.75  137  76.13  0.56  96.55  141  79.33  0.56  96.72  Fsc  0.0336        −0.00917        −0.00737        Fst  0.0624        0.0345        0.03278        Fct  0.0299        0.04325        0.0399        Numbers identified by bold font represent statistically significant values (0.05). View Large Figure 3. View largeDownload slide Map of North American illustrating haplotype distributions when samples are partitioned according to Migratory Flyways, Bird Conservation Regions and Genetic Regions. Pie charts represent the haplotype frequency for the given Migratory Flyway, Bird Conservation Region and Genetic Region. The colours in each map represent the different proposed management units for each model. Figure 3. View largeDownload slide Map of North American illustrating haplotype distributions when samples are partitioned according to Migratory Flyways, Bird Conservation Regions and Genetic Regions. Pie charts represent the haplotype frequency for the given Migratory Flyway, Bird Conservation Region and Genetic Region. The colours in each map represent the different proposed management units for each model. Population fluctuations A unimodal mismatch distribution was detected for the Mediterranean, Palearctic and Nearctic lineages, indicating a more recent demographic expansion, although the only lineage that had a significant P-value for raggedness was the Nearctic (Fig. 2, Table 2). Within the Nearctic groupings, raggedness and SSD were both statistically significant when the samples were grouped into Migratory Flyways. When grouped into BCRs, raggedness was significant for the grouping of individuals from Alaska, Alberta, British Columbia and Yukon as well as the grouping of individuals from California, Idaho and Oregon. When samples were partitioned in accordance with the Genetic Regions, raggedness was only significant for the grouping of our most northern individuals from Alaska, Alberta, British Columbia and Yukon. SSD was not significant for either BCRs or Genetic Regions, adding support to the recent expansion hypothesis. Figure 2. View largeDownload slide Mismatch distribution for each of the lineages: Mediterranean (M), Palearctic (P) and Nearctic (N). The vertical bars represent the observed frequencies while the line represents the expected frequencies under a model of sudden expansion. Figure 2. View largeDownload slide Mismatch distribution for each of the lineages: Mediterranean (M), Palearctic (P) and Nearctic (N). The vertical bars represent the observed frequencies while the line represents the expected frequencies under a model of sudden expansion. DISCUSSION When the range of a species occupies the Holarctic, one of two phylogeographical patterns typically emerge. North American and Eurasian lineages occur in sympatry with a distinct North American haplogroup resulting from recent gene flow across Beringia or the Bering Strait, with examples such as the raven (Corvus corax) (Omland, Baker & Peters, 2006) and the peregrine falcon (Falco peregrinus) (Bell et al., 2014). Alternatively, due to a lack of gene flow across Beringia or the Bering Strait, speciation occurs and distinct species are found in both North America and Eurasia with examples such as white-tailed eagle and bald eagle (Hailer et al., 2006) or three-toed woodpecker (Picoides tridactylus) and American three-toed woodpecker (Picoides dorsalis) (Zink et al., 2000). Our results support the conclusions of Nebel et al. (2015) that there are two distinct mtDNA lineages: a Mediterranean lineage and a Holarctic lineage (Fig. 1). As also determined by Nebel et al. (2015), the Mediterranean lineage represents a monophyletic group, with individuals characterized by these haplotypes being found only in Mediterranean locations. The Holarctic lineage represents haplotypes found throughout the Palearctic and the Nearctic, with only one haplotype (H1) being truly Holarctic. In North America, H1 is the most numerous and most widespread haplotype, making up 45% of individuals sequenced, as opposed to the Palearctic where it represents few (<1%) individuals with a strong location bias in Europe. The widespread range of H1 together with its central location in the haplotype network suggests that this is an ancestral haplotype. The pattern of having a single shared haplotype between hemispheres was also found when using the mitochondrial control region in ravens (Omland et al., 2006). The other haplotypes defined as Holarctic by Nebel et al. (2015) are only found on a single continent, but our data do not recover reciprocal monophyly between the Palearctic and Nearctic haplotypes (Figs 1, 2). Despite not detecting reciprocal monophyly between the Nearctic and Palearctic haplotypes, pairwise genetic differentiation indicates little to no current gene flow between the Nearctic and Palearctic linages (ΦST = 0.517). Additionally, pairwise comparison of ΦST between the Nearctic and Mediterranean indicates almost complete lack of gene flow (ΦST = 0.9139). Among these three large geographical regions, there does appear to be some gene flow occurring between the Palearctic and Mediterranean lineages with a pairwise ΦST of 0.3. Our data also suggest that the Nearctic population of golden eagles has undergone a recent population expansion. When comparing the Nearctic to the Palearctic and Mediterranean lineages, both Tajima’s D and Fu’s F were statistically significant. Together with the unimodal mismatch analysis, which had significant values for raggedness, this supports a recent population expansion. The first mitochondrial analysis of North American golden eagles examined individuals representing the California Channel Islands and the California mainland (Sonsthagen et al., 2012a). Among the 42 individuals for which they were able to obtain mitochondrial sequence data, Sonsthagen et al. (2012a) detected five haplotypes. More recently, Craig et al. (2016) examined the same mitochondrial region of golden eagles from Alaska, California and Idaho. Among the 49 individuals they examined, they detected eight haplotypes, three unique to birds in their study and the five haplotypes previously detected in southern California golden eagles by Sonsthagen et al. (2012a). Additionally, Craig et al. (2016) reported the possibility of two region-specific haplotypes, one that appeared to be restricted to California and a second that appeared to be restricted to Idaho. To these previous mitochondrial studies of North American golden eagles, we add data from 115 additional individuals. Overall, we detected 19 haplotypes with 11 being newly detected variants. Our measures of haplotype diversity are higher for the Nearctic samples than detected by Nebel et al. (2015) for Old World samples, but our measures of nucleotide and haplotype diversity are similar to the results obtained by Craig et al. (2016). Regarding the conclusion of Craig et al. (2016), our data do not support the possibility of their haplotypes being region-specific because N14, which Craig et al. (2016) suggested was restricted to California, was also detected in Oregon in our analyses. Similarly, haplotype N10, which Craig et al. (2016) suggested was restricted to Idaho, was detected in samples from Oregon and New Mexico. Of the 11 new haplotypes that we identified, eight were region-specific. Haplotypes N5 and N17 were only found in Oregon, N6, N7 and N9 were only found in Alaska, N13 was only found in California, N16 was only found in Idaho and N18 was only found in Colorado (Table 1). Because these sequences were only single occurrences, we suspect that with more thorough sampling, they may be found to be in other areas. The only Holarctic haplotype (H1) was found in every sampling location in the Nearctic except Oklahoma, where only one individual was sampled (Table 1). Due to the wide geographical distribution of this haplotype and that locations bordering Oklahoma had this haplotype, we suspect that if sampling efforts were increased in future studies, this haplotype would be found in Oklahoma as well. Due to the uncertainty in how best to manage North American golden eagles (U.S. Fish and Wildlife, 2016), we partitioned our samples so that the proposed Migratory Flyways and BCRs could be tested as well as the results indicated in studies utilizing nuclear DNA (Genetic Regions). Although the hierarchical AMOVAs failed to reveal statistically significant levels of overall genetic differentiation for any of the potential groupings, the analysis of pairwise ΦST did reveal several patterns of significant genetic differentiation (Table 4). Most notably, for both the BCRs and the Genetic Regions, statistically significant pairwise genetic differentiation was detected between the grouping of northern golden eagles (Alaska and Canada) and golden eagles from California, Idaho and Oregon. Moreover, we detected significant genetic differentiation between the California, Idaho and Oregon birds and those from Montana and Wyoming. Both of these groupings are found in different Conservation Regions and Flyways. Table 4. Pairwise ΦST values for Migratory Flyways, Bird Conservation Regions and Genetic Regions   Migratory Flyways`    Bird Conservation Regions    Genetic Regions    West  East    AK-AB-BC-YT  CA-ID-OR  MO-WY  NE-OK  AZ-CO-NM    AK-AB-BC-YT  CA-ID-OR  MO-WY  NE-OK-AZ- CO-NM  Central  0    AK-AB-BC-YT  0          AK-AB-BC-YT  0        Pacific  0.018  0  CA-ID-OR  0.090  0        CA-ID-OR  0.090  0            MT-WY  −0.026  0.117  0      MT-WY  −0.0262  0.117  0          NE-OK  −0.062  0.068  −0.0532  0    NE-OK-AZ- CO-NM  −0.0002  0.015  0.0178  0        AZ-CO-NM  0.020  −0.003  0.044  0.0072  0              Migratory Flyways`    Bird Conservation Regions    Genetic Regions    West  East    AK-AB-BC-YT  CA-ID-OR  MO-WY  NE-OK  AZ-CO-NM    AK-AB-BC-YT  CA-ID-OR  MO-WY  NE-OK-AZ- CO-NM  Central  0    AK-AB-BC-YT  0          AK-AB-BC-YT  0        Pacific  0.018  0  CA-ID-OR  0.090  0        CA-ID-OR  0.090  0            MT-WY  −0.026  0.117  0      MT-WY  −0.0262  0.117  0          NE-OK  −0.062  0.068  −0.0532  0    NE-OK-AZ- CO-NM  −0.0002  0.015  0.0178  0        AZ-CO-NM  0.020  −0.003  0.044  0.0072  0            Numbers in bold are statistically significant at 0.05. View Large Based on the samples that we analysed, the primary difference between the BCRs and the Genetic Regions is that for the BCRs our samples representing Arizona, Colorado and New Mexico are considered to belong to a group genetically differentiated from our samples from Nebraska and Oklahoma, whereas for the Genetic Regions approach the samples from these five states are combined into a single unit. The additional partitioning of samples in the BCR detected significant genetic differentiation between samples from the Northern Plains (Montana and Wyoming) and the southern Rocky Mountain regions (Arizona, Colorado and New Mexico). These findings based on BCRs and Genetic Regions are similar to results reported by Doyle et al. (2016) in which the Alaska and California eagles were significantly different from the other western eagle population. This was also reported by Van Den Bussche et al. (2017), in which the Alaska and British Columbia eagles created a significant unit and the Idaho, California and Oregon eagles created a significant unit. Finally, when partitioning the Nearctic samples into Migratory Flyways, BCRs or Genetic Regions, the flyway subgroupings as well as the most northern golden eagles (Alaska and Canada) and the California, Idaho and Oregon golden eagles appear to have experienced a recent population expansion (Table 2). CONCLUSION Our study provides another example of intermediate divergence and the importance of studying the continuum of divergence for understanding geographical speciation, species limits and conservation priorities (Omland et al., 2006). This study also provides the first insight into the phylogeographical distribution and partitioning of mitochondrial genetic variation for North American golden eagles. To aid in determining more accurate historical events and to aid in determining proper management units, a more thorough sampling scheme should be carried out utilizing nuclear loci. Currently, single nucleotide polymorphisms have been used to begin to elucidate the genetic structure of North American golden eagles, but with limited numbers of samples from throughout the range (Doyle et al., 2016; Van Den Bussche et al., 2017). Future studies that couple mitochondrial and nuclear loci can provide greater insight into the demography of golden eagles by revealing patterns of genetic variation at both maternally and biparentally inherited loci. SUPPORTING INFORMATION Additional Supporting Information may be found in the online version of this article at the publisher’s web-site Table S1. Locality, age, provider and haplotype information for the 115 samples sequenced for the Nearctic portion of the study. ACKNOWLEDGMENTS We thank Bryan Bedrosian, Rob Domenech, Garth Herring, Gail Kratz, Brian Millsap, Gary Roemer, Victor Roubidoux and Brian Smith for contributing golden eagle blood samples from their own research and rehabilitation efforts. Without their kindness and contribution, we would not have had the geographical coverage represented in this study. We further thank the Iowa Tribe of Oklahoma for funding this study in collaboration with the Grey Snow Eagle House and the Ford Foundation for funding Megan Judkins’ research assistantship while the work for this project was performed. We extend our appreciation to the other members of Megan Judkins’ dissertation committee, Drs Meredith Hamilton, Michi Tobler, Andrew Doust and Jim Lish, whose review and feedback of an earlier version of the manuscript was critical. Finally, we thank the four anonymous reviewers of a previous version of the manuscript for their helpful comments. REFERENCES Bayard de Volo S. 2008. Genetic studies of Northern Goshawks (Accipiter Ggentilis): Genetic tagging and individual identification from feathers, and determining phylogeography, gene flow and population history for goshawks in North America . Unpublished D. Phil. Dissertation, Colorado State University. Bell DA, Griffiths CS, Caballero IC, Hartley RR, Lawson RH. 2014. Genetic evidence for global dispersal in the peregrine falcon (Falco peregrinus) and affinity with the Taita falcon (Falco fasciinucha). Journal of Raptor Research  48: 44– 53. Google Scholar CrossRef Search ADS   Bielikova M, Ficek A, Valkova D, Turna J. 2010. Multiplex PCR amplification of 13 microsatellite loci for Aquila chrysaetos in forensic applications. Biologia  65: 1081– 1088. Google Scholar CrossRef Search ADS   Bourke BP, Frantz AC, Lavers CP, Davison A, Dawson DA, Burke TA. 2010. Genetic signatures of population change in the British golden eagle (Aquila chrysaetos). Conservation Genetics  11: 1837– 1846. Google Scholar CrossRef Search ADS   Clement M, Posada D, Crandall KA. 2000. TCS: A computer program to estimate gene genealogies. Molecular Ecology  9: 1657– 1659. Google Scholar CrossRef Search ADS PubMed  Craig EH, Adams JR, Waits LP, Fuller MR, Whittington DM. 2016. Nuclear and mitochondrial DNA analyses of golden eagles (Aquila chrysaetos canadensis) from three areas in western North America; initial results and conservation implications. PLoS ONE  11: 1– 15. Google Scholar CrossRef Search ADS   Doyle JM, Katzner TE, Roemer GW, CainIII JW, Millsap BA, McIntyre CL, Sonsthagen SA, Fernandez NB, Wheeler M, Bulut Z, Bloom PH, DeWoody A. 2016. Genetic structure and viability selection in the golden eagle (Aquila chrysaetos): a vagile raptor with a Holarctic distribution. Conservation Genetics  17: 1– 16. Google Scholar CrossRef Search ADS   Excoffier L, Lischer HEL. 2010. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources  10: 564– 567. Google Scholar CrossRef Search ADS PubMed  Fielding AH, Whitfield DP, McLeod DRA. 2006. Spatial association as an indicator of the potential for future interactions between wind energy developments and golden eagles Aquila chrysaetos in Scotland. Biological Conservation  131: 359– 369. Google Scholar CrossRef Search ADS   Godoy JA, Negro JJ, Hiraldo F, Donazar JA. 2004. Phylogeography, genetic structure and diversity in the endangered bearded vulture (Gypaetus barbatus, L.) as revealed by mitochondrial DNA. Molecular Ecology  13: 371– 390. Google Scholar CrossRef Search ADS PubMed  Hailer F, Helander B, Folkestad AO, Ganusevich SA, Garstad S, Hauff P, Koren C, Nygard T, Volke V, Vila C, Ellegren H. 2006. Bottlenecked but long-lived: high genetic diversity retained in white-tailed eagles upon recovery from population decline. Biology Letters  2: 316– 319. Google Scholar CrossRef Search ADS PubMed  Hewitt G. 2000. The genetic legacy of the Quaternary ice ages. Nature  405: 907– 913. Google Scholar CrossRef Search ADS PubMed  Hewitt GM. 2004. Genetic consequences of climatic oscillations in the Quaternary. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences  359: 183– 195. Google Scholar CrossRef Search ADS PubMed  Hull JM, Girman DJ. 2005. Effects of Holocene climate change on the historical demography of migrating sharp-shinned hawks (Accipiter striatus velox) in North America. Molecular Ecology  14: 159– 170. Google Scholar CrossRef Search ADS PubMed  Longmire JL, Maltbie M, Baker RJ. 1997. Use of “lysis buffer” in DNA isolation and its implication for museum collections. Occasional Papers, The Museum, Texas Tech University  163: 1– 3. Maddison DR, Maddison WP. 2005. MacClade 4: Analysis of phylogeny and character evolution . Sunderland, MA: Sinauer Associates, Inc., Publishers. Masuda R, Noro M, Kurose N, Nishida-Umehara C, Takechi H, Yamazaki T, Kosuge M, Yoshida MC. 1998. Genetic characteristics of endangered Japanese golden eagles (Aquila chrysaetos japonica) based on mitochondrial DNA D-loop sequences and karyotypes. Zoo Biology  17: 111– 121. Google Scholar CrossRef Search ADS   Millsap BA, Harmata AR, Stahlecker DW, Mikesic DG. 2014. Natal dispersal distance of bald and golden eagles originating in the coterminous United States as inferred from band encounters. Journal of Raptor Research  48: 13– 23. Google Scholar CrossRef Search ADS   Morneau F, Tremblay JA, Todd C, Chubbs TE, Maisonneuve C, Lemaitre J, Katzner T. 2015. Known breeding distribution and abundance of golden eagles in Eastern North America. Northeastern Naturalist  22: 236– 247. Google Scholar CrossRef Search ADS   Nebel C, Gamauf A, Haring E, Segelbacher G, Villers A, Zachos FE. 2015. Mitochondrial DNA analysis reveals Holarctic homogeneity and a distinct Mediterranean lineage in the golden eagle (Aquila chrysaetos). Biological Journal of the Linnean Society  116: 328– 340. Google Scholar CrossRef Search ADS   Omland KE, Baker JM, Peters JL. 2006. Genetic signatures of intermediate divergence: Population history of Old and New World Holarctic ravens (Corvus corax). Molecular Ecology  15: 795– 808. Google Scholar CrossRef Search ADS PubMed  Pagel JE, Kritz KJ, Millsap BA, Robert K. 2013. Bald eagle and golden eagle mortalities at wind energy facilities in the contiguous United States. Journal of Raptor Research  47: 311– 315. Google Scholar CrossRef Search ADS   Petit RJ, El Mousadik A, Pons O. 1998. Identifying populations for conservation on the basis of genetic markers. Conservation Biology  12: 844– 855. Google Scholar CrossRef Search ADS   Poulakakis N, Antoniou A, Mantziou G, Parmakelis A, Skartsi T, Vasilakis D, Elorriaga J, De La Puente J, Gavashelishvili A, Ghasabyan M, Katzner T, McGrady M, Batbayar N, Fuller M, Natsagdorj T. 2008. Population structure, diversity, and phylogeography in the near-threatened Eurasian black vultures Aegypius monachus (Falconiformes; Accipitridae) in Europe: Insights from microsatellite and mitochondrial DNA variation. Biological Journal of the Linnean Society  95: 859– 872. Google Scholar CrossRef Search ADS   Sonsthagen SA, Coonan TJ, Latta BC, Sage GK, Talbot S. 2012a. Genetic diversity of a newly established population of golden eagles on the Channel Islands, California. Biological Conservation  146: 116– 122. Google Scholar CrossRef Search ADS   Sonsthagen SA, Rosenfield RN, Bielefeldt J, Murphy RK, Steward AC, Stout WE, Driscoll TG, Bozek MA, Sloss BL, Talbot SL. 2012b. Genetic and morphological divergence among Cooper’s Hawk (Accipiter cooperii) populations breeding in North-Central and Western North America. The Auk  129: 427– 437. Google Scholar CrossRef Search ADS   Stauber E, Finch N, Talcott PA, Gay JM. 2010. Lead poisoning of bald (Haliaeetus leucocephalus) and golden (Aquila chrysaetos) eagles in the U.S. inland Pacific Northwest region – an 18-year retrospective study: 1991–2008. Journal of Avian Medicine and Surgery  24: 279– 287. Google Scholar CrossRef Search ADS PubMed  Suchentrunk F, Haller H, Ratti P. 1999. Gene pool variability of a golden eagle (Aquila chrysaetos) population from the Swiss Alps. Biological Conservation  90: 151– 155. Google Scholar CrossRef Search ADS   U.S. Fish and Wildlife Service. 2016. Bald and Golden Eagles: Population demographics and estimation of sustainable take in the United States, 2016 update . Washington: Division of Migratory Bird Management. Van Den Bussche RA, Judkins ME, Montague MJ, Warren WC. 2017. A resource of genome-wide single nucleotide polymorphisms (SNPs) for the conservation and management of golden eagles (Aquila chrysaetos). Journal of Raptor Research  51: 368– 377. Google Scholar CrossRef Search ADS   Watson J, Brockie K, Watson D. 2010. The golden eagle, 2nd edn . New Haven: Yale University Press. Wink M, Clouet M, Goar LJ, Barrau C. 2004. Sequence variation in the cytochrome b gene of subspecies of Golden Eagles Aquila chrysaetos. Alauda  72: 153– 157. Zink RM, Barrowclough GF, Atwood JL, Blackwell-Rago RC. 2000. Genetics, taxonomy, and conservation of the threatened California Gnatcatcher. Conservation Biology  14: 1394– 1405. Google Scholar CrossRef Search ADS   © 2017 The Linnean Society of London, Biological Journal of the Linnean Society

Journal

Biological Journal of the Linnean SocietyOxford University Press

Published: Feb 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve Freelancer

DeepDyve Pro

Price
FREE
$49/month

$360/year
Save searches from Google Scholar, PubMed
Create lists to organize your research
Export lists, citations
Access to DeepDyve database
Abstract access only
Unlimited access to over
18 million full-text articles
Print
20 pages/month
PDF Discount
20% off