Phylogenetic Implication of Large Intergenic Spacers: Insights from a Mitogenomic Comparison of Prosopocoilus Stag Beetles (Coleoptera: Lucanidae)

Simple Summary Insect mitochondrial genomes (mitogenomes) show high diversity in some lineages. In the mitogenome of some Coleoptera species, a large intergenic spacer (IGS) has been identified. However, very little is known about mitogenomes of lucanid beetles. In this work, to provide further insight into the phylogenic relationships among species in lucanid beetles (genus Prosopocoilus), two Prosopocoilus species (Prosopocoilus castaneus and Prosopocoilus laterotarsus) were newly sequenced and comparatively analyzed. Significantly, the two newly sequenced Prosopocoilus species contained a large IGS located between trnI and trnQ. Our phylogenomic analyses showed that P. castaneus and P. laterotarsus were clustered in a clade with typical Prosopocoilus species (Prosopocoilus confucius, Prosopocoilus blanchardi, and Prosopocoilus astacoides). These results provide valuable data for the future study of the phylogenetic relationships in this genus. Abstract To explore the characteristics of mitogenomes and discuss the phylogenetic relationships within the genus Prosopocoilus, the mitogenomes of two species (P. castaneus and P. laterotarsus) were newly sequenced and comparatively analyzed. The arrangement of the mitogenome in these two lucanid beetles was the same as that in the inferred ancestral insect, and the nucleotide composition was highly biased towards A + T as in other lucanids. The evolutionary rates of 13 protein-coding genes (PCGs) suggested that their evolution was based on purifying selection. Notably, we found evidence of the presence of a large IGS between trnI and trnQ genes, whose length varied from 375 bp (in P. castaneus) to 158 bp (in P. laterotarsus). Within the large IGS region, a short sequence (TAAAA) was found to be unique among these two species, providing insights into phylogenomic reconstruction. Phylogenetic analyses were performed using the maximum likelihood (IQ-TREE) and Bayesian (PhyloBayes) methods based on 13 protein-coding genes (PCGs) in nucleotides and amino acids (AA) from published mitogenomes (n = 29). The genus Prosopocoilus was found to constitute a distinct clade with high nodal support. Overall, our findings suggested that analysis of the characteristics of the large IGS (presence or absence, size, and location) in mitogenomes of the genus Prosopocoilus may be informative for the phylogenetic and taxonomic analyses and for evaluation of the genus Prosopocoilus, despite the dense sampling materials needed.


