Complete Mitochondrial Genome Characterization and Phylogenomics of the Stingless Bee, Heterotrigona itama (Apidae: Meliponini)

TL;DR

This paper analyzes the complete mitochondrial genome of Heterotrigona itama, a vulnerable stingless bee species, to understand its genetic makeup and evolutionary relationships for conservation and research purposes.

Contribution

The study provides the first complete mitochondrial genome characterization of Heterotrigona itama, revealing its gene structure and phylogenetic placement.

Findings

The mitogenome is 15,318 bp with 13 protein-coding genes, 22 tRNA genes, and 2 rRNA genes.

The genome has a strong A+T bias (75.41%) and conserved gene arrangement similar to Tetragonula species.

Phylogenetic analysis shows H. itama is more closely related to bumblebees than honeybees.

Abstract

Stingless bees serve an important part in tropical ecosystems by pollinating and producing medicinal honey. Heterotrigona itama is a commercially significant stingless bee species in Southeast Asia. However, overexploitation has raised questions about its conservation status in some areas, notably Thailand. We examined H. itama’s entire mitochondrial genome to gain a better understanding of its genetic traits and evolutionary ties. Our findings shed light on gene organization, nucleotide composition, codon usage, and evolutionary relationships among other corbiculate bees. This knowledge will help to further taxonomy, population genetics, and stingless bee conservation efforts.

With increasing demand for stingless bee honey, meliponiculture has gained widespread attention. Heterotrigona itama is one of the most economically important species. However, excessive exploitation for commercial purposes has led to population declines, and the species is now considered vulnerable in Thailand. Despite its ecological and economic significance, genomic and taxonomic information on H. itama remains limited. In this study, we sequenced and characterized the complete mitochondrial genome (mitogenome) of H. itama to explore its genome structure and phylogenetic position. The circular mitogenome is 15,318 bp in length and consists of 13 protein-coding genes (PCGs), 22 transfer RNA (tRNA) genes, and 2 ribosomal RNA (rRNA) genes. The genome exhibits a strong A+T bias (75.41%), which affects codon usage and amino acid composition. Isoleucine, methionine, and phenylalanine were the most commonly encoded amino acids. Gene arrangement was highly conserved and closely resembled that of Tetragonula species. Phylogenetic analyses confirmed that H. itama clusters with other stingless bees and is more closely related to bumblebees than to honeybees. Several gene rearrangements suggest a high degree of mitogenomic plasticity. This study provides essential genomic resources for future studies in systematics, phylogenetics, population genetics, and conservation of stingless bees in the Meliponini tribe.

Figures

Full text

1. Introduction

The highly eusocial stingless bees are among the most diverse corbiculate bees, classified under the tribe Meliponini of the Apinae subfamily. At least 600 species have been described and are distributed globally between the Tropic of Cancer and the Tropic of Capricorn [1,2]. The beneficial traits of stingless bee pollinators include their loyalty and constancy to specific flowers, ability to sustain perennial populations and colonies, harmless stings, ease of beekeeping, and strong worker recruitment behavior [3,4,5]. Additionally, stingless bee honey production is increasing because of its rich medicinal value [6,7,8]. Meliponine bees play essential roles in tropical ecosystems and produce valuable hives. They are also an essential example in understanding the evolution, systems, and diversity of social insects [9].

In Thailand, Heterotrigona itama, Cockerell, 1918 is one of the most important extensively used species in meliponiculture. This species is distributed throughout Thailand, Malaysia, Singapore, and Indonesia [10]. Meliponiculture of H. itama includes propagation in apiaries to sell colonies, produce honey products, and pollinate crops [11,12]. Due to the rising demand for stingless bee honey, interest in meliponiculture has also increased. Therefore, beekeeping H. itama can potentially increase household income in many rural areas throughout Thailand. The excessive exploitation of wild stingless bees for commercial purposes has reduced their population size, resulting in H. itama being allocated an at risk conservation status in Thailand.

