Morphological identification primarily followed the criteria in Fujita (2010) and Kawakami (2023). Specimens of H. davisoni, H. nigritus, H. sanuchi, H. howdeni, H. forsteri, H. mniszechi, H. bowringi, Weinreichius perroti Lacroix, 1978 and Yumikoi makii Arnaud and Miyashita, 2006 were derived from the second to fifth generations of artificial breeding, while the remaining specimens were collected from wild populations between 2020 and 2023 near their respective type localities. We confirm that the specimens used in this study were obtained ethically following all applicable international and local permitting requirements.

The specimens analyzed in this study were dried and preserved in the insect collection of Nanjing Normal University, Nanjing, China. Genomic DNA was extracted from the prothoracic muscle tissue of each species using the TIANamp Genomic DNA Kit (TIANGEN, Beijing, China). DNA concentrations were measured using a Nanodrop 2000 spectrophotometer.

All species’ genomic DNA was sequenced utilizing next-generation sequencing (NGS) on the Illumina NovaSeq platform. Each sample’s library was prepared using the TruSeq™ DNA Sample Prep Kit (insert size 400 bp), and all libraries were sequenced in PE150 mode (paired-end, 2 × 150 bp). Approximately 4 GB of raw data were generated for each sample. The raw data underwent quality filtering with Fastp (Chen et al. 2018). MitoZ was used for the assembly of mitogenomes (Meng et al. 2019), with the kmers-megahit = 79, 99, 119, 141. To verify the accuracy of NGS sequencing results, COI gene fragments of all samples were sequenced using Sanger sequencing with a pair of universal primers (LCO1490 as a forward primer, 5’-GGTCAACAAATCATAAAGATATTGG-3’; HCO2198 as a reverse primer, 5’-TAAACTTCAGGGTGACCAAAAAATCA-3’) (Folmer et al. 1994). The r-Taq polymerase (Takara, Beijing, China) was used for the polymerase chain reaction (PCR) based on the method of Yanai et al. (2017).

The annotation of mitogenomes was performed using the MITOS WebServer (Bernt et al. 2013) and MitoZ (Meng et al. 2019). These annotations were subsequently verified through alignments with available Hexarthrius sequences from GenBank, utilizing ClustalW within MEGA v11 (Tamura et al. 2021). The secondary structures of tRNAs were analyzed using tRNAscanSE v2.0 (Lowe and Chan 2016) in conjunction with the MITOS WebServer. Circular maps of the mitogenomes were generated using the visualization module of MitoZ (Meng et al. 2019).

The nucleotide composition was computed using MEGA v11, while the AT-skew and GC-skew were derived from the formulas: AT-skew = (A–T)/(A+T) and GC-skew = (G–C)/(G+C) (Perna and Kocher 1995). A sliding window analysis with a window size of 200 bp and a step size of 20 bp was conducted to estimate nucleotide diversity (Pi) across the 13 protein-coding genes (PCGs) using DnaSP v6 (Rozas et al. 2017). Genetic distances of COI genes are frequently utilized for insect species identification (Cox et al. 2013; Zheng et al. 2022), which were calculated based on the Kimura 2-parameter model (Kimura 1980) using MEGA v11. The ratio of nonsynonymous to synonymous rates (Ka/Ks) for the 13 PCGs was also determined using DnaSP v6.

A total of 28 mitogenomes from eight genera were included in the phylogenetic reconstruction, comprising 18 sequences from 14 Hexarthrius and 9 sequences from six related genera (Huang and Chen 2013; Kim and Farrell 2015). The mitogenome of Dorcus hopei (Saunders, 1854) served as the outgroup (Table S1). Of the 28 mitogenomes analyzed, 23 were sequenced in this study, while the remaining five were sourced from GenBank (Table S1).