Mitogenome Assembly, Annotation, and Analysis
High-quality reads were applied to de novo assemble using IDBA-UD [38] with minimum and maximum k values of 80 and 240 bp, respectively. We then evaluated the accuracy of the assembly by re-mapping clean reads to the mitogenome assemblies in Geneious v6.1.7 (Biomatters Ltd., Auckland, New Zealand), allowing mismatches of up to 2%, a minimum overlap of 100 bp, and a maximum gap size of 3 bp. The preliminary annotations were made on the MITOS Web Server with the invertebrate mitochondrial code (http://mitos.bioinf.unileipzig.de/inxdex.py, accessed on 12 September 2021). Proteincoding regions were identified by aligning published mitochondrial sequences in Geneious v6.1.7 (Biomatters Ltd.), whereas tRNA genes and their secondary structures were inferred using tRNA scan-SE v2.0 [39]. Owing to they could not be identified using tRNA scan-SE, 16S ribosomal RNA (rrnL) and 12S ribosomal RNA (rrnS) were determined according to sequence similarity with related species. Relative synonymous codon usage (RSCU) and Nucleotide composition were computed using MEGA-X [40]. Composition skew analysis was performed using the following formulas: AT skew = (A -T)/(A + T) and GC skew = (G -C)/(G + C) [41]. The evolutionary rates (Ka/Ks ratios) for each PCG were calculated using DnaSP v5.0 [42]. The mitogenomic map was drawn using the CG View server V 1.0 [43].

Phylogenetic Analyses
Phylogenetic analyses were performed using five Prosopocoilus stag beetles, twentyone other stag beetles as an ingroup, and three scarab beetles (Cheirotonus jansoni, Protaetia brevitarsis, and Rhopaea magnicornis) as an outgroup (Table 2). Protein coding genes (PCGs) of these species were extracted based on GenBank annotations using GenScalpel [44]. The nucleotide and amino acid sequences of the 13 PCGs were used as the dataset to construct the Bayesian inference (BI) and maximum likelihood (ML) phylogenetic trees, respectively. All these nucleotide and amino acid sequences were aligned using translation align and geneious align, respectively, in the programme Geneious 9.0.5. Gaps and ambiguous sites were removed from the protein alignment to generate a 10,878 bp nucleotide dataset and a corresponding amino acid dataset (3626 amino acids). We then set the model selections as Akaike information criterion (AIC), greedy search algorithm, and unlinked branch lengths to estimate the best fitting schemes in the programme PartitionFinder 2.1.1 [45]. The best-fit partitioning schemes and evolutionary models for the nucleotide and amino acid datasets are presented for ML analyses (Table S1).
BI and ML analyses were conducted using PhyloBayes MPI 1.5a [46] and IQ-TREE web server [47], respectively. Then, BI analysis was performed on the CIPRES Science Gateway [48] and the site-heterogeneous mixture model (GTR + CAT) was chosen [49]. Two independent chains starting from a random tree were run for 20,000 cycles, with trees being sampled every 10 cycles. The initial 25% of trees for each MCMC run were discarded as burn-in. A consensus tree was computed from the remaining 1500 trees combined from two runs, and the two runs converged at a maxdiff of less than 0.1. In ML analyses, the "Auto" option was set under optimal evolutionary models, and the phylogenomic trees were constructed using an ultrafast bootstrap approximation approach with 10,000 replicates. Phylogenomic trees were viewed and rooted with the three species in Scarabaeidae as the outgroup in Figtree v1.4.3 (http://beast.bio.ed.ac.uk/figtree, accessed on 1 June 2022).

Genome Organization and Base Composition
In this study, we obtained one complete mitogenome of P. castaneus and one nearly complete mitogenome of P. laterotarsus. The total lengths of the two mitogenomes ranged from 17,523 bp (P. castaneus) to 17,333 bp (P. laterotarsus). Both new mitogenomes shared the ancestral type for insects [55,56]; thus, these mitogenomes were composed of 13 PCGs, 22 tRNAs, two rRNAs, and a control region. We recovered only a partial control region of the mitogenome of P. laterotarsus (Table 3). Of the 37 genes, nine PCGs and 14 tRNA genes were located on the majority strand (J-strand), with the remaining four PCGs, two rRNA genes, and eight tRNAs genes on the minority strand (N-strand) ( Figure 1). The mitochondrial genome map of Prosopocoilus species is shown in Figure 1. Both mitochondrial genomic arrangements were similar to other published stag beetle mitogenomes [4,20,23,50]. All 22 tRNA genes were detected in these two newly sequenced mitogenomes (Table 3, Figure 1). The lengths of the tRNA genes ranged from 60 to 70 bp ( Table 3). All tRNA genes displayed the representative clover-leaf secondary structure (Figures S1 and S2), whereas trnS1 (AGN) lacked the dihydrouridine (DHU) arm, and it was replaced with a simple loop, which was common in other lucanid mitogenomes [4,23,[50][51][52][53][54]. There is a similar nucleotide proportion among the two newly sequenced mitogenomes, with a high content of AT ( Table 4). The AT skews of PCGs, tRNA genes, and rRNA genes in the P. castaneus and P. laterotarsus mitogenomes were −0.16/−0.18, 0.05/0.05 and −0.11/−0.08, respectively ( Table 4). The AT-skew values were negative in PCGs but positive in the tRNA genes within the two mitogenomes, consistent with those found in other lucanid mitogenomes [4,20,23,50]. Additionally, the nucleotide compositions of lucanid mitochondrial genomes corresponded well to the AT skew generally observed in other beetles [8,[11][12][13].

Protein-Coding Genes and Codon Usage
The lengths, nucleotide proportion, and codon usages of 13 PCGs in these two new mitogenomes were nearly the same as those in the ancestral type of insects (Tables 3 and 4; Figures 1 and 2). All 13 PCGs were identified in these two new mitogenomes. Twelve PCGs of these two new mitogenomes used ATN (where N represents A, C, G, or T) as initiation codons, with the exception of cox1, which was initiated with AAN. Notably, AAN is an accepted conventional start codon in many beetle mitogenomes [8,11,15,23,57]. In the P. castaneus and P. laterotarsus mitogenomes, six PCGs shared the typical stop codons TAA and TAG, whereas in the remaining genes, an incomplete stop codon consisting of T or TA was inferred ( Table 3). All of the new mitogenomes had incomplete stop codons, as described in other stag beetles [4,23,51] and other insects [58,59], which have been demonstrated to produce functional stop codons in polycistronic transcription cleavage and polyadenylation processes [8,13,14]. Comparisons of five mitogenomes of Prosopocoilus stag beetles showed that the cox1 genes of P. castaneus, P. laterotarsus, P. blanchardi, P. astacoides, and P. confucius shared the same incomplete stop codon (T). By contrast, complete terminators (TAA) are utilized in Kirchnerius guangxii, Epidorcus gracilis, and the related species Serrognathus platymelus. Thus, we assumed that these species shared a similar preference for stop codon adoption and may have a closer relationship, as confirmed by the phylogenetic analysis (see below). The average AT contents of the 13 PCGs were 67.32% and 67.63% in P. castaneus and P. laterotarsus, respectively. The PCGs encoded by the J-strand displayed T-skews (T > A) and G-skews (G > C), whereas others encoded by the N-strand displayed T-skews and C-skews (C > G). The characteristics of the relative synonymous codon usages (RSCU) in these two new mitogenomes showed that codons including A or T at the third position were overrepresented compared with the other synonymous codons (Figure 2), reflecting the nucleotide bias of insect mitogenomes [8]. There was also a high AT content at the third codon site, indicating a high background mutational pressure towards AT nucleotides [60].

Evolutionary Rates of PCGs
Five available mitogenomes from the genus Prosopocoilus (from P. castaneus, P. laterotarsus, P. confucius, P. astacoides, and P. blanchardi) were used to calculate the evolutionary rate of its PCGs. For each PCG, the ratio of non-synonymous substitution (Ka) to synonymous substitution (Ks) was calculated (Figure 3). The NADH dehydrogenase subunits (nad1-6 and 4l) and ATP synthase subunits (atp8 and atp6) had higher Ka/Ks ratios than the cytochrome oxidase subunits (cox1, cox2, and cox3) and cytochrome b (cytb). This phenomenon suggested that various functional genes in the mitochondria of Prosopocoilus have undergone different selection pressures during the evolutionary process. The order of Ka/Ks for 13 protein genes was as follows: atp8 > nad6 > nad5 > nad4 > nad4l > nad2 > nad3 > nad1 > atp6 > cox2 > cytb > cox3 > cox1. Notably, the Ka/Ks ratios of the 13 PCGs among the five mitogenomes of Prosopocoilus were all less than 1.0, suggesting that these functional genes were all under strong purifying selection, as was reported in other insects [8,14,60]. The slowest and fastest evolution rates were observed for cox1 and atp8 genes, respectively  Figure 3). Furthermore, the cox1 had the smallest evolutionary rate, indicating that positive selection was less powerful for this gene than functional constraints, as found in other insects or stag beetles [8,61,62].

Evolutionary Rates of PCGs
Five available mitogenomes from the genus Prosopocoilus (from P. castaneus, P. laterotarsus, P. confucius, P. astacoides, and P. blanchardi) were used to calculate the evolutionary rate of its PCGs. For each PCG, the ratio of non-synonymous substitution (Ka) to synonymous substitution (Ks) was calculated ( Figure 3). The NADH dehydrogenase subunits (nad1-6 and 4l) and ATP synthase subunits (atp8 and atp6) had higher Ka/Ks ratios than the cytochrome oxidase subunits (cox1, cox2, and cox3) and cytochrome b (cytb). This phenomenon suggested that various functional genes in the mitochondria of Prosopocoilus have undergone different selection pressures during the evolutionary process. The order of Ka/Ks for 13 protein genes was as follows: Notably, the Ka/Ks ratios of the 13 PCGs among the five mitogenomes of Prosopocoilus were all less than 1.0, suggesting that these functional genes were all under strong purifying selection, as was reported in other insects [8,14,60]. The slowest and fastest evolution rates were observed for cox1 and atp8 genes, respectively ( Figure 3). Furthermore, the cox1 had the smallest evolutionary rate, indicating that positive selection was less powerful for this gene than functional constraints, as found in other insects or stag beetles [8,61,62].

Intergenic Spacers
In the two newly sequenced mitogenomes, the sizes of intergenic spacers were examined, varying from 1 bp to 375 bp in P. castaneus, and 1 bp to 158 bp in P. laterotarsus ( Table 3). The large intergenic spacers in these Prosopocoilus mitogenomes were located between trnI and trnQ, consistent with those in the previously studied species of P. blanchardi (4051 bp), P. astacoides (375 bp), and P. confucius (580 bp) [4,23]. In Lucanidae, the notably large IGSs between trnI and trnQ were only found at this particular position in mitogenomes of the genus Prosopocoilus at present. Moreover, a short repetitive sequence with the sequence TAAAA was identified within the large IGS in both of the two newly sequenced mitogenomes. Compared with three previously reported Prosopocoilus species, the sites and number of times this short repetitive sequence (TAAAA) appeared and the lengths of the intergenic region when it was repeated were different among the five sequenced Prosopocoilus mitogenomes (Figure 4). In these five Prosopocoilus mitogenomes, the short sequence (TAAAA) appeared four times (P. castaneus and P. laterotarsus), five times (P. blanchardi), and three times (P. confucius and P. astacoides), and the length of the intergenic region between the repetitive sequence (TAAAA) differed among these five species (Figure 4). Additionally, in these five species of Prosopocoilus, we also detected an intergenic spacer of 18 bp between trnS2 (UCN) and nad1 with the conserved motif (TACTAAA), similar to other reported lucanid beetles [4,20,[50][51][52][53][54]. The IGS between trnS2 (UCN) and nad1 is common in Coleoptera mitogenomes, but varies in length [4,20,23,[50][51][52][53][54]. The characteristics of the large IGS between trnI and trnQ are only present in the representatives of this genus among the family Lucanidae; therefore, we propose that this feature may be synapomorphic for the members of the genus Prosopocoilus. This large IGS between trnI and trnQ could support the taxonomic positions of the other 195 species currently within Prosopocoilus. For example, Epidorcus gracilis, which was previously considered a species in the genus Prosopocoilus, does not contain this large IGS [4]. Previous studies have supposed that large intergenic spacers have possible roles in insect mitogenomic evolution owing to their existence in most insect species, despite their irregular appearance [55,56]. It is possible that this difference may have been caused by environmental selection pressures at the time during the evolutionary process and thus could serve as a useful phylogenetic signal [20]. The IGS between trnI and trnQ may be a useful marker for distinguishing Prosopocoilus from its closely related and indistinguishable genera with respect to this intergenic spacer that may exist in all studied Prosopocoilus species while it is absent in other genera. The present study has provided meaningful implications into the roles of these IGSs in the phylogenetic analysis of Prosopocoilus stag beetles. From our phylogeny is that the species that follow the genus Prosopocoilus have in common the presence of this large IGS that the implication could be in the role of this non-coding region in the evolution of the mitogenome of these species. Figure 3. Evolution rate of mitochondrial protein-coding genes among five Prosopocoilus stag beetle species in the present study (P. laterotarsus, P. castaneus, P. confucius, P. astacoides, and P. blanchardi). Ka represents non-synonymous substitution rate, Ks represents synonymous substitution rate, and Ka/Ks represents the evolution rate of each PCG.

Intergenic Spacers
In the two newly sequenced mitogenomes, the sizes of intergenic spacers were examined, varying from 1 bp to 375 bp in P. castaneus, and 1 bp to 158 bp in P. laterotarsus ( Table 3). The large intergenic spacers in these Prosopocoilus mitogenomes were located between trnI and trnQ, consistent with those in the previously studied species of P. blanchardi (4051 bp), P. astacoides (375 bp), and P. confucius (580 bp) [4,23]. In Lucanidae, Figure 3. Evolution rate of mitochondrial protein-coding genes among five Prosopocoilus stag beetle species in the present study (P. laterotarsus, P. castaneus, P. confucius, P. astacoides, and P. blanchardi). Ka represents non-synonymous substitution rate, Ks represents synonymous substitution rate, and Ka/Ks represents the evolution rate of each PCG. marker for distinguishing Prosopocoilus from its closely related and indistinguishable genera with respect to this intergenic spacer that may exist in all studied Prosopocoilus species while it is absent in other genera. The present study has provided meaningful implications into the roles of these IGSs in the phylogenetic analysis of Prosopocoilus stag beetles. From our phylogeny is that the species that follow the genus Prosopocoilus have in common the presence of this large IGS that the implication could be in the role of this non-coding region in the evolution of the mitogenome of these species. Figure 4. Composition of the large intergenic spacer between trnI and trnQ among the five sequenced mitochondrial genomes of Prosopocoilus in the present study. The grey-shaded region represents the short sequence repeat (TAAAA). The green-shaded region represents the size of the spacers between the short sequence (TAAAA) repeats.

Phylogenetic Analyses
In this study, the nucleotide and amino acid sequences of the 13 PCGs were used to reconstruct the phylogenetic relationship through BI and ML inference methods and generate four nearly identical topologies (Figures 5 and 6). Monophyly of the family Lucanidae was strongly supported (BPP = 1, MLP = 100), consistent with the phylogeny inferred from the multi-gene fragments in previous works [1,35]. Within the family Lucanidae ( Figures 5 and 6, blue shaded), close relationships were observed among Prosopocoilus, Dorcus, and Serrognathus. According to the topologies, the current genus Prosopocoilus was monophyletic and included the five representatives discussed in this study (BPP = 1, MLP = 100). P. castaneus, P. astacoides, and P. blanchardi formed the sister clade of the P. confucius and P. laterotarsus clade (Figures 5 and 6, yellow shaded), whereas the subclade of K. guangxii was a sister group to the subclade of E. gracilis + S. platymelus. These two subclades were both included in Dorcus sensu lato. Additionally, from our phylogenomic tree, Prosopocoilus was a sister group to the clade ((Hexarthrius + Rhaetus) + Pseudorhaetus) consistent with the results of a study by Lin et al. [4]. The three species of the D. velutinus complex, Dorcus curvidens hopei and Dorcus parapllelipipedus, were clustered in a monophyletic clade, as reported by Chen et al. [20]. These results were also consistent with previous morphological comparisons [26,27]. Five species (P. laterotarsus, P. castaneus, P. confucius, P. astacoides, and P. blanchardi) shared high similarities with the type species of Prosopocoilus (Prosopocoilus cavifrons Hope & Westwood, 1845) in their external characteristics and the genitalia traits of both sexes [26,27], suggesting that these five representatives may follow the genus Prosopocoilus sensu stricto. The three others K. guangxii, E. gracilis, and S. platymelus showed partial morphological overlap with typical members of Prosopocoilus, Dorcus, and Serrognathus, and their taxonomic positions have been frequently adjusted among different genera by different taxonomists and amateurs [27,33,[63][64][65]. Our phylogenetic analyses based on the mitogenomic data supported that the five species exhibited Figure 4. Composition of the large intergenic spacer between trnI and trnQ among the five sequenced mitochondrial genomes of Prosopocoilus in the present study. The grey-shaded region represents the short sequence repeat (TAAAA). The green-shaded region represents the size of the spacers between the short sequence (TAAAA) repeats.

Phylogenetic Analyses
In this study, the nucleotide and amino acid sequences of the 13 PCGs were used to reconstruct the phylogenetic relationship through BI and ML inference methods and generate four nearly identical topologies ( Figures 5 and 6). Monophyly of the family Lucanidae was strongly supported (BPP = 1, MLP = 100), consistent with the phylogeny inferred from the multi-gene fragments in previous works [1,35]. Within the family Lucanidae ( Figures 5 and 6, blue shaded), close relationships were observed among Prosopocoilus, Dorcus, and Serrognathus. According to the topologies, the current genus Prosopocoilus was monophyletic and included the five representatives discussed in this study (BPP = 1, MLP = 100). P. castaneus, P. astacoides, and P. blanchardi formed the sister clade of the P. confucius and P. laterotarsus clade (Figures 5 and 6, yellow shaded), whereas the subclade of K. guangxii was a sister group to the subclade of E. gracilis + S. platymelus. These two subclades were both included in Dorcus sensu lato. Additionally, from our phylogenomic tree, Prosopocoilus was a sister group to the clade ((Hexarthrius + Rhaetus) + Pseudorhaetus) consistent with the results of a study by Lin et al. [4]. The three species of the D. velutinus complex, Dorcus curvidens hopei and Dorcus parapllelipipedus, were clustered in a monophyletic clade, as reported by Chen et al. [20]. These results were also consistent with previous morphological comparisons [26,27]. Five species (P. laterotarsus, P. castaneus, P. confucius, P. astacoides, and P. blanchardi) shared high similarities with the type species of Prosopocoilus (Prosopocoilus cavifrons Hope & Westwood, 1845) in their external characteristics and the genitalia traits of both sexes [26,27], suggesting that these five representatives may follow the genus Prosopocoilus sensu stricto. The three others K. guangxii, E. gracilis, and S. platymelus showed partial morphological overlap with typical members of Prosopocoilus, Dorcus, and Serrognathus, and their taxonomic positions have been frequently adjusted among different genera by different taxonomists and amateurs [27,33,[63][64][65]. Our phylogenetic analyses based on the mitogenomic data supported that the five species exhibited typical characteristics of Prosopocoilus, whereas the other three, which had partial morphological features, formed a different clade from the genus Prosopocoilus. We also found that Dorcus was a sister of the clade ((Epidorcus + Serrognathus) + Kirchnerius). The relationships among these lucanid representatives have been described in previous studies [4,62]. The clade Prosopocoilus shares the presence of the large IGS, providing robust information for determining the taxonomic positions of the species following the genus Prosopocoilus. Additionally, our results continue to support that mitogenomes could be very useful molecular tools for improving our understanding of the phylogeny of problematic taxa in Lucanidae, despite the dense sampling materials needed. typical characteristics of Prosopocoilus, whereas the other three, which had partial morphological features, formed a different clade from the genus Prosopocoilus. We also found that Dorcus was a sister of the clade ((Epidorcus + Serrognathus) + Kirchnerius). The relationships among these lucanid representatives have been described in previous studies [4,62]. The clade Prosopocoilus shares the presence of the large IGS, providing robust information for determining the taxonomic positions of the species following the genus Prosopocoilus. Additionally, our results continue to support that mitogenomes could be very useful molecular tools for improving our understanding of the phylogeny of problematic taxa in Lucanidae, despite the dense sampling materials needed.

Conclusions
Our study presents and describes the mitochondrial genomes of two stag beetles from the genus Prosopocoilus. The two mitogenomes of the genus Prosopocoilus maintained the typical gene content and organization of the ancestor mitogenome organization as proposed. The evolutionary rates of 13 PCGs among our studied Prosopocoilus (including three previously reported species) indicated that their evolution occurred under purifying selection. The large intergenic spacer was identified in each of the five Prosopocoilus mitogenomes, and comparisons suggested that the characteristics of large intergenic spacers (presence or absence, size, and location) may be phylogenetically meaningful for evaluat- Figure 6. Phylogenomic relationships of Lucanidae inferred from mitogenomes. Three other scarab species (Cheirotonus jansoni, Rhopaea magnicornis, and Protaetia brevitarsis) were used as the outgroup. The phylogenetic topology was obtained based on the data set AA. Blue shading shows Lucanidae and yellow shading shows Prosopocoilus. The numbers at the nodes are BI posterior probabilities (left) and ML bootstrap values (right).

Conclusions
Our study presents and describes the mitochondrial genomes of two stag beetles from the genus Prosopocoilus. The two mitogenomes of the genus Prosopocoilus maintained the typical gene content and organization of the ancestor mitogenome organization as proposed. The evolutionary rates of 13 PCGs among our studied Prosopocoilus (including three previously reported species) indicated that their evolution occurred under purifying selection. The large intergenic spacer was identified in each of the five Prosopocoilus mitogenomes, and comparisons suggested that the characteristics of large intergenic spacers (presence or absence, size, and location) may be phylogenetically meaningful for evaluating the genus Prosopocoilus. Our phylogenomic analyses including two newly sequenced species, further supported that P. castaneus and P. laterotarsus were clustered in a clade with typical Prosopocoilus species (P. confucius, P. astacoides, and P. blanchardi). Although we were unable to fully confirm that this large IGS was present in all Prosopocoilus species owing to limited sample availability, our findings established a new potential candidate for molecular identification of this genus. Moreover, our findings suggest that variations in the quantity, sequence, and location of IGSs may also be important signals for phylogenomic and evolutionary studies at lower taxonomic levels if these details become available for more taxa in the future. Finally, our results also indicated that mitogenomes could provide useful molecular evidence for improving our understanding of the evolution of stag beetles.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/ani12131595/s1. Table S1: The best-fit schemes and evolutionary models for two datasets. Figure S1: The secondary structure of 22 tRNAs for P. castaneus. Figure S2: The secondary structure of 22 tRNAs for P. laterotarsus.  Data Availability Statement: All the mitochondrial genome sequences were submitted to GenBank (accessions ON401054 and ON401055, as in Table 2), and they will be accessible when the article is published.