Comparison of the mitochondrial genomes of three geographical strains of Apis laboriosa indicates high genetic diversity in the black giant honeybee (Hymenoptera: Apidae)

Abstract Apis laboriosa is the largest honeybee that lives mainly on cliff faces, with strong migratory ability. In this study, we firstly sequenced and assembled two complete mitochondrial genomes of A. laboriosa isolated from two distant locations in China (Chongqing and Shangri‐La regions). Combined with the published mitochondrial genome of A. laboriosa from Nepal, comparative genomic analyses were conducted to gain insight into the genetic diversity of giant honeybees from different geographical distributions. The mitochondrial genomes of A. laboriosa from Chongqing and Shangri‐La regions were 15,579 and 15,683 bp in length, respectively, both larger than that from Nepal with the length of 15,510 bp. Three mitochondrial genomes all harbor 37 common genes and present the same AT bias and the frequency of codon usage. However, the fragments including COX1, SSUrRNA, LSUrRNA, and the AT‐rich region of the mitochondrial genome from Shangri‐La region demonstrate distinctive insertions and deletions compared to those from Chongqing and Nepal regions. Phylogenetic trees of mitochondrial genomes show that A. laboriosa from Chongqing is most closely related to that from Nepal, rather than to Shangri‐La. Genetic distance between Shangri‐La and Chongqing or Nepal was even larger than that between the various subspecies of Apis mellifera. Overall, these results unmark that A. laboriosa in different geographical distributions can exhibit high genetic diversity at the mitochondrial genomic level, and therein, A. laboriosa from Shangri‐La may be the subspecies. All these studies will contribute to our understanding of the geographical distribution and genetic differentiation of black giant honeybee in Asian region.


| INTRODUC TI ON
The mitochondria have been widely applied in studying the organism's evolutionary origin and genetic diversity due to its fast evolutionary rate, simplified structure, and effective genetic information (Sun et al., 2010;Wang et al., 2013;Wolstenholme, 1992). As an important pollinating insect, worker bees in honeybee colony descend from one exclusive queen and matches the mitochondrial matrilineal pattern, which were considered as one of the ideal model organisms for studying genetic diversity based on mitochondrial DNA.
Previous researches have found that mitochondrial DNA was differentiated among honeybees of species, subspecies, and geographical groups and can be used to explore the genetic diversity and evolution of honeybee Franck et al., 2000;Hall & Muralidharan, 1989;Smith et al., 1991;Yu et al., 2019).
Apis laboriosa is an important pollinator species usually discovered in the Pan-Himalaya Mountains in southeastern Asia (Smith, 1871). This insect has a high resistance to adversity and strong collection ability, which are features of great economic and scientific value (Batra, 1996;Woyke et al., 2003). However, due to live on inaccessible cliff faces and strong migratory abilities, it is hard for people to collect samples and implement the scientific research.
Up to date, a few studies involved in biology behaviors and morphological comparisons have been conducted on black giant honeybees, which will contribute to the taxonomic classification of A. laboriosa comparing with other honeybee species Kuang et al., 1983;Lo et al., 2010;Underwood, 1990;Woyke et al., 2008).
However, there is no information about morphological comparisons and genetic diversity about intraspecies of this honeybee. Recent study has disclosed the complete mitochondrial genome and whole genome of A. laboriosa (Lin et al., 2021;Takahashi et al., 2018), which offers the reference to deeply studying the genetic diversity of A.

laboriosa.
Previous reports have shown that populations of A. laboriosa have declined in some countries and regions due to habitat destruction, human activities, and bee diseases (Allen et al., 1990;Joshi et al., 2004). Therefore, in-depth research on the genetic diversity of A. laboriosa is conducive to the protection and utilization of the species' scarce resources. In this study, we collected A. laboriosa samples from Jinyun Mountain in Chongqing and Shangri-La in Yunnan province to sequence their mitochondrial genome. After the acquisition of the complete mitochondrial genomes of the honeybees from Chongqing and Shangri-La using Sanger sequencing and assembly. comparative genome analysis was performed to systematically explore the genome variation of A. laboriosa from three geographical origins, and then might be helpful to give insight into geographical distribution and genetic differentiation of this extraordinary giant honeybee.

| Sample collection
Two giant honeybee colonies of A. laboriosa was discovered in Jinyun Mountain, Chongqing (29.62°N, 106.28°E), and Shangri-La, Yunnan (28.56°N, 100.35°E), as shown in Figure 4c, respectively. The 10 adult workers for each colony were collected as samples and then stored in 75% ethanol at 4°C.