The nucleotide sequences of the 13 PCGs were aligned with MAFFT L-INS-i strategy (Katoh and Standley 2013) in PhyloSuite v1.2.2 (Zhang et al. 2020). Ambiguous alignment sites were removed with Gblocks v0.91b (Talavera and Castresana 2007). Individual genes were then concatenated within PhyloSuite v1.2.2. Two datasets were created for phylogenetic reconstruction: (1) a PCGs matrix including all codon positions of the 13 PCGs, and (2) a PR matrix comprising the 13 PCGs and two rRNA genes. Phylogenetic reconstructions were performed using Bayesian Inference (BI) and Maximum Likelihood (ML) methods based on both datasets. The optimal partitioning model for each dataset was selected via ModelFinder using Bayesian Information Criterion (BIC) and Akaike Information Criterion corrected (AICc) (Kalyaanamoorthy et al. 2017). BI trees were constructed using MrBayes v3.2.6 through the online CIPRES Science Gateway (Miller et al. 2011; Ronquist et al. 2012), with settings for two parallel runs of four Markov chains over 10 million generations (sampling every 1,000 generations), discarding 25% as burn-in. ML analyses were conducted with RAxML v8.2.0, employing the GTRGAMMAI model and 1,000 bootstrap replicates (Stamatakis 2014). Phylogenetic trees were edited using FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree).

Divergence times among species were estimated using BEAST v2.6.3 (Bouckaert et al. 2014), applying the PCG matrix and BI phylogenetic topology derived from previous analyses. In BEAUti, the same partitioning and site models established in the phylogenetic analysis were utilized. The root of the tree was calibrated at the divergence point between D. hopei and the other species, following the findings of Kim and Farrell (2015). The molecular dating indicated a root divergence time of 62.4 Ma and a split between Hexarthrius and Prosopocoilus Hope & Westwood, 1845 at 40.2 Ma, based on fossil evidence and prior phylogenetic analyses (Kim and Farrell 2015). A birth-death tree prior and a fixed substitution rate of 0.0115 (corresponding to an arthropod mitochondrial molecular clock of 2.3% per million years) (Brower 1994) were implemented. Two independent MCMC runs were executed, each consisting of 200 million generations with sampling every 1,000 generations. Tracer v1.5 (Drummond and Rambaut 2007) was employed to examine posterior estimates and ensure that the effective sample size (ESS) exceeded 200 for critical parameters. The tree files from both runs were combined using LogCombiner v1.5.3, with the first 10% of generations discarded as burn-in. Based on the combined tree file, TreeAnnotator v1.5.3 (Drummond and Rambaut 2007) was used to calculate the consensus tree and annotate divergence times. The final topology was visualized and refined using FigTree v1.4.4.

Ancestral area reconstruction was conducted using RASP (Reconstruct Ancestral State in Phylogenies) v4.0 (Yu et al. 2020) with the BioGeoBEARS package (Matzke 2014). Each Hexarthrius species was considered endemic, given their narrow and allopatric distributions. Six areas of endemism were defined: (A) the Western Ghats of South India, represented by H. davisoni; (B) a triangular region defined by the eastern half of the Himalayas to the north, the Hengduan Mountains of China to the east, and the Myanmar Shan Plateau to the west, encompassing seven Hexarthrius species (H. parryi, H. forsteri, H. vitalisi, H. melchioritis, H. mniszechi, H. aduncus, and H. bowringi); (C) the eastern region of area B, comprising southern China and northern Indo-China (north of 15° N), which includes H. parryi and H. vitalisi; (D) the southern part of the Indo-China Peninsula (south of 15° N), the Malay Peninsula, and the islands of Sumatra and Java, featuring six species (H. parryi, H. sanuchi, H. nigritus, H. buquetti, H. rhinoceros, and H. mandibularis); (E) Borneo, with two species, H. parryi and H. mandibularis; and (F) Mindanao and Luzon Islands of the Philippines, represented by H. howdeni. The specific distributions of ingroup species were inputted, and the BI-PCG tree topology from earlier analyses was applied in the reconstruction. The DIVALIKE model was identified as the most suitable, with a maximum limit of six ancestral areas.

Ancestral state reconstruction (ASR) was performed for eight morphologically significant characters: (1) male clypeolabrum; (2) granulation on the male mandible; (3) basal teeth of the male mandible; (4) denticles between the male mandibular base and major tooth; (5) male antennal club; (6) male dorsal head protuberance; (7) male elytra coloration; and (8) female head structure. These reconstructions were mapped onto the BI tree derived from the PCG matrix, employing parsimony methods in Mesquite v3.81 (Maddison and Maddison 2023). Character histories were traced utilizing unordered state transformations.

