The Javan Hawk-Eagle (Nisaetus bartelsi), endemic to the stunning landscapes of Java, Indonesia, stands as a national symbol embodying both cultural significance and ecological importance. It is designated as the national bird of Indonesia, and its representation within the national emblem, GARUDA, symbolizes the strength and identity of the Indonesian state^1,2. The presence of this hawk-eagle in the GARUDA emblem elevates its status, making it a critical species for conservation efforts as it embodies national pride and cultural identity^3–6. The Javan Hawk-Eagle (N. bartelsi) plays a critical role in Indonesia's biodiversity as both a top predator and a symbol of national pride, being the country's national bird. As an endemic species, it is uniquely adapted to the specific ecological conditions found on the island of Java ^3,7. Despite its emblematic status, the Javan Hawk-Eagle is currently listed as Endangered, facing serious threats from rapid habitat loss, fragmentation, and human disturbance^3,89. The once expansive forests of Java, which provided crucial habitat for this species, have become fragmented, resulting in population declines that warrant immediate conservation action¹⁰. The Javan hawk eagle is also threatened by wildlife trade^2,11,12, particularly the illegal pet trade, despite being declared Indonesia's National Rare/Precious Animal in 1993³. These pressures underline the urgency for science-based conservation strategies that not only protect habitat, but also incorporate information on the species’ genetic status and evolutionary history.

Molecular techniques have become essential for conservation, enabling the identification of distinct populations and the evaluation of their genetic diversity, thereby enhancing conservation efforts. ^13–16. Genetic studies based on mitochondrial DNA (mtDNA) have become a cornerstone for understanding population dynamics and evolutionary relationships of many species ^17–20. However, despite the critical need for such information, a complete mitochondrial genome sequence for Nisaetus bartelsi has not yet been made available. Current genetic research on Nisaetus bartelsi is limited to identification through COI ²¹ and molecular sexing techniques ⁵. The absence of this genomic information hampers efforts to resolve phylogenetic relationships within the Accipitridae family and constrains our ability to implement informed conservation strategies. The necessity of this research cannot be overstated, as the complete mitochondrial genome will serve as an essential resource for various studies, including phylogenetics^22,23, population genetics^24,25, and evolutionary biology^22,26,27. Mitochondrial genomes are particularly useful in providing insight into evolutionary histories due to their rapid mutation rates and maternal inheritance patterns ^15,17,28,29. Furthermore, understanding the complete mtDNA sequence of the Javan Hawk-Eagle will enhance our knowledge of its genetic diversity and population structure, which are crucial for effective conservation planning.

This study aims to generate and characterize the complete mitochondrial genome of Nisaetus bartelsi, providing the first mitochondrial genomic resource for this endemic and endangered species. Specifically, we sequenced the full mitogenome using high-coverage Illumina reads, annotated all gene features, analyzed nucleotide composition and codon usage patterns, and reconstructed the phylogenetic position of the Javan hawk-eagle within Accipitridae using time-calibrated methods. By establishing this mitogenomic baseline, we anticipate that future researchers will be able to conduct high-resolution population-genetic studies, apply forensic markers to detect illegal trade, and better understand the species' evolutionary history and genetic distinctiveness. The availability of a complete and well-annotated mitochondrial genome reference will strengthen the scientific foundation for conservation planning and management of one of Indonesia's most threatened and culturally significant raptors.

Sample Collection and DNA Extraction

A naturally shed feather of Nisaetus bartelsi was collected from the Natural Resources Conservation Agency of East Java (BKSDA Jawa Timur). Sampling authorization and genetic material access were granted by the Ministry of Forestry under permit no. 254/2024. Total genomic DNA was isolated using the DNeasy Blood & Tissue Kit (Qiagen, Germany) following the manufacturer's instructions. DNA quality and quantity were assessed using the Qubit dsDNA High Sensitivity assay (Thermo Fisher Scientific) and Agilent TapeStation 4150 (Agilent Technologies). The extracted DNA yielded a concentration of 281.7 ng/µL with an A260/A280 absorbance ratio of 1.81, and TapeStation analysis indicated high-molecular-weight genomic DNA with a mean fragment size of approximately 58,347 bp and a DNA Integrity Number (DIN) of 8.1, confirming suitability for downstream sequencing applications.

Library Preparation and High-Throughput Sequencing

