# Genetic Diversity ofApis spp. in Thailand Inferred from 28S rRNA Nuclear and Cytochrome b Mitochondrial Gene Sequences

Authors: Thitipan Meemongkolkiat; Atsalek Rattanawannee; Chanpen Chanchao
Journal: Psyche
 (2019)

---

## Abstract

Knowledge of the genetic diversity ofApisspp. is important in order to provide a better understanding of breeding strategies that relate to the conservation of wild species and colony survival of farmed species. Here, honeybees of fiveApisspecies were collected from 12 provinces throughout Thailand. After DNA extraction,28S rRNA nuclear (710 bp) and cytochrome b (cytb) mitochondrial (520 bp) gene fragments were sequenced. Homologous sequences (nucleotide identity of over 95%) were obtained from GeneBank using the BLASTn algorithm, aligned, and analysed by maximum likelihood and Bayesian inference phylogenetics. For28S rRNA, a low genetic variation was detected within species (haplotype diversity ranging from 0.212 to 0.394), while 19 polymorphic sites were detected between species. Although the relative haplotype diversity was high, a low nucleotide divergence was found (0.7%), with migratory species. Forcytb, the sequence divergence ranged from 0.24 to 3.88% within species and 7.35 to 13.07% between species. The divergence ofcytb was higher than that of28S rRNA.A. cerana showed two distinct clades between Southern Thailand and the other regions. Groups ofA. cerana (Asian cavity-nesting),A. mellifera (European cavity-nesting),A. dorsata (giant open-nesting), andA. florea andA. andreniformis (dwarf bees) were defined in the28S rRNA andcytb tree topologies.


---

## Body

## 1. Introduction

The diversity of honeybees (Apis spp.) is important in agricultural areas because they play a crucial role in pollinating a wide variety of crop flowers [1]. Large forest areas have been increasingly fragmented for agriculture and other land management, with an increase in larger monoculture and industrial areas. Furthermore, climate change has become more serious. These changes have likely affected the diversity ofApis spp. since each species has different food plant and habitat preferences.The diversity ofApis spp. in Thailand has been reported on over the past few decades. Morphological variation was reported inA. florea [2],A. andreniformis [3], andA. dorsata [4], as well as the genetic variation ofA. dorsata andA. mellifera [5, 6]. Studies on the genetic variation ofApis have used mitochondrial DNA gene fragments [7, 8], single nucleotide polymorphisms [9], and allozymes [10] to differentiate the breeds and provide a more powerful approach to determine the degree of genetic variation among and between honeybee populations.Mitochondrial genes are the most popular molecular markers used to determine the degree of genetic diversity in insects [11] since the mitochondrial genome typically shows higher average mutation rates (10x) than the nuclear genome [12] due to nucleotide imbalance [13] and the small effective population size associated with effective haploid inheritance across eukaryotic taxa [14]. For eusocial insects, the polyandrous mating system seems to insert the higher substitution rates in mitochondrial genes than in nuclear genes since the mitochondrial DNA is inherited exclusively through the maternal lineage [15]. Thus, the inheritance of mitochondrial DNA is transmitted via only one queen, while the nuclear genome is derived from the parental line. This leads to a bias between the effective size of the nuclear and mitochondrial DNA pools [16].However, the use of only mitochondrial DNA as a genetic marker to resolve phylogenetic studies can be insufficient and possibly lead to erroneous results. Accordingly, the phylogenetic relationships between nuclear and mitochondrial genomes have been compared and combined [17]. Nuclear gene introns appear to be acceptable for this analysis. Using both maximum likelihood (ML) and neighbour-joining analyses of the mitochondrialNADH dehydrogenase subunit 2(ND2) and nuclearelongation factor 1-α (EF1-α) intron ofApis spp. produced the same basic topologies ofApis spp. [18].Although the resolution of such a circumscribed taxonomic group asApis spp. has been reported, uncertainty remains concerning important aspects of the phylogenetic relationships within honeybees, and more data are required to reconfirm the status of honeybee diversity in Thailand. Thus, this work aimed to investigate the level of genetic variation and phylogenetic relationships within and amongApis spp. populations in Thailand using partial sequence data from the nuclear28S rRNA region compared with the mitochondrialcytochrome b gene (cytb). These data will provide complementary information and a more satisfactory comprehension of the phylogenetic variation and the origin within and amongApis spp. on an evolutionary time scale.


## 2. Materials and Methods

### 2.1. Sample Collection

A total of 76 colonies of honeybees, comprisingApis cerana (26),A. mellifera (11),A. dorsata (15),A. florea (22) andA. andreniformis (2) were collected from 12 provinces in eastern, central, northern, western, and southern Thailand, ranging from 8° 46′ 04.4′′ N to 18° 45′ 34.1′′ N latitude and 98°57′ 56.3′′ E to 102° 27′ 09.0′′ E longitude during 2016 and 2017 (Figure 1 and Supplement 1). The samples were frozen with dry ice and stored at -20°C in the laboratory until used.Figure 1
Sampling sites ofApis spp. in Thailand.


### 2.2. Partial Nucleotide Sequences of the 28S rRNA and Cytb Genes