Owing to their small size and diversity in terms of external morphology [13], as well as overlapping identification traits among species, many stingless bee species have failed to be identified, leaving the systematic and evolutionary histories of many species uncertain [14]. These factors influence thorough research and the efficient utilization of stingless bee resources. Some previous studies have presented important references for the positions of high-level taxa, whereas sufficient taxonomic data remain rare for some smaller genera, including Heterotrigona [12,15,16,17].

Although morphology-based classification and identification are critical, researchers typically require considerable practice to precisely identify the subtle differences among some stingless bee species [18]. On the other hand, DNA data are less reliant on expertise and can be effectively used to identify species-level taxa and clarify conflicting evolutionary histories [19]. High mutation rate and non-recombination, mitochondrial genome data are most valuable for identifying and classifying cryptic stingless bee taxa [14]. Due to their advantages as a biological marker and easy availability, mitochondrial genome data have been extensively used in identification of cryptic species, population genetics studies, evaluation of microevolutionary processes, molecular phylogenetic and systematics investigations, and genome characterization analysis [14,20].

Several previous studies indicated that complete mitochondrial genome (mitogenome) sequences are particularly effective for examining species biodiversity and conservation [18,21,22]. The molecule of bee mitogenome is a circular double-stranded one that contains 13 protein-coding genes (PCGs), 22 transfer RNAs (tRNAs), and 2 ribosomal RNAs (rRNAs) [23,24]. Additionally, the mitochondrial genome also consists of the noncoding and D-loop regions, which controls both replication and transcription of the whole mitogenome [18,25].

The objectives of this study were (1) to analyze the whole mitogenome of H. itama; (2) examine the structural features, base composition, bias of codon usage, and tRNA gene alteration; and (3) construct a phylogenetic tree with other hymenopterans of the family Apidae, gene order, and gene arrangement. These mitochondrial genome sequences will be useful for studying the systematic analysis and evolutionary history of this species, thereby providing insights for taxonomic and conservation biology.

2. Materials and Methods

2.1. Sample Collection and Species Identification

The adult worker specimens of H. itama were obtained from managed colonies belonging to two commercial meliponaries located in the Narathiwat (three samples: 6°24′11″ N; 101°42′22″ E) and Nakhon Sri Thammarat (two samples: 8°05′09″ N; 99°52′46″ E) provinces. Bee samples were collected from the entrance of the hive. All bee samples were initially preserved in absolute ethanol, kept at −20 °C, and deposited in the insect specimen room of the Department of Entomology, Faculty of Agriculture, Kasetsart University (ENTO‒AGR‒KU), Thailand, and the specimen’s voucher was archived (ENTOKU Voucher HI01‒HI05). Species confirmation was carried out based on the taxonomic literature guidelines provided by Samsudin et al. [16] and Trianto et al. [17]. Species identification was confirmed by DNA barcoding at the Ento‒AGR‒KU.

2.2. Ethics Statement

In this study, written consent was obtained from stingless beekeepers, who approved and agreed to participate in the research by allowing bee sample collection from their farms. Farm owners assisted in preparing colonies for sample collection. No special permits were required, as the study did not involve endangered or protected species. The number of stingless bee samples collected was minimal, and ethical research standards were followed. All animal experiments adhered to the guidelines set by the Animal Experiment Committee of Kasetsart University (approval no. ACKU66−AGR−015).

2.3. DNA Extraction

The thorax of a single worker of H. itama per colony was used to extract genomic DNA using a DNeasy Blood and Tissue Kit (Qiagen, Germantown, MD, USA). The quality of extracted genomic DNA was quantified using 0.8% agarose gel electrophoresis. The quantity of the DNA was confirmed by spectrophotometric measurements at 260 and 280 nm using a Tecan Microplate Reader Infinite 200 Pro (Thermo Fisher Scientific, Zurich, Switzerland). Then, a REPLI-g mitochondrial DNA kit (Qiagen, Germantown, MD, USA) was used to isolate mitochondrial DNA from the genomic DNA according to the production manual.

2.4. Sequencing, Assembly, Annotation, and Analysis