The 18 complete mitogenomes of Hexarthrius range from 17,435 bp in H. mandibularis to 19,176 bp in H. rhinoceros rhinoceros (Figs 1A–F, S1–S3). Among the Hexarthrius clade, which includes Hexarthrius and five related genera (Rhaetulus, Pseudorhaetus, Weinreichius, Yumikoi, and Rhaetus), Rhaetulus crenatus rubrifemoratus Nagai, 2000 possesses the largest mitogenome (19,236 bp) (Fig. 1C). Twenty-one of the 23 new mitogenomes exhibit the ancestral insect mitochondrial organization of 13 PCGs, 22 tRNAs, two rRNAs and one control region, with 23 genes encoded on the J-strand and 14 on the N-strand. Two previously unknown gene rearrangements were identified: an extra trnS (UGA) in H. mandibularis and a complex trnI-NCR-trnQ-NCR-trnQ-NCR-trnM cluster in R. crenatus rubrifemoratus (Fig. 1A, C).

The A+T content of the Hexarthrius clade mitogenomes varies from 65.5% in H. mandibularis to 69.5% in R. crenatus, with H. mandibularis and H. howdeni showing relatively reduced A+T bias (Fig. S4A). The 23 new mitogenomes displayed either positive AT-skews and negative GC-skews, or the reverse. Closely related subspecies exhibited similar base compositions, except H. parryi deyrollei Parry, 1864 and H. parryi paradoxus Möllenkamp, 1898, which differed in AT-skew direction.

The 13 PCGs showed the ratio of non-synonymous (Ka) to synonymous (Ks) substitution rates (Ka/Ks) ranging from 0.015 (COI) to 0.180 (ATP8), confirming strong purifying selection across all genes (Fig. S4B). COI was the most conservative gene, showing utility for species identification. Nucleotide diversity (Pi) varied from 0.116 (ND1) to 0.195 (ND6), with ND6 showing the highest variability (Fig. S4C).

We calculated the genetic distances of the complete COI genes from 28 mitogenomes (Supporting Information, Table S2). Within Hexarthrius, pairwise distances varied from 0.005 (H. nigritus–H. sanuchi) to 0.176 (H. forsteri kiyotamii–H. buquetti), averaging 0.139 (Supporting Information, Table S2). The unusually low value between H. nigritus and H. sanuchi and the high intersubspecific divergence (0.102) between H. vitalisi vitalisi Didier, 1925 and H. vitalisi tsukamotoi Nagai, 1998 warrant taxonomic reassessment.

Both BI and ML analyses based on both datasets (PCGs and PR) generated congruent topologies with strong support for the monophyly of Hexarthrius (Fig. 2). Two well-supported clades were recovered (Figs 3, 4): a Himalayan clade (H. forsteri, H. vitalisi, H. melchioritis, H. mniszechi, H. aduncus, H. bowringi) and a Tropical clade (H. davisoni, H. parryi, H. sanuchi, H. nigritus, H. buquetti, H. rhinoceros, H. mandibularis, H. howdeni).

Rhaetulus formed the earliest branch within a larger “Hexarthrius clade”, comprising Rhaetulus + ((Pseudorhaetus + (Weinreichius+Yumikoi)) + (Rhaetus+Hexarthrius)). Prosopocoilus was supported as the sister group to this clade (Fig. 2).

Ancestral area reconstruction indicated that the ancestral range of Hexarthrius+Rhaetus was centered in the eastern Himalayas and Hengduan Mountains of China (Endemism B) (Figs 3, 4). This region also represents the ancestral range for the broader clade comprising (Hexarthrius+Rhaetus) and (Pseudorhaetus + (Weinreichius+Yumikoi)), as well as all nodes within the Himalayan clade (Fig. 4). In contrast, Endemism D was inferred as the ancestral area for Pseudorhaetus, Weinreichius, Yumikoi, and most members of the Tropical clade, except for H. mandibularis and H. howdeni, whose ancestral range was estimated as D+F (Fig. 4).

Across the Hexarthrius clade, BioGeoBEARS identified five vicariance, eight dispersal, and one extinction event (the latter at the divergence of H. davisoni and its sister group) (Fig. 4).

