Introduction

The Golden Eagle, Aquila chrysaetos, is an accipitrid of boreal distribution that inhabits a variety of open and semi-open biotopes from sea level to 3,630 m in altitude in biomes as tundra, chaparral, temperate grassland, temperate deciduous forest, and coniferous forest (^{De León-Girón et al., 2016}; ^{Flesch et al., 2020}; ^{Kochert et al., 2002}; ^{Watson et al., 2011}). This emblematic species is known to occur in Mexico in the arid and semiarid environments of the northern and central regions, including some locations as south as Oaxaca (^{Bolger et al., 2014}; ^{De León-Girón et al., 2016}; ^{Howell & Webb, 1995}; ^{Rodríguez-Estrella, 2002}; ^{Rodríguez-Estrella et al., 1991}, ²⁰²⁰). Although the demography of the Golden Eagle has been determined with different genetic methods that confirm the structuring of its populations (^{Craig et al., 2016}; ^{Doyle et al., 2016}), little is known about the genetic identity of populations of this eagle in its distribution range. Northwestern Baja California is considered the greatest nesting and conservation potential region for the Golden Eagle in Mexico (^{De León-Girón et al., 2016}; ^{Rodríguez-Estrella, 2002}; ^{Rodríguez-Estrella et al., 1991}; ^{Tracey et al., 2017}), with vast areas of habitats still pristine for the population conservation of this species shared with the United States of America (^{Craig et al., 2016}; ^{De León-Girón et al., 2016}, ²⁰²⁴; ^{Doyle et al., 2014}; ^{Katzner et al., 2023}).

Studies based on mitochondrial DNA (mtDNA) allow determining the current state of conservation of species, as well as the identification of significant evolutionary and management units with their demographic aspects (^{Moritz, 1994}). With the advent of next-generation sequencing (NGS) technologies (especially on complete mitochondrial genomes), new studies with non-model species became more common. As was the case for various species of the family Accipitridae, including the Golden Eagle. ^{Doyle et al. (2014)} were the first to describe the nuclear and mitochondrial genome of the Golden Eagle from the central California region, in the USA.

In this study we focus on the synteny of mitochondrial genes, the phylogenetic relationships, and the mitochondrial genetic diversity among Golden Eagle individuals in northwestern Baja California, Mexico. Given the high migration potential of this species, individuals collected in northwestern Baja California are predicted to be phylogenetically related to the westernmost Golden Eagle population in the USA. Therefore, studies focused on examining the genetic variation at the mitochondrial level will allow better decision-making for management and conservation programs for the species in a binational spectrum. In the same way, our study will be a valuable reference for analyzing other Golden Eagle populations in Mexico and other regions of its distribution range.

Materials and methods

Tissue samples of 6 Golden Eagle specimens from northwestern Baja California, Mexico were obtained. These specimens were found dead in different agriculture valleys of the Mediterranean region in northwestern Baja California, Mexico, between the municipalities of Rosarito and San Quintín, from 1995 to 2014, and deposited as vouchers in the Bird Collection of the Science Faculty, at the Autonomous University of Baja California (UABC), campus Ensenada (Fig. 1, Supplementary material: T1).

Map by Rafael Hernández Guzmán.

Figure 1 Geographical location of the Golden Eagle voucher specimens examined for mitogenome in the northwestern Baja California, Mexico. Geographic coordinates for each individual can be found in Supplementary material T1.  

For the DNA extraction, tissue samples of the pads of the feet of 5 Golden Eagle specimens referred to here as Ach 1 to Ach 5 were obtained. DNA using the nucleic acid purification method by differential saline precipitation (^{Aljanabi & Martinez, 1997}) was extracted. Both feathers and blood for the Ach 6 sample (Supplementary material: T1) were used and extracted DNA with the Qiagen Blood & Tissue kit. The quality and quantity of the extracted DNA in both cases was evaluated using Nanodrop, Qubit and agarose gel electrophoresis.