Double-stranded DNA fragments were processed using the Nextera XT DNA Library Preparation Kit (Illumina, USA) according to manufacturer protocols. Prepared libraries were purified with the Agencourt AMPure XO magnetic bead system (Beckman Coulter, USA) to remove contaminants and small fragments. Sequencing was performed on an Illumina NextSeq 500 platform (Illumina, USA) generating paired-end reads of 150 base pairs each (2 × 150 bp), using the TruSeq DNA Library Preparation Kit for library construction. Multiple samples were pooled and cleaned using Highprep magnetic beads (Magbio, USA) before sequencing. Raw base calls in BCL format were converted to FASTQ files using bcl2fastq v.2.20 (Illumina).

Quality Control and Sequence Processing

Raw sequencing reads were assessed and processed using FASTP v.0.23.1³⁰, which trimmed low-quality bases (Phred score < 15), removed adapter sequences, and filtered reads shorter than 50 bp. Following quality control, the processed read dataset was de novo assembled using MEGAHIT v.1.2.9³¹ with default k-mer parameters (21, 29, 39, 59, and 79 bp). Assembly coverage depth was evaluated by aligning processed reads to assembled contigs using BWA-Aligner v.0.7.13³², and contig graphs were visualized in Bandage v.0.8.1³³ to assess assembly quality and identify circular or near-circular sequences characteristic of organellar genomes.

Mitochondrial Genome Identification