Divergence between the Himalayan and Tropical clades was dated to approximately 8.9 million years ago (Ma) (95% highest posterior density [HPD]: 3.50–19.30 Ma) (Fig. 4). Within the Himalayan clade, H. forsteri diverged first (~5.23 Ma; 95% HPD: 1.28–12.62 Ma), followed by H. vitalisi (~3.64 Ma; 95% HPD: 1.09–7.52 Ma). The H. melchioritis + H. mniszechi and H. aduncus + H. bowringi clades separated around 2.48 Ma (95% HPD: 0.68–5.77 Ma).

In the Tropical clade, H. davisoni diverged earliest (~5.23 Ma; 95% HPD: 1.92–11.47 Ma). The H. mandibularis + H. howdeni clade split around 3.47 Ma (95% HPD: 1.09–6.81 Ma), while the H. parryi + (H. sanuchi + H. nigritus) and H. buquetti + H. rhinoceros clades diverged at ~2.56 Ma (95% HPD: 0.64–5.25 Ma).

The Hexarthrius clade diverged from Prosopocoilus approximately 39.67 Ma (95% HPD: 33.55–45.93 Ma). Within the Hexarthrius clade, Rhaetulus branched off at ~27.53 Ma (95% HPD: 19.74–37.05 Ma); Pseudorhaetus diverged from the (Weinreichius+Yumikoi) clade at ~7.77 Ma (95% HPD: 0.72–19.22 Ma); Weinreichius and Yumikoi separated around 1.05 Ma (95% HPD: 0.11–5.83 Ma). The sister taxa Rhaetus and Hexarthrius diverged around 23.77 Ma (95% HPD: 18.66–28.78 Ma).

Prior to this study, gene rearrangements in Lucanidae were sparsely reported, with only one known case of the translocation of trnL (UUR) into the control region in Sinodendron yunnanense Král, 1994 (Lin et al. 2017). Our findings expand that limited record by revealing two previously undescribed rearrangements in Lucanidae: an extra copy of trnS (UGA) inserted in the control region of H. mandibularis and a complex trnI-NCR-trnQ-NCR-trnQ-NCR-trnM cluster in R. crenatus rubrifemoratus (Fig. 1A, C). Gene rearrangements of this type are uncommon in Coleoptera (Timmermans and Vogler 2012; Li et al. 2016; Lin et al. 2017; Jeong et al. 2020; Ge et al. 2022; Yang et al. 2023) and the duplications reported here are unique among beetles.

Among the major models proposed for mitochondrial gene rearrangements, i.e., recombination, tandem replication with non-random loss (TDNL), illicit priming by tRNAs, and tandem duplication-random loss (TDRL), the TDRL mechanism best fits the patterns observed in R. crenatus rubrifemoratus (Yang et al. 2023). The presence of NCR between trnI and trnQ in Prosopocoilus (noted previously as synapomorphic for that genus) and Rhaetulus and its absence in other genera further suggests a sequence of duplication and differential loss across the clade (Kim and Farrell 2015). Our phylogenetic analysis confirmed Prosopocoilus as the sister group to the Hexarthrius clade and Rhaetulus as the earliest offshoot within this clade, which aligns with the scenario of TDRL (Fig. 2). We therefore propose that an ancestral insertion of an NCR between trnI and trnQ preceded subsequent lineage-specific TDRL events: retention of a simplified trnI–NCR–trnQ arrangement in Prosopocoilus, triplication of an NCR–trnQ block followed by partial loss in the ancestor of Rhaetulus (yielding the observed trnI–NCR–trnQ–NCR–trnQ–NCR–trnM cluster), and independent loss of redundant NCRs/trnQ copies in the remaining genera, ultimately restoring the canonical trnI–trnQ–trnM order in many lineages (Fig. 2). This sequence of events explains both the distribution of NCRs in the group and the phylogenetic placement of genera recovered by our mitogenomic analyses.

An integrated reading of the phylogeny, divergence times and ancestral-area reconstructions indicates that the core Hexarthrius lineage and its close relatives originated in the eastern Himalaya–Hengduan region (Endemism B). We infer that the early divergence of Rhaetulus from the remainder of the clade was driven largely by vicariance associated with the complex topography of the Hengduan–Himalayan orogen; uplift and attendant habitat fragmentation would have promoted long-term isolation between north-south faunal elements (Xing and Ree 2017) (Fig. 4). The split between the (Pseudorhaetus + (Weinreichius+Yumikoi)) assemblage from the (Hexarthrius+Rhaetus) lineage is best interpreted as an initial dispersal event across the mountain complex followed by vicariance between eastern and western slopes.