Extracted DNA (from 77 to 210 ng/µL per sample quantified in Nanodrop for Ach1-Ach5 individuals and 0.581 ng/µL for Ach6 quantified in Qubit) was purified with SpeedBeads magnetic beads (Thermo-Scientific, Waltham, MA, USA) at a 1.2:1 beads:DNA ratio and resuspended in 50 µl of TLE1X buffer. Libraries for shotgun sequencing were prepared, for which purified DNA from each individual was fragmented by sonication in a Bioruptor® (Diagenode, Liege, Belgium) using 2 rounds, each consisting of 5 cycles of 30 sec of sonication and 30 sec without sonication at the highest setting. DNA fragments were then prepared using the Kapa Biosystems® Hyper Prep Kit (KR0961-v4.15, Roche, Basel, Switzerland), with which end repair and A-tailing were performed, adapters were ligated, and PCR was performed with indexed primers (^{Glenn et al., 2020}) for 14 cycles. Following amplification, fragment size selection was performed by double- dip with SpeedBeads magnetic beads (Thermo-Scientific, Waltham, MA, USA) that allow to preserve fragments between ~250 and ~700 base pairs (bp) for each of the samples, which are the appropriate insert size for the Illumina platform. Sequencing was performed paired end on 2 different platforms. We sequenced the Ach1-Ach5 samples on the Illumina MiSeq V.3 at the Georgia Genomics and Bioinformatics Core (GGBC) to generate 300- bp paired- end fragments, and we sequenced Ach6 on an Illumina HiSeq 4000 at the Oklahoma Medical Research Foundation Clinical Genomics Center to generate 150-bp paired-end fragments.

Bioinformatic analysis. The quality of the raw sequences was evaluated using the FastQC program (Andrews, 2010). Subsequently, the sequences using a standard treatment were filtered using the Trimmomatic software (^{Bolger et al., 2014}) with 4 steps. In the first step, the sequencing adapters and over-represented sequences were removed by means of the ILLUMINACLIP function. The second step discarded sequences with a quality value below 30 QS with the AVGQUAL function. The third step removed fragments below 25 QS with the SLIDINGWINDOW function. The fourth step did an even stricter cleanup using MAXINFO. This step conserved sequences of at least 100 base pairs and was configured to preferentially conserve longer sequences with a value of 0.3. For the Ach 5 specimen, an extended and modified version of the standard treatment was used due to the low quality of sequences. This involved applying the SLIDINGWINDOW function, which modifies the quality value.

