Comparative Analyses Reveal Mitogenome Characteristics of Halictidae and Novel Rearrangement (Hymenoptera: Apoidea: Anthophila)

TL;DR

This study sequenced four new bee mitogenomes, revealing gene rearrangements and phylogenetic relationships in Halictidae, a key pollinator group.

Contribution

The discovery of a novel gene rearrangement pattern in bees and the sequencing of four new Halictidae mitogenomes.

Findings

Five gene rearrangement patterns were identified, including one novel pattern in bees (ND4-trnP-ND4L-trnT).

Phylogenetic analysis showed Rophitinae as the basal group within Halictidae.

Mitogenomes had high AT content (82.55–86.44%) and specific codon usage preferences like UUA and UCA.

Abstract

Species of Halictidae are renowned for their outstanding pollination function and variable social behavior, while our understanding of the evolutionary history of this group has been hindered by the scarce molecular data. In this study, we newly obtained four mitogenomes of Halictidae. Using comprehensive comparative genomic methods, we analyzed the mitogenome base composition and codon usage of Halictidae which are similar to published bee species. We found five gene rearrangement patterns, including one novel pattern of bees. Furthermore, we combined published data to reconstruct phylogenetic relationships based on Bayesian Inference and Maximum Likelihood methods. Our study enriches molecular databases for research on bees and provides resources for further study on the evolutionary biology of Halictidae.

Halictidae, as a major pollinator family in bees, has significant ecological value. However, the insufficient molecular data for this group has limited our understanding of the evolutionary history of this group. Herein, we newly sequenced and assembled four mitogenomes of Halictidae, including three species of Nomiinae and one species of Rophitinae. We analyzed the characters of the newly obtained mitogenomes, including nucleotide composition, sequence length, and gene rearrangements. The length of the newly sequenced mitogenomes ranged from 16,492 to 21,192 bp, and all newly obtained mitogenomes contained 22 tRNAs, 13 protein-coding genes, two rRNAs, and one control region. Their AT content (%) ranged from 82.55 to 86.44. Relative synonymous codon usage analysis showed that UUU, UUA, and AUU were the preferred codons. The relative synonymous codon usage > 2 of mostly newly sequenced species was as follows: UUA > UCA > CGA. All newly obtained mitogenomes show gene rearrangement; we found five gene rearrangement patterns in total. Notably, ND4-trnP-ND4L-trnT was the first reported gene rearrangement pattern in bees. In addition, we reconstructed the phylogenetic relationships of Halictidae based on 10 species (eight ingroups and two outgroups), using Bayesian Inference and Maximum Likelihood approaches. Phylogenetic analysis showed that Rophitinae was the basal group within Halictidae.

Figures

Full text

1. Introduction

Among bees, the family Halictidae contains major pollinators, with more than 4400 described species [1,2,3]. Species of this group can be found in nearly all terrestrial environments; in many temperate areas of the world, halictids dominate other bees in number of individuals, except the species Apidae [3,4,5]. Halictidae contains four subfamilies, including Rophitinae, Nomiinae, Halictinae, and Nomioidinae [1,4].

Nomiinae, as one of the subfamilies within the family Halictidae, contains more than 620 described species; Asia and Africa are the hotspots of richness centers [2,3,4]. Nomiinae species hold significant value in agricultural economic development. Certain species within this group are actively managed for crop pollination [6]. For instance, Nomia melanderi serves as an effective pollinator of alfalfa and has been extensively utilized in the United States since the 1990s [6,7]. In addition, this species is the only solitary ground-nesting bee that is artificially managed [6,8,9]. Nevertheless, the studies of this group have predominantly concentrated on the description of new species or the review of a single genus or subgenus. Consequently, there is a significant deficiency in molecular data, and the phylogenetic relationships within Nomiinae remain inadequately understood. Michener, based on morphological characters, divided Nomiinae into 11 genera and 19 subgenera, which were supported by most bee researchers [10]. As research progressed, Baker, Pesenko, and Pauly suggested that the subgenera Pseudapis, Nomiapis, Maynenoiay, Maculonomia, Gnathonomia, and Austronomia need to be treated at the genus level [5,11,12]. However, Pauly mainly focused on the taxa of Africa and the Western Pacific, and both studies only used morphological characters without molecular data, but most bee researchers disagree with those points. Actually, there are significant differences among these subgenera in terms of morphological characteristics. For instance, certain species within Gnathonomia and Maculonomia lack colored bands on the apical margin of the metasomal terga [2,13], which does not conform to the diagnostic criteria of the genus Nomia. It is insufficient to rely on morphological characteristics; additional molecular data should be incorporated. However, only one published mitochondrial genome exists for the >600 species of Nomiinae (accessed on 20 May from NCBI database). Thus, the deficient molecular data of Nomiinae have seriously limited our understanding of the evolutionary history and phylogenetic relationships of this group, and also of Halictidae in general.