In this study, the Illumina NextSeq 500 platform was used to sequence the mitogenome library of H. itama. Primary FASTQ files were generated. Then, the total number of nitrogenous bases, number of reads, and quality score of 30 (Q30) were determined. The fastp software (version 0.12.4) [26] was utilized to remove the adapter, junk sequences, and ambiguous reads from the high-quality data. The filtered sequences were then aligned with reference sequences using BWA software (version 0.7.17) [27]. Consensus sequences were extracted and analyzed for protein-coding genes (PCGs), transfer RNAs (tRNAs), and ribosomal RNAs (rRNAs) using SAM software (http://samtools.sourceforge.net) [28]. Finally, gene prediction for the invertebrate mitochondrial genome was performed using MITOS web server (http://mitos.bioinf.uni-leipzig.de/index.py) [29].

To construct larger stretches of DNA contigs, de novo Bruijn graph-based assembly was performed using short-read data files. The resulting contigs were utilized for annotation of the mitogenome. The Megahit assembler software (version 1.2.9) [30] implemented in the MitoZ package [31] was used to generate the mitogenomes. The quality of the genome assembly was quantified using QUAST [32]. The whole mitochondrial genome of H. itama was further annotated and subsequently substantiated with three mitogenome sequences, including Tetragonula iridipenis (OQ139639), Tetragonula pagdeni (OK336459.1), and Lepidotrigona flavibasis (MN747147), obtained from the GenBank. The secondary structure of tRNA genes of the H. itama mitogenome was also predicted using MITOS. Then, the GCview [33] was performed to examine the entire mitogenome map, GC content, and GC skewness.

To examine the relative synonymous codon usage (RSCU), codon usage, and the A+T bias of the 13 PCGs, the MEGA X software (version 10.2.6) [34] was performed. Then, the AT and GC skews were examined following the report of Perna and Kocher [35]. The intergenic spacers and overlapping regions between the genes of H. itama mitogenome were manually calculated. Finally, the complete sequences of mitochondrial genome of H. itama were deposited in the NCBI database (GenBank accession no: PQ759010 and PQ766639).

2.5. Phylogenomic Reconstruction

The taxonomic and phylogenetic relationships of H. itama with other corbiculate bees (Meliponini, Apini, and Bombini) were examined using 19 complete mitochondrial genome sequences obtained from GenBank.

To estimate the best-fit model for phylogeny construction of each partitioning scheme, the software Partition Finder [36] was used with a corrected Akaike Information Criterion (AICc). The GTR + I model was selected using AIC for nucleotide alignment analysis. All complete mitogenome sequences from 19 different corbiculate bee species, along with an outgroup comprising Colletes gigas (KM978210), Euaspis polynesia (MT909816.1), and Megachile sculpturalis (KT223644), were aligned using MEGA X software. Maximum likelihood (ML) and Bayesian inference (BI) methods were used to determine phylogenetic relationship using the online CIPRES Science Gateway [37]. The ML analysis was conducted using IQ-TREE [38], followed by 1000 ultrafast bootstrap replicates [39] to assess the bootstrap support (BS) values of the tree topology. The BI analysis was executed with MrBayes (version 3.2.7) software [40]. A clade in the reconstructed phylogenetic tree was considered strongly supported if the ultrafast bootstrap support was at least 95% and the Bayesian bipartition posterior probability was 0.95 or higher [39,41].

2.6. Gene Rearrangement Assessment

The complete mitochondrial genome of H. itama, along with those of 22 related bee species downloaded from GenBank, was analyzed for gene order and potential gene rearrangements using PhyloSuite software (http://phylosuite.jushengwu.com/) [42]. The results were then visualized on the iTOL web server [43]. To quantitatively compare the rearrangements in each mitogenome, the rearrangement score (RS) and the rearrangement frequency (RF) of each individual gene were computed using qMGR software (http://qmgr.hnnu.edu.cn/) [44].

3. Results

3.1. Organization and Nucleotide Composition

The complete circular mitochondrial genome of the Thai commercial stingless bee H. itama was 15,318 bp (Figure 1). Thirty-four sequence components were identified, including 22 tRNA genes, 2 rRNA genes, and 13 PCGs (Table 1 and Figure 1). The strand localization demonstrated that the H- (+) and L-strands (−) contained 28 and 9 genes, respectively (Table 1). The maximum and minimum lengths of the mitochondrion genome of H. itama were 1653 bp (nad5) and 57 bp (tRNA-Ser), respectively. The nitrogenous base composition in the H. itama mitogenome sequence was found to be 35.06% adenine (A), 42.42% thymine (T), 11.30% guanine (G), and 12.22% cytosine (C) (Table 2). The A+T bias varied between 65.22% and 85.24% across the 13 protein-coding genes (PCGs) and was 75.41% for the entire mitogenome. The AT skew ranged from −0.214 to −0.014 across the thirteen protein-coding genes (PCGs), showing a relatively stronger bias compared to the two rRNA genes, rrnS (0.070) and rrnL (0.014) (Table 2).

3.2. PCGs and Codon Usage Bias

Of the 15,318 bp in the complete mitogenome sequence of H. itama, 13 PCGs comprising 11,045 bp (72.10%) were identified. These PCGs varied in length from 168 bp (atp8) to 1653 bp (nad5). Only the nad1 gene was found on the L-strand (−), whereas the other 12 PCGs were positioned on the H-strand (+). Five start codons (ATT, ATA, ATC, ATG, and TTG) were identified in the mitogenome of H. itama. The ATT codon was used by six PCGs (cox1, cox2, atp8, nad4l, nad5, and nad1) as their start codon, ATG was adopted by three genes (cox3, nad4, and cytb), ATA was used by two genes (nad3 and nad6), and ATC and TTG were used by the nad2 and atp6 genes, respectively. The TAA was characterized by eight PCGs (atp8, nad1, nad2, nad3, nad4L, nad5, nad6, and cytb) as their stop codon, while TAG was used by the rest of the four genes, and the uncompleted stop codon “T” was used by nad4. The RSCU and amino acid usage in the PCGs are presented in Table 3. The most commonly occurring amino acids were arginine (Arg), leucine (Leu), valine (Val), serine (Ser), proline (Pro), threonine (Thr), alanine (Ala), and glycine (Gly), while tryptophan (Trp) was the least frequent (Figure 2). For the codon usage in the mitochondrial genome of H. itama, the three most commonly used were UUU, AUU, and AUA (Table 3).

3.3. Transfer RNAs and Ribosomal RNAs

All 22 tRNA genes were distributed over the mitogenome of H. itama, varying between 57 (trnS1^Ser^) and 70 (trnY^Tyr^) bp in length. The H-strand (+) of the mitogenome contained 14 tRNA genes, while the remaining 8 tRNA genes were situated on the L-strand (−) (Table 1 and Figure 3). The 22 tRNAs in the H. itama mitogenome sequence had a combined length of 1469 bp, accounting for 9.59% of the entire mitogenome (Table 1). These 22 tRNA genes exhibited characteristic cloverleaf secondary structures, except for trnS1^Ser^ (UCU), which lacked a dihydrouridine (DHU) stem and loop (Figure 3).

Both rrnS and rrnL genes were located on the H-strand (+) of the mitochondrial genome. Together, they constituted 13.42% (2056 bp in length) of the entire H. itama mitogenome (Table 1). Positioned between tRNA^Gln^ and tRNA^Val^ genes, the rrnS gene is 761 bp long and has an A+T content of 75.82%. The rrnL gene, measuring 1295 bp with an A+T content of 77.52%, was situated between tRNA^Val^ and tRNA^Ser^ (Figure 1).

3.4. Intergenic Spacer and Overlapping Region

In the mitochondrial genome of H. itama, six overlapping regions between genes, totaling 8 bp in length, were identified. The lengths of these regions ranged from 1 bp, as observed between tRNA^Trp^ and tRNA^Tyr^, tRNA^Leu^ and cox2, nad6 and cytb, and cytb and tRNA^Ser^, to 2 bp, found between rrnS and tRNA^Val^, as well as rrnL and tRNA^Leu^. The intergenic spacers in the H. itama mitogenome ranged from 1 to 100 bp in length, comprising 770 bp, and were distributed across 24 regions throughout the mitogenome (Table 1). The three largest intergenic spacers were located between tRNA^Pro^ and nad4l, tRNA^Arg^ and tRNA^Gln^, and tRNA^Arg^ and nad6, with lengths of 100, 100, and 87 bp, respectively.

3.5. Phylogenomic Relationship

To gain a comprehensive understanding of genome-level evolution in H. itama, we conducted phylogenomic analyses using concatenated nucleotide sequences from 13 PCGs, 2 rRNA genes, and 22 tRNA genes, along with those from 22 representative hymenopteran insects with complete mitogenomes available in GenBank.

The phylogenetic analysis clearly grouped three major clades corresponding to the tribes Meliponini, Bombini, and Apini, with Colletes gigas (KM978210), Euaspis polynesia (MT909816), and Megachile sculpturalis (KT223644) serving as the outgroup. All nine meliponine-producing bees formed a monophyletic clade (Figure 4). H. itama from the present study was first grouped with two other species of Tetragonula (T. iridipennis: OQ139639; T. pagdeni: OK336459) and then combined with three other genera of stingless bees: Lepidotrigona (L. terminate: MN737481; L. flavibasis: MN747147), Tetragonisca (T. angustula: OR030859), and Melipona (M. fusciculata: MH680930, M*. scutellaris*: KP202303, and M. bicolor: AF466146).

3.6. Gene Rearrangement

The gene order of the H. itama mitochondrial genome closely resembled that of two Tetragonula species but differed from the gene arrangements found in six others previously reported stingless bee mitogenomes, which exhibited several gene rearrangements (Figure 5). The gene rearrangement events aligned with the phylogenetic tree derived from mitogenome sequences (Figure 4 and Figure 5). Additionally, variations in the position of tRNA^Lys^ in the stingless bee mitogenome were observed. In Heterotrigona and Tetragonula, a tRNA^Lys^‒tRNA^Ala^‒tRNA^Ile^‒tRNA^Met^ block was found, whereas in Lepidotrigona and Melipona, the block was altered to tRNA^Met^‒tRNA^Lys^‒tRNA^Ala^‒tRNA^Ile^ and tRNA^Ile^‒tRNA^Ala^‒tRNA^Lys^‒tRNA^Met^, respectively (Figure 5).

In this study, we calculated the RS of each mitogenome by summing the RS values of all genes in the mitogenome. The RF of each gene was also computed using the mitogenome. Our results indicate that H. itama was a highly active group, comparable to stingless bees of the genus Tetragonula, and exhibited the highest RS among corbiculate bees (Figure 6). Among the individual genes, tRNA genes exhibited higher rearrangement frequencies, as observed in the PCGs (Figure 7). Among tRNA genes, tRNA^Leu^ showed the highest rearrangement frequency. Apart from tRNA genes, PCGs and rRNA genes also underwent significant rearrangements in H. itama. The conserved nad1‒tRNA^Leu^‒rrnL‒tRNA^Val^‒rrnS block in Apini and Bombini was transposed into Heterotrigona and Tetreagonula (Figure 5). The conserved block of nad1‒tRNALeu‒rrnL‒tRNAVal‒rrnS found in Apini and Bombini was transposed in Heterotrigona and Tetragonula (Figure 5). Additionally, the tRNA^Phe^‒nad5‒nad4‒nad4l block in Apis species was transposed and inverted into H. itama. These observations indicate that H. itama has a relatively active mitochondrial genome.

4. Discussion

The mitochondrial genome of insects, including hymenopteran insects, is typically shown as a circular double-stranded molecule of approximately 14–19 kb in length. It generally contains 37 genes, comprising 13 PCGs, 2 ribosomal RNA (rRNA) genes, and 22 transfer RNA (tRNA) genes [24]. In this study, the total length of the complete mitochondrial genome of H. itama is 15,318 bp. This is consistent with the mitochondrial genome sizes of other stingless bees (T. iridipennis: 15,045 bp, T. pagdeni: 16,061 bp, L. flavibasis: 15,408 bp, L. terminata: 15,431 bp, and M. bicolor: 15,001 bp) and honeybees (A. mellifera: 16,336, A. cerana: 15,895 bp, A. florea: 15,993 bp, and A. dorsata: 15,892 bp). However, the mitochondrial genome of H. itama is slightly shorter than that of bumblebees, such as B. terrestris (17,232 bp) and B. canariensis (17,300 bp).

An apparent bias towards the nucleotide bases A and T (AT bias) was observed in the H. itama mitogenome sequence. Among the four nitrogenous bases, base T was the most abundant, whereas base G was the least abundant in the entire mitogenome. Furthermore, our findings confirm that the nucleotide composition of stingless bees has a high A+T content [14,24,45,46]. In general, both AT and GC skewness can be used to describe the base composition bias in nucleotide sequences [18]. Clary and Wolstenholme [47] reported that bias in nucleotide usage can be attributed to the active use of specific bases by DNA polymerases during mitochondrial DNA replication. Moreover, Xia [48] reported that subsequent transcription generates an AT bias in species relying on mitochondrial efficiency to sustain a high metabolic rate, as breaking A-T bonds requires less energy. This study revealed a similar pattern, with a negative AT skewness value (−0.083), suggesting a preference for base T over base A in the H. itama mitogenome.

As observed in other stingless bees, such as M. bicolor [49], M. scutellaris [50], L. terminata [51], L. flavibasis [14], T. pagdeni [46], and T. iridipennis [18], all 13 PCGs were allocated throughout the mitogenome of H. itama, with slight variations in size and position. The start codon (ATN) and the incomplete stop codon (T), which are commonly found in hymenopteran bees, were also observed in H. itama. No anomalous initiation codons were observed in any of the 13 PCGs of the H. itama mitogenome, unlike L. flavibasis (GTG start codon for nad2) [14] and Apis cerana (TTG start codon for cox1) [52]. In the universal genetic code, the TGA codon serves as a stop codon; however, in most mitochondrial genomes‒except those of higher plants–it encodes tryptophan [53]. The majority of stop codons use TAA for termination, though an incomplete stop codon (nad4) was identified in the H. itama mitogenome. This suggests that polyadenylation likely occurs post-transcription to finalize the termination codon [52]. Moreover, incomplete stop codons are commonly found in protein-coding genes (PCGs) of the Hymenoptera mitochondrial DNA sequenced so far [18,52,54,55,56]. Additionally, RSCU analysis indicated a preference for AT-rich codons in the H. itama mitogenome, primarily attributed to an AT mutational bias [46].

All 22 tRNA genes exhibited the typical cloverleaf structure, except for tRNASer, which lacked a dihydrouridine (DHU) arm and instead formed a simple loop (Figure 3). In most insect and metazoan mitochondria, tRNA^Ser^ (UCU) has an unpaired DHU arm, tRNA^Ser^ (UCU) generally has an unpaired DHU arm, and tRNA^Ser^ (UGA) retains the standard cloverleaf structure [24,57]. The H. itama tRNA^Ser^ (UCU) and tRNA^Ser^ (UGA) structures confirmed this hypothesis. The putative secondary structures of tRNAs are similar to that of T. iridipenis [18], L. flavibasis [14], and T. pagdeni [46], indicating their similar function. The anticodons of tRNAs were similar to that of meliponine bees reported by various studies [14,18,46,49]. In the tRNAs of the H. itama mitogenome, the amino acid acceptor stem (15 bp) and the anticodon loop (7–9 bp) are highly conserved, resembling those observed in T. pagdeni [46], L. flavibasis [14], and A. cerana [52].

Mismatched base pairs in mitochondrial tRNAs have been documented in various hymenopteran insects, including T. iridipenis [18], T. pagdeni [46], L. flavibasis [14], M. bicolor [49], A. cerana [52], A. mellifera [45], Bombus ignitus [55], Diadegma semiclausum [56], and Evania appendigaster [58]. In this study, 10 mismatched base pairs were detected in the tRNAs of H. itama. Specifically, three tRNA genes (tRNA^Cys^, tRNA^Gln^, and tRNA^Thr^) exhibited single T-T mismatches in the acceptor stem, while seven others (tRNA^Lys^, tRNA^Ala^, tRNA^Leu^, tRNA^Asp^, tRNA^Arg^, tRNA^Val^, and tRNA^Pro^) had single G-T mismatches in the same region. Additionally, G-T mismatches were identified in the DHU arm of tRNA^Gln^, tRNA^Pro^, and tRNA^Phe^ (Figure 3). The arrangement of two ribosomal RNA genes in the H. itama mitogenome was consistent with that of T. iridipennis and T. pagdeni but differed from those found in other stingless bees (Lepidotrigona and Melipona species), bumblebees, and honeybees. These variations resulted in their placement in distinct clades of the phylogenetic tree [14,18,49].

Phylogenetic analyses based on complete mitochondrial sequences are more reliable than those using partial or combined partial mitochondrial and nuclear gene sequences, as they produce topologies with strong support values [14,18]. In this study, the complete mitogenome of H. itama, along with those of 22 common hymenopteran species obtained from GenBank, was utilized for phylogenetic analysis. Our results indicate that H. itama and other stingless bee species are closely related with strong support (Figure 4), supporting taxonomic classification using morphological and molecular data. Three major clades were constructed within the Apidae family: Apini, Bombini, and Meliponini, which were highly supported in all analyses. These findings align with those of previous studies [18,46,58,59]. Our results strongly support a relationship among hymenopteran insects.

Most hymenopteran mitogenomes exhibit similar gene order [18,60,61]. Generally, gene rearrangements, including tandem duplication, inverse transposition, genetic transposition, and gene inversion, can cause the loss of gene function [24,26]. Consistent with that in other stingless bees, the absence of a conserved tRNA^Asp^–tRNA^Lys^ block suggests that tRNA^Lys^ movement may have occurred early in the progenitors of the stingless bee lineage, thus separating stingless bees from other hymenopterans (Apini and Bombini tribes) [46,52]. Our results also reveal that the mitogenome of H. itama shows a unique translocation that disrupted the relatively frequent intact nad1–tRNA^Leu^–rrnL–tRNA^Val^–rrnS block observed in hymenopteran mitogenomes. This finding implies that H. itama has experienced a high degree of mitogenomic plasticity. The mitogenomes of most invertebrate species largely retain their ancestral gene order. Nonetheless, natural selection may not entirely remove mitogenome rearrangements, reflecting the evolutionary dynamics and rates within a particular species or genus [62].

In the tropics and subtropics, numerous agricultural crops and wild plants depend mainly on insect pollinators for cross-pollination. Among insect pollinators, honeybees, stingless bees, and bumblebees, all belonging to the family Apidae, are the most important pollination vectors. The analysis of the H. itama mitogenome has enhanced our knowledge of hymenopteran mitochondrial genomes, offering new genetic markers for research in population genetics, systematics, phylogenomics, and biogeography of stingless bees.

5. Conclusions

The complete mitogenome of H. itama was sequenced using next-generation Illumina technology. The circular mitogenome spans 15,318 bp and comprises 37 sequence elements, including 13 protein-coding genes (PCGs), 2 rRNA genes, and 22 tRNA genes. Phylogenetic and gene order analyses revealed that H. itama is closely related to stingless bees of the genus Tetragonula but remains relatively distinct from the six previously reported stingless bee mitogenomes. The genetic data obtained from this study will aid future research in comparative genomics, molecular systematics, and population genetics of various hymenopteran bee species.