To selectively identify mitochondrial contigs, all assembled sequences were compared against the RefSeq Mitochondrial database (https://ftp.ncbi.nlm.nih.gov/refseq/release/mitochondrion/) using DIAMOND v.2.0.9.147 in BlastX mode with an e-value threshold of 10⁻⁵. Candidate mitochondrial contigs were then manually inspected within Geneious R11.1.5 (Biomatters Ltd.) to confirm identity and assess genome completeness and circularity.

Mitochondrial Genome Annotation

Protein-coding genes, ribosomal RNA genes, and transfer RNA genes were annotated using MITOS2 ³⁴, which integrates hidden Markov models and comparative genomics. Transfer RNA genes were independently verified using tRNAscan-SEv.1.21³⁵ with the vertebrate mitochondrial genetic code. Secondary structure predictions for tRNA genes were inferred using standard cloverleaf models. A circular map of the mitogenome was generated using Proksee³⁶, an online visualization tool for mitochondrial genomic data.

Phylogenetic and Divergence Time Analysis

Phylogenetic trees were constructed from 13 concatenated protein-coding genes (CYTB, ND1-ND6, COX1-COX3, ATP8, ATP6) from 30 accipitrid and related species. The optimal model of sequence evolution was identified as GTR + Γ + I using jModelTest 2³⁷. Maximum-Likelihood (ML) inference was performed using IQ-TREE v.2³⁸ with model selection and ultrafast bootstrap approximation (1,000 replicates). Bayesian analysis and divergence-time estimation were conducted using BEAST v.2.7.6³⁹ with two calibration points: Aquilinae crown age is estimated at 9.8–15.4 Ma and Aegypiinae at 7.4–12.6 Ma (95% HPD), based on Nagy and Tökölyi (2014). An uncorrelated lognormal relaxed molecular clock and Yule speciation process were employed with 50 million generations. MCMC convergence was assessed using Tracer v.1.7, confirming effective sample sizes (ESS) greater than 200 for all parameters. Final time-calibrated trees were summarized by discarding the initial 25% of trees (burn-in) and computing maximum clade credibility trees with node ages and 95% HPD intervals using TreeAnnotator⁴⁰. Trees were visualized in FigTree v.1.4.4⁴¹.

The 13 protein-coding genes collectively span 10,711 bp, representing approximately 59.6% of the total mitochondrial genome. This proportion is consistent with the average PCG length of 10,709.5 bp observed across 27 raptor species in the comparative analysis, indicating remarkable conservation of coding sequence length within Accipitridae.

The protein-coding genes exhibit a distinct nucleotide composition pattern that differs from the overall genome composition. The PCG nucleotide distribution is: T (24.1%), C (34.3%), A (27.6%), and G (14.0%), resulting in a GC content of 48.3% and AT content of 51.7%. Notably, the GC content in PCGs (48.3%) is slightly higher than the whole genome GC content (46.5%), reflecting the coding constraint and codon usage bias in functional genes^46,47. This elevated GC content in protein-coding regions is attributable to coding constraints that limit nucleotide variability at functionally critical positions, particularly at first and second codon positions where substitutions more frequently result in non-synonymous changes⁴⁸. Codon usage bias analysis in mitochondrial genomes has demonstrated that natural selection, rather than mutational pressure alone, predominantly shapes codon preference to optimize translational efficiency and maintain protein function⁴⁹. The balance between mutational pressure and purifying selection acts to preserve optimal codon usage patterns that facilitate efficient translation by the mitochondrial ribosomal machinery⁵⁰.​ The protein-coding genes display significant compositional asymmetry. The GC skew value of -0.40 indicates a strong bias toward cytosine over guanine on the coding strand, which is a characteristic feature of vertebrate mitochondrial genomes and reflects strand-specific mutational pressures during replication. This compositional asymmetry arises from the asymmetric replication mechanism of mitochondrial DNA, wherein the heavy (H) strand remains in a single-stranded state for extended periods during replication⁵¹⁵². During this single-stranded exposure, cytosine residues undergo spontaneous deamination at rates over 100-fold higher than in double-stranded DNA, converting C to uracil (subsequently repaired as thymine), thereby generating the characteristic excess of thymine on the H-strand and complementary guanine depletion⁵³. The magnitude of GC skew has been shown to correlate with the distance from replication origins, with maximal skew occurring at sites of longest single-strand exposure. The conserved negative GC skew observed across Accipitridae species reflects the shared ancestral replication mechanism and origin positioning typical of avian mitochondrial genomes^54,55.​ The AT skew of 0.06 is slightly positive, showing a modest bias toward adenine over thymine. These skew values are highly consistent with those observed in related raptor species, where GC skew ranges from − 0.39 to -0.44 and AT skew ranges from 0.05 to 0.09. This narrow range of variation in compositional skewness across Accipitridae indicates that the mutational and selective forces shaping mitochondrial genome composition have remained relatively constant throughout the evolutionary history of this family.​

When compared to 27 related Accipitridae species from GenBank, Nisaetus bartelsi shows PCG characteristics that fall well within the typical range for this taxonomic group (Table 2). The total PCG length of 10,711 bp is identical to the family average, and both nucleotide composition and skewness metrics align closely with congeneric and confamilial species, suggesting conservation of mitochondrial genome architecture across Accipitridae. This conservation of protein-coding gene characteristics across Accipitridae reflects strong purifying selection maintaining mitochondrial genome structure and function⁵⁴. The similarity in PCG length, nucleotide composition, and compositional skewness among accipitrid species indicates that the mitochondrial genome organization has remained stable throughout the diversification of this family, despite variation in body size, ecological niche, and geographic distribution. Such conservation suggests that the mitochondrial genome architecture represents an optimal configuration for the energetic demands and life-history characteristics shared by Accipitridae raptors⁵⁶.

The relative synonymous codon usage (RSCU) analysis of the 13 mitochondrial protein-coding genes in Nisaetus bartelsi reveals structured pattern of codon bias (Fig. 2). Codons with RSCU > 1 are preferentially used, whereas codons with RSCU < 1 are under-represented⁵⁷. In the Javan hawk-eagle mitogenome, most preferred codons end with A or U at the third position, and several codons show particularly strong over-representation, such as UUA and CUA for leucine, UCC for serine, CCC for proline, ACA for threonine, GCA for alanine, CAA for glutamine, AAC for asparagine, AAA for lysine, GAC for aspartate, CUG/GGC for glycine, and CGC/CGA for arginine (RSCU often > 1.5–2.5). In contrast, many G- or C-ending synonymous codons such as UUG (Leu), UCG (Ser), ACG (Thr), GAG (Glu), AAG (Lys), AGG/AGA (Arg), and GGG (Gly) are strongly suppressed (RSCU near 0–0.5), indicating a directional bias in synonymous codon usage. This preference for A/U-ending codons is consistent with the overall slight AT bias of the mitochondrial protein-coding region and mirrors patterns reported for many animal mitochondrial genomes, where RSCU analyses commonly show enrichment of codons ending in A or U at the third position⁵⁸. Such codon usage patterns are thought to arise from a combination of mutation pressure (favoring A/T over G/C at third codon positions) and natural selection on translational efficiency and accuracy, including adaptation to the cellular tRNA pool^59,60. Comparative studies across diverse taxa have demonstrated that mitochondrial codon usage bias is typically weak to moderate but non-random, with optimal codons often corresponding to those decoded by more abundant tRNA species and to genomes with similar base composition⁶¹, suggesting that both nucleotide composition constraints and selection shape RSCU profiles.​

In N. bartelsi, the strong over-representation of specific A/U-ending codons (e.g., UUA for Leu, UCC for Ser, AAA for Lys, CAA for Gln) and under-representation of their G/C-ending counterparts indicates that the codon bias is not purely neutral but likely reflects a balance between mutational tendencies and selective forces optimizing mitochondrial protein synthesis. Similar A/U-ending preferences and RSCU distributions have been observed in other avian mitogenomes, where codons such as CTA/ATA/ACA (Leu, Met, Thr) or CGA/CUA display high RSCU values, further supporting that the Javan hawk-eagle shares the general codon usage signatures of bird mitochondrial genomes. Overall, the RSCU pattern of N. bartelsi is typical for vertebrate mitochondria and provides additional evidence that its mitochondrial genome evolution is governed by compositional bias and selection acting on synonymous sites, in line with broader observations across animal mitogenomes

The mitochondrial genome of Nisaetus bartelsi contains 22 tRNA (Fig. 3) genes with a combined length of 1,519 bp, accounting for 8.45% of the entire mitogenome. These tRNA genes are discontinuously distributed throughout the whole mitogenome, interspersed between protein-coding genes and ribosomal RNA genes. The tRNAs show an asymmetric strand distribution, with 14 tRNAs encoded on the forward strand (+) and 8 tRNAs encoded on the reverse strand (-).

The size of individual tRNA genes ranges from 3 bp to 74 bp, with most tRNAs falling within the typical range of 66–73 bp. The longest tRNA genes are trnL2 and trnS2, each spanning 74 bp. All 22 tRNA genes exhibit the typical cloverleaf secondary structure characteristic of functional transfer RNAs in vertebrate mitochondrial genomes^62,63, as shown in the secondary structure prediction (Fig. 3). Each tRNA can be folded into the canonical structure consisting of the acceptor stem, D-arm, anticodon arm, variable loop, and TψC arm, which are essential for their aminoacylation and translation functions⁶².​

Two ribosomal RNA genes are present in the mitochondrial genome. The 12S rRNA (rrnS) gene is 971 bp in length and is located at positions 4391–5361 on the forward strand, flanked by trnF upstream and trnV downstream. The 16S rRNA (rrnL) gene is substantially larger at 1,567 bp and is positioned at locations 5455–7021, also on the forward strand, situated between trnV and trnL2. The total length of the two rRNA genes is 2,538 bp, representing 14.12% of the entire mitogenome. This proportion is consistent with the typical rRNA content observed in avian mitochondrial genomes, where ribosomal RNA genes constitute a substantial fraction of the non-coding functional elements essential for mitochondrial protein synthesis^64,65.​

The Nisaetus bartelsi mitochondrial genome exhibits the characteristic dual control region arrangement that is conserved in most Accipitridae mitogenomes. This arrangement consists of a functional control region (CR1) and a pseudo-control region (CR2), which originated through tandem duplication followed by gradual degeneration of the duplicated copy.​ CR1 (the functional control region) is located between trnT and trnP, with an estimated length of approximately 1,163 bp. This region contains the conserved regulatory elements essential for mitochondrial DNA replication and transcription.

Notably, CR2 harbors a 17× tandem repeat sequence (Fig. 4), representing a highly repetitive motif that may play roles in replication initiation and regulation of mtDNA copy number. Such tandem repeats are commonly found in accipitrid control regions and show high variation in repeat number among species^66–70.​ CR2 (the pseudo-control region) is positioned between trnE and trnF, spanning approximately 1,265 bp. This region represents a degenerate copy of the control region that has lost most conserved regulatory elements and is undergoing evolutionary decay. CR2 typically accumulates more nucleotide substitutions and contains variable numbers of tandem repeats (VNTRs) but lacks the functional domains present in CR1⁵⁴. The presence of CR2 is characteristic of the gene order rearrangement (GO-IV type) found in Accipitridae, distinguishing this family from other raptors with different levels of control region duplication or reduction^54,71,72.​ Together, the two control regions comprise approximately 2,428 bp, accounting for roughly 13.5% of the total mitochondrial genome. The dual control region structure has been hypothesized to be associated with longevity in birds⁷³, potentially by protecting cells from age-related mitochondrial deletions or by increasing the flexibility of mitochondrial responses to environmental changes.

Phylogenetic of Javan Hawk-Eagle ( Nisaetus bartelsi )

The time-calibrated phylogeny inferred from complete mitochondrial genomes of 30 taxa robustly resolves the position of the Javan hawk-eagle within Accipitridae. Nisaetus bartelsi is recovered as part of a well-supported Aquilinae clade together with Spizaetus nipalensis and Spizaetus alboniger, with both Bayesian posterior probabilities and maximum-likelihood (ML) bootstrap values at the relevant nodes equal to 1.0/100. The topology shows N. bartelsi as sister to S. alboniger, and this pair in turn sister to S. nipalensis, indicating that the three “hawk-eagle” taxa form a tightly clustered monophyletic group that has diverged relatively recently compared with deeper splits among major accipitrid subfamilies.

The time-calibrated phylogenetic tree provides divergence time estimates at each node, where the height parameter (with 95% highest posterior density; HPD- interval) represents the estimated age of that node in millions of years ago (Mya). The root height of approximately 30.1 (HPD:25.7–34.2) Mya, spanning the late Oligocene to early Miocene, indicates that the crown group comprising Pandion haliaetus, Elanus caeruleus, and Accipitridae originated near the Oligocene–Miocene boundary. This estimate is congruent with previous time-scaled phylogenies of Accipitriformes^74–76, which place the initial radiation of raptorial birds in the Oligocene followed by extensive diversification throughout the Miocene and Pliocene, and it overlaps the oldest fossil accipitrids documented from late Oligocene–early Miocene deposits in Europe, North America, and Australia^77–81.​

Within subfamilies, divergence times inferred from node heights and HPD intervals reveal a largely Miocene origin of the major accipitrid lineages, with subsequent Pliocene–Pleistocene species-level radiations. In Aquilinae, the crown group diverged mainly in the mid- to late Miocene (height 12.45;HPD 10.85–13.98 Mya), in line with multilocus analyses that date the origin of “true eagles” to this interval. The node uniting N. bartelsi with S. alboniger and S. nipalensis has a relatively young (3.36;2.83–3.92 Mya), implying a Pliocene–Pleistocene diversification of the hawk-eagle clade, whereas the split between this group and Hieraaetus/Aquila occurs in the middle Miocene (12.45; 10.85–13.98 Mya), suggesting a deeper separation between forest hawk-eagles and open-country eagles. This pattern is consistent with the fossil record of large aquiline eagles from the early–middle Miocene and with scenarios of Miocene habitat differentiation followed by Quaternary diversification in forested Southeast Asia^82–85.​

In Aegypiinae, the crown group diverged within the Miocene, with very short internal branches (11.1;9.51–12.68 Mya) among Gyps species indicative of late Miocene–Pliocene radiations of Old World vultures; these ages match late Miocene and Pliocene fossils attributed to gypaetine and aegypiine vultures from Europe and Africa. The Accipitrinae clade diverged during the Miocene (17.18;14.66–19.4 Mya) and exhibits short terminal branches among Accipiter species, suggesting rapid, recent diversification during the late Miocene–Pliocene that parallels their ecological and morphological radiation as forest hawks, as reported in previous multilocus and ultraconserved-element studies^86–88. Buteoninae and related kite lineages (Haliastur, Milvus) diverged in the mid-Miocene, with tip ages extending into the Pliocene (16.5;14.2–18.8 Mya), consistent with fossil buteonines from the middle–late Miocene and molecular estimates of a long buteonine history followed by more recent species-level radiations^79,82,89,90. Finally, the positions of E. caeruleus and P. haliaetus at or near the base of the tree, with large node heights in the late Oligocene–early Miocene (24.2;20.8–27. 5 Mya), confirm their status as early-diverging lineages relative to core Accipitridae and align with independent genomic studies^91–93 that treat them as separate subfamilies (Elaninae and Pandionidae) branching near the base of Accipitriformes.

Overall, the high posterior probabilities and ML bootstrap values across the tree, together with divergence times that are consistent with both previous molecular estimates and the fossil record, indicate that complete mitochondrial genomes provide a coherent and well-supported framework for understanding the temporal diversification of Accipitridae. Within this framework, Nisaetus bartelsi emerges as a very recent derivative within a young Pliocene–Pleistocene hawk-eagle clade nested in an Aquilinae lineage that originated during the Miocene radiation of large raptorial birds.​

The complete mitochondrial genome and divergence-time framework of Nisaetus bartelsi generated here provide a critical genetic baseline for conservation of this endemic and critically endangered Javan raptor. High-coverage mitogenomic data offer reliable markers for population monitoring, identification of illegal trade specimens, and assessment of genetic diversity and connectivity among fragmented populations, which are essential parameters for designing effective management units and translocation or captive-breeding programs. The placement of N. bartelsi as a recently diverged lineage within a young hawk-eagle clade further underscores its irreplaceable evolutionary distinctiveness and supports prioritizing its habitats for protection in regional conservation planning, in line with growing evidence that preserving phylogenetic diversity helps maintain ecosystem resilience and options for future adaptation under rapid environmental change