Rophitinae is a relatively small subfamily within Halictidae, with 261 described species [1]. Rophitinae is the only subfamily of Halictidae in which all members are solitary [4]. Species within Rophitinae exhibited restricted floral specialization, particularly the genus Dufourea [14,15]. Dufourea is a major genus of Rophitinae; species of this genus are nonmetallic black or dull greenish or bluish metallic, sometimes with the metasoma red [4]. Alves-dos-Santos (2003) observed that mouthparts of bees specializing in Pontederiaceae, such as Dufourea novaeangliae, have evolved hairs to remove pollen from the concealed anthers [16]. Despite the existence of some studies on this group, research on the mechanisms of molecular evolution remains limited. Notably, the available molecular data for this group are insufficient, constraining evolutionary research.

Mitochondrial genomes, as important molecular markers, have typically been employed in investigations pertaining to phylogeny, evolutionary history, speciation, and phylogeography within the taxonomic group of insects [17,18,19]. The maternal mode of inheritance, which is characterized by a high degree of substitution and the ease of obtaining these genomes, renders them particularly advantageous for the aforementioned research [20,21]. Insects' mitochondrial genomes typically span between 14,000 and 20,000 bp in length, encompassing two ribosomal RNAs (rRNAs), 22 transfer RNAs (tRNAs), 13 protein-coding genes (PCGs), and one non-coding control region (CR) [21]. Insects’ mitochondrial genomes have highly conserved composition and organization [21]. Furthermore, the structural characteristics of mitogenomes can offer additional insights and corroborating evidence for taxonomic classification [21,22]. In recent years, the number of complete mitochondrial genomes of insects that have been sequenced has increased, advancing the rapid development of next-generation sequencing technologies to resolve the structure comparison and evolutionary history for different groups [20,21,22,23,24].

Herein, we newly sequenced, assembled, and annotated four mitogenomes of the family Halictidae; all data have been submitted to the NCBI database. Combined with one published mitochondrial genome, we analyzed the main structure, evolutionary rate, and substitution of all newly obtained mitogenomes. In addition, we combined them with six published mitochondrial genomes of Halictidae to reconstruct the phylogenetic relationships of Halictidae based on Bayesian Inference (BI) and Maximum Likelihood (ML) approaches.

2. Materials and Methods

2.1. Taxon Sampling and Sequencing

Here, four bee species within Halictidae were newly sequenced and assembled (detailed information are shown in Table 1). All specimens were stored in 99% ethanol at −20 °C before species identification and DNA extraction. To compare the mitochondrial characters of the four newly sequenced species and also for the phylogenetic analyses, we downloaded one mitochondrial genome of Nomia (Nomia chalybeata) from GenBank. Two taxonomists (Dan Zhang and Zeqing Niu) identified all species based on morphological characteristics [2,3,4].

Legs of each sample were used to extract whole genome DNA with a Qiagen DNeasy Blood & Tissue Kit (Qiagen, Venlo, The Netherlands) and Qubit^®^ DNA Assay Kit in Qubit^®^ 2.0. A flurometer (ThermoFisher, Waltham, MA, USA) was used to measure the concentration of the whole genome DNA. The Illumina NovaSeq 6000 platform was used to generate sequencing libraries of 150 bp paired-end reads with an insert size of 350 bp. Short, adapter-contaminated, and low-quality reads of raw data were removed via Trimmomatic v0.32 [25].