Total genomic DNA was extracted from the thorax using a QIAamp® DNA Mini Kit (catalog # 51304, Qiagen). The28S rRNA andcytb gene fragments were amplified using the polymerase chain reaction (PCR) with the specific28S rRNA (forward: 5′-AGAGAGAGTTCAAGAGTACGTG-3′ and reverse: 5′-TAGTTCACCATCTTTCGGGTCCC-3′) andcytb (forward: 5′-TGAAATTTTGGATCAATTCTTGG -3′ and reverse: 5′- TCCAAGAGGATTAGATGATCCAG -3′) primer pairs [3, 19]. Each PCR amplification was performed in a 25 μl final reaction volume comprising 12.5 μl 2X EmeraldAmp® PCR master mix (catalog # RR300A, Takara), 1 μl of each of the 10 μM forward and reverse primer, at least 30 ng of genomic DNA template and nuclease-free distilled H2O. The thermocycling condition was 94°C for 2 min 30 sec followed by 35 cycles of 94°C for 1 min, x  °C for 1 min (x = 65°C for28S rRNA and 50°C forcytb) and 72°C for 1 min and then a final 72°C for 10 min. After electrophoresis on a 2% (w/v) agarose gel at 80 V for 30 min, PCR products were visualized by UV transillumination after ethidium bromide staining. The expected product sizes of the amplified28S rRNA andcytbfragments were 710 and 520 bp, respectively. The amplified DNA fragments were purified using a QIAquick Gel Purification Kit (catalog # 28704, Qiagen) and then commercially direct sequenced at Macrogen Inc. (Korea). The obtained sequences were used to search for homologous sequences in the National Centre for Biotechnology Information (NCBI) GenBank database using the BLASTn algorithm and aligned.


### 2.3. Sequences Analysis and Diversity Indices

The MEGA7 program version 7.0 [20] was used to generate the multiple sequence alignments and to obtain the number of substitutions between sequences, nucleotide compositions, and the transition or transversion ratios for the28S rRNA andcytb sequences. The genetic variation within and between species was determined as the number of haplotypes (No.), number of polymorphic sites (s), haplotype diversity (h) [21], average number of nucleotide differences (k), and nucleotide diversity (π) [22], using the DNAsp V.5.0 program [23]. The sequence divergence (%) with Jukes and Cantor [24] between and within species and the average number of nucleotide differences were estimated using the MEGA7 V7.0 and ARLEQUIN ver 3.5.2.2 software [25].


### 2.4. Phylogenetic Analysis

The ML and Bayesian inference (BI) methods were used to reconstruct the phylogenetic relationships in the28S rRNA andcytb sequences. The best-fit models of nucleotide substitution were determined using the Akaike information criterion (AIC) for ML [26] and Bayesian information criterion (BIC) for BI [27] as implemented in the Kakusan 4 program [28]. The ML analysis was performed in Treefinder with 1000 bootstrap replicates to estimate branch confidence values [29]. Tree topologies with bootstrap values (BV) of 70% or greater were regarded as sufficiently resolved. The BI was performed using the MrBayes V3.1 program [30] using the selected model with four simultaneous chains run for at least 1,000,000 generations starting from a random tree, with tree sampling every 100 generations. Both ML and BI gave a consensus of the phylogenetic tree.


## 2.1. Sample Collection

A total of 76 colonies of honeybees, comprisingApis cerana (26),A. mellifera (11),A. dorsata (15),A. florea (22) andA. andreniformis (2) were collected from 12 provinces in eastern, central, northern, western, and southern Thailand, ranging from 8° 46′ 04.4′′ N to 18° 45′ 34.1′′ N latitude and 98°57′ 56.3′′ E to 102° 27′ 09.0′′ E longitude during 2016 and 2017 (Figure 1 and Supplement 1). The samples were frozen with dry ice and stored at -20°C in the laboratory until used.Figure 1
Sampling sites ofApis spp. in Thailand.


## 2.2. Partial Nucleotide Sequences of the 28S rRNA and Cytb Genes

Total genomic DNA was extracted from the thorax using a QIAamp® DNA Mini Kit (catalog # 51304, Qiagen). The28S rRNA andcytb gene fragments were amplified using the polymerase chain reaction (PCR) with the specific28S rRNA (forward: 5′-AGAGAGAGTTCAAGAGTACGTG-3′ and reverse: 5′-TAGTTCACCATCTTTCGGGTCCC-3′) andcytb (forward: 5′-TGAAATTTTGGATCAATTCTTGG -3′ and reverse: 5′- TCCAAGAGGATTAGATGATCCAG -3′) primer pairs [3, 19]. Each PCR amplification was performed in a 25 μl final reaction volume comprising 12.5 μl 2X EmeraldAmp® PCR master mix (catalog # RR300A, Takara), 1 μl of each of the 10 μM forward and reverse primer, at least 30 ng of genomic DNA template and nuclease-free distilled H2O. The thermocycling condition was 94°C for 2 min 30 sec followed by 35 cycles of 94°C for 1 min, x  °C for 1 min (x = 65°C for28S rRNA and 50°C forcytb) and 72°C for 1 min and then a final 72°C for 10 min. After electrophoresis on a 2% (w/v) agarose gel at 80 V for 30 min, PCR products were visualized by UV transillumination after ethidium bromide staining. The expected product sizes of the amplified28S rRNA andcytbfragments were 710 and 520 bp, respectively. The amplified DNA fragments were purified using a QIAquick Gel Purification Kit (catalog # 28704, Qiagen) and then commercially direct sequenced at Macrogen Inc. (Korea). The obtained sequences were used to search for homologous sequences in the National Centre for Biotechnology Information (NCBI) GenBank database using the BLASTn algorithm and aligned.


## 2.3. Sequences Analysis and Diversity Indices

The MEGA7 program version 7.0 [20] was used to generate the multiple sequence alignments and to obtain the number of substitutions between sequences, nucleotide compositions, and the transition or transversion ratios for the28S rRNA andcytb sequences. The genetic variation within and between species was determined as the number of haplotypes (No.), number of polymorphic sites (s), haplotype diversity (h) [21], average number of nucleotide differences (k), and nucleotide diversity (π) [22], using the DNAsp V.5.0 program [23]. The sequence divergence (%) with Jukes and Cantor [24] between and within species and the average number of nucleotide differences were estimated using the MEGA7 V7.0 and ARLEQUIN ver 3.5.2.2 software [25].


## 2.4. Phylogenetic Analysis

The ML and Bayesian inference (BI) methods were used to reconstruct the phylogenetic relationships in the28S rRNA andcytb sequences. The best-fit models of nucleotide substitution were determined using the Akaike information criterion (AIC) for ML [26] and Bayesian information criterion (BIC) for BI [27] as implemented in the Kakusan 4 program [28]. The ML analysis was performed in Treefinder with 1000 bootstrap replicates to estimate branch confidence values [29]. Tree topologies with bootstrap values (BV) of 70% or greater were regarded as sufficiently resolved. The BI was performed using the MrBayes V3.1 program [30] using the selected model with four simultaneous chains run for at least 1,000,000 generations starting from a random tree, with tree sampling every 100 generations. Both ML and BI gave a consensus of the phylogenetic tree.


## 3. Results

### 3.1. Molecular Data and Sequence Analysis

A total of 106 sequences (73 sequences from28S rRNA and 33 sequences fromcytb) were obtained and have been deposited in GenBank with the accession codes as listed in Supplement 1. The average nucleotide composition of the28S rRNA sequence did not vary much between theApisspecies, with an across all species average of A (17.33%), T (21.05%), G (31.84%), and C (29.78%), giving a high GC content of 61.6% (Supplement 2). The averagecytb sequence composition varied slightly between species, with an average between all species of A (34.2%), T (42.5%), G (11.8%), and C (11.4%), with a high A+T content of 77.2% (Supplements 2 and 3).For the28S rRNA, 75 individuals were analysed and a total of 612 positions were analysed in the final dataset. Transversion sites were less frequent than transitions (Table 1). Forcytb, 33 sequences were analysed, including the 1st, 2nd, and 3rd codons and noncoding positions separately, giving a total of 325 positions in the final dataset with 76 polymorphic sites, of which 27 sites showed transitions, 42 sites showed transversions, and seven sites carried both types of base substitution (Table 1).Table 1
Frequency (%) of transitions and transversions of28S rRNA and cytb.
Species(No. of sequences) Transitions Transversions A/G T/C A/T A/C G/T G/C  28S rRNA A. cerana (26)∗ 0 0 19.34 23.94 26.07 30.68 A. florea (22)∗ 34.21 0 12.74 15.41 17.48 20.15 A. dorsata (15)∗ 52.79 0 9.13 11.12 12.49 14.48 A. mellifera (11)∗ 42.72 0 10.94 13.68 14.93 17.67 All species (75)∗∗(includingA. andreniformis) 32.69 29.73 7.22 8.85 9.94 11.57  cytb A. cerana (10)∗ 8.65 79.41 4.53 2.68 3.29 1.44 A. florea (11)∗ 0 0 38.30 23.58 26.42 11.70 A. mellifera (9)∗ 0.08 99.54 0.14 0.08 0.10 0.04 All species (33)∗∗(including A. andreniformis and A. dorsata) 25.98 28.56 17.28 10.35 12.39 5.46 Remark: “∗” and “∗∗” show the nucleotide compositions within and among, respectively,Apis species from the different regions in Thailand.Based on the sequence alignment within species, the28S rRNA revealed low genetic diversity indices and low genetic variation among AsianApisspp. (Tables 2 and 3). The sequence divergence (pairwise comparison) over the28S rRNA withinA. cerana,A. florea,andA. dorsata was derived from 612, 638, and 627 bp, respectively, using 26, 22, and 15 sequences from across Thailand (Table 3). The sequence alignment between species revealed 19 polymorphic sites comprising seven A/G transitions, six T/C transitions, two A/T transversions, three G/T transversions, and one G/C transversion. The between species transition and transversion mutations in the28S rRNA sequence were 62.42% and 37.58%, respectively. Our results also demonstrated that the intraspecific sequence divergence was lower than the interspecific sequence divergence (Table 4). Forcytb, only one and two sequences were obtained forA. dorsata andA. andreniformis, respectively. The number of base differences per sequence from 385 positions is shown in Table 4, revealing 33 different sequence variants in five species. The sequence alignment between species revealed 76 polymorphic sites, of which there were 11 A/G transitions, 20 T/C transitions, 37 A/T transversions, one G/T transition, one G/C transition, and seven sites that carried both types of base substitution. The transition and transversion mutations incytb sequences between species were 54.54% and 55.46%, respectively. The mean sequence divergence for all the variants was 8.9%. In general, the sequence divergence was high among species, ranging from 6.5 to 25.2%, and lower between individuals of the same species.Table 2
Molecular diversity indices within the28S rRNA and cytb sequences in each Apis species in Thailand.
Species(No. of sequences) No. S h(± SD) k π ± SD)  28S rRNA A. cerana (26) 2 1 0.212 ± 0.097 0.212 0.00033 ± 0.00041 A. florea (22) 3 2 0.394 ± 0.119 0.541 0.00085 ± 0.00086 A. dorsata (15) 4 10 0.371 ± 0.153 1.333 0.00213 ± 0.00491  Cytb A. cerana (10) 8 16 0.933 ± 0.077 5.956 0.0154 ± 0.0147 A. florea (11) 3 3 0.473 ± 0.162 1.200 0.00201 ± 0.00250 Number of haplotypes (No.), number of polymorphic sites (s), haplotype diversity (h) [21], average number of nucleotide differences (k), and nucleotide diversity (π) [22] with standard deviation (SD).Table 3
Estimated sequence divergence over the28S rRNA sequence pairs within each species. Sequence divergence (%) within species is above the diagonal, while the average number of nucleotide differences is below the diagonal.
A. cerana Central (Samutsongkram) West (Kanchanaburi, Ratchaburi) North (Chiangmai, Nan) East (Chanthaburi) South (Chumphon, Suratthani, Nakhon si thammarat)  Central - 0.03 0.07 0.04 0.04 West 0.25 - 0.06 0 0 North 0.45 0.40 - 0.06 0.06 East 0.25 0 0.40 - 0 South 0.25 0 0.40 0 -  A. florea Central (Samutsongkram, Nakhonsawan) West (Kanchanaburi, Ratchaburi) East (Chanthaburi) South (Chumphon) North (Pichit)  Central - 0.08 0.03 0.03 0.01 West 0.57 - 0.07 0.07 0.01 East 0.20 0.50 - 0.00 0.09 South 0.20 0.50 0.00 - 0.09 North 0.64 0.77 0.60 0.60 -  A. dorsata Central (Samutsongkram, Nakhonsawan) West (Kanchanaburi) North(Nan) South (Suratthani)  Central - 0.24 0.29 0.24 West 1.50 - 0.05 0.00 North 1.83 0.33 - 0.05 South 1.50 0.00 0.33 -Table 4
Estimated sequence divergence in the28S rRNA and cytb sequence pairs between species. The pairwise comparison uses 612 and 383 bp of 28SrRNA and cytb sequence, respectively.
28S rRNA A. cerana A. mellifera A. dorsata A. florea A. andreniformis  A. cerana - 0.55 0.61 1.20 1.17 A. mellifera 3.39 - 0.81 1.05 1.05 A. dorsata 3.76 4.98 - 1.13 1.09 A. florea 7.29 6.42 6.86 - 0.03 A. andreniformis 7.11 6.36 6.67 0.23 -  cytb A. mellifera A. andreniformis A. dorsata A. cerana A. florea  A. mellifera - 11.48 10.58 10.22 13.07 A. andreniformis 34.61 - 12.89 12.49 7.35 A. dorsata 34.00 38.50 - 12.59 13.45 A. cerana 32.50 37.40 38.30 - 12.19 A. florea 39.00 22.77 42.00 36.57 - Sequence divergence (%) within species is above the diagonal, while the average number of nucleotide differences is below the diagonal.Based on thecytb sequence alignment, thecytb analysis ofA. cerana (n = 10) revealed eight haplotypes from 16 polymorphic sites that contained four A/G transitions, eight T/C transitions, and four A/T transversions (Supplement 4). TheA. cerana populations from the south of Thailand showed a clear sequence divergence from the other regions of Thailand. The sequence divergence (%) within species and the average number of nucleotide differences observed in pairwise comparison amongA. cerana are shown in Table 5. The highest sequence divergence was detected when HAP1–HAP6 (includingA. cerana from east, central, and west regions) was compared with HAP7 and HAP8 (from south regions) (Supplement 4). These results suggest thatA. cerana in southern Thailand were from different populations than those in the other regions of Thailand.Table 5
Estimated sequence divergence overcytb sequence pairs (383 bp) within A. cerana (10 samples) and A. florea (11 samples).
A. cerana Chantaburi (East) Kanchanaburi (West) Samutsongkram (Central) Suratthani (South)  Chantaburi (East) - 0.84 1.05 3.88 Kanchanaburi (West) 3.22 - 0.26 3.06 Samutsongkram (Central) 4.00 1.00 - 2.79 Suratthani (South) 14.50 11.50 10.5 -  A. florea Chantaburi (East) Chumphon (South) Pichit (North) Ratchaburi (West) Samutsongkram (Central)  Chantaburi (East) - 0 0.24 0 0.74 Chumphon (South) 0 - 0.24 0 0.74 Pichit(North) 1.00 1.00 - 0.24 0.49 Ratchaburi (West) 0 0 1.00 - 0.74 Samutsongkram (Central) 3.00 3.00 2.00 3.00 - Sequence divergence (%) with Jukes and Cantor between species is above the diagonal and the average number of nucleotide differences is below the diagonal.InA. florea, the pairwise comparison of thecytbsequences suggested thatA. florea from two colonies from SA (central regions) were slightly different from the other regions (Table 5). Interestingly, the sequence divergence amongA. cerana andA. mellifera (10.22%) was lower than that amongA. dorsata,A. florea, andA. andreniformis (Table 4). Meanwhile,A. florea showed less sequence divergence fromA. andreniformis (7.35%) than fromA. cerana andA. mellifera (12.19% and 13.07%, respectively).


### 3.2. Phylogenetic Analysis

The best-fit model of evolution for constructing the phylogenetic tree by ML analysis was found by the AIC to be TVM + gamma, while that for the BI tree under the BIC was HKY85 + homogeneous. The ML trees constructed for the28S rRNA andcytb gene sequences were essentially the same as the corresponding BI tree and so only the ML tree topologies are illustrated here (Figures 2 and 3).Figure 2
Phylogenetic relationship of honeybee species in Thailand based upon the 28S rRNA sequences by ML analysis.Trigona collina (accession # FJ042138) was utilized as the outgroup. Node support inferred from Bayesian PP and ML BV is shown. Sample (species/geographic location) codes are as in Supplement 1.Figure 3
Phylogenetic relationship of honeybee species in Thailand based upon cytb sequences by BI and includingApisspp. from GenBank (accession # KC294229, EU040176, EU040168, EU040165, JX982136, KC170303, AP017941, AP017942, and GQ162109).Trigona collina (accession # AY575081) was utilized as the outgroup. Node support as Bayesian PP and ML BV is shown. Sample (species/geographic location) codes are as in Supplement 1.For the28S rRNA sequences (612–638 bp), the phylogenetic relationships withinA. cerana,A. florea, andA. dorsata from all geographical regions showed no evidence of geographical clustering in the phylogenetic tree, although the nonnativeA. mellifera showed some divergence within the species. Overall, a low genetic variation within all AsianApisspp. was noticed based on the 75 sequences of this study plus two sequences ofA. andreniformis from GenBank and usingTrigona collina as an outgroup (Figure 2). TheApisspp. were subdivided into four major clades (BV range from 50 to 98%, PP range from 0.60 to 0.95), strongly supporting the monophyly ofApisspecies in Thailand (BV = 100%, PP = 1.0). Interestingly, all fourApis species were clearly grouped into their respective species clade, except for the two dwarf honeybee species,A. florea andA. andreniformis, which were not completely resolved (Figure 2). This reflects the low sequence divergence between these dwarf species, which were basal in theApis groups (Figure 2).For thecytb region (383 bp in length), the ML phylogenetic dendrograms based on 33 sequences from this study from allApis species including only one and two sequences ofA. dorsata andA. andreniformis, respectively, plus five sequences ofApisspp. from GenBank were constructed withT. collina as the outgroup. The obtained tree revealed four major clades (Figure 3), where eachApis species was clearly separated from the others, includingA. florea andA. andreniformis. AlthoughA. florea in Thailand was divided into small subgroups with three haplotypes in the population from Samutsongkram (central region), there was a low level of genetic variation incytb in theA. florea populations. The basal group wasA. dorsata. ForA. cerana (clade 1), the 10A. cerana samples in Thailand of this study plus threeA. cerana sequences from GenBank were subdivided into the two large groups of Southern Thailand and the other regions, plus one subgroup in the Eastern region (Figure 3). TheA. cerana populations from Suratthani were clearly separated from the other regions (BV = 100%; PP = 0.95), while the samples from Eastern Thailand showed a slight difference from the other regions although this had only weak support (BV = 64%). This low support might reflect the nucleotide differences among the haplotypes (Supplement 4) that are still placed in the same clade. These results suggested thatA. cerana in Southern Thailand are from a different population to those in the other regions of Thailand.


## 3.1. Molecular Data and Sequence Analysis

A total of 106 sequences (73 sequences from28S rRNA and 33 sequences fromcytb) were obtained and have been deposited in GenBank with the accession codes as listed in Supplement 1. The average nucleotide composition of the28S rRNA sequence did not vary much between theApisspecies, with an across all species average of A (17.33%), T (21.05%), G (31.84%), and C (29.78%), giving a high GC content of 61.6% (Supplement 2). The averagecytb sequence composition varied slightly between species, with an average between all species of A (34.2%), T (42.5%), G (11.8%), and C (11.4%), with a high A+T content of 77.2% (Supplements 2 and 3).For the28S rRNA, 75 individuals were analysed and a total of 612 positions were analysed in the final dataset. Transversion sites were less frequent than transitions (Table 1). Forcytb, 33 sequences were analysed, including the 1st, 2nd, and 3rd codons and noncoding positions separately, giving a total of 325 positions in the final dataset with 76 polymorphic sites, of which 27 sites showed transitions, 42 sites showed transversions, and seven sites carried both types of base substitution (Table 1).Table 1
Frequency (%) of transitions and transversions of28S rRNA and cytb.
Species(No. of sequences) Transitions Transversions A/G T/C A/T A/C G/T G/C  28S rRNA A. cerana (26)∗ 0 0 19.34 23.94 26.07 30.68 A. florea (22)∗ 34.21 0 12.74 15.41 17.48 20.15 A. dorsata (15)∗ 52.79 0 9.13 11.12 12.49 14.48 A. mellifera (11)∗ 42.72 0 10.94 13.68 14.93 17.67 All species (75)∗∗(includingA. andreniformis) 32.69 29.73 7.22 8.85 9.94 11.57  cytb A. cerana (10)∗ 8.65 79.41 4.53 2.68 3.29 1.44 A. florea (11)∗ 0 0 38.30 23.58 26.42 11.70 A. mellifera (9)∗ 0.08 99.54 0.14 0.08 0.10 0.04 All species (33)∗∗(including A. andreniformis and A. dorsata) 25.98 28.56 17.28 10.35 12.39 5.46 Remark: “∗” and “∗∗” show the nucleotide compositions within and among, respectively,Apis species from the different regions in Thailand.Based on the sequence alignment within species, the28S rRNA revealed low genetic diversity indices and low genetic variation among AsianApisspp. (Tables 2 and 3). The sequence divergence (pairwise comparison) over the28S rRNA withinA. cerana,A. florea,andA. dorsata was derived from 612, 638, and 627 bp, respectively, using 26, 22, and 15 sequences from across Thailand (Table 3). The sequence alignment between species revealed 19 polymorphic sites comprising seven A/G transitions, six T/C transitions, two A/T transversions, three G/T transversions, and one G/C transversion. The between species transition and transversion mutations in the28S rRNA sequence were 62.42% and 37.58%, respectively. Our results also demonstrated that the intraspecific sequence divergence was lower than the interspecific sequence divergence (Table 4). Forcytb, only one and two sequences were obtained forA. dorsata andA. andreniformis, respectively. The number of base differences per sequence from 385 positions is shown in Table 4, revealing 33 different sequence variants in five species. The sequence alignment between species revealed 76 polymorphic sites, of which there were 11 A/G transitions, 20 T/C transitions, 37 A/T transversions, one G/T transition, one G/C transition, and seven sites that carried both types of base substitution. The transition and transversion mutations incytb sequences between species were 54.54% and 55.46%, respectively. The mean sequence divergence for all the variants was 8.9%. In general, the sequence divergence was high among species, ranging from 6.5 to 25.2%, and lower between individuals of the same species.Table 2
Molecular diversity indices within the28S rRNA and cytb sequences in each Apis species in Thailand.
Species(No. of sequences) No. S h(± SD) k π ± SD)  28S rRNA A. cerana (26) 2 1 0.212 ± 0.097 0.212 0.00033 ± 0.00041 A. florea (22) 3 2 0.394 ± 0.119 0.541 0.00085 ± 0.00086 A. dorsata (15) 4 10 0.371 ± 0.153 1.333 0.00213 ± 0.00491  Cytb A. cerana (10) 8 16 0.933 ± 0.077 5.956 0.0154 ± 0.0147 A. florea (11) 3 3 0.473 ± 0.162 1.200 0.00201 ± 0.00250 Number of haplotypes (No.), number of polymorphic sites (s), haplotype diversity (h) [21], average number of nucleotide differences (k), and nucleotide diversity (π) [22] with standard deviation (SD).Table 3
Estimated sequence divergence over the28S rRNA sequence pairs within each species. Sequence divergence (%) within species is above the diagonal, while the average number of nucleotide differences is below the diagonal.
A. cerana Central (Samutsongkram) West (Kanchanaburi, Ratchaburi) North (Chiangmai, Nan) East (Chanthaburi) South (Chumphon, Suratthani, Nakhon si thammarat)  Central - 0.03 0.07 0.04 0.04 West 0.25 - 0.06 0 0 North 0.45 0.40 - 0.06 0.06 East 0.25 0 0.40 - 0 South 0.25 0 0.40 0 -  A. florea Central (Samutsongkram, Nakhonsawan) West (Kanchanaburi, Ratchaburi) East (Chanthaburi) South (Chumphon) North (Pichit)  Central - 0.08 0.03 0.03 0.01 West 0.57 - 0.07 0.07 0.01 East 0.20 0.50 - 0.00 0.09 South 0.20 0.50 0.00 - 0.09 North 0.64 0.77 0.60 0.60 -  A. dorsata Central (Samutsongkram, Nakhonsawan) West (Kanchanaburi) North(Nan) South (Suratthani)  Central - 0.24 0.29 0.24 West 1.50 - 0.05 0.00 North 1.83 0.33 - 0.05 South 1.50 0.00 0.33 -Table 4
Estimated sequence divergence in the28S rRNA and cytb sequence pairs between species. The pairwise comparison uses 612 and 383 bp of 28SrRNA and cytb sequence, respectively.
28S rRNA A. cerana A. mellifera A. dorsata A. florea A. andreniformis  A. cerana - 0.55 0.61 1.20 1.17 A. mellifera 3.39 - 0.81 1.05 1.05 A. dorsata 3.76 4.98 - 1.13 1.09 A. florea 7.29 6.42 6.86 - 0.03 A. andreniformis 7.11 6.36 6.67 0.23 -  cytb A. mellifera A. andreniformis A. dorsata A. cerana A. florea  A. mellifera - 11.48 10.58 10.22 13.07 A. andreniformis 34.61 - 12.89 12.49 7.35 A. dorsata 34.00 38.50 - 12.59 13.45 A. cerana 32.50 37.40 38.30 - 12.19 A. florea 39.00 22.77 42.00 36.57 - Sequence divergence (%) within species is above the diagonal, while the average number of nucleotide differences is below the diagonal.Based on thecytb sequence alignment, thecytb analysis ofA. cerana (n = 10) revealed eight haplotypes from 16 polymorphic sites that contained four A/G transitions, eight T/C transitions, and four A/T transversions (Supplement 4). TheA. cerana populations from the south of Thailand showed a clear sequence divergence from the other regions of Thailand. The sequence divergence (%) within species and the average number of nucleotide differences observed in pairwise comparison amongA. cerana are shown in Table 5. The highest sequence divergence was detected when HAP1–HAP6 (includingA. cerana from east, central, and west regions) was compared with HAP7 and HAP8 (from south regions) (Supplement 4). These results suggest thatA. cerana in southern Thailand were from different populations than those in the other regions of Thailand.Table 5
Estimated sequence divergence overcytb sequence pairs (383 bp) within A. cerana (10 samples) and A. florea (11 samples).
A. cerana Chantaburi (East) Kanchanaburi (West) Samutsongkram (Central) Suratthani (South)  Chantaburi (East) - 0.84 1.05 3.88 Kanchanaburi (West) 3.22 - 0.26 3.06 Samutsongkram (Central) 4.00 1.00 - 2.79 Suratthani (South) 14.50 11.50 10.5 -  A. florea Chantaburi (East) Chumphon (South) Pichit (North) Ratchaburi (West) Samutsongkram (Central)  Chantaburi (East) - 0 0.24 0 0.74 Chumphon (South) 0 - 0.24 0 0.74 Pichit(North) 1.00 1.00 - 0.24 0.49 Ratchaburi (West) 0 0 1.00 - 0.74 Samutsongkram (Central) 3.00 3.00 2.00 3.00 - Sequence divergence (%) with Jukes and Cantor between species is above the diagonal and the average number of nucleotide differences is below the diagonal.InA. florea, the pairwise comparison of thecytbsequences suggested thatA. florea from two colonies from SA (central regions) were slightly different from the other regions (Table 5). Interestingly, the sequence divergence amongA. cerana andA. mellifera (10.22%) was lower than that amongA. dorsata,A. florea, andA. andreniformis (Table 4). Meanwhile,A. florea showed less sequence divergence fromA. andreniformis (7.35%) than fromA. cerana andA. mellifera (12.19% and 13.07%, respectively).


## 3.2. Phylogenetic Analysis

The best-fit model of evolution for constructing the phylogenetic tree by ML analysis was found by the AIC to be TVM + gamma, while that for the BI tree under the BIC was HKY85 + homogeneous. The ML trees constructed for the28S rRNA andcytb gene sequences were essentially the same as the corresponding BI tree and so only the ML tree topologies are illustrated here (Figures 2 and 3).Figure 2
Phylogenetic relationship of honeybee species in Thailand based upon the 28S rRNA sequences by ML analysis.Trigona collina (accession # FJ042138) was utilized as the outgroup. Node support inferred from Bayesian PP and ML BV is shown. Sample (species/geographic location) codes are as in Supplement 1.Figure 3
Phylogenetic relationship of honeybee species in Thailand based upon cytb sequences by BI and includingApisspp. from GenBank (accession # KC294229, EU040176, EU040168, EU040165, JX982136, KC170303, AP017941, AP017942, and GQ162109).Trigona collina (accession # AY575081) was utilized as the outgroup. Node support as Bayesian PP and ML BV is shown. Sample (species/geographic location) codes are as in Supplement 1.For the28S rRNA sequences (612–638 bp), the phylogenetic relationships withinA. cerana,A. florea, andA. dorsata from all geographical regions showed no evidence of geographical clustering in the phylogenetic tree, although the nonnativeA. mellifera showed some divergence within the species. Overall, a low genetic variation within all AsianApisspp. was noticed based on the 75 sequences of this study plus two sequences ofA. andreniformis from GenBank and usingTrigona collina as an outgroup (Figure 2). TheApisspp. were subdivided into four major clades (BV range from 50 to 98%, PP range from 0.60 to 0.95), strongly supporting the monophyly ofApisspecies in Thailand (BV = 100%, PP = 1.0). Interestingly, all fourApis species were clearly grouped into their respective species clade, except for the two dwarf honeybee species,A. florea andA. andreniformis, which were not completely resolved (Figure 2). This reflects the low sequence divergence between these dwarf species, which were basal in theApis groups (Figure 2).For thecytb region (383 bp in length), the ML phylogenetic dendrograms based on 33 sequences from this study from allApis species including only one and two sequences ofA. dorsata andA. andreniformis, respectively, plus five sequences ofApisspp. from GenBank were constructed withT. collina as the outgroup. The obtained tree revealed four major clades (Figure 3), where eachApis species was clearly separated from the others, includingA. florea andA. andreniformis. AlthoughA. florea in Thailand was divided into small subgroups with three haplotypes in the population from Samutsongkram (central region), there was a low level of genetic variation incytb in theA. florea populations. The basal group wasA. dorsata. ForA. cerana (clade 1), the 10A. cerana samples in Thailand of this study plus threeA. cerana sequences from GenBank were subdivided into the two large groups of Southern Thailand and the other regions, plus one subgroup in the Eastern region (Figure 3). TheA. cerana populations from Suratthani were clearly separated from the other regions (BV = 100%; PP = 0.95), while the samples from Eastern Thailand showed a slight difference from the other regions although this had only weak support (BV = 64%). This low support might reflect the nucleotide differences among the haplotypes (Supplement 4) that are still placed in the same clade. These results suggested thatA. cerana in Southern Thailand are from a different population to those in the other regions of Thailand.


## 4. Discussion

A total of 7328S rRNA and 33cytb sequences were obtained. The overall sequence diversity of both gene fragments within and among species was low, especially for28S rRNA that showed only a slight sequence divergence for all haplotypes (0.7%). This suggests that there is little variation in28S rRNA withinApisspp. in Thailand. The28S rRNA gene is often used as a representative nuclear gene at the family-species level and above, where, for example, it was reported that the most conserved28S-D3 region could be used to defineEncarsia spp. [31].However, using28S rRNA alone may be insufficient, especially at the species level and below. In addition to28S rRNA, other nuclear genes have been used, likeitpr [32] and theEF1-α intron [18], in order to resolve the family level phylogenies in bees and other insects. Nevertheless, the nuclearEF1-α intron has been reported to be useful for classifying honeybees between species, but not for resolving within species [18]. Likewise, it was reported that the28S rRNA was unlikely to recover deep divergences at the family level in bees [33]. In contrast, in combination with other genes, such asopsin,cytb, and16S rRNA,28S rRNA appeared to contribute to resolving basal bee divergence.In this study, the number of transversions incytb was higher than the number of transitions in eachApisspecies (Table 1). This is in agreement with the transitional bias reported previously [18]. WithinA. cerana, thecytb sequence divergence ranged from 0.26 to 3.88% but, interestingly, was higher in the southern populations than in the other regions of Thailand (Table 5). The phylogenetic results are consistent with the molecular diversity indices ofA. cerana, whereA. cerana populations in southern Thailand were separated from those in the East, West, and Central regions (Figure 3). The result could be explained by the biogeography ofA. cerana. In addition, the data concur with the previous PCR-RFLP analysis [34].ForA. florea, only a very lowcytb sequence divergence was detected (0–0.74%), with only one haplotype. Only two colonies from Samutsongkram province were slightly different from the others (Figure 3).Overall, the sequence divergence ofcytb (8.9 ± 5.0%) was higher than that for28S rRNA (0.7 ± 0.5%), consistent with a higher evolutionary rate for mitochondrial than nuclear genomes [11]. Although the mitochondrialcytbshowed a higher rate of variation, it still showed a low variation in theseApis spp. among the different geographical regions. This may reflect thatApisspp. are highly eusocial and polyandrous, which reduces the substitution rate and counteracts the effect of genetic drift by increasing the effective population size [35]. The ancestral genome in diverging populations will be maintained and so the genetic distances between populations are reduced [36]. The other possible factor leading to a low genetic variation is migration behaviour. The migration distance ofA. dorsata has been reported to be 50 km wide and up to 200 km [37]. InA. cerana, although they have migratory swarms, flying dozens of km [38], they have a strong tendency to swarm 5–6 times a year and are well adapted to new areas since it is more resistant to some diseases and pest problems in each area [39].The ML and BI phylogenetic analyses of the28S rRNA andcytb sequences showed the same topologies with four major clades amongApis spp. (for ML: Figures 2 and 3). UsingT. collina as the outgroup, the28S rRNA tree placed the dwarf bees in a basal position within the genus, whileA. dorsata was located in the basal position in thecytb tree, consistent with previous analysis [40]. Although both thecytb and28S rRNA phylogenetic trees in this study could not clearly resolve the relationship between theseApisspecies, they could separateA. dorsata from others (Figures 2 and 3). A close relationship has been reported betweenA. dorsata andA. cerana in the28S rRNA tree and molecular indices [32]. AlthoughA. cerana andA. mellifera have been reported to be closely related to each other [18], they have ecological and behavioural differences [41]. In this study, the28S rRNAandcytb topologies clearly separatedA. cerana andA. mellifera in Thailand from each other, although they were still closer to each other compared to the otherApis species.Despite the fact thatA. dorsata andA. florea are open-nesting bees,A. dorsata was closer toA. cerana andA. mellifera than toA. florea andA. andreniformis in this study. This can be explained by the evolution of behaviour inApis, where cavity-nesting is secondarily derived [42]. Moreover,A. dorsata has a buzzing dance in daylight and at night, whileA. florea has only a silent waggle dance and cavity-nesting bees adapted this buzzing waggle dance in a dark condition. Thus, a loss of visual cues was compensated for with sound.


## 5. Conclusion

Overall, the genetic variation ofApisspp. from 64 populations within 12 provinces in Thailand showed low genetic diversity indices and variation within and between species. This may indicate that habitat fragmentation, monocrop agriculture, industrial management, and climate changes have not affected the diversity ofApis spp. in Thailand yet, but clearly further work is of vital importance in establishing the actual genetic diversity and relationship withinApis subspecies, including the analysis of other genes for better resolution.



---

*Source: 5823219-2019-02-20.xml*