Genetic Relationships and Signatures of Adaptation to the Climatic Conditions in Populations of Apis cerana Based on the Polymorphism of the Gene Vitellogenin

Simple Summary The Oriental honey bee Apis cerana, similar to the Western honey bee Apis mellifera, is distributed in different climatic conditions, and each of them is subdivided into more than thirty subspecies and ecotypes. Their sustainability depends on adaptations to the local climate conditions. The VG gene is involved in the development, reproduction, labor division, and caste differentiation of honey bees. We found the nucleotide sequences of the VG gene reflect the adaptation of honey bees to the local climate conditions. The VG gene sequences are acceptable tools to study the sustainability, genetic structure, and adaptation of A. cerana populations and can be applied in conservation genetics of local honey bee subspecies. Abstract Apis cerana and Apis mellifera are important honey bee species in Asia. A. cerana populations are distributed from a cold, sharply continental climate in the north to a hot, subtropical climate in the south. Due to the Sacbrood virus, almost all A. cerana populations in Asia have declined significantly in recent decades and have recovered over the past five years. This could lead to a shift in the gene pool of local A. cerana populations that could affect their sustainability and adaptation. It was assumed that adaptation of honey bees could be observed by comparative analysis of the sequences of genes involved in development, labor division, and caste differentiation, such as the gene Vitellogenin VG. The VG gene nucleotide sequences were used to assess the genetic structure and signatures of adaptation of local populations of A. cerana from Korea, Russia, Japan, Nepal, and China. A. mellifera samples from India and Poland were used as the outgroup. The signatures of adaptive selection were found in the local population of A. cerana using VG gene sequence analysis based on Jukes–Cantor genetic distances, cluster analysis, dN/dS ratio evaluation, and Tajima’s D neutrality test. Based on analysis of the VG gene sequences, Apis cerana koreana subspecies in the Korean Peninsula were subdivided into three groups in accordance with their geographic localization from north to south. The VG gene sequences are acceptable tools to study the sustainability and adaptation of A. cerana populations.


Introduction
Apis cerana is the second most widely distributed honey bee species in Asia after Apis mellifera [1]. A. cerana is the basic pollinator and producer of apiculture products in Asia [2]. Similar to A. mellifera, A. cerana has high potential for further selection based on genetic markers [3]. The native range of A. cerana embraces northern to southern Asia from Russia to Indonesia and eastern and western Asia from Japan to Afghanistan [4]. A. cerana occupies a large range of climatic conditions at different altitudes and latitudes. Similar to A. mellifera, A. cerana has high genetic and phenotypic variations across a range and is subdivided into several ecotypes and subspecies [3,5,6]. The most northern distribution of A. cerana reaches 47 • 54 N latitude in Khabarovsky Krai of Russia [7]. Each subspecies of A. cerana in Japan (A. c. japonica), in China (A. c. cerana), in Korea (A. c. koreana), in Russia (A. c. ussuriensis), and in Vietnam (A. c. indica) are adapted to their native climate and contain genes linked with adaptation [7,8]. There are no published genetic markers connected with adaptation of A. cerana. All populations of A. cerana are understudied using genetic markers. A. cerana subspecies subdivision and their distribution range are studied insufficiently.
The gene pool of A. cerana in all populations significantly changed after the Sacbrood virus eliminated 95% of colonies [4,9,10]. Development of agriculture and increased treatment with pesticides have led to additional population declines [3,8,9]. The decline in the populations of A. cerana entails a reduction in the biodiversity of local ecosystems. A. cerana and A. mellifera are not equivalent pollinators of native plants; they cannot replace each other [3,4,11].
Within its natural range, the local subspecies of A. cerana has undergone hybridization that violates the adaptive complexes of its gene pool. Additionally, the genetic pool of local populations of subspecies A. cerana may be completely replaced by the genetic pool of introduced subspecies A. cerana. Conservation of the subspecies A. cerana cannot take place without molecular genetic methods. Genetic markers related to adaptation should be found in genes involved in the development and caste differentiation of A. cerana [12,13]. Molecular genetic markers make it possible to develop basic strategies for preserving A. cerana and study its phylogeography.
Vitellogenin, VG, is a basic reproductive protein and endocrine development factor, which, together with juvenile hormone JH and insulin-like growth factor IGF, is responsible for the caste differentiation and development of A. cerana. VG is a source of the royal jelly produced by workers and fed to larvae, so just as it affects queen egg laying, it may also affect the rearing of larvae hatching from those eggs. Moreover, VG plays an important role in regulating the labor division that can impact colony fitness by affecting colony foraging rates, foraging efficiency, and growth rates. The VG gene has high rates of molecular evolution in the honey bee, as evident by high rates of within-A. mellifera polymorphisms and high rates of divergence in the genus Apis. The VG gene is located on the fourth chromosome of A. mellifera, which has the highest recombination rate (27.6 cM ⁄ Mb) among all 16 chromosomes. The recombination rate in the VG gene is 5.9 times higher than in other genes [12][13][14][15][16]. The experiment with RNAi interference of VG gene expression paces feeding behavior and induces shifts in the caste differentiation of A. mellifera [17]. The VG gene is connected to the spatial distribution and population structure of A. mellifera [12,17]. Comparative analysis of VG gene showed adaptation signatures in A. mellifera populations [13]. The VG gene has not been previously studied in A. cerana in the aspects of population genetics.
The VG gene could be useful for the study of the evolutionary adaptation and population structure of A. cerana. Comparative analysis of VG gene sequences promises to yield interesting results regarding the relationship and distribution of local A. cerana populations in Korea, Russia, Japan, Nepal, and China. The signatures of adaptive selection could be found in the local population of A. cerana based on Jukes-Cantor genetic distances, cluster analysis, and Tajima's D neutrality test. The aim of the paper is to study the genetic signatures of climate adaptation in Apis cerana populations using nucleotide sequences of the VG gene. The objectives of the paper are to sequence the VG gene exons and perform comparative cluster analysis with reference sequences from GenBank.
In this paper, the signatures of adaptive evolution were found in the local population of A. cerana using VG gene sequence analysis based on Jukes-Cantor genetic distances, cluster analysis, and Tajima's D neutrality test. The subdivision into three groups in accordance with climate adaptation was found in Apis cerana koreana subspecies population in the Korean Peninsula. Thus, the VG gene sequences are acceptable tools to study the sustainability and adaptation of A. cerana populations.