Within Hexarthrius itself, the phylogeny resolves a Himalayan clade that remained largely restricted to montane areas and a Tropical clade that expanded southward into Indochina and Sundaland. The divergence between these clades (~8.9 Ma) and subsequent species-level splits (late Miocene–Pleistocene) coincide with episodes of regional orogeny and climatic change that alternately connected and fragmented forest habitats. We infer a pattern in which ancestral populations dispersed from the Himalayan–Hengduan core into lower-elevation corridors, with some lineages colonizing southern India and Sundaland; local extinctions in the ancestral area (Endemism B) and repeated Pleistocene sea-level fluctuations then promoted insular isolation and rapid differentiation in the Tropical clade (Hanebuth et al. 2011; He et al. 2023). The distribution of H. howdeni (Philippines) and its sister taxon H. mandibularis (Borneo) is consistent with cross-sea dispersal and subsequent island isolation; the present geological configuration of the Philippines predates these species’ divergence, indicating that vicariance through recent tectonics is unlikely to be the sole driver of their separation (Morrison 2014) (Fig. 4). Overall, our results support a two-phase scenario of mountain-origin Vicariance followed by lowland/insular dispersal and island-driven speciation.

Our phylogenetic results provide a framework to interpret the evolution of adult morphology in Hexarthrius (Figs S5, S6). Morphological distinctions between males of the Himalayan and Tropical clades are notable, including: (1) mandibular granules are small or absent in the Himalayan clade but enlarged in the Tropical clade (except in H. davisoni); (2) denticle number along the inner mandibular margin is usually < 5 in the Himalayan clade and > 5 in the Tropical clade (excluding H. davisoni, H. howdeni, and H. forsteri); and (3) the clypeolabrum is strongly protruding in the Tropical clade but smooth in the Himalayan clade (with exceptions in H. mandibularis and H. howdeni). However, these characters frequently conflict with phylogenetic relationships. For example, H. davisoni, although belonging to the Tropical clade, resembles Himalayan species in mandible shape and denticle patterns, and H. aduncus is difficult to distinguish morphologically from H. vitalisi, despite being sister to H. bowringi, a lineage whose mandibular and cephalic morphology differs substantially (Figs 3, S5, S6).

The evolution of antennal club segmentation also reflects phylogenetic history. A five-segmented antennal club is inferred as the ancestral state (Figs S5, S6). The unique development of a sixth segment in H. forsteri, caused by dense setation on the basal segment, represents an autapomorphy. Expansion of the antennal club is known to evolve independently in other lucanid genera (e.g., Prosopocoilus) and appears unrelated to phylogenetic structure within Hexarthrius.

Mandibular dentition provides additional insight (Figs S5, S6). The ancestral state for the genus consists of two distinct basal teeth (one dorsal and one ventral). This morphology is retained in Rhaetulus, Pseudorhaetus, Rhaetus, and in H. forsteri and H. parryi + (H. sanuchi + H. nigritus). In other lineages, the distal basal tooth has been modified. In the lineage comprising (H. aduncus + H. bowringi) + (H. melchioritis + H. mniszechi) + H. vitalisi, the teeth evolved into a single, medially directed large tooth, sometimes bifurcated. In H. melchioritis and H. mniszechi, subsequent reduction led to the presence or absence of only a ventral tooth. In the Tropical clade, H. davisoni independently evolved a similar medially pointed tooth, whereas other Tropical clade species exhibit a continuous row of small denticles replacing the dorsal tooth.

The denticles between the basal teeth and the main tooth represent another primitive character shared across Hexarthrius and related genera (Figs S5, S6). Denticle reduction occurred convergently in the Himalayan clade and in H. davisoni. In the ancestors of the clade (H. aduncus + H. bowringi) + (H. melchioritis + H. mniszechi) + H. vitalisi, denticles were reduced to a single remnant or eliminated, but occasional individuals of H. aduncus and H. bowringi exhibit partial reappearance of multiple denticles, suggesting reversibility. The absence of denticles in H. howdeni is interpreted as a consequence of miniaturization, as small males across species exhibit reduced dentition.