| Determination of mitogenome sequences
One worker from each colony was randomly selected and total DNA was extracted from the head and chest muscles by CTAB (cetyltrimethylammonium bromide) method (Shi et al., 2020). Six pairs of specific primers covering the entire sequence of the reference mitochondrial genome were designed for synthesis using NCBI Primer-BLAST tools (https://www.ncbi.nlm.nih.gov/tools/ prime r-blast). The details of the primer sequences used in this study were summarized in Table 1. PCR amplification was carried out using a total reaction volume of 20 μL, which consisted of 1 μL    (Minh et al., 2020;Ronquist et al., 2012), with the bootstrap value set to 1000 replicates and 1000,000 iterations. Also, the p-distance model of MEGA was used for genetic distance analysis.

| Analysis of mitochondrial genome assembly and annotation
As shown in Figure 1, the complete mitochondrial genome lengths of the A. laboriosa from Chongqing and Shangri-La are 15,579 and 15,683 bp, respectively. Both of the genomic lengths were higher than that from Nepal, which was reported as 15,510 bp in length.

| Codon bias analysis of coding genes
We analyzed the frequency and relative synonymous codon usage

TA B L E 2
The annotation information of three complete mitochondrial genomes from three different regions (Chongqing/Shangri-La/ Nepal).

| Analysis of mitochondrial genome structure and sequence variation
We

| Phylogenetic and genetic distance analysis
Two phylogenetic trees based on the mitochondrial genomes and 13 protein-coding genes were constructed using three methods including IQTREE, MEGA, and MrBayes ( Figure 4). Both trees showed that the A. laboriosa from Chongqing, Shangri-La, and Nepal can be assembled into one group and were sistered to another giant honeybee, Apis dorsata. Also, consistent results from two trees showed that Chongqing and Nepal grouped together and form a sister taxon to Shangri-La, suggesting that the A. laboriosa from Chongqing was more closely related to that from Nepal. Pairwise comparison of

F I G U R E 2
The statistics of codon bias of 13 protein-coding genes of A. laboriosa from Chongqing (a), Shangri-La (b), and Nepal (c).

F I G U R E 3
The multiple sequence alignment of A. laboriosa genes from three regions. a: COX1, b: Small subunits of rRNA, c: Large subunits of rRNA, d: AT-rich region.
genetic distance also showed that the genetic distance of Chongqing and Nepal (0.0039) was smaller than Shangri-La and Chongqing (0.0201), or Shangri-La and Nepal (0.0200), supporting the same conclusion that the A. laboriosa from Chongqing and Nepal have a closer relationship (Figure 4d). Significantly, the genetic distances between Shangri-La and Chongqing or between Shangri-La and Nepal, were significantly greater than the genetic distance between various pairwise subspecies of Apis mellifera (0.0027-0.0178), as shown in Figure 4d. Based on the results above, we propose that A.
laboriosa from Shangri-La may be a geographical subspecies.

| DISCUSS ION
Black giant honeybee, A. laboriosa, is usually found in the Himalayas and Hengduan Mountains at an altitude of about 1200-3500 m (Underwood, 1990). Due to specific geographical distribution and rocky seclusion (Roubik et al., 1985;Woyke et al., 2012), field observation and sample collection become considerable difficulty compared with other semidomesticated bees within the Apis genus.
Thus, few studies involved in morphological character, behavior, and genetic difference have been conducted on this unusual honeybee.
In this study, two complete mitochondrial genomes of A. laboriosa from Chongqing and Shangri-La regions were obtained, with the addition of the third published mitochondrial genome of A. laboriosa from Nepal (Takahashi et al., 2018), it offers a great opportunity to illustrate the genetic diversity among different regions of Previous studies show that the AT-rich region is usually responsible for regulating gene duplication and transcription in the mitochondrial genomes of insects, and its sequence length varies considerably among different insect groups and even within genera (Boore, 1999;Cook et al., 2005). Moreover, rRNA is also commonly used for the taxonomic identification of species (Gao et al., 2011;Li et al., 1995;Rougemont et al., 2004). This seems to indicate the special characteristics of the black giant honeybee in Shangri-La. In addition, both full-length mitochondria and 13 PCGs trees showed that Chongqing and Nepal grouped together and form a sister taxon to Shangri-La, and the genetic distances between Shangri-La and Chongqing or between Shangri-La and Nepal, were significantly greater than the genetic distance between various pairwise subspecies of A. mellifera.
Therefore, we speculate that the A. laboriosa from Shangri-La may be a geographical subspecies. Of course, morphological classification will be required so as to provide further evidence for distinguishing subspecies diversification, according to the classic anatomy methods (Mi & Tan, 2007;Ruttner et al., 1978). Due to the limitations in the number of collected bee samples from Shangri-La and Nepal, morphological analysis cannot be carried forward in this study. Previous study (Cao, 2012) (Smith, 1871;Underwood, 1990). However, the location of the black giant honeybee found in Chongqing is only 224 m.
Hence, we presumed that the bee colony in Chongqing may be those honeybees that were margined from high altitudinal regions. More samples need to be collected in future to prove our presumptions.
In conclusion, our results unmark that A. laboriosa in different geographical distributions exhibits high genetic diversity at the mitochondrial genomic level, and therein A. laboriosa from Shangri-La may approach the subspecies level of differentiation. Although more evidence is still needed to support these inferences. Project administration (equal).  Technical System (No. 2022 [12]).

CO N FLI C T O F I NTER E S T S TATEM ENT
We certify that there is no conflict of interest with any financial organization regarding the material discussed in the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
All the newly sequenced genomic data in this study are deposited in GenBank database with accession OP764682 and OP764683.