Materials and Methods
The worker bees of A. cerana were collected from the following five Asian countries: (1) South Korea (provinces Gyeonggido, Gangwondo, Chungcheongbukdo, Gyeongsangbukdo, Jeollabukdo, Gyeongsangnamdo), (2) Japan (prefectures Hiroshima, Hokkaido), (3) Taiwan (Central Taiwan region), (4) Nepal (province Bagmati Pradesh), and (5) Russia (Primorsky Territory) ( Figure 1, Table S1). Worker bees were captured and directly displaced in liquid nitrogen before being stored at −80 • C. The genomic DNA was isolated from the thoraxes of honey bees using the G-spin Total DNA Extraction Kit (iNtRON Biotechnology, South Korea, Cat No. 17045), according to the manufacturer's instructions at the Division of Life Sciences of Incheon National University. The species of all collected honey bee samples were determined morphologically [5]. All PCR products were purified using the QIAquick PCR Purification Kit (250) (QI-AGEN, Hilden, Germany) following the instructions of the manufacturer. The PCR products were sequenced on both strands in Macrogen Inc. (Seoul, Korea) using the ABI 3730xl 96-capillary DNA analyzer (Applied Biosystems, Foster City, CA, USA) and the ABI BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA). All sequences were analyzed and submitted to GenBank (National Center for Biotechnology Information, NCBI) (Table S1). In total, 210 sequences of six exons of 35 honey bee samples (A. cerana and A. mellifera) were deposited to GenBank (Table S1). Comparative analysis of the VG gene sequences of A. cerana and A. mellifera was performed using Clustal W methods of alignment of nucleotide sequences of honey bee samples retrieved from GenBank in MEGA v. 10.0.5 at the Division of Life Sciences of Incheon National University [18]. The sequences of the VG gene exons of A. cerana and A. mellifera were retrieved from GenBank (Table S1). In order to reduce the number of samples on the dendrogram, the samples from the same regions were combined into regional groups.
Pairwise nucleotide sequence divergences were estimated using Unipro UGENE 1.28 (UNIPRO, Moscow, Russia) and CLC Genomics Workbench 11 (CLCbio, Aarhus, Denmark) based on the VG gene sequences with the Jukes-Cantor model [19]. Pairwise Euclidean distances between A. cerana populations based on morphology data were calcu-  (Table 1), 0.5 µL of template DNA, and 8.6 µL of nuclease free water under following conditions: All PCR products were purified using the QIAquick PCR Purification Kit (250) (QIA-GEN, Hilden, Germany) following the instructions of the manufacturer. The PCR products were sequenced on both strands in Macrogen Inc. (Seoul, Korea) using the ABI 3730xl 96-capillary DNA analyzer (Applied Biosystems, Foster City, CA, USA) and the ABI BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA). All sequences were analyzed and submitted to GenBank (National Center for Biotechnology Information, NCBI) (Table S1). In total, 210 sequences of six exons of 35 honey bee samples (A. cerana and A. mellifera) were deposited to GenBank (Table S1).
Comparative analysis of the VG gene sequences of A. cerana and A. mellifera was performed using Clustal W methods of alignment of nucleotide sequences of honey bee samples retrieved from GenBank in MEGA v. 10.0.5 at the Division of Life Sciences of Incheon National University [18]. The sequences of the VG gene exons of A. cerana and A. mellifera were retrieved from GenBank (Table S1). In order to reduce the number of samples on the dendrogram, the samples from the same regions were combined into regional groups.
Pairwise nucleotide sequence divergences were estimated using Unipro UGENE 1.28 (UNIPRO, Moscow, Russia) and CLC Genomics Workbench 11 (CLCbio, Aarhus, Denmark) based on the VG gene sequences with the Jukes-Cantor model [19]. Pairwise Euclidean distances between A. cerana populations based on morphology data were calculated using Statistica 8.0. Based on pairwise alignments, amino acid identity (%) was calculated for homologous genes. The significance of Tajima's D test was determined using coalescent simulations with 10,000 runs as implemented in DNAsp [20].
A phylogenetic tree was constructed using the neighbor-joining method, the Jukes-Cantor model, and 1000 bootstrap replications [21]. The percent of evolutionary divergence, number of nucleotide differences, number of transitions and transversions, and dN/dS ratio of nonsynonymous (dN) to synonymous (dS) substitutions were determined in MEGA X based on the VG gene sequences using the Nei-Gojobori method [22].
The dN/dS ratio between all A. cerana samples was evaluated. The dN/dS ratio remains one of the most popular and reliable measures of evolutionary pressures in coding regions. The dN/dS ratio provides information about the evolutionary forces operating on a particular gene [23]. The dN/dS ratio quantifies the mode and strength of selection by comparing synonymous substitution rates (dS)-assumed to be neutral-with nonsynonymous substitution rates (dN), which are exposed to selection as they change the amino acid composition of a protein. Under neutrality, dN/dS = 1. For genes that are subject to functional constraint such that nonsynonymous amino acid substitutions are deleterious and purged from the population, i.e., genes under negative selection, dN/dS < 1. For the genes under positive selection, dN/dS > 1 [23].
We performed Tajima's D neutrality test for A. cerana populations from different local populations based on the VG gene sequences. Tajima's D test is a statistical method for testing the neutral mutation hypothesis by DNA polymorphism [24]. Tajima's D is a statistical method for testing the neutral mutation hypothesis by DNA polymorphism [24]. The purpose of Tajima's D test is to distinguish between a DNA sequence evolving neutrally without selection and non-neutrally under selection. A neutrally evolving DNA sequence contains mutations with no effect on the fitness and survival of an organism. Tajima's D test computes a standardized measure of the total number of segregating sites in the sampled DNA and the average number of mutations between pairs in the sample. A negative Tajima's D signifies an excess of low-frequency polymorphisms relative to expectation, indicating population size expansion after a bottleneck or a selective sweep. A positive Tajima's D signifies low levels of both low-and high-frequency polymorphisms, indicating a decrease in population size and/or balancing selection [24].

