Abstract

Knowledge of the genetic diversity of Apis spp. 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 five Apis species 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. For 28S 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. For cytb, the sequence divergence ranged from 0.24 to 3.88% within species and 7.35 to 13.07% between species. The divergence of cytb was higher than that of 28S rRNA. A. cerana showed two distinct clades between Southern Thailand and the other regions. Groups of A. cerana (Asian cavity-nesting), A. mellifera (European cavity-nesting), A. dorsata (giant open-nesting), and A. florea and A. andreniformis (dwarf bees) were defined in the 28S rRNA and cytb tree topologies.

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 of Apis spp. since each species has different food plant and habitat preferences.

The diversity of Apis spp. in Thailand has been reported on over the past few decades. Morphological variation was reported in A. florea [2], A. andreniformis [3], and A. dorsata [4], as well as the genetic variation of A. dorsata and A. mellifera [5, 6]. Studies on the genetic variation of Apis 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 mitochondrial NADH dehydrogenase subunit 2 (ND2) and nuclear elongation factor 1-α (EF1-α) intron of Apis spp. produced the same basic topologies of Apis spp. [18].

Although the resolution of such a circumscribed taxonomic group as Apis 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 among Apis spp. populations in Thailand using partial sequence data from the nuclear 28S rRNA region compared with the mitochondrial cytochrome b gene (cytb). These data will provide complementary information and a more satisfactory comprehension of the phylogenetic variation and the origin within and among Apis spp. on an evolutionary time scale.

2. Materials and Methods

2.1. Sample Collection

A total of 76 colonies of honeybees, comprising Apis cerana (26), A. mellifera (11), A. dorsata (15), A. florea (22) and A. 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.

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). The 28S rRNA and cytb gene fragments were amplified using the polymerase chain reaction (PCR) with the specific 28S rRNA (forward: 5′-AGAGAGAGTTCAAGAGTACGTG-3′ and reverse: 5′-TAGTTCACCATCTTTCGGGTCCC-3′) and cytb (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 for 28S rRNA and 50°C for cytb) 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 amplified 28S rRNA and cytb fragments 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 the 28S rRNA and cytb 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 the 28S rRNA and cytb 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 from 28S rRNA and 33 sequences from cytb) were obtained and have been deposited in GenBank with the accession codes as listed in Supplement 1. The average nucleotide composition of the 28S rRNA sequence did not vary much between the Apis species, 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 average cytb 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 the 28S 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). For cytb, 33 sequences were analysed, including the , , and 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).

Based on the sequence alignment within species, the 28S rRNA revealed low genetic diversity indices and low genetic variation among Asian Apis spp. (Tables 2 and 3). The sequence divergence (pairwise comparison) over the 28S rRNA within A. cerana, A. florea, and A. 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 the 28S 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). For cytb, only one and two sequences were obtained for A. dorsata and A. 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 in cytb 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.

Based on the cytb sequence alignment, the cytb analysis of A. 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). The A. 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 among A. cerana are shown in Table 5. The highest sequence divergence was detected when HAP1–HAP6 (including A. cerana from east, central, and west regions) was compared with HAP7 and HAP8 (from south regions) (Supplement 4). These results suggest that A. cerana in southern Thailand were from different populations than those in the other regions of Thailand.

In A. florea, the pairwise comparison of the cytb sequences suggested that A. florea from two colonies from SA (central regions) were slightly different from the other regions (Table 5). Interestingly, the sequence divergence among A. cerana and A. mellifera (10.22%) was lower than that among A. dorsata, A. florea, and A. andreniformis (Table 4). Meanwhile, A. florea showed less sequence divergence from A. andreniformis (7.35%) than from A. cerana and A. 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 the 28S rRNA and cytb 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).

For the 28S rRNA sequences (612–638 bp), the phylogenetic relationships within A. cerana, A. florea, and A. dorsata from all geographical regions showed no evidence of geographical clustering in the phylogenetic tree, although the nonnative A. mellifera showed some divergence within the species. Overall, a low genetic variation within all Asian Apis spp. was noticed based on the 75 sequences of this study plus two sequences of A. andreniformis from GenBank and using Trigona collina as an outgroup (Figure 2). The Apis spp. were subdivided into four major clades (BV range from 50 to 98%, PP range from 0.60 to 0.95), strongly supporting the monophyly of Apis species in Thailand (BV = 100%, PP = 1.0). Interestingly, all four Apis species were clearly grouped into their respective species clade, except for the two dwarf honeybee species, A. florea and A. andreniformis, which were not completely resolved (Figure 2). This reflects the low sequence divergence between these dwarf species, which were basal in the Apis groups (Figure 2).