Other male secondary sexual traits show similar evolutionary trajectories. Both large mandibular granules and dorsal protuberances of the head originated in the clade comprising (H. mandibularis + H. howdeni) + ((H. buquetti + H. rhinoceros) + (H. parryi + (H. sanuchi + H. nigritus))). These characters were subsequently lost in H. howdeni, again likely as a result of miniaturization. The protruding clypeolabrum occurs in ((H. buquetti + H. rhinoceros) + (H. parryi + (H. sanuchi + H. nigritus))) and in the monotypic genus Weinreichius, but ancestral state reconstruction indicates independent origins.

Male elytral coloration exhibits evolutionary lability (Figs S5, S6). Orange elytra occur sporadically within Hexarthrius, being fixed in H. parryi and Y. makii, but geographically variable in H. sanuchi, H. vitalisi, and R. crenatus (Kawakami 2023; Kawakami 2024). Within H. sanuchi, multiple color forms occur within the same locality, indicating polymorphism rather than geographic structuring (Kawakami 2023).

Body size evolution in the H. mandibularis + H. howdeni lineage is striking. Hexarthrius mandibularis represents the largest Hexarthrius species (up to 119.5 mm), whereas H. howdeni is the smallest (max. 63 mm) (Kawakami 2023). Despite this extreme difference, genetic distance (0.064) and divergence times indicate a recent shared origin (Figs 3, 4; Supporting Information, Table S2). This rapid size divergence likely reflects island-driven evolutionary mechanisms and highlights a broader issue: reliance on characters of large males in previous lucanid taxonomy may have obscured true evolutionary relationships, particularly when large males are absent or rare.

In contrast to male traits, only female head punctation provides phylogenetically consistent signal. Females of the Tropical clade possess coarse, dense punctures, whereas Himalayan clade females have smooth to weakly punctate heads (Fig. S6). Our analysis suggests that the coarse head morphology likely represents the ancestral state for the Hexarthrius clade, while the smooth head evolved independently in Rhaetus and the Himalayan clade of Hexarthrius (Fig. S5H). The intermediate condition in H. vitalisi and H. forsteri reflects transitional morphological evolution, matching their phylogenetic placement. This result indicates that female morphology, despite being overlooked historically, may be more informative for phylogeny than highly variable male characters.

Species boundaries in stag beetles remain problematic (Huang and Chen 2013). Our COI K2P distance matrix (Supporting Information, Table S2) identified two anomalies: an unusually low distance of 0.005 between H. sanuchi and H. nigritus, and a deep divergence of 0.102 between H. vitalisi vitalisi (collected from Yunnan) and H. vitalisi tsukamotoi (Guangxi). The near-identity of COI sequences in H. sanuchi and H. nigritus, coupled with their overlapping morphological variation and unstable diagnostic characters (Araya 2014; Kawakami 2023), strongly suggests conspecificity or recent gene flow. However, because our samples include multi-generation, captive-bred material and we lacked direct examination of type specimens, we refrain from proposing formal synonymy pending broader sampling and examination of types.

Conversely, the deep COI divergence among H. vitalisi subspecies indicates cryptic diversity within what is currently treated as a single species (Supporting Information, Table S2). Morphological variation within H. vitalisi is geographically structured: some populations (e.g. H. vitalisi tsukamotoi and an unnamed central Vietnam group) share orange elytra and similarly shaped mandibles, whereas others (e.g. H. vitalisi vitalisi and an unnamed southern China group) display dark brown elytra and distinct pronotal and mandibular morphologies. These patterns, together with geographically disjunct distributions, imply multiple independently evolving lineages that deserve taxonomic reassessment. Until additional specimens, type-material comparisons and nuclear loci are examined, we conservatively retain current subspecific assignments but highlight these cases as priorities for taxonomic revision.

Across the genus Hexarthrius, the average interspecific K2P distance is 0.14, with most interspecific values above 0.05 and intersubspecific values below 0.05 (Supporting Information, Table S2). While this empirical threshold is useful as a preliminary guide, exceptions documented here (notably H. sanuchi / H. nigritus and the H. vitalisi complex) underscore the need for integrative approaches combining mitogenomes, nuclear markers, morphology (including females), and geographic sampling to delimit species robustly in Hexarthrius and allied genera such as Prosopocoilus.