2.2. Assembly, Annotation, and Composition Analyses

We used two different methods to assemble the mitogenome: (1) Illumina reads were assembled by NOVOPlasty v3.8.3 (Brussel, Belgium) [26] and the k-mer sizes of 23–39 bp; (2) High-quality reads were used for assembly via IDBA-UD v1.1.3 (Boston, MA, USA) [27], while the minimum and maximum k-values were 40 and 120 bp, respectively. Geneious 2020.2.1 [28] was performed to compare mitogenome sequences, which were obtained by the above methods and then merged into a single sequence. tRNAscan SE [29] was used for the secondary structure analysis of tRNAs. Clustal Omega in Geneious (2020.2.1) was applied to annotate the rRNAs and PCGs. MEGA X [30] was performed to check the boundaries of PCGs and rRNAs. SeqKit v0.16.0 (Chongqing, China) [31] was implemented to check the nucleotide composition of each gene and the bias of the nucleotide composition. The AT-skew and GC-skew were calculated as follows: AT-skew = (A − T)/(A + T), and GC-skew = (G − C)/(G + C). The relative synonymous codon usage (RSCU) of the newly sequenced species was performed by MEGA X. Mitogenome maps of each newly sequenced species were generated by an online server CGview (https://cgview.ca/, accessed on 30 April 2025).

2.3. Phylogenetic Analysis

To reconstruct the phylogenetic relationships of Halictidae, we selected 10 samples, including eight ingroups and two Colletidae species as outgroups based on prior research of bees [32]. Thirteen PCGs and two rRNA genes of all 10 species were used for phylogenetic analyses. The nucleotide and protein sequences of all samples were aligned by MAFFT v7.450 (Osaka, Japan) with the L-INS-I method [33]. Sequence trimming was conducted via Trimal v1.4.1 (Barcelona, Spain) [34] with “-automated1” strategy. Finally, five matrices were generated to reconstruct the phylogenetic relationship, and FASconCAT-G v1.04 (Santa Cruz, CA, USA) [35] was performed to concatenate each matrix: (1) cds_faa matrix, including all PCG amino acid reads; (2) cds_fna matrix, containing all PCG nucleotide reads; (3) cds_rrna matrix, including all PCG and two rRNA nucleotide reads; (4) cds12_fna matrix, containing all PCG nucleotide reads except the third codon positions; (5) cds12_rrna matrix, including all PCG nucleotide reads, which removed the third codon positions and two rRNA genes.

For all matrices, we used ML and BI methods to reconstruct the phylogenetic relationship within Halictidae. For ML analyses, we used ModelFinder [36] in IQ-TREE 2 (Canberra, ACT, Australia) [37] to select the best substitution model for each matrix. In addition, we used the posterior mean site frequency (PMSF) model [38] (‘−m − mtART + C60 + FO + R’) to reconstruct the phylogenetic relationship of Halictidae in IQ-TREE for matric cds_faa to minimize long-branch attraction artifacts. Phylobayes-MPI (Montréal, QC, Canada) [39] was applied to generate the BI tree using the site-heterogeneous mixture model (−m CAT + GTR). We performed two independent Markov Chain Monte Carlo (MCMC) chains, each with 10,000,000 generations, stopping once we achieved satisfactory convergence (maxdiff < 0.3). Finally, iTOL was used to beautify the phylogenetic tree (https://itol.embl.de/upload.cgi, accessed on 20 April 2025).

3. Results

3.1. Mitogenomic Organization

We sequenced about 6 Gb of raw data for each species. Four species (Nomia thoracica, Lipotriches guihongi, Lipotriches capitata, Dufoure subclavicrus) of Halictidae were obtained, all of which were linear mitogenomes. All newly obtained mitogenomes have been submitted to GenBank under accession numbers PQ050612–PQ050615 (Table 1). For newly obtained mitogenomes, the whole length ranged from 16,492 (L. capitata) to 21,192 bp (D. subclavicrus, Table 2). We also identified one CR, two rRNAs, 22 tRNAs, and 13 PCGs (Figure 1). All newly sequenced mitogenomes showed positive AT skew, except L. capitata, ranging from −0.031 (L. capitata) to 0.047 (D. subclavicrus). The GC skew for all newly reported mitogenomes was negative, ranging from −0.364 (L. guihongi) to −0.112 (L. capitata, Table 2). The GC content (%) ranged from 13.55 (L. capitata) to 17.38 (D. subclavicrus), and the AT content (%) ranged from 82.55 (D. subclavicrus) to 86.41 (L. capitata, Table 2).

3.2. Protein-Coding Genes, Codon Usage, and Evolutionary Rates

The sizes of PCGs, tRNAs, and rRNAs of newly sequenced species are similar. The PCG lengths of the newly obtained species ranged from 11,075 (D. subclavicrus) to 11,142 (L. guihongi). PCGs of all newly obtained mitogenomes showed negative AT skew, ranging from −0.122 (L. capitata) to −0.069 (D. subclavicrus). The GC skew for PCGs of all newly reported mitogenomes were negative, except L. capitata, ranging from −0.153 (D. subclavicrus) to 0.026 (L. capitata, Table 2). The GC content (%) of PCGs ranged from 13.39 (L. capitata) to 18.18 (D. subclavicrus), and the AT content (%) ranged from 81.82 (D. subclavicrus) to 86.61 (L. capitata, Table 2).

CR sizes of L. capitata were the shortest (1576 bp), while D. subclavicrus exhibited the longest CR sizes (6572 bp). All newly sequenced and assembled mitogenomes were started with ATN. However, varying start codons were found: ATP6 and CYTB used ATG as the start codon in two species; COX1, COX2, ND1, ND2, ND5, and ND6 used ATA as the start codon in one species; ATP8, COX2, ND3, ND4, ND6, and COX1 used ATT as the start codon in two species; ND4L, COX1, and COX2 used ATC as the start codon in one species; COX1 used TTG as the start codon in one species and so on (Figure S1). TAA or TAG are the common stop codons of this group, whereas COX1, ND4, and ND5 in one species have TA or T as the stop codons, which were incomplete termination codons (Figure S2). The termination codons of PCGs are frequently found to be incomplete in insects, and they are commonly completed through a process known as polyadenylation, which occurs after the excision of the downstream tRNA gene [40,41,42,43].

Twenty-two tRNAs were identified, ranging from 50 to 77 bp in length. The tRNA characters of newly obtained mitogenomes are as follows: AT content (%) ranged from 84.11 (D. subclavicrus) to 88.20 (L. capitata); the GC content (%) ranged from 11.80 (L. capitata) to 15.89 (D. subclavicrus); and the AT skew and GC skew exhibited positive values (Table 2).

Four newly sequenced mitochondrial genomes exhibited congruent relative synonymous codon usage (RSCU) patterns (Figure 2), which showed the RSCU values for all synonymous codons corresponding to the 22 amino acids, utilizing the 62 available codons in the 13 PCGs of Nomiinae species. Compared to the newly reported species, we found that the preferred codons were UUU, UUA, and AUU. Meanwhile, the frequently utilized amino acids were Leu2, Phe, and IIe, indicating that A/T was the preferred nucleotide composition. The codons with RSCU > 2 values for most newly sequenced species were as follows: UUA > UCA > CGA, while for L. capitata, they were UUA > UCU > GGA (Figure 3, Tables S1–S4). For amino acid, NNA and NNU were the most frequently used codons, reflecting the AT preference in nucleotide composition (Tables S1–S4).

3.3. Gene Rearrangement

In this study, we used the gene arrangements of Drosophila yakuba as the ancestral state [44]. Gene rearrangements existed in all newly sequenced species, and five gene rearrangement patterns were found in total (Figure 3). The gene cluster trnI-trnQ-trnM was frequently rearranged to trnM-trnI-trnQ in three species (D. subclavicrus, L. capitata, N. thoracica, Figure 3); the trnK-trnD tRNA block was rearranged in all newly sequenced species (COX2-trnK-trnD to COX2-trnD-trnK); the trnR was inverted in all newly sequenced Nomiinae species (ND3-trnA-trnR-trnN transformed into ND3-trnR-trnA-trnN); trnP was inverted in L. capitata and N. thoracica (ND4-ND4L-trnT-trnP to ND4 -ND4L-trnP-trnT, Figure 3). The ND4-ND4L-trnT-trnP gene cluster was rearranged to ND4-trnP-ND4L-trnT in L. guihongi (Figure 3), which was the first instance of mitochondrial gene rearrangement found in Apoidea.

3.4. Phylogenetic Relationships

Five matrices were used to reconstruct the phylogenetic relationships of the Halictidae: cds_faa (5375 sites), cds_fna (7148 sites), cds_rrna (12,580 sites), cds12_fna (7148 sites), and cds12_rrna (9006 sites). We used five matrices to generate four phylogenetic topologies based on BI and ML approaches (Figure 4 and Figures S3–S9). These results showed that all of them supporting Rophitinae were the basal group within the Halictidae (Figure 4 and Figures S3–S9). For the Nomiinae, there were four different topologies, all of which show that Austronomia is not a subgenus within the genus Lipotriches, and the subgenus Gnathonomia may be treated as a genus within Nomiinae. Topology 1 (T1) was generated from cds_faa, cds_rrna using BI methods, cds_faa, cds_fna, cds_rrna using the partition model, and cds_faa using the PMSF model: (((Gnathonomia, Lipotriches), Nomia), Austronomia) (Figure 4A, Figures S3, S4 and S7–S9). This topology suggested that the subgenus Gnathonomia is relatively close to the genus Lipotriches. Topology 2 (T2) was generated using cds_fna, cds12_fna based on the BI method, and cds12_fna based on the partition model: ((Nomia, Austronomia), Lipotriches, Gnathonomia) (Figure 4B, Figures S5 and S6). Here, the genus Nomia is close to Austronomia, and the Gnathonomia is the basal lineage. Topology 3 (T3) was generated using the matric cds12_rrna based on the partition model: (((Gnathonomia, Lipotriches), Austronomia), Nomia) (Figure 4C). In T3, the subgenus Gnathonomia is close to the genus Lipotriches. Topology 4 (T4) was generated by the matric cds12_rrna based on BI methods: (((Gnathonomia, Lipotriches), Nomia), Austronomia) (Figure 4D). This topology also supported that the subgenus Gnathonomia is close to the genus Lipotriches. Overall, the phylogenetic relationships of the subfamily Nomiinae were highly unstable with different models.

4. Discussion

4.1. Mitochondrial Genome Structure

We newly sequenced and assembled four mitochondrial genomes within Halictidae. All newly obtained mitochondrial genomes had typical circular double-stranded structures, and the genes exhibited high conservation in length. The arrangement order and nucleotide composition of newly reported data are similar to those of other published insect mitochondrial genomes [18,19,22,23,24]. The length of the mitogenomes was significantly different among species (Figure 2, Table 2), mainly because of variations in control region size. The nucleotide composition of newly sequenced species exhibited a considerable bias toward A + T, and the A + T content was higher than 80% in all species (81.82–86.61%, Table 2), which was commonly observed in published Hymenoptera mitochondrial genomes [45,46]. There were no remarkable differences in base composition discovered among the newly sequenced mitochondrial genomes (Table 2). For insects, AT skew in mitogenomes exhibits positive values in general, while our results showed a negative AT skew for L capitata, and all GC skew values were negative (Table 2), indicating a preference for thymine (T) and guanine (G) bases in Hymenoptera, and similar results have been published for other Hymenoptera and bee species [45,46].

4.2. Mitogenomic Gene Rearrangements

Gene rearrangement in insect mitochondrial genomes was uncommon but has been reported in Thysanoptera, Psocodea, Hemiptera, Hymenoptera, and Trichoptera [47,48,49,50,51]. To date, four main types of gene rearrangements have been observed in insects: local inversion, remote inversion, gene shuffling, and translocation [52]. In Apoidea, inversion and inverse transposition events occurred at a relatively high frequency, similarly to observations in other hymenopterans [50,53,54]. Indeed, insects in the order Hymenoptera are known to exhibit an exceptionally high level of mitochondrial gene rearrangement, which often carries important phylogenetic signals [17,46,47,51]. However, our understanding of the mechanism underlying gene rearrangement within Halictidae remains limited, largely due to the limited samples.

Contrasted with the putative ancestral gene rearrangement of insects, all four newly sequenced mitochondrial genomes presented gene rearrangement (Figure 3). We found five gene rearrangements in the newly sequenced species, and all rearranged genes were tRNA genes (Figure 3). In this study, COX2-trnD-trnK and ND3-trnR-trnA-trnN were observed in all newly sequenced species. trnM-trnI-trnQ was also found in most species of our study, except L. guihongi. Gene shuffling (trnQ/trnM, trnW/trnC-trnY, and trnK/trnD) was found to be a significant rearrangement pattern in bees, which had been proposed as an important type of gene rearrangement in Hymenoptera [52,53,54]. The gene rearrangement CR-trnI-trnQ-trnM-nad2-trnW-trnC-trnY was the main hot spot in bees, and our study showed that the gene rearrangement hot spot of Halictidae also exists in this region (Figure 3). In addition, this study described a special rearrangement of the ND4-trnP-ND4L-trnT in L. guihongi (Figure 3), which was the first instance of mitochondrial gene rearrangement found in bees, providing a foundation for research on Halictidae biological characteristics. In theory, 37 genes of the mitochondrial genome have considerable potential for rearrangement. Therefore, different species are unlikely to adopt the same gene sequence. The shared gene sequence is more likely to indicate a common evolutionary history [17,46,47]. Most gene rearrangement patters in this study were observed in all newly sequenced species, showing that those gene rearrangements in tRNAs may be a common feature within the Halictidae.

4.3. Phylogenetic Relationship Analyses

Herein, we used mitochondrial genomes to reconstruct the relationships within Halictidae. Five matrices based on BI and ML approaches yielded four different topologies; all topologies supported that the Rophitinae was the basal linage of Halictidae (Figure 4 and Figures S3–S9), which was similar to the results of Peskon and Danforth [55,56,57]. The former used morphological characteristics, while the latter used nuclear genes to reconstruct the phylogenetic relationships within Halictidae 15 years ago [55,56,57].

Halictidae is an important bee family, not only because its members contribute to pollination, but also because of their social behavior [4]. However, there have been relatively few studies on the phylogeny of this group in recent years, especially for Nomiinae. Previous studies of Nomiinae mainly focused on morphology. The molecular data are seriously inadequate, and the relationships in this group have led to great controversy based on morphological characters. Michener divided Nomiinae into 11 genera and 19 subgenera, while Pauly treated almost all subgenera as independent genera [3,4,10,58,59]. Pauly’s classification perspectives have not been widely supported by other taxonomists. In this study, all topologies showed that Austronomia was not a subgenus within the genus Lipotriches, and Gnathonomia may also be treated as a genus rather than as a subgenus of Nomia (Figure 4 and Figures S3–S9), which closely resembles the generic system of Pauly [5]. Unfortunately, due to limited samples and phylogenetic signal markers, we have not yet resolved the phylogenetic relationships of Nomiinae. In the future, more samples of genera and subgenera will be included to explore the phylogenetic relationships within this group, and we will use more efficient data, including low-coverage whole-genome genome and ultra-conserved element data to resolve the relationships within Nomiinae, which show great promise in phylogenetic research [60,61,62].

5. Conclusions

In this study, we newly sequenced, assembled, and annotated four mitogenomes of Halictidae. Compared to previously published data, all newly obtained mitogenomes had similar structural characters and nucleotide compositions. In addition, gene rearrangement analysis found one novel model: ND4-trnP-ND4L-trnT. Thus, our study provided more available mitogenomes for the study of Nomiinae and Halictidae. In adding two Halictinae species, we reconstructed the phylogenetic relationships within Halictidae based on BI and ML methods. All topologies showed that Rophitinae is the basal group within Halictidae. Our research provides a molecular database for research on Halictidae. Further study, including the addition of more samples and using more phylogenetic signal markers (e.g. UCEs, whole-genome sequencing), are needed to define the phylogenetic relationships within the subfamily and the evolutionary biology of Halictidae.