For the cytb region (383 bp in length), the ML phylogenetic dendrograms based on 33 sequences from this study from all Apis species including only one and two sequences of A. dorsata and A. andreniformis, respectively, plus five sequences of Apis spp. from GenBank were constructed with T. collina as the outgroup. The obtained tree revealed four major clades (Figure 3), where each Apis species was clearly separated from the others, including A. florea and A. andreniformis. Although A. 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 in cytb in the A. florea populations. The basal group was A. dorsata. For A. cerana (clade 1), the 10 A. cerana samples in Thailand of this study plus three A. 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). The A. 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 that A. cerana in Southern Thailand are from a different population to those in the other regions of Thailand.

4. Discussion

A total of 73 28S rRNA and 33 cytb sequences were obtained. The overall sequence diversity of both gene fragments within and among species was low, especially for 28S rRNA that showed only a slight sequence divergence for all haplotypes (0.7%). This suggests that there is little variation in 28S rRNA within Apis spp. in Thailand. The 28S 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 conserved 28S-D3 region could be used to define Encarsia spp. [31].

However, using 28S rRNA alone may be insufficient, especially at the species level and below. In addition to 28S rRNA, other nuclear genes have been used, like itpr [32] and the EF1-α intron [18], in order to resolve the family level phylogenies in bees and other insects. Nevertheless, the nuclear EF1-α intron has been reported to be useful for classifying honeybees between species, but not for resolving within species [18]. Likewise, it was reported that the 28S rRNA was unlikely to recover deep divergences at the family level in bees [33]. In contrast, in combination with other genes, such as opsin, cytb, and 16S rRNA, 28S rRNA appeared to contribute to resolving basal bee divergence.

In this study, the number of transversions in cytb was higher than the number of transitions in each Apis species (Table 1). This is in agreement with the transitional bias reported previously [18]. Within A. cerana, the cytb 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 of A. cerana, where A. 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 of A. cerana. In addition, the data concur with the previous PCR-RFLP analysis [34].

For A. florea, only a very low cytb 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 of cytb (8.9 ± 5.0%) was higher than that for 28S rRNA (0.7 ± 0.5%), consistent with a higher evolutionary rate for mitochondrial than nuclear genomes [11]. Although the mitochondrial cytb showed a higher rate of variation, it still showed a low variation in these Apis spp. among the different geographical regions. This may reflect that Apis spp. 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 of A. dorsata has been reported to be 50 km wide and up to 200 km [37]. In A. 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 the 28S rRNA and cytb sequences showed the same topologies with four major clades among Apis spp. (for ML: Figures 2 and 3). Using T. collina as the outgroup, the 28S rRNA tree placed the dwarf bees in a basal position within the genus, while A. dorsata was located in the basal position in the cytb tree, consistent with previous analysis [40]. Although both the cytb and 28S rRNA phylogenetic trees in this study could not clearly resolve the relationship between these Apis species, they could separate A. dorsata from others (Figures 2 and 3). A close relationship has been reported between A. dorsata and A. cerana in the 28S rRNA tree and molecular indices [32]. Although A. cerana and A. mellifera have been reported to be closely related to each other [18], they have ecological and behavioural differences [41]. In this study, the 28S rRNA and cytb topologies clearly separated A. cerana and A. mellifera in Thailand from each other, although they were still closer to each other compared to the other Apis species.

Despite the fact that A. dorsata and A. florea are open-nesting bees, A. dorsata was closer to A. cerana and A. mellifera than to A. florea and A. andreniformis in this study. This can be explained by the evolution of behaviour in Apis, where cavity-nesting is secondarily derived [42]. Moreover, A. dorsata has a buzzing dance in daylight and at night, while A. 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 of Apis spp. 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 of Apis spp. in Thailand yet, but clearly further work is of vital importance in establishing the actual genetic diversity and relationship within Apis subspecies, including the analysis of other genes for better resolution.

Data Availability

The authors confirm that all data in the manuscript are freely available.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was financially supported by the Science Achievement Scholarship of Thailand and the Anniversary of Chulalongkorn University Fund (Ratchadaphiseksomphot Endowment Fund).

Supplementary Materials

Supplement 1: detail of sample collection and GenBank accession numbers. Supplement 2: frequency (%) of the nucleotide compositions of 28S rRNA and cytb regions. Supplement 3: mean frequency (%) for base compositions at different codon positions in the cytb gene fragment. Supplement 4: the 16 polymorphic sites among eight cytb haplotypes of A. cerana in Thailand. Only positions that are different from haplotype HAP1 are indicated. (Supplementary Materials)