Assembly and annotation. Two strategies were followed to assemble the mitochondrial genomes: the reference based and the de novo strategy (^{Machado et al., 2015}, ²⁰¹⁸). Mapping against a reference mtDNA genome was performed using the Bowtie2 program version 2.3.4.2 with clean reads (^{Langmead & Salzberg, 2012}). The reference mitochondrial genome was from the same species A. chrysaetos, with the GenBank accession number of KF905228.1. De novo assembly only for samples Ach 2 to Ach 6 was performed using the A5-miseq pipeline (^{Coil et al., 2015}). Due to the large number of reads after the filter (Supplementary material: T2) sample Ach 1 was analyzed with the Velvet software version 1.1 (^{Zerbino & Birney, 2008}). The mitogenome was annotated from the scaffolds obtained from the de novo assemblies. We determined scaffolds longer than 4 Kpb by applying a search with the BLAST tool (Basic Local Alignment Search Tool) of the NCBI (National Center for Biotechnology Information) website. The scaffolds were aligned to achieve greater coverage considering the length of the reference genome used previously. Afterwards, RNAweasel (http://megasun.bch.umontreal.ca/RNAweasel/) was implemented to the identification of tRNA’s (transfer RNA), rRNA’s (ribosomal RNA) and introns and, in turn, MFannot (http://megasun.bch.umontreal.ca/RNAweasel/) for the identification of proteins and open reading frames (^{Sieber et al., 2018}).

In addition to the previous annotation, MITOS web server (^{Bernt et al., 2013}) was used to assist in the annotation of de novo mitochondrial genomes, allowing gene names, tRNA and rRNA secondary structures, and codon usage to be obtained. Finally, a manual curation of the annotations was used to review the 6 existing reading frames with the UGENE software (^{Okonechnikov et al., 2012}). The control region was identified based on 99% similarity with an A. chrysaetos partial control region sequence (EF459579.1). Genome Vx software (^{Conant and Wolfe, 2008}) was used to map the Golden Eagle’s mitogenome, using individual Ach 6 to achieve this analysis.

In the annotation of the de novo assembly, the ND3 protein appeared divided into 2 fractions (see Results), a trait that was not reported by ^{Doyle et al. (2014)}. Therefore, an experimental verification was required through PCR amplification and its subsequent Sanger sequencing. We designed primers from de novo assembly using the NCBI Primer-BLAST program: ND3Ach-1 (5’GCCTGATACTGGCACTTCGT3’ and 5’CCCTATCA-ATCTGACCCACCG3’) which generates a 716 bp fragment and ND3Ach-2 (5’CTTCTTCGTCGCTACAGGCT3’ and 5’CCTTCCACCGAACCCACTTAA3’) which generates a 774 bp fragment. Eight Golden Eagle samples from the Ornithological Collection of the Faculty of Sciences of the Autonomous University of Baja California were used for PCR amplification. These samples correspond to the 6 individuals used for sequencing and 2 extra samples that are not part of the mitogenome assembly study. The conditions of the PCR reaction were as follows: 6 min at 94°C, 35 cycles of denaturation, annealing, and extension at 94°C for 30 sec, 51°C for 30 sec and at 72°C for 1.5 min, respectively. A final elongation stage at 70°C for 7 min, and a final conservation stage at 15°C. We sent the amplified fragments for sequencing to the SeqXcel.inc Company in San Diego, California, and analyzed with the ABI PRISM® 3130xl Genetic Analyzer DNA kit.

Synteny and search for polymorphisms. A posteriori synteny analysis was performed between the Golden Eagle mitogenome reported by ^{Doyle et al. (2014)} and the newly assembled mitochondrial genome of this study. This analysis was carried out with the MAUVE software (^{Darling et al., 2004}). The same reference mitogenome was used to call SNP’s and INDEL’s for the 6 genomes assembled in this work using the SAMtools program (^{Li, 2011}). Filtering the call quality of the variants involved discarding those with a value less than 20 Qs and keeping those with at least 5X depth (^{Li, 2011}). Each variant was manually evaluated, retaining SNPs and INDELs where at least half of the aligned reads showed the change and verified its existence in at least one other individual (^{Fridjonsson et al., 2011}).

Phylogenetic relationships. Phylogenetic relationships were inferred using the individual gene dataset and the mitogenome. The individual gene fragments analyzed consisted of the ND2, ND3, coxI and Cytb recovered from Aquila chrysaetos individuals from GenBank (Supplementary material: T3) and those produced in this study. The mitogenome analysis for phylogenetic relationships consisted of 10 complete mitogenomes: 6 obtained in this study, 2 reference samples of Aquila chrysaetos (LR822062 from the United Kingdom and NC_024087.1 from California, USA), 1 reference sample classified as Aquila heliaca (NC_035806.1), but has been reported to actually consist of a A. chrsyasetos individual (^{Sangster & Lukesenburg, 2021}) and 1 outgroup (Aquila nipalensis, GenBank accession number NC _ 045042.1). Both datasets were aligned using the MAFFT 7 algorithm (^{Katoh & Standley, 2013}). Removing the non -conserved regions between the sequences of the multiple alignments was done using the online program Gblocks with the default parameters (^{Castresana, 2000}). The nucleotide substitution model and the mutation rate were determined using JModelTest (^{Posada, 2008}), considering the Bayesian information criterion (BIC).

Two phylogenetic reconstruction methods were used: Maximum Likelihood (ML) and Bayesian analysis. In IQ-TREE (^{Nguyen et al., 2015}) the ML method (^{Felsenstein, 1981}) was run with the DNA substitution model selected by ModelFinder (^{Kalyaanamoorthy et al., 2017}) and 1,000 bootstrap pseudo replicates were performed (^{Hoang et al., 2017}). Bayesian inference was used with MrBayes program (^{Ronquist et al., 2012}); for this, 4 Markov Monte Carlo chains (MCMC) were implemented with a total of 10 million generations, and a sampling every 1,000 generations. To construct the consensus phylogenetic tree, a 25% of burn-in was applied. The support of the nodes was evaluated via the posterior probability values. The results of both approaches were visualized with FigTree 1.4 (^{Rambaut, 2009}). For the mitogenome phylogenetic tree, the dataset on genes that were annotated across the entire mitogenome were partitioned (Supplementary material: T4). To determine the most appropriate DNA substitution model for each gene, the Akaike Information Criterion test implemented in jModelTest2 (^{Darriba et al., 2012}; ^{Guindon & Gascuel, 2003}) was used (see selected models for each gene in Supplementary material: T4). The partitioned mitogenome dataset was used to estimate a Bayesian phylogenetic tree in MrBayes (^{Ronquist et al., 2012}). The Bayesian analysis was carried out by running 4 MCMC chains for 5 million generations and saving the trees every 1,000 generations. The consensus phylogenetic tree was constructed using a 10% burn-in.

Genetic diversity. Individual DNA segments from the coxI, Cytb, ND2, and -ND3 genes were used to estimate the genetic diversity parameters for the A. chrysaetos species complex since these were the genes that had a relative amount of A. chrysaetos DNA sequences available in GenBank. Using DNA sequences from the Baja California individuals produced in this study and reference sequences from the species complex (Supplementary material: T3) the following indices were estimated with DnaSP 6 (^{Rozas et al., 2017}): the number of segregating sites (S), number of haplotypes (Nh), haplotype diversity (Hd), and nucleotide diversity (π). Genetic diversity indices were calculated among all sequences included in the respective datasets and within specific groups of sequences according to phylogenetic clades or geographic origin of samples. In addition, the genetic distance (uncorrected p-distance) between and within the phylogenetic clades observed within the mitochondrial genetic lineages in PAUP software was estimated (Swofford, 1993). The clade OTU1 included sequences from a diverse geographical background, including all the DNA sequences of the Baja California golden eagles. In contrast, the A. chrysaetos DNA sequences grouped in the other closely related clade were classified as OTU2. Additionally, the DNA sequence of an individual of A. chrysaetos from California (Southern Sierra Nevada) was analyzed separately to compare its genetic distance with those from Baja California, given the geographical proximity between both populations.

Results

Raw sequencing reads ranged from 3,324,579 (individual Ach 1) to 539,825 (individual Ach 4). The number of total sequences after the filter ranged between 3,251,746 (Ach 1) and 526,875 (Ach4) (Supplementary material: T2A).

Assembly and annotation. There were large differences in assembly between individuals; for example, the average depth in the reference mapping for individual Ach 1 is 30.5X, while for Ach 6 is 168.3X, even though Ach 6 had fewer total reads than Ach 1. Assemblies for individuals Ach 1-5 showed fewer aligned sequences (Supplementary material: T2A). For de novo assembly, the average depth for individual Ach 6 is 8.0X. The longest scaffold generated (17,472 bp) allowed us to assemble the entire mitogenome of the Ach 6 individual, eliminating the need to align multiple scaffolds to recover the full length (Supplementary material: T2B). The data yield was inferior when doing de novo assembly for Ach 1-5 individuals, and none of these individuals allowed the recovery of the complete mitochondrial genome (Supplementary material: T3B).

Based on the above results, we annotated and curated the mitochondrial genome using only the complete mitogenome obtained by de novo assembly of the Ach 6 individual. All 3 annotation tools identified tRNAs, and none of the tools identified introns within the sequence. The MITOS program (Fig. 2, Supplementary material: T4) was the only software that identified all the genetic elements, such as the 22 tRNAs, 2 rRNAs, and the 13 protein-coding genes that make up the mitogenome. This last annotation was the most complete, so the curation was made from it. We found a region with 7 nucleotide bases that is composed of an AGA stop codon of the ND3 protein, followed by an adenine (A) and finally an ATC start codon that codes for isoleucine. This region divides the ND3 fragment into the 2 regions proposed here (ND3a and ND3b, Fig. 2).

Figure 2 Mitogenome annotated by de novo assembly for Golden Eagle (Ach6) from Baja California, Mexico. In the central part of the diagram in a ring form is shown a scale of the length of mitochondrial genome clockwise. Enlarged view of the section showing the division into 2 fractions of the mitochondrial gene for NADH dehydrogenase subunit 3 (ND3). The stop codon of the ND3a fraction is highlighted in red, the insertion of nucleotide A in yellow, and the start codon of the ND3b fraction in green. The ND3 gene is composed entirely of 352 bp; the ND3a fraction is 207 bp long, while the ND3b fraction is 144 bp long. 

Synteny and search for polymorphisms. The linear order of the mitogenome of Ach 6 coincides with the reference genome (Accession number MT319112.1, Supplementary material: T5). We consider the SNPs and INDELs identified in the Ach 6 reliable, given that the average coverage is 168.3X, and their call quality was high (Supplementary material: T2A). Most SNPs and INDELs are found within genes that code for some protein, transfer RNA and ribosomal genes. The coxI gene presented the greatest number of SNPs with respect to the reference genome in the 6 individuals, while the Cytb and tRNA F genes presented the greatest length changes in the INDELs. In addition, in the ND3 gene, there is a consistent change in the sequences of the 6 individuals with respect to the reference genome (Table 1).

Phylogenetic relationships. With the use of the mitogenome we constructed a phylogenetic tree that identified 2 clades with strong support (1): OTU1 that includes all the individuals analyzed in this study, plus 2 GenBank sequences (Aquila chryseatos, LR82062.1 from the United Kingdom and A. heliatica, NC_ 035806.1), and OTU2 that corresponds to a GenBank sequence of a specimen collected in California (Fig. 3). For their part, the phylogenetic trees obtained with the fragments of individual genes (ND2, ND3, Cytb and coxI) did not resolve the divergence between individuals from Baja California and California obtained with the mitogenome, but all except the ND3 gene showed a divergence (OTUs 1 and 2), mainly with respect to eagles from the European continent (when information on the collection location is reported, Supplementary material: F1).

Figure 3 Bayesian phylogenetic tree based on the mitogenome for Golden Eagle voucher specimens from northwestern Baja California, Mexico. Posterior probabilities are shown above internal nodes. GenBank accession numbers are positioned next to reference sequences. Aquila nipalensis was used as an outgroup for rooting the tree. Two Operational Taxonomic Units are identified (OTU 1, 2). 

Genetic diversity. We observed low genetic diversity, for instance, ND3, ND2, and coxI each had a single segregating site (Table 2), resulting in low nucleotide diversity values (0.0005, 0.0008, and 0.0002, respectively). In contrast, Cytb was the most polymorphic, with a total of 10 segregating sites and a nucleotide diversity an order of magnitude larger (0.00297) (Table 2). The genetic diversity found within OTU1 was higher than that found in OTU2 (Table 2), although OTU1 was not a geographically homogeneous clade among the trees of the different mitochondrial genes, for example in the case of ND2 it was made up of North American eagles, including the individuals of this study, but coxI apart from the previous ones was constituted with eagles from Sweden, Norway and Japan (Table 2).

Discussion

Despite using 2 sequencing platforms for the mitogenome of Golden Eagle from northwestern Baja California, Mexico, we successfully assembled the mitochondrial genome of all 6 individuals by the reference method. This included some individuals (Ach 1-5) with non-consensus regions. ^{Bolger et al. (2014)} recommend a pre-processing step of the readings before any analysis, be it assembly by reference or de novo, and mention that if the library and identification adapters are not removed, they can be incorporated in the final assembly. In our case, the pre-processing of the sequences positively influenced the performance of both assemblies and the SNP’s call since the quality of the sequence of the 6 individuals was considerably improved after pre-processing in all aspects of quality reporting.

The most reliable mitogenome assembly was that of Ach 6 individual due to its higher alignment rate, its greater depth, and its 100% coverage of the reference genome. This ensemble can be considered true and not an artifact due to the procedures used, such as sampling every 500 generations and a burn-in value of 600,000 generations (^{Lerner & Mindell, 2005}). In addition, we analyzed the data with Bowtie2, which implements an alignment strategy based on the FM-Index and Burrows-Wheeler. These programs have been shown to work well for applications such as INDEL discovery (^{Lindner & Friedel, 2012}) and for aligning long sections of the reference genome (^{Thankaswamy-Kosalai et al., 2017}). ^{Hunt et al. (2014)} pointed out that the best result of assembling a genome de novo is when it is contained in a single scaffold for the entire mitochondrial DNA molecule or for each chromosome in the case of the nuclear genome. The Golden Eagle mitochondrial genome assembly of individual Ach 6 was the best because it was found within a single scaffold (17,472 bp) and its depth was at least 8.0X (Supplementary material: T2B; ^{Baker, 2012}).

The extension and position of the genes in the mitogenome of the Golden Eagle described in this study (Fig. 2, Supplementary material: T4) showed differences compared to those reported by ^{Doyle et al. (2014)} from an individual sampled in central California. For example, the size of the mitogenome described here is larger (17,472 bp) than the one reported by ^{Doyle et al. (2014)} (17,332 bp), with a difference of 140 bp. When experimentally analyzing the size of the ND3 gene amplified fragments, the size was between 700 and 800 bp, that is, the theoretically expected size. Another difference consisted in the fraction of the gene encoding the ND3 protein reported here as split into 2 segments (ND3a and ND3b). The distinction of the 2 subunits is, within the reading frame, a region made up of 7 nucleotides between the ATC start codon that starts the protein, and the TAA stop codon that marks the end of the protein. The first 3 of these 7, code for a new AGA stop codon, followed by an adenine (A) nucleotide and 3 nucleotides that code for a new ATC start codon. This causes 2 coding portions to coexist, annotated as ND3, separated by a single nucleotide that changes the reading frame, causing the second portion to appear with its start and stop codons. ^{Mindell et al. (1998)} report a nucleotide that is not translated within the ND3 protein sequence. ^{Eberhard and Wright (2016)} mentioned that such a trait is observed in the entire order Psittaciformes (parrots, parakeets, and allies).

Specifically, a variety of distinct indels have been found within several of the mitochondrial protein-coding genes in Psittaciformes (^{Eberhard & Wright, 2016}), with variations in terms of their evolutionary origin since some are present in all Psittaciformes and some are more recent in origin and only found within specific taxa. ^{Slack et al. (2003)} reported variation in the length of the mitochondrial ND6 gene in other avian taxa. Overall, the biological implications of these variations remain unclear until more studies on the proteins are performed to determine if the protein compositions become altered by these variations. However, the experimental verification carried out in the present study strengthens our annotation as 2 fractions for the ND3 region of the A. chrysaetos mitogenome. Also, we noted the regulatory non-coding region (the pseudo control region, ψCR) to be highly conserved among the A. chrysaetos individuals analyzed, which is a relatively conserved region within the order Accipitriformes (cf. ^{Song et al., 2015}) and is also found in the same position between the tRNA-E and tRNA-F, as shown by ^{Liu et al. (2017)} for the mitogenome of Accipiter gularis, another congeneric species.

Synteny and search for polymorphisms. The assembled Golden Eagle mitogenome did not present changes in gene order with respect to the reference (KF905228.1, a male Golden Eagle from southern Sierra Nevada, California). As expected, neither did it present changes in the order of the genes with respect to other species of the same group of birds (Accipitridae), such as Buteo buteo (^{Haring et al., 2001}), Accipiter virgatus (^{Song et al., 2015}), Aquila fasciata (^{Jiang et al., 2015}) and Accipiter gularis (^{Liu et al., 2017}). In addition, it is consistent with one of the general models of the birds’ mitogenome that is characterized by a duplication of the control region and other adjacent genes, which were subsequently degraded, giving rise to the ψCR (^{Eberhard & Wright, 2016}). SNP’s and INDEL’s observed in the mitogenome of the Ach 6 individual meet the needed quality (probability of not being an error) and sequence depth (^{Li, 2011}), so they can be considered as potential markers that should be tested in future population studies, considering some technical aspects, and thus avoid the call of false positives (^{The 1000 Genomes Project Consortium, 2015}).

Reconstruction of phylogenetic relationships. The phylogenetic tree classified all Golden Eagle individuals into 2 mitogenome clades (OTU1 and OTU2). All individuals from Baja California, an individual from the United Kingdom, and an A. heliaca individual clustered within OTU1 (Fig. 3), while OTU2 was composed of a single Golden Eagle individual from California. The genetic distance between the 2 OTUs was 0.002 nucleotide differences per site (Supplementary material: T6). This pattern is not surprising given the poor sampling of taxa for both the mitogenomes and individual genes of A. chrysaetos. Once additional mitogenomes are available for this species, it will become possible to increase the phylogenetic resolution that could confirm or reject the potential presence of multiple mitochondrial lineages circulating within Golden Eagle populations. Nevertheless, given the broad geographic area included in the mitogenome dataset (i.e., an A. chrysaetos individual from the United Kingdom, as well an A. heliaca individual), and the low genetic diversity observed within these mitochondrial sequences (see below), suggests that at the mitochondrial level, the genetic diversity found within A. chrysaetos populations is low.

The use of genes or fragments of the mitochondrial genome allows for the analysis to include a larger data set in the number of locations and individuals. The topology of 3 of the 4 individual gene trees (ND2, Cytb and coxI) did support the presence of at least 2 mitochondrial clades (Supplementary material: F1). However, these gene trees were not consistent when grouping the North American Golden Eagle individuals on the same OTUs.

Future studies targeting molecular markers with a higher mutational rate (i.e., microsatellites) or genomics approach could potentially uncover the presence of genetic differences that are not well defined at the mitochondrial level but might be biologically important. One such difference is seen in the species’ bivalent behavior and its migration pattern in the west (^{Mcintyre & Lewis, 2016}) and northwest of the USA, and through Canada (Bedrosian et al., 2018). The mitogenomic lineages that potentially infer the difference between the eagles of Baja California and central California could correspond to the area where 2 ecological populations are coexisting (^{De León-Girón et al., 2016}). The first population is migratory, originating in western North America (Oregon), while the second lineage is resident in Baja California. According to ^{Craig et al. (2016)}, this lineage has 1 of the less frequent haplotypes of the mitochondrial control region and is restricted to California. The ability of the Golden Eagle to adapt to different environments (^{Judkins & Van Den Bussche, 2018}), the large dispersal distances of reproductive or floating between Mexico and the USA (^{De León-Girón et al., 2016}, ²⁰²⁴; ^{Rodríguez-Estrella et al., 2020}; ^{Tracey et al., 2017}), and its extensive home range of reproductive pairs (^{D’Addario et al., 2019}), would support these hypotheses. Besides, the presence of these groups of individuals (reproductive and non-reproductive) in Baja California would be part of the reproductive behavior of the species (^{Watson et al., 2011}), that is, the process of succession and substitution of reproductive pairs in the region (^{De León-Girón et al., 2016}).

Genetic diversity. The values of mitochondrial genetic diversity calculated for most genes (except Cytb) were near zero (Table 2). We expected this result, considering that they are genes that code for proteins with highly conserved regions (^{Dawnay et al., 2007}) and are subject to biochemical limitations that cause high levels of homoplasy (^{Faria et al., 2007}). However, ^{Bates et al. (2003)} mentioned that ND2 is genetically more diverse than Cytb, while ^{Faria et al. (2007)} stated that ND2 is one of the most variable mitochondrial genes within birds and, therefore, regularly implemented in population genetics. A pattern we did not observe in A. chrysaetos since Cytb was the gene with the highest genetic diversity observed, with a haplotype diversity of 0.154 and 0.667 for OTU1 and OTU2, respectively, while both OTUs had a haplotype diversity of zero for the ND2 gene (Table 2).

Implications for conservation. Golden eagles are recognized for presenting large expanses of territory, and their resident, migratory and floating populations (reproductive adults without territories), promote the population gene flow (^{Craig et al., 2016}; ^{Poessel et al., 2022}). The mitochondrial DNA data produced in this study confirm a very close genetic relationship between the northwestern Baja California individuals and those from the USA. Therefore, the execution of binational conservation programs by the agencies of both countries (SEMARNAT and USFWS) are priorities for the “Californian-Baja Californian” metapopulation of Golden Eagle. It is necessary to increase the number of young individuals with satellite tracking in the southern Baja California peninsula, to continue the genetic characterization (with different types of markers, mitochondrial DNA, microsatellites and SNPs obtained by Next Generation Sequencing) and population monitoring to evaluate the structure and connectivity of the species in both countries.

Our findings warrant developing joint conservation efforts between the governments of Mexico and the USA to monitor and preserve the North American Golden Eagle.