Results
The VG gene sequences of A. cerana collected from Korea, Russia, Japan, Nepal, and China were first sequenced (with an average size of 4125 bp) and then aligned with the VG gene sequences of A. cerana and outgroup sequences of A. mellifera (with an average size of 4128bp) retrieved from GenBank (Table S1, Figure S1). The VG gene sequences of A. mellifera were used as an outgroup for A. cerana samples. Both honey bee species, A. cerana and A. mellifera, have a similar number and size of VG gene exons, and differ only in nucleotide sequences ( Figure S1).
All samples of A. cerana from Korea, Russia, Japan, Nepal, and China were pooled into twenty populations according to their geographical origin. All three samples of A. mellifera from India and Poland were pooled into one outgroup. The grouping of honey bee samples according to their spatial distribution was applied in order to reduce the number of groups so that the phylogenetic tree and tables were not cumbersome and were easy to understand. The population genetic parameters, which included the number of nucleotide differences (transitions + transversions), the ratio of nonsynonymous to synonymous substitutions (dN/dS), and Jukes-Cantor genetic distances, were assessed on twenty-one populations (twenty populations of A. cerana and one outgroup population of A. mellifera) (Tables S2-S4 The phylogenetic tree was constructed using VG gene sequences based on Jukes-Cantor genetic distances using the neighbor-joining method of clustering and 1000 bootstrap replications ( Figure 2). Korean A. cerana populations were subdivided into three clusters. Japanese A. cerana populations were distantly joined to the Korean A. cerana group. All Chinese and Nepalese A. cerana populations were clustered together. The Russian A. cerana population was located apart from all A. cerana populations (Figure 2). The Jukes-Cantor genetic distances varied from 0.002 to 0.015 between A. cerana populations and from 0.065 to 0.070 between A. cerana and A. mellifera populations. The minimal Jukes-Cantor genetic distances varied from 0.002 to 0.003 between Korean (Goseong-Hongcheon, Goseong-Samcheok, Goseong-Sancheong, Goseong-Okcheon, Goseong-Uisung, Gangneung-Samcheok, Hongcheon-Samcheok, Hongcheon-Uisung, Samcheok-Okcheon, Gangneung-Sancheong, Hongcheon-Okcheon, Gangneung-Wanju, Sancheong-Wanju, Seoul-Wanju) A. cerana populations. The maximal Jukes-Cantor genetic distances varied from 0.012 to 0.015 between Russian and Korean (Vladivostok-Uisung, Vladivostok-Goseong), Russian and Japanese (Vladivostok-Kyoto), Russian and Chinese (Vladivostok-Jiangxi, Vladivostok-Yunnan), and Russian and Nepalese (Vladivostok-Kathmandu) A. cerana populations.
The phylogenetic tree was constructed using VG gene sequences based on Jukes-Cantor genetic distances using the neighbor-joining method of clustering and 1000 bootstrap replications (Figure 2). Korean A. cerana populations were subdivided into three clusters. Japanese A. cerana populations were distantly joined to the Korean A. cerana group. All Chinese and Nepalese A. cerana populations were clustered together. The Russian A. cerana population was located apart from all A. cerana populations (Figure 2).  In addition, the forty samples of A. cerana from Korea, Russia, Japan, Nepal, and China were grouped into five major populations according to the location of the country. All three samples of A. mellifera from India and Poland were pooled into one outgroup as representatives of one sister species. The aim of the paper is to study populations of A. cerana in comparison with outgroup sister species A. mellifera. The grouping of A. cerana samples into five major populations and A. mellifera samples into one major outgroup population provides insight into the global population structure of A. cerana in Asia. The population genetic parameters, which included the number of nucleotide differences (transitions plus transversions), the ratio of nonsynonymous to synonymous substitutions (dN/dS), and Jukes-Cantor genetic distances, were assessed on six populations (five populations of A. cerana and one outgroup population of A. mellifera) ( Table 2).
The average number of nucleotide differences, including transitions and transversions, varied from 20.6 to 54.3 between A. cerana populations and from 255.1 to 269.8 between A. cerana and A. mellifera populations. The minimal number of nucleotide differences, including transitions and transversions, varied from 20.6 to 28.8 between Japanese-Korean and Japanese-Nepalese A. cerana populations. The maximal number of nucleotide differences, including transitions and transversions, varied from 51.5 to 54.3 between Russian-Chinese and Russian-Nepalese A. cerana populations ( Table 2).   (Tables 2 and S4).
We constructed the neighbor-joining phylogenetic tree and the three-dimensional PCA plot based on Jukes-Cantor genetic distances between VG gene sequences of five A. cerana and one outgroup A. mellifera populations ( Figure 3A,B). A. cerana populations from Korea and Japan were clustered together and separately from populations from Russia, Nepal, and China ( Figure 3A).
The Tajima's neutrality test was assessed for all A. cerana populations based on the VG gene sequences (Table 3). Additionally, the number of segregating sites, nucleotide diversity, number of synonymous and nonsynonymous sites, and the number of nucleotide differences were counted on the basis of the VG gene sequences (Table 3). diversity, number of synonymous and nonsynonymous sites, and the number of nucleotide differences were counted on the basis of the VG gene sequences (Table 3).

Discussion
The genetic structure of the A. cerana populations in the Korean Peninsula was insufficiently studied previously using nuclear genetic markers. The Korean Peninsula is a predominantly mountainous landscape with various altitudes, ranging from a continental climate in the north to a subtropical climate in the south. South Korea, part of the East Asian Monsoon region, has a temperate climate with four distinct seasons. The movement of air masses from the Asian continent exerts a greater influence on South Korea's weather than air movement from the Pacific Ocean. Winters are usually long, cold, and dry, whereas summers are short, hot, and humid. Spring and autumn are rather short [25].
In this project, we analyzed A. cerana from a wide range of the Korean Peninsula, from south to north and from east to west in a comparison with A. cerana from northern (Russia, Japan) and southern (China, Nepal) countries and outgroup samples of A. mellifera, to find genetic signatures of adaptations to the local environments based on the nucleotide polymorphism of the nuclear VG gene. Both A. cerana and A. mellifera species have similar exon-intron structures and differ only by nucleotide sequences (Figure 2). This makes it possible to provide comparative analysis samples from different regions with adaptation to cold and hot climates. Due to the VG gene's involvement in development and behavior, its sequences can reflect local adaptations to various climates. Samples of A. cerana adapted to similar environmental conditions should have fewer genetic differences in VG gene sequences than samples adapted to different environmental conditions.

Discussion
The genetic structure of the A. cerana populations in the Korean Peninsula was insufficiently studied previously using nuclear genetic markers. The Korean Peninsula is a predominantly mountainous landscape with various altitudes, ranging from a continental climate in the north to a subtropical climate in the south. South Korea, part of the East Asian Monsoon region, has a temperate climate with four distinct seasons. The movement of air masses from the Asian continent exerts a greater influence on South Korea's weather than air movement from the Pacific Ocean. Winters are usually long, cold, and dry, whereas summers are short, hot, and humid. Spring and autumn are rather short [25].
In this project, we analyzed A. cerana from a wide range of the Korean Peninsula, from south to north and from east to west in a comparison with A. cerana from northern (Russia, Japan) and southern (China, Nepal) countries and outgroup samples of A. mellifera, to find genetic signatures of adaptations to the local environments based on the nucleotide polymorphism of the nuclear VG gene. Both A. cerana and A. mellifera species have similar exon-intron structures and differ only by nucleotide sequences (Figure 2). This makes it possible to provide comparative analysis samples from different regions with adaptation to cold and hot climates. Due to the VG gene's involvement in development and behavior, its sequences can reflect local adaptations to various climates. Samples of A. cerana adapted to similar environmental conditions should have fewer genetic differences in VG gene sequences than samples adapted to different environmental conditions. According to the VG gene analysis, populations of A. cerana from northern Asian countries, such as China, Japan, Korea, and Russia, grouped together and separately from populations of A. cerana from southern Asian countries, such as Nepal.
The genetic analysis was provided with 40 samples of A. cerana from Korea, Russia, Japan, Nepal, and China from 20 geographical populations and outgroup A. mellifera samples from India and Poland. The average number of nucleotide differences between the VG gene sequences, including transitions and transversions, varied from 20.6 to 28.8 between A. cerana populations from different countries, which is an indication of the genetic divergence level. The minimal number of nucleotide differences (20.6) was observed between Korean and Japanese populations of A. cerana (Table S2), which proves their close genetic relationship due to recent active movement of bee colonies between countries. The nucleotide differences between the A. cerana and A. mellifera outgroup populations were highest (255.1-269.8), which demonstrates their significant divergence. We counted the number of nucleotide differences between the VG gene sequences within the northern populations of A. cerana, which varied from 10 to 52 (the average was 21.9). Additionally, we counted the number of nucleotide differences in the VG gene sequences between the northern and southern populations of A. cerana, which varied from 22 to 61 (the average was 32.4) (Table S2). Here, the majority of the nucleotide differences of the VG gene sequences between the northern and southern populations of A. cerana were higher than those within the northern populations of A. cerana. The number of nucleotide substitutions likely reflects the magnitude of local adaptations, which are minimal between populations from similar climatic conditions. Additionally, the similarity in VG gene sequences can be caused by human-assisted movement of A. cerana colonies between countries.
The dN/dS ratio between all A. cerana samples varied from 0.158 to 0.999, which demonstrates the functional significance of the VG gene such that nonsynonymous amino acid substitutions are deleterious and purged from the population. Indeed, the VG gene affects multiple aspects of social organization of honey bee colonies and plays a key role in the caste differentiation and maintenance of honey bees. The dN/dS ratio between the northern populations of A. cerana varied from 0.158 to 0.999 (the average was 0.46), which indicates strong negative selection for the VG gene, which characterizes the high significance of the VG gene's conservativeness in the adaptation of A. cerana to the cold northern climate. Positive selection leads to fixation changes in genes, and negative, or purifying, selection leads to rejection changes in genes. The dN/dS ratio between the northern and southern populations of A. cerana varied from 0.289 to 0.947 (the average was 0.49) (Table S3), which indicates less strong negative selection of the VG gene in evolution and the high significance of the VG gene polymorphism in the adaptation of A. cerana to the different climates of southern and northern Asia. In comparison, the average dN/dS ratio for 3256 nuclear genes was 0.1068 between 3 species of Formicidae, 0.1033 between 3 species of Polistinae and Vespinae wasps, and 0.1086 between 10 species of Apis, Bombus, and Tetragonula. It was 0.1025 between three species of Siricoidea and 0.1005 between five species of Cynipoidea, which characterized the negative selection of these genes in evolution and their role in adaptation to the local climate [26]. It has been reported that genes involved in the labor division, caste differentiation, and development are affected by strong negative selection in social insects. Authors estimated the strength of negative selection at 5 UTRs, 3 UTRs, introns, exons, intergenic sites, and overall genomes in A. mellifera, B. impatiens, P. dominula, and S. invicta and found that 0-fold sites with nonsynonymous changes and overall genomes experienced the strongest negative selection in all four eusocial species [27]. On the contrary, another study showed evidence for positive selection on the VG gene in A. mellifera and A. cerana, which leads to fixation of the nonsynonymous mutations, bias of the allele frequency spectrum, growth of the genetic differentiation at neutral sites, and magnification of the linkage disequilibrium relative to other genes [13]. Hence, the VG gene is very important for the adaptation of A. cerana; most nucleotide variations reduce its sustainability and are eliminated from populations.
The Jukes-Cantor genetic distance is the genetic relatedness between two nucleotide sequences calculated as the sum of the substitutions that have accumulated between them since they diverged from their common ancestor during evolution, assessed using the Jukes-Cantor model. The Jukes-Cantor genetic distances based on the VG gene sequences between A. cerana populations varied from 0.002 to 0.015, and between A. cerana and A. mellifera populations varied from 0.065 to 0.070. The Jukes-Cantor genetic distances based on the VG gene sequences between the northern populations of A. cerana varied from 0.002 to 0.013 (the average was 0.005). The Jukes-Cantor genetic distances based on the VG gene sequences between the northern and southern populations of A. cerana varied from 0.005 to 0.015 (the average was 0.010) (Table S4). In comparison, the Jukes-Cantor genetic distance was 0.018 between subspecies of A. mellifera (A. m. carnica, A. m. ligustica, A. m. sicula, A. m. iberica, A. m. adami, A. m. macedonica, A. m. anatoliaca, A. m. syriaca, A. m. intermissa) based on the COX1 gene of mtDNA [28]; was 0.17 between 2 honey bee species (A. mellifera, A. cerana) based on complete mtDNA [29]; varied from 0.001 to 0.014 between subspecies of A. cerana; varied from 0.075 to 0.081 between species A. mellifera and A. cerana based on the VG gene of nDNA [8]; was 0.033 between 15 species of Anastrepha based on the COX1 gene of mtDNA [30]; and varied from 0.002 to 0.039 within species and from 0.074 to 0.131 between 5 species of Apis based on the CYTB gene of mtDNA [31]. Thus, the Jukes-Cantor genetic distances based on the VG gene sequences between A. cerana populations from similar climatic conditions is half that between A. cerana populations from different climatic conditions. This can explain the important role of the VG gene in the adaptation of A. cerana populations to different climatic conditions from northern to southern Asia. The phylogenetic tree constructed using the neighbor-joining method with the Jukes-Cantor genetic distances based on the VG gene sequences (Figure 3) shows the informativeness of the gene in phylogenetic and ecological studies of A. cerana. Additionally, we can see the genetic relationship between A. cerana populations on the phylogenetic tree and three-dimensional PCA plot ( Figure 3). Moreover, there was a genetic differentiation based on the VG gene sequences between Korean populations of A. c. koreana that is related to the geographic localization and climatic conditions of southern, central, and northern provinces of the Republic of Korea.
The results of Tajima's D test were counted for each population of A. cerana. Almost all of the estimated D values in A. cerana bees were not significantly different from zero, implying that A. cerana honey bees are in a state of mutation drift equilibrium [17]. The number of segregating sites (S) was lowest in the Nepalese population of A. cerana and highest in the Chinese population of A. cerana. This may be due to the high genetic diversity that should be observed in the largest Chinese population of A. cerana. The nucleotide diversity was highest in the Russian population of A. cerana and lowest in the Korean population of A. cerana. This can be explained by the fact that the A. cerana population in Russia is feral under natural selection. There are no managed A. cerana populations in Russia ( Table 3).
The values of Tajima's D test in all A. cerana populations were statistically significant. The Korean population of A. cerana had the biggest negative value (−0.685), which means that this population could be in a state of bottleneck, selective sweep, purifying selection, and negative selection. Indeed, the Korean population of A. c. cerana experienced an extreme population decline with the kSBV virus, which killed 95% of the local population [4,9,10,32]. All other populations of A. cerana also had the signatures of the population size expansion after a selective sweep in the past (Table 3). For comparison, Tajima's D test values based on mitochondrial COX1 gene sequences varied from −1.369 to 2.383 in the Guangdong Province population of A. c. cerana, from −2.293 to −0.048 in the Guangxi Autonomous Region population of A. c. cerana, from −1.670 to 1.276 in the Hainan Island population of A. c. cerana, and from −1.622 to 2.383 in the mainland China population of A. c. cerana [33]. The values of Tajima's D test in Russian, Korean, Japanese, Chinese, and Nepalese populations of A. cerana based on VG gene sequences were similar to those in China based on the mitochondrial COX1 gene. Thus, both nuclear and mitochondrial genes showed similar results for A. cerana populations, which were characterized by an excess of low-frequency polymorphisms relative to expectations, indicating population size expansion after a bottleneck due to the Sacbrood virus. The Tajima's D test values and dN/dS ratio can be used as indicators for negative selection and the role of some genes in adaptation of local A. cerana populations [33].

Conclusions
We performed a thorough genetic analysis on A. cerana populations from South Korea in comparison with A. cerana populations from Russia, Japan, Nepal, and China and the outgroup population of A. mellifera. Based on the VG gene sequences, we assessed the negative selection that characterizes adaptation in A. cerana populations to the local climatic conditions. The cluster and principal component analysis based on the sequences of the VG gene divided populations of A. cerana according to their climatic and geographic distribution into southern and northern groups. The local populations of A. c. koreana were subdivided according to their geographical distribution into southern, northern, and central Korean clusters. The Jukes-Cantor genetic distances based on the VG gene sequences showed the A. cerana populations from northern countries such as Korea, Russia, and Japan are more closely related to each other than the southern populations of A. cerana. The signatures of adaptation were found in the local population of A. cerana. All A. cerana populations showed signs of population size expansion following recent declines. The VG gene sequences can be used as informative markers for monitoring the genetic structure and adaptation processes in A. cerana populations.
in molecular genetics and experimental biotechnology work, for collaboration in honey bee research and for sharing their ideas on honey bee adaptation.