Skip to main content
Download PDF

Geometric morphology and population genomics provide insights into the adaptive evolution of Apis cerana in Changbai Mountain

BMC Genomics volume 23, Article number: 64 (2022) Cite this article

Abstract

Background

Exploration of adaptive evolution of organisms in response to environmental change can help to understand the evolutionary history of species and the underlying mechanisms of adaptation to local environments, thus guiding future conservation programmes. Before the introduction of Apis mellifera in China, eastern honey bees (Apis cerana) were the only species used for beekeeping in this region. In the mountains of Changbai, populations of A. cerana are considered a distinct ecotype of the species which formed through the distinct selective pressures in this area over time.

Result

We performed a measure of 300 wing specimens of eastern honey bees and obtained the geometric morphological variation in the wing of A. cerana in Changbai Mountain. A total of 3,859,573 high-quality SNP loci were yielded via the whole-genome resequencing of 130 individuals in 5 geographic regions.

Conclusion

Corresponding geometric morphology and population genomics confirmed the particularity of the A. cerana in Changbai Mountain. Genetic differentiation at the subspecies level exists between populations in Changbai Mountain and remaining geographic regions, and a significant reduction in the effective population size and an excessive degree of inbreeding may be responsible for a substantial loss of population genetic diversity. Candidate genes potentially associated with cold environmental adaptations in populations under natural selection were identified, which may represent local adaptations in populations. Our study provided insights into the evolutionary history and adaptation of A. cerana in Changbai Mountain, as well as its conservation.

Peer Review reports

Background

Nowadays, the number of pollinators has been declining considerably; in particular, the invertebrate pollinators represented by honeybees and butterflies are on the verge of extinction [1, 2]. Apis cerana is native to Asia with a long history of managed beekeeping, the populations have also experienced severe declines in recent decades [3, 4].In Northeast China, Korean Peninsula, and even the Russian Far East in a broader sense, populations of A. cerana are considered a distinct ecotype of the species [5]. This ecotype wields a crucial role in maintaining the biodiversity of Changbai Mountains. It has adapted to the long winters and short summers of Changbai Mountain. As a crucial part of the ecological food chain in Changbai Mountain, a reduction in the number of existing A. cerana bees is bound to affect the survival of upstream predators in the future. In terms of economic value, A. cerana in Changbai Mountain bears excellent traits such as rapid reproduction, strong maintenance of colonies, collection capability, resistance to adversity, and cold tolerance [6]. Thus, it is a popular breeding species among beekeepers, and colony breeding and wild hive nectar harvesting remain as an important source of income for people living in economically underdeveloped mountain areas. Unfortunately, the biotic and abiotic factors, such as ecological invasion of A. mellifera, habitat loss and fragmentation, overexploitation of biological resources, intensive agricultural development, pathogenic microbial infection, overuse of pesticides, global warming, and atmospheric pollution, have caused a massive decline in the population size and distribution density of A. cerana in Changbai Mountain [7, 8]; its population size has decreased from more than 40,000 colonies in 1983 to 19,000 colonies in 2008 [9], and then to less than 10,000 colonies at present, with an alarming decline rate of 75%, resulting in an endangered–sustained state. Understanding the adaptive evolution of organisms towards environmental changes not only helps clarify the evolutionary history of species and the potential mechanism of adapting to the local environment but also guides for future protection plans[10,11,12]. With the dramatic changes occurring in the genetic structure and diversity of pollinating insects worldwide, considering the changes in climate, ecology, and other environmental factors, introduced or translocated exotic colonies may alter the germplasm gene pool of local A. cerana in Changbai Mountain. Collectively, the diversity of A. cerana in Changbai Mountain should be further clarified.

Shape and veins of insect wings contain a wealth of ecological and behavioural information [13,14,15]. Geometric morphology analyses of bee wings have been used for the identification of species and subspecies [16,17,18]. The wing veins of hymenopteran insects have tended to decrease during the evolutionary history [19]. It has been reported that the capped broods of Apis cerana cerana were incubated at constant lower temperature until emergence, resulting in the increase of wing veins of adults, which may be the wing veins that disappeared in the process of evolution [20]. Thus, changes in wing size and morphology under environmental stress can advance our understanding on the adaptive evolution of A. cerana in Changbai Mountain. In addition, there are differences in adaptation to various climates between eastern honey bees and western honey bees. It is generally believed that the adaptation of western honey bees to temperate climate conditions is a key step in their evolution [21]. Native bees show higher vitality and environmental tolerance than non-native bees [22,23,24], therefore A. cerana populations inhabiting the mountains of Changbai may have a more effective social defense system than A. mellifera through the combination of complex molecular and behavioral mechanisms [25]. A. cerana in Changbai Mountain has developed distinct traits in response to mountain conditions, particularly strong resistance to the cold (able to naturally overwinter in extreme cold conditions at − 40 °C) [9]. While population genomics, morphometrics and adaptation of western honey bees to cold environments have been extensively studied [18, 26], a similar assessment on A. cerana in Changbai Mountain is lacking.

In this study, we collected samples from 64 sampling sites in the Changbai Mountain, northern, southern, Qinghai–Tibetan, and northwestern regions. Moreover, resequencing data from 66 A. cerana samples of each ecotype were downloaded from the literature [27]. A comprehensive analysis of geometric morphometry and population genomics revealed the variation patterns of wing morphology, population genetic diversity, and population history of A. cerana in Changbai Mountain, and helped in identifying candidate genes related to their adaptation to the local cold environment. This paper offers an in-depth investigation of the genetic basis and uniqueness of the environmental adaptation of A. cerana in Changbai Mountain, thus rendering directions for the scientific and effective conservation of the population.

Results

Population genetic structure analysis

The samples were collected from North, Northwest, Qinghai-Tibet, South and Changbai Mountain, among which samples from Changbai Mountain were procured covering the entire natural range of the populations (Fig. 1, Additional file 1: Table S1). Whole-genome resequencing yielded a total of 243.37 Gb of data. Genome alignment resulted in an average depth of 13.32 × (Additional file 2: Table S2). For the population genetics analysis, we also downloaded 165.44 Gb of resequenced data of 66 A. cerana samples from the published literature (Additional file 3: Table S3). After quality control, a total of 3.86 million SNPs were identified, of which 253,656 were in exonic regions and 52,161 were in non-homologous regions (Additional file 4: Table S4), and a total of 2,237,786 Indels were detected, of which 15,588 were in exonic regions (Additional file 4: Table S5). Population genetic structure was first evaluated using ADMIXTURE [28] (Fig. 2A). The result showed that when K = 2–5, most Jilin individuals and all Heilongjiang individuals in Changbai Mountain were from the same ancestor, and all Liaoning individuals were of mixed lineage and had different ancestral origins. In particular, when K = 4 and K = 5, the origin of the Liaoning individuals was more complex, with individual samples containing a mixture of 3 or more ancestral subgroups, and few Jilin samples showed different degrees of admixture. According to the cross-validation errors, when K = 4, the CV error value (0.22113) was the lowest, which was the best K value. All regions were divided into 4 genetic components, the individuals from Tibetan (XZ01-XZ08) and Hainan (HI01-HI10) formed an ancestral cluster (Fig. 2A), whereas individuals from Southern region showed a more highly mixed state than the individuals from Northern and Northwestern due to their richness in geographic ecotypes. Afterwards, PCA was performed using GCTA [29] (Fig. 2B), and the results of all three PCAs separated most of the individuals of Changbai Mountain and clustered them together separately (Additional file 5: Fig. S1). The Liaoning individuals (LN01-LN09) were mixed with those from the northern, northwestern, Qinghai–Tibetan, and southern regions, further supporting the results of genetic structure analysis.

Fig. 1
figure 1

The positions of the 130 individuals in 5 geographic regions in China. The sampling location of the Changbai Mountain area is highlighted as a yellow enlarged area map. Triangles represent the sampling locations of the samples downloaded; circles represent the sampling locations of the samples from northern, southern, northwestern, and Qinghai-Tibetan regions; and stars represent the sampling locations of the samples from Changbai Mountain

Full size image
Fig. 2
figure 2

Results of population genetic structure analysis. (A) Genetic structure of A. cerana populations from 5 regions. Groupings of samples from two to five ancestral clusters are shown. (B) Scatter plot of PC1 versus PC2 for the A. cerana populations

Full size image

Genetic differentiation and genetic diversity

The nucleotide diversity (π) and Fst values of A. cerana for five geographic regions and Changbai Mountain were analysed. Among A. cerana individuals in each geographical region, the highest π was found in the individuals from the North region, and the lowest was in the Changbai Mountain (Fig. 3A). The highest π was found in the Liaoning and the lowest in the Heilongjiang among individuals from Changbai Mountain (Fig. 3B). Paired Fst ranged from 0.0094 to 0.2609, with an average of 0.1133, and A. cerana in various regions exhibited moderate genetic differentiation [27, 30, 31]. The mean Fst was 0.2294 between the individuals from Changbai Mountain and the A. cerana from the northern, northwestern, southern, and Qinghai–Tibetan regions, showing a high degree of genetic differentiation (Table 1). Paired Fst for samples from each region within the Changbai Mountain individuals ranged from − 0.0031 to 0.2030 (Table 2). Fst for Jilin individuals from JL-BS, JL-JL, JL-TH, and JL-YB were all below 0.02, indicating an extremely low degree of genetic differentiation. The average Fst between the individuals from Heilongjiang and Jilin was 0.1002, which indicated moderate to low differentiation. Notably, the individuals from Liaoning reached a high degree of genetic differentiation, with an average Fst of 0.1506.

Fig. 3
figure 3

The result of nucleotide diversity(π). (A) A. cerana from five regions. (B) A. cerana in Changbai Mountain

Full size image
Table 1 Paired Fst genetic distances between individuals from 5 regions of A. cerana
Full size table
Table 2 Paired Fst genetic distances between individuals from Changbai Mountain region
Full size table

We evaluated expected heterozygosity (He), observed heterozygosity (Ho), inbreeding coefficient (F), and linkage disequilibrium (LD) for each geographic regions of A. cerana, and A. cerana in Changbai Mountain had the lowest He, highest Ho, and highest F (Table 3), indicating that they had a low level of genetic diversity, a high level of genetic differentiation from other geographic regions of A. cerana, and a high degree of inbreeding. LD analyses showed that the Changbai Mountain region exhibited the slowest LD decay (Fig. 4). Presumably, inbreeding decline may be responsible for the lowest genetic diversity, which is consistent with the results of the analysis of genetic diversity indicators.

Table 3 Genetic diversity information of individuals from 5 regions of A. cerana
Full size table
Fig. 4
figure 4

The analysis of LD for A. cerana

Full size image

Demographic history

The PSMC model is a reliable method to estimate the changes in effective population size [32]. The quantitative variation in the effective population size of A. cerana in various regions, including Changbai Mountain populations, coincided with the cyclical variation in historical global temperature, which peaked during the most recent major interglacial period of the Pleistocene, approximately 13,000–11,000 years ago, when the global climate was warm, followed by a sharp decline, consistent with the findings of Chen et al. [27] who reported that the effective population size of A. cerana peaked at MIS5. The difference lied in the fact that the A. cerana population in Changbai Mountain exhibited a more moderate expansion and contraction of the effective population and the population size was at a lower level than A. cerana populations in other regions (Fig. 5A).

Fig. 5
figure 5

Effective population size of A. cerana. Red line represents the effective population size of an A. cerana in Changbai Mountain. (A) PSMC; (B) SMC +  + 

Full size image

The SMC +  + method is an extension of the PSMC method for estimating more precisely the effective population size within 10,000 years [33]. From the MIS4 to the Holocene period, the A. cerana population in Changbai Mountain experienced two significant population contractions and one significant population expansion compared with populations in other regions (Fig. 5B).

Selective sweep analysis

We performed a selective sweep analysis between Changbai Mountain individuals and Hainan individuals (Additional file 2: Table S2&Additional file 3: Table S3). Genomic regions with the top 5% values for the Fst and θπ ratios (log2(θπ ratio)) in each pairwise comparison were considered to be potentially selective-sweep regions (Fig. 6A). Genes were categorised according to the reference genome for functional gene annotation and KEGG and GO functional categories. A total of 273 genes were selected in the Changbai Mountain population (Additional file 6: Table S6), and they were involved in 1621 GO and 40 KEGG pathways, with 99 GO and 5 KEGG pathways significantly enriched at a threshold of P ≤ 0.05 (Additional file 7: Table S7; Table S8).

Fig. 6
figure 6

(A) Plot of the selective sweep analysis between the individuals from Changbai Mountain and the individuals from Hainan based on Fst and θπ. The horizontal coordinates are the θπ ratio values and the vertical coordinates are the Fst values, which correspond to the frequency distribution plot above and the frequency distribution plot on the right, respectively, and the dot plots in the middle represent the corresponding Fst and θπ ratio values in different windows, where the blue and green areas are the top 5% areas selected by θπ and the red areas are the top 5% areas selected by Fst; (B) A selection elimination analysis plot based on Fst and θπ for individuals from Changbai Mountain and Tibet

Full size image

In this work, we paid considerable attention to the genes related to cold adaptation in the Changbai Mountain region. According to the decreasing order of Fst values, we focused on the genes involved in the top 20 bins. Among these, four structural constituent of cuticle protein and four lipid storage droplets surface-binding protein were included (Additional file 8: Table S9). Based on previous studies on cold tolerance of honeybees, we found many of these genes involved in the GO terms and KEGG pathways to be associated with cold adaptation in insects [34, 35]. In the GO classification, the selected signals were mainly enriched in the metabolism of energy substances, synthesis and production processes, transmembrane transporter activity, ion transport, channel activity, homeostasis, localisation, regulation, binding, enzyme activity, response, signal transduction, and sensory perception. Ether lipid metabolism, glycerophospholipid metabolism, and sphingolipid metabolism are the KEGG pathways in which higher Fst values are jointly involved by selective genes.

The pathway enrichment of the selected genes combined with the functional annotation of the genes themselves resulted in the identification of several highly enriched functional categories, including structural constituent of cuticle (4), zinc ion binding (4), serine/threonine-protein kinase (5), lipid storage droplets surface-binding protein (4), and G-protein-coupled receptor 1 family (4) (Additional file 6: Table S6). Through the selective sweep analysis of population evolution, we attempted to determine the candidate genes that may be related to local adaptation. Because the Tibetan and Changbai Mountain populations belong to the cold-tolerant honeybee ecotype, the same method was used to perform the selective sweep analysis between the Tibetan and Hainan individuals (Fig. 6B). In total, 163 selected genes were obtained for the Tibetan individuals, and 22 genes were identical to the selected region of the Changbai Mountain individuals (Additional file 9: Table S10), of which two enriched genes belonged to the G-protein-coupled receptor 1 family.We speculated that this gene family might wield a role in the genetic basis of cold adaptation in the Changbai Mountain population.

Morphometric analyses of wings across populations

For the morphological analysis, 30 out of 64 individuals were selected. They were obtained from 6 northern regions, 6 northwestern regions, 3 Qinghai–Tibetan regions, 8 southern regions, and 7 Changbai Mountain regions. Among these, the samples from Changbai Mountain were JL01 (belongs to JL-YB), JL06 (belongs to JL-JL), JL07 (belongs to JL-TH), and JL13 (belongs to JL-BS) in Jilin Province; LN01 and LN09 in Liaoning Province; and HL03 in Heilongjiang Province (Fig. 1; Additional file 1: Table S1).

The wing size is expressed in terms of the centriod size (CS) of the wing. One-way ANOVA revealed a significant difference in CS values of forewings of A. cerana distributed in the five geographic regions (F = 54.38, P < 0.001). The results of multiple tests by the LSD method showed that the individuals in the northwestern region, Qinghai–Tibetan region, and Changbai Mountain had larger forewings, and no notable differences were observed between the individuals. Individuals from Changbai Mountain differed significantly from the individuals in the northern and southern regions in terms of forewing size (P < 0.05) (Fig. 7A). Hindwing data demonstrated a highly significant difference in the hindwing CS values among individuals from each region (F = 36.08, P < 0.001). A. cerana from the Qinghai–Tibetan region, northwestern regions, and Changbai Mountain had larger hindwings, and non-significant differences were observed between these individuals, whereas the individuals from Changbai Mountain differed notably (P < 0.05) from the northern and southern region individuals (Fig. 7B). Stepwise regression analyses indicated that the mean CS values of forewings and hindwings at each sampling site increased with the latitude and altitude (P < 0.05), whereas no significant correlation with longitude was observed (P > 0.05) (Additional file 10: Table S11).

Fig. 7
figure 7

Boxplot of CS of A. cerana in wing size across geographical populations. (A) Centroid size of forewing of A. cerana. (B) Centroid size of hindwing of A. cerana. Note:Data with different letters indicated significant variations among individuals from each regions at 0.05 level

Full size image

The first two principal components of forewings and hindwings occupied 77.74% and 83.29% of the total variance, respectively, both of which were able to illustrate main differences in forewings and hindwings among the geographical regions. Results of forewing PCA showed that individuals from Jilin province and Heilongjiang province were concentrated in the positive half-axis of PC1 and PC2, whereas an individual from Liaoning province (LN09) was in the negative half-axis of PC1 and PC2 (Additional file 11: Fig. S2). Results of hindwing PCA showed that individuals from Jilin province and Heilongjiang province were concentrated in the positive half-axis of PC1, whereas an individual from Liaoning province (LN01) was in the negative half-axis of PC1 (Additional file 11: Fig. S2). Taken together, a certain degree of wing morphological differentiation was observed within A. cerana in Changbai Mountain, with no significant differentiation observed between the individuals from Jilin and Heilongjiang, and significant differentiation observed between Liaoning individuals.

Discussion

A certain number of wild bee colonies are present in Changbai Mountain, mostly in the hollow caves of trees, rock caves in steep mountain cliffs, rock crevices, underground caves, and other hidden places. Eastern honey bees have the germplasm characteristic of flight escape, and a certain number of wild colonies are obtained from various reared colonies and are constantly replenished. In recent decades, human capture and the outbreak of sacbrood virus have greatly changed the population dynamics of eastern honey bees in Changbai Mountain [25, 36]. The wild population suffered a devastating loss, resulting in an increase in the proportion of the wild colonies obtained from the managed colonies. A large reduction in the number of wild bees leads to a decrease in the effective population, which may directly lead to a decrease in the genetic diversity of the Apis cerana in Changbai Mountain, and a low level of genetic diversity tends to increase inbreeding. The previous analysis of genetic diversity of A. cerana in Changbai Mountain, which were based on SSR and mtDNA markers, also reflected the current situation of a single genetic structure and low genetic diversity [37,38,39]. It has the highest degree of LD, which is likely to have experienced a serious evolutionary bottleneck event, while showing a very low genetic polymorphism and a slow level of decay. Therefore, we believe that the declining effective population size and loss of genetic diversity explained the main intrinsic factors for the endangerment of A. cerana in Changbai Mountain.

The results of population history analysis show that biogeography and climate change have seriously affected the Changbai Mountain population, and glaciation and local volcanic events have likely brought about the bottleneck effect of population history. Available evidence demonstrated that the A. cerana was affected by cyclical climate changes during the Quaternary ice age [40], and that volcanic eruption activity at Tianchi, Changbai Mountain, continued throughout the Quaternary [41, 42], inevitably impacting the population evolution and distribution patterns of the existing A. cerana in Changbai Mountain. According to the Quaternary glacial remains around the Tianchi volcanic cone, the basic geomorphic features of the Tianchi caldera were formed in the MIS4 and were maintained to the present day after basic stereotypes [42], and the fossil glacier scale of the last glacial period was larger in the early (MIS4) than in the late (MIS2) [43], which implied that the paleoclimate of Changbai Mountain was greatly influenced by glaciation at that time. We therefore focused on labelling the changes in the number of effective population size in the MIS4. MIS4 belongs to the early stage of the last glacial period. The glacier development scale is large. The contraction of the effective population size in the cold period indicated that the decline of global temperature may be detrimental to the development of A. cerana.

The population decline during the MIS4, which lasted until the last glacial maximum of approximately 15,000 years ago, is presumed to be mainly due to the paleoglaciation, grand development of the ice sheet, low snowline, and small dust flow, indicative of a cold and wet climate [44], and the harsh natural ecological climate brought about a population decline, which corroborated with the analysis results obtained by the PSMC method. A massive explosive eruptive activity occurred at Tianchi volcano during the Tianwenfeng period at 50,600 years ago, which formed the first breakout crater [45], and climate change and ecological degradation caused by strong volcanism might have aggravated the decrease in A. cerana. The effective size of the A. cerana in Changbai Mountain increased significantly for once from the last glacial maximum of 15,000 years ago to the Holocene, indicating that global warming was beneficial to the Changbai Mountain population. The period of maximum extent of continental icesheet was 18,000 years ago, and since then, glaciers began to retreat. The Holocene warm period in the northeast was the best period for vegetation development, and the enhanced summer monsoon precipitation due to warming led to ecological recovery [46]. Global warming accelerated the rapid expansion of honeybee populations. The effective population size declined significantly again at approximately 5000–2000 years ago. During the middle Holocene period, approximately 8000–2500 years ago, the climate was warm and humid, the temperature at its peak may have been more than 2 °C higher than modern times [47], the ecological conditions were favourable for species reproduction, while the effective population of honeybees drastically reduced. According to volcanic geology and historical records, there have been several eruptions of the Changbai Mountain 5000 years ago [48, 49], which produced ash drifts up to thousands of kilometers away and deposition thicknesses of several centimeters [50, 51], causing large-scale biological extinctions and serious impacts on the global climate, environment, and ecology [52,53,54]. In addition, a small peak in the number of hardy plants such as pine, spruce, fir, and birch occurred in the Changbai Mountain around 5000 years ago, which may indicate a climatic cooling event [55], and a similar European Ulmus decline occurred in the pollen fossil belt 5000 years ago in eastern China [56]. Therefore, it is reasonable to speculate that the combination of volcanic heat and abnormal cooling led to a significant decline in the effective population size of A. cerana in Changbai Mountain during 5000 years ago.

We explored the selected genes in the evolution of the Changbai Mountain population from multiple perspectives by conducting the selective sweep analysis. The KEGG pathways involved in the selected genes with a high FST value are all related to the metabolism of lipid energy substances, and energy metabolism is crucial for the cold resistance of honeybees [26, 35]. The overwintering period of the Changbai Mountain population lasts for six months, which is a longer period of cold domestication than the colonies living in the warm environment experience in the south, and its cold tolerance is exercised for a long time to develop superior resistance to cold. Cold tolerance exercise in insects is normally an energy-consuming process, and cold tolerance adaptive traits are genetically and evolutionarily stable [57], indicating that the population has evolved a complete adaptive capacity for storage and metabolism of cold-resistant substances under long-term low temperature stress. The result of wing geometric morphometric showed that the CS values of A. cerana populations in each geographical region increased with latitude and altitude. Body size in insects may be related to cold resistance, highland populations of A. cerana are larger, darker and have longer body hair [12]. The superiority of the Changbai Mountain population for wing and body sizes plays a role in the survival competition of A. cerana. Larger wings imply larger insects who adapt to cold conditions at high altitudes and long overwintering periods through an increase in body size [26, 61]. The cold environment of Changbai Mountain challenges the thermoregulation, flight, and reproductive abilities of honeybees, which have larger CS for its forewings and hindwings, implying that individuals can carry more weight in flight [62], have a greater foraging range [63], and have an increased nectar collection ability [64, 65].

The subspecies-level genetic differences among the individuals from Changbai Mountain and Northern, Northwestern, Southern, and Qinghai-Tibetan regions were observed [30, 31, 66], reflecting the special status of populations in the evolutionary process. Some researchers have suggested that A. cerana in Changbai Mountain were morphologically similar to the Korean and Japanese populations, and may be classified as the same taxon [67, 68]. Worker bees usually have a bifurcated protrusion on the radiomedial crossvein and a high cubital index of the forewing, which is a major feature that distinguishes the Changbai Mountain population from other ecotype populations (Fig. 8), and this feature is often used as a bee morphological identification index for genetic classification studies [9]. Zhou et al. [69] confirmed significant morphological differences and similar degrees of differentiation in samples collected from Jilin, Hainan, and Shaanxi, with the cubital index of Jilin A. cerana being significantly different from that of the other two honeybees, and suggested that the Jilin population may be a novel subspecies. Gene flow blocking is responsible for population differentiation, prompting differentiation into independent populations that produce genetic traits adapted to the local environment [37]. A recent study revealed the population structure and rapid radiation history of eastern honey bees in mainland Asia, where northeastern populations radiated from a central ancestral group with 6 peripheral subspecies, including Taiwan, Hainan, Bomi, Aba, and Chinghai, and adapted independently to different habitats [70]. On the other hand, A. cerana in Changbai Mountain bears varying degrees of genetic differentiation within its population. The genetic differences were negligible among the local individuals of Jilin province, which exhibited a low genetic differentiation compared with Heilongjiang individuals, whereas a significant and major genetic differentiation was noted between the Jilin and Liaoning individuals. It is also proved by the results of wing geometric morphology.

Fig. 8
figure 8

Location of the geometric landmarks on the wing of A. cerana in Changbai Mountain. (A) 20 landmarks on the right forewing; the red circle is the small protrusion vein on the radiomedial crossvein. (B) 10 landmarks on the right hindwing

Full size image

In summary, our research shows that eastern honey bees from Changbai Mountain have genetic and morphological singularity. However, populations in the core area of Changbai Mountain are under threat from the outer reserves, it may lead to irreversible loss of local adaptations in native populations, genetic homogenisation, and pose a risk of genetic contamination [71, 72]. Therefore, local subspecies must be effectively protected. First, it is imperative to strengthen the management of core germplasm reserves and restrain introductions of colonies from non-germplasm core areas to reduce the risk of genetic infiltration of native bees and introduction of pathogens and/or parasites. Second, the feasibility of controlling queen mating should be considered, special areas in the core area where local bees can gather and where mating can be controlled should be designated. Third, conducting basic research on population reproduction, foraging, mating behaviour, and production traits is necessary for successful population conservation, and using appropriate molecular biology techniques to assess the native honeybees population size and regular genetic monitoring of populations will help maintain their purity. Finally, wild colonies experience different selection pressures relative to managed colonies, which can provide insights into the process of adaptation to local environmental conditions, and these honeybees can serve as a genetic reserve gene pool for future managed populations.

Conclusions

The body and wing sizes of A. cerana in Changbai Mountain showed a gradual change with altitude and latitude. A certain degree of wing morphological differentiation was observed within A. cerana in Changbai Mountain. Genetic differences were observed between Changbai Mountain and other geographic regions at subspecies level, which has a lower effective population size and genetic diversity, and a higher degree of inbreeding. Selected genomic regions associated with adaptation of populations to cold environments were identified, and these genes may represent local adaptations of population, which serve as excellent candidates for further research. Geometric morphological and colony genomic evidence confirmed the special evolutionary status of A. cerana in Changbai Mountain, which will contribute to future conservation programmes.

Methods

Honeybee collections

A. cerana samples were collected from 64 sampling sites in 5 geographical regions of China, including northern, northwestern, Qinghai–Tibetan, southern regions, and Changbai Mountain. All the samples were taken from germplasm reserves or original rearing colonies, with priority given to samples from wild hives (e.g. tree cave hives) and traditional-type hives (e.g. tree barrel hives). The samples were transported in anhydrous ethanol and stored in an ultra-low temperature refrigerator at − 80 °C.

Sequencing, quality control, and alignment reference genome

One workerbee was randomly selected from each colony in the sampling site for whole-genome resequencing. First, genomic DNA was extracted and quantified precisely, and DNA samples with a total volume of ≥ 1.5 μg were used for library construction. The inspected DNA samples were randomly fragmented into 350-bp-long fragments through Covaris, followed by high-quality library construction using the TruSeq Library Construction Kit, the library was prepared through various steps such as end repair, polyA tailing, sequencing junction, purification, and PCR amplification. The constructed library was sequenced by Illumina. Qubit3.0 was then used for preliminary quantification, and the insert size of the library was measured by Agilent 2100 after diluting the library to 1 ng/μL. Then, the effective concentration of the library was accurately quantified by Q-PCR. Finally, according to the effective concentration of library and data output requirements, the Illumina platform was used for PE150 sequencing.

Effective high-quality sequencing data were mapped to the reference genome [73] by using BWA [74] software (set parameter mem -t 4 -k 32 -M), and then the potential PCR duplicates were removed by employing SAMtools [75]. The reference genome version used was GCA_002290385.1_ApisCC1.0_genomic.fa with a genome size of 228,791,026 bp.

Variation detection and annotation

A total of 130 samples were collected for variant detection. SNP detection was performed using SAMtools [75], and high-quality SNPs were filtered and screened by (1) SNPs with a support number above 3, (2) missing ratio < 10%, and (3) minor allele frequency > 1%. The SNP detection results were then functionally annotated using ANNOVAR software[76].

Population structures

Based on the results of the population structure analysis, the SNP data obtained from the Jilin and Heilongjiang samples were selected to perform genetic diversity index calculation, LD analysis, population history analysis, and selection elimination analysis. GCTA64 software was used for PCA, and the first three principal components were plotted. The distance matrix was calculated using Treebest-1.9.2 software, and based on this data, a phylogenetic tree was constructed using the neighbour-joining method. The bootstrap values of the tree were obtained after up to 1000 calculations.

Population genetics analyses

VCFtools was used to calculate π, Fst, He, Ho, and F. For calculating π, the parameter was set to “–window-pi 40,000 –window-pi-step 20,000”; and for calculating Fst, the parameter was set to “–fst-window-size 40,000 –fst-window-step 20,000”.

Linkage disequilibrium analysis

LD analysis was performed using Haploview software. The SNP spacing on the genome was fitted to the linkage disequilibrium coefficient (r2), and data with r2 greater than 0.1 were used as the decay value of the LD.

Effective population size

PSMC and SMC +  + methods were used. The parameters were set as -N30 -t15 -r5 -p 4 + 25*2 + 4 + 6, the mutation rate was set as 5.3 × 10−9, and the generation time was set as an estimate of 1 year [26].

Selective sweep

The A. cerana samples from Changbai Mountain (JL01-JL07, JL09-JL14 and HL01-HL05), Hainan (HI01-HI10), and Tibet (XZ01-XZ08) were selected for selection sweep analysis (Additional file 2: Table S2&Additional file 3: Table S3). We compared the Changbai Mountain individuals with the Hainan individuals and the Tibet individuals with the Hainan individuals. We calculated the genome-wide distribution of Fst values and θπ ratios for each sliding window (in 20-kb windows with 10-kb step size). The Fst values were Z-transformed, and the θπ ratios were log2-transformed. The putative selection targets were extracted based on being in the top 5% of log-odds ratios for both Fst and θπ. The GO enrichment analysis result was compared with the Pfam database to obtain GO terms, and enrichment analysis was performed using goseq. KOBAS was adopted for KEGG enrichment, the species used was A. mellifera, and the BlastX parameter was set to -evalue 1e-5 -m8.

Dissection and collection of morphometric data

10 worker bees were randomly selected from 30 sampling site, and 300 honeybee samples were collected, which basically covered all climatic regions in China (Additional file 1: Table S1). The right forewings and hindwings of the worker bees were removed with forceps, and the wings were fixed by pressing them in the middle by using two slides to make specimens. The specimens were then photographed using a German Leica LEICA-M165FC fluorescence microscopy imaging system. During image acquisition, all specimens were photographed under the same standard and scale.

The forewings and hindwings of the samples were digitally marked using TPS Util and TPS Dig software, respectively [77, 78], with 20 points for the forewing and 10 points for the hindwing (Fig. 8). All the landmarks were the intersection points of the stable wing veins [18, 79]. The obtained data of the wings were input into Coordgen 6.0 [79], with the 1st and 7th coordinate points as the baseline for the forewings and the 2nd and 10th coordinate points as the baseline for the hindwings. The coordinate values of the landmark were overprinted using the Generalized Procrustes Analysis to remove the influence of non-morphological factors, while the average profile information of each population was calculated on the basis of the GLS method [79].

Geometric morphometric analysis of wings

One-way ANOVA and LSD multiple tests were performed using SPSS software [79], and the results were presented as box plots. The correlation between wing size variation in elevation, longitude, and latitude was determined by calculating the CS of the wing and regressing the geographical information of the sample sites where they were located. The PCA analyses were performed using PAST software [80].

Availability of data and materials

All sequencing data generated during the current study have been submitted to the NCBI Short Read Archive (SRA) under the BioProject accession number PRJNA730587. All other data analysed in the present study are included in this published article (and its Supplementary Information files).

Abbreviations

A cerana:

Apis cerana

A mellifera:

Apis mellifera

JL-BS:

Jilin-Baishan

CS:

Centroid size

DNA:

Deoxyribonucleic acid

F:

Inbreeding coefficient

Fst:

F-statistics Genetic differences among population

GO:

Gene ontology

He:

Expected heterozygosity

HL:

Heilongjiang

Ho:

Observed Heterozygosity

JL-JL:

Jilin-Jilin

KEGG:

Kyoto Encyclopedia of Genes and Genomes

KOBAS:

KEGG Orthology Based Annotation System

LN:

Liaoning

LD:

Linkage disequilibrium

LM:

Landmark

mtDNA:

Mitochondrial deoxyribonucleic acid

MIS:

Marine isotope stages

PC:

Principal component

PCA:

Principal component analysis

π:

Nucleic acid diversity

PSMC:

Pairwise sequentially markovian coalescent

RH:

Raohe

SMC +  + :

Sequentially markovian coalescent

SSR:

Simple sequence repeats

JL-TH:

Jilin-Tonghua

TPS:

Thin plate spline

JL-YB:

Jilin-Yanbian

References

  1. Potts SG, Imperatriz-Fonseca V, Ngo HT, Aizen MA, Biesmeijer JC, Breeze TD, Dicks LV, Garibaldi LA, Hill R, Settele J, et al. Safeguarding pollinators and their values to human well-being. Nature. 2016;540(7632):220–9. https://doi.org/10.1038/nature20588.

    Article  CAS  PubMed  Google Scholar 

  2. Aslan CE, Zavaleta ES, Tershy B, Croll D. Mutualism disruption threatens global plant biodiversity: a systematic review. PLoS ONE. 2013;8(6): e66993. https://doi.org/10.1371/journal.pone.0066993.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Theisen-Jones H, Bienefeld K. The Asian honey bee (Apis cerana) is significantly in decline. Bee World. 2017;93(4):90–7. https://doi.org/10.1080/0005772X.2017.1284973.

    Article  Google Scholar 

  4. Chen C, Liu ZG, Luo YX, Xu Z, Wang SH, Zhang XW, Dai RG, Gao JL, Chen X, Guo HK, et al. Managed honeybee colony losses of the Eastern honeybee (Apis cerana) in China (2011–2014). Apidologie. 2017;48(5):692–702. https://doi.org/10.1007/s13592-017-0514-6.

    Article  CAS  Google Scholar 

  5. Radloff SE, Hepburn C, Hepburn HR, Fuchs S, Hadisoesilo S, Tan K, Engel MS, Kuznetsov V. Population structure and classification of Apis cerana. Apidologie. 2010;41(6):589–601. https://doi.org/10.1051/apido/2010008.

    Article  Google Scholar 

  6. Liu NN, Wang XM, Li JL, He JM. Production performance and biological characteristics of Changbai mountain Chinese bee during the blooming periods of Tilia. Spec Wild Econ Animal and Plant Res. 2020;42(6):5–9. https://doi.org/10.16720/j.cnki.tcyj.2020.06.002.

    Article  Google Scholar 

  7. Brown MJF, Paxton RJ. The conservation of bees: a global perspective. Apidologie. 2009;40(3):410–6. https://doi.org/10.1051/apido/2009019.hal-00892052.

    Article  Google Scholar 

  8. Zayed A. Bee genetics and conservation. Apidologie. 2009;40(3):237–62. https://doi.org/10.1051/apido/2009026.

    Article  Google Scholar 

  9. The National Animal Genetic Resources Committee. Animal genetic resources in China bees. Beijing: China Agriculture Press; 2011.

    Google Scholar 

  10. Balanya J, Huey RB, Gilchrist GW, Serra L. The chromosomal polymorphism of Drosophila subobscura: a microevolutionary weapon to monitor global change. Heredity (Edinb). 2009;103(5):364–7. https://doi.org/10.1038/hdy.2009.86.

    Article  CAS  Google Scholar 

  11. Liu Z, Tan X, Orozco-terWengel P, Zhou X, Zhang L, Tian S, Yan Z, Xu H, Ren B, Zhang P et al: Population genomics of wild Chinese rhesus macaques reveals a dynamic demographic history and local adaptation, with implications for biomedical research. Gigascience 2018, 7(9). http://dx.doi.org/https://doi.org/10.1093/gigascience/giy106

  12. Montero-Mendieta S, Tan K, Christmas MJ, Olsson A, Vilà C, Wallberg A, Webster MT. The genomic basis of adaptation to high-altitude habitats in the eastern honey bee (Apis cerana). Mol Ecol. 2019;28(4):746–60. https://doi.org/10.1111/mec.14986.

    Article  PubMed  Google Scholar 

  13. Bai Y, Dong JJ, Guan DL, Xie JY, Xu SQ. Geographic variation in wing size and shape of the grasshopper Trilophidia annulata (Orthoptera: Oedipodidae): morphological trait variations follow an ecogeographical rule. Sci Rep. 2016;6:32680. https://doi.org/10.1038/srep32680.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Ning Y, Duan C, Chang QQ, Hou XH. Geometric morphometric analysis of wing shape variation and the genetic relationships of six species of the subgenus Culicoides (Diptera: Ceratopogonidae). Journal of Environmental Entomology. 2020;42(1):147–51. https://doi.org/10.1603/me10110.

    Article  Google Scholar 

  15. Shevtsova E, Hansson C, Janzen DH, Kjaerandsen J. Stable structural color patterns displayed on transparent insect wings. Proc Natl Acad Sci USA. 2011;108(2):668–73. https://doi.org/10.1073/pnas.1017393108.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Francoy TM, Wittmann D, Drauschke M, Müller S, Steinhage V, Bezerra-Laure MAF, De Jong D, Gonçalves LS. Identification of Africanized honey bees through wing morphometrics: two fast and efficient procedures. Apidologie. 2008;39(5):488–94. https://doi.org/10.1051/apido:2008028.

    Article  Google Scholar 

  17. Tofilski A. Using geometric morphometrics and standard morphometry to discriminate three honeybee subspecies. Apidologie. 2008;39(5):558–63. https://doi.org/10.1051/apido:2008037.

    Article  Google Scholar 

  18. Eimanifar A, Brooks SA, Bustamante T, Ellis JD. Population genomics and morphometric assignment of western honey bees (Apis mellifera L.) in the Republic of South Africa. BMC Genomics. 2018;19(1):615. https://doi.org/10.1186/s12864-018-4998-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Gauld I, Bolton B. The Hymenoptera. Oxford: Oxford University Press; 1988.

    Google Scholar 

  20. BF Zhou XJ Zhu Y Li 2011 New mutations in hind wing vein of Apis cerana cerana (Hymenoptera: Apidae) induced by lower developmental temperature Acta Ecologica Sinica 31 5 1387 92 CNKI:SUN:STXB.0.2011-05-023

  21. Han F, Wallberg A, Webster MT. From where did the Western honeybee (Apis mellifera) originate? Ecol Evol. 2012;2(8):1949–57. https://doi.org/10.1002/ece3.312.Epub2012Jul12.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Büchler R, Costa C, Hatjina F, Andonov S, Meixner MD, Conte YL, Uzunov A, Berg S, Bienkowska M, Bouga M, et al. The influence of genetic origin and its interaction with environmental effects on the survival of Apis mellifera L. colonies in Europe. Journal of Apicultural Research. 2014;53(2):205–14. https://doi.org/10.3896/IBRA.1.53.2.03.

    Article  Google Scholar 

  23. Meixner MD, Büchler R, Costa C, Francis RM, Hatjina F, Kryger P, Uzunov A, Carreck NL. Honey bee genotypes and the environment. J Apic Res. 2015;53(2):183–7. https://doi.org/10.3896/IBRA.1.53.2.07.

    Article  Google Scholar 

  24. Zammit-Mangion M, Meixner M, Mifsud D, Sammut S, Camilleri L. Thorough morphological and genetic evidence confirm the existence of the endemic honey bee of the Maltese Islands Apis mellifera ruttneri: recommendations for conservation. J Apic Res. 2017;56(5):514–22. https://doi.org/10.1080/00218839.2017.1371522.

    Article  Google Scholar 

  25. Park D, Jung JW, Choi BS, Jayakodi M, Lee J, Lim J, Yu Y, Choi YS, Lee ML, Park Y, et al. Uncovering the novel characteristics of Asian honey bee, Apis cerana, by whole genome sequencing. BMC Genomics. 2015;16:1–16. https://doi.org/10.1186/1471-2164-16-1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Chen C, Liu ZG, Pan Q, Chen X, Wang HH, Guo HK, Liu SD, Lu HF, Tian SL, Li RQ, et al. Genomic analyses reveal demographic history and temperate adaptation of the newly discovered honey bee subspecies Apis mellifera sinisxinyuan n. ssp. Mol Biol Evol. 2016;33(5):1337–48. https://doi.org/10.1093/molbev/msw017.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Chen C, Wang HH, Liu ZG, Chen X, Tang J, Meng F, Shi W. Population genomics provide insights into the evolution and adaptation of the eastern honey bee (Apis cerana). Mol Biol Evol. 2018;35(9):2260–71. https://doi.org/10.1093/molbev/msy130.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;19(9):1655–64. https://doi.org/10.1101/gr.094052.109.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: a tool for genome-wide complex trait analysis. Am J Hum Genet. 2011;88(1):76–82. https://doi.org/10.1016/j.ajhg.2010.11.011.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Wright S. Evolution and the genetics of populations. Chicago: University of Chicago Press; 1978.

    Google Scholar 

  31. Li YC, Chao TL, Fan YH, Lou DL, Wang GZ. Population genomics and morphological features underlying the adaptive evolution of the eastern honey bee (Apis cerana). BMC Genomics. 2019;20(1):869. https://doi.org/10.1186/s12864-019-6246-4.

    Article  CAS  Google Scholar 

  32. Li H, Durbin R. Inference of human population history from individual whole-genome sequences. Nature. 2011;475(7357):493–6. https://doi.org/10.1038/nature10231.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Terhorst J, Kamm JA, Song YS. Robust and scalable inference of population history from hundreds of unphased whole genomes. Nat Genet. 2017;49(2):303–9. https://doi.org/10.1038/ng.3748.

    Article  CAS  PubMed  Google Scholar 

  34. Yu SH, Yang P, Sun T, Qi Q, Wang XQ, Chen XM, Feng Y, Liu BW. Transcriptomic and proteomic analyses on the supercooling ability and mining of antifreeze proteins of the Chinese white wax scale insect. Insect Sci. 2016;23(3):430–7. https://doi.org/10.1111/1744-7917.12320.

    Article  CAS  PubMed  Google Scholar 

  35. Xu K, Niu QS, Zhao HT, Du YL, Jiang YS. Transcriptomic analysis to uncover genes affecting cold resistance in the Chinese honey bee (Apis cerana cerana). PLoS ONE. 2017;12(6): e0179922. https://doi.org/10.1371/journal.pone.0179922.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Choe SE, Nguyen TT, Hyun BH, Noh JH, Lee HS, Lee CH, Kang SW. Genetic and phylogenetic analysis of South Korean sacbrood virus isolates from infected honey bees (Apis cerana). Vet Microbiol. 2012;157:32–40. https://doi.org/10.1016/j.vetmic.2011.12.007.

    Article  CAS  PubMed  Google Scholar 

  37. Hao LN, Zhou SJ, Zhu XJ, Xu XJ, Cai ZB, Niu QS, Yuan CY, Chen DH, Zhou BF. Genetic differentiation and genetic diversity analysis of Apis cerana in northeast China. J of Northeast Agric Un. 2019;50(9):35–43. https://doi.org/10.19720/j.cnki.issn.1005-9369.2019.09.005.

    Article  Google Scholar 

  38. Yu YL, Zhou SJ, Xu XJ, Zhu XJ, Zhou BF. Analysis on genetic diversity of Apis cerana cerana in Changbai Mountains. J of Fujian Apiculture and Forestry Univ (Nature Science Edition). 2013;42(6):643–7. https://doi.org/10.13323/j.cnki.j.fafu(nat.sci.).2013.06.021.

    Article  Google Scholar 

  39. FC Gao 2015 Study on genetic diversity of Apis cerana cerana populations from Changbai Mountain by microsatellite makers Apiculture of China. 66 7 20 2 CNKI:SUN:ZGYF.0.2015-07-008

  40. Whitfield C W, Behura S K, Berlocher S H, Clark A G, Johnston J S, Sheppard W S, Smith D R, Suarez A V, Weaver D, Tsutsui N D. Thrice out of Africa: ancient and recent expansions of the honey bee. Apis mellifera Science. 2006;314(5799):642–5. https://doi.org/10.1126/science.1132772.

    Article  CAS  PubMed  Google Scholar 

  41. Fan QC, SUi JL, Wang TH, Li N, Sun Q. Eurption history and magma evolution of the trachybasalt in the Tianchi volcano Changbaishan. Acta Petrologica Sinica. 2006;22(6):1449–57. https://doi.org/10.3321/j.issn:1000-0569.2006.06.001.

    Article  CAS  Google Scholar 

  42. Wang ZB, Wang JY, Li BJ. Copupling relation between cyclic eruption of Quaternary volcano and ancient glacier in Tianchi. Changbai Mts East China Geology. 2020;41(2):108–15. https://doi.org/10.16788/j.hddz.32-1865/P.2020.02.002.

    Article  CAS  Google Scholar 

  43. JJ Li 1992 Glacial relics of monsoonal Asia in the last glaciation Quaternary Sci. 4 332 40 CNKI:SUN:DSJJ.0.1992-04-004

  44. Ono Y, Aoki T, Hasegawa H, Dali L. Mountain glaciation in Japan and Taiwan at the global Last Glacial Maximum. Quatern Int. 2005;138–139:79–92. https://doi.org/10.1016/j.quaint.2005.02.007.

    Article  Google Scholar 

  45. Pan B. Eruptive history of the Changbaishan Tianchi Volcano, China/DPRK, since the Last Pleistocene. Recent Dev in World Seismology. 2017;11:39–41. https://doi.org/10.1016/j.jvolgeores.2020.106870.

    Article  CAS  Google Scholar 

  46. XQ Li C Zhao XY Zhou 2019 Vegetation pattern of Northeast China during the special periods since the Last Glacial Maximum Scientia Sinica Terrae. 49 8 1213 30 CNKI:SUN:JDXK.0.2019-08-004

  47. Wang SW, Cai JN, Zhu JH, Gong DY. Studies on climate change in China. Climatic and Environ Res. 2002;7(2):137–45. https://doi.org/10.3969/j.issn.1006-9585.2002.02.002.

    Article  CAS  Google Scholar 

  48. Cui ZX. The study on volcanic eruption and forest conflagration since holocene quaternary in Changbai MT. Geographical Research. 1997;16(1):92–7. https://doi.org/10.11821/yj1997010012.

  49. Li JT. Synopsis of research of volcanic disaster concerning Changbai Mountain. Seismological Res of Northeast China. 1998;14(2):59–66. https://doi.org/10.6038/cjg20160314.

    Article  Google Scholar 

  50. Machida H, Arai F. Extensive ash falls in and around the Sea of Japan from large late Quaternary eruptions. J Volcanol Geoth Res. 1983;18(1–4):151–64. https://doi.org/10.1016/0377-02738390007-0.

    Article  Google Scholar 

  51. Sun CQ, Plunkett G, Liu JQ, Zhao HL, Sigl M, McConnell JR, Pilcher JR, Vinther B, Steffensen JP, Hall V. Ash from Changbaishan Millennium eruption recorded in Greenland ice: implications for determining the eruption’s timing and impact. Geophysical Res Lett. 2014;41(2):694–701. https://doi.org/10.17521/cjpe.2015.0216.

    Article  Google Scholar 

  52. Machida H, Moriwaki H, Zhao DC. The recent major eruption of Changbai Volcano and its environmental effects. Geogr rep of Tokyo Metropolitan Univ. 1990;25:1–20.

    Google Scholar 

  53. QJ Liu Z Wang SX Wang 1993 Effects of recent volcanic eruptions on alpine and subalpine vegetation in the Changbai Moutain Scientia Geographica Sinica. 13 1 57 6 110.13249/j.cnki.sgs.1993.01.011

  54. Liu JQ, Chen SS, Guo WF, Sun CQ, Zhang ML, Guo ZF. Research advances in the Mt. Changbai Volcano. Bull of Mineralogy Petrology and Geochem. 2015;34(4):710–23. https://doi.org/10.3969/j.issn.l007-2802.2015.04.005.

    Article  CAS  Google Scholar 

  55. Mu KH, Zhang W, Niu YB, Yan L, Li CC. Research on the climatic change during the Holocene in Northeastern China. Territory&Natural Resour Study. 2008;2:48–50. https://doi.org/10.16202/j.cnki.tnrs.2008.02.001.

    Article  Google Scholar 

  56. WY Li YL Liang 1985 Vegetation and environment of the hypsithermal interval of holocene in Eastern Hebei plain Acta Botanica Sinica. 27 6 640 51 CNKI:SUN:ZWXB.0.1985-06-013

  57. Jing XH, Kang L. Geographical variation in egg cold hardiness: a study on the adaptation strategies of the migratory locust Locusta migratoria L. Ecological Entomol. 2003;28:151–8. https://doi.org/10.1046/j.1365-2311.2003.00497.x.

    Article  Google Scholar 

  58. Zhao XM, Jia P, Gou X, Liu WM, Ma EB, Zhang JZ. Expression and functional analysis of endocuticle structural glycoprotein gene LmAbd-5 in Locusta migratoria. Scientia Agricultura Sinica. 2017;50(10):1817–26. https://doi.org/10.3864/j.issn.0578-1752.2017.10.007.

    Article  Google Scholar 

  59. Francis RM, Nielsen SL, Kryger P. Varroa-virus interaction in collapsing honey bee colonies. PLoS ONE. 2013;8(3): e57540. https://doi.org/10.1371/journal.pone.0057540.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Liu L, Gong Z, Guo X, Xu B. Cloning, structural characterization and expression analysis of a novel lipid storage droplet protein-1 (LSD-1) gene in Chinese honeybee (Apis cerana cerana). Mol Biol Rep. 2012;39(3):2665–75. https://doi.org/10.1007/s11033-011-1020-y.

    Article  CAS  PubMed  Google Scholar 

  61. Underwood BA. Thermoregulation and energetic decision-making by the honeybees Apis cerana, Apis dorsata and Apis laboriosa. J Exp Biol. 1991;157:19–34. https://doi.org/10.1242/jeb.157.1.19.

    Article  Google Scholar 

  62. Bai Y, Dong FD, Bao KO, Qiu AN, Liang RJ, Wang H. Using geometric morphometrics to quantify the sexual dimorphism of Pantala flavescens. Chin J of Appl Entomol. 2015;52(2):363–9. https://doi.org/10.7679/j.issn.2095-1353.2015.040.

    Article  Google Scholar 

  63. Greenleaf SS, Williams NM, Winfree R, Kremen C. Bee foraging ranges and their relationship to body size. Oecologia. 2007;153(3):589–96. https://doi.org/10.1007/s00442-007-0752-9.

    Article  PubMed  Google Scholar 

  64. Abou-Shaara HF, Ahmed ME. Characterisation and tracking changes of morphological characteristics in honey bee, Apis mellifera, Journal of Entomological and Acarological Research. 2015;47(3)https://doi.org/10.1016/j.sjbs.2021.02.050.

  65. Zhu X, Zhou S, Xu X, Wang J, Yu Y, Yang K-c, Luo Q, Xu Y, Wang S, Zhou B. Morphological differentiation in Asian honey bee (Apis cerana) populations in the basin and highlands of southwestern China. J of Apicultural Res. 2017;56(3):203–9. https://doi.org/10.1080/00218839.2017.1306374.

    Article  Google Scholar 

  66. Gloag R, Ding G, Christie JR, Buchmann G, Beekman M, Oldroyd BP. An invasive social insect overcomes genetic load at the sex locus. Nat Ecol Evol. 2016;1(1):11. https://doi.org/10.1038/s41559-016-0011.

    Article  PubMed  Google Scholar 

  67. Tan K, Ge FC, Zhao R, He SY. Morphometric characters of Changbaishan A. cerana. Journal of Bee. 2004;869(6):8–9. https://doi.org/10.3969/j.issn.1003-9139.2004.06.002.

    Article  Google Scholar 

  68. Zhang LJ, Fan X, Tan K, He SY. Morphometric study of Apis cerana in Changbai Moutain. J of Yunnan Agric Univ. 2006;21(4):511–6. https://doi.org/10.16211/j.issn.1004-390x(n).2006.04.024.

    Article  Google Scholar 

  69. Zhou SJ, Zhu XJ, Xu XJ, Zheng XJ, Zhou BF. Assessing of Geometric Morphometrics Analyses in Microtaxonomy of the Apis cerana Fabricius (Hymenoptera: Apidae) within China. J Kansas Entomol Soc. 2016;89(4):297–305. https://doi.org/10.2317/0022-8567-89.4.297.

    Article  Google Scholar 

  70. Ji YK, Li XA, Ji T, Tang JB, Qiu LF, Hu JH, Dong JX, LuoLiuFrandsen SQSLPB, et al. Gene reuse facilitates rapid radiation and independent adaptation to diverse habitats in the Asian honeybee. Sci Adv. 2020;6(51):eabd3590. https://doi.org/10.1126/sciadv.abd3590.

    Article  CAS  PubMed  Google Scholar 

  71. De la Rúa P, Jaffé R, Dall’Olio R, Muñoz I, Serrano J. Biodiversity conservation and current threats to European honeybees. Apidologie. 2009;40(3):263–84. https://doi.org/10.1186/s12711-019-0518-y.

    Article  Google Scholar 

  72. Muñoz I, Henriques D, Jara L, Johnston JS, Chavez-Galarza J, De La Rua P, Pinto MA. SNPs selected by information content outperform randomly selected microsatellite loci for delineating genetic identification and introgression in the endangered dark European honeybee (Apis mellifera mellifera). Mol Ecol Resour. 2017;17(4):783–95. https://doi.org/10.1111/1755-0998.12637.

    Article  CAS  PubMed  Google Scholar 

  73. Diao Q, Sun L, Zheng H, Zeng Z, Wang S, Xu S, Zheng H, Chen Y, Shi Y, Wang Y, et al. Genomic and transcriptomic analysis of the Asian honeybee Apis cerana provides novel insights into honeybee biology. Sci Rep. 2018;8(1):822. https://doi.org/10.1038/s41598-017-17338-6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60. https://doi.org/10.1093/bioinformatics/btp324.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. Genome project data processing S: The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9. https://doi.org/10.1093/bioinformatics/btp352.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010;38(16): e164. https://doi.org/10.1093/nar/gkq603.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Rohlf, FJ. tpsDIG, version 2.11. Department of Ecology and Evolution, State University of New York; Stony Brook, New York; 2008b.

  78. Rohlf, FJ. tpsUtil, version 1.40. Department of Ecology and Evolution, State University of New York; Stony Brook, New York; 2008a.

  79. Qiu LF, Wei ZM, Wang JJ, Liu HX, Zhao BX, Zhang Y, Lian ZM. Geometric analysis of morphological variation in the wing of Apis cerana cerana from the Qinling-Daba Mountain Areas. Chin J of Appl Entomol. 2018;55(3):503–13. https://doi.org/10.7679/j.issn.2095-1353.2018.064.

    Article  Google Scholar 

  80. Bai Yi, Dai DF, Bao KO, Qiu AN, Liang RJ, Wang HR. Using geometric morphometrics to quantify the sexual dimorphism of Pantala flavescens. Chin J of Appl Entomol. 2015;52(2):363–9.

    Google Scholar 

Download references

Acknowledgements

The authors are especially grateful to all beekeepers, friends of apiculture in the sampling area, and other contributors who provided the valuable samples for this work. We also wish to give special thanks to all technicians in the project, involved State Key Laboratory for Molecular Biology of Special Economic Animals and Key Laboratory for Bee Genetics and Breeding of Jilin Province. Finally, we would also thank Novegene Technology Company Limited for their help.

Funding

This work was financially supported by Funds of Central government guided local science and technology development funds for Jilin Province basic research project (No. 202002060JC), and Jilin Province science and technology development Project (No. 20200703002ZP). The funders provided the financial support to the research, but had no role in the design of the study, analysis, interpretations of data and in writing the manuscript.

Author information

Authors and Affiliations

  1. State Key Laboratory for Molecular Biology of Special Economic Animals, Institute of Special Animal and Plant Sciences of Chinese Academy of Agricultural Sciences, Changchun, People’s Republic of China

    Liu Nannan, Liu Huamiao, Ju Yan, Li Yang, Wang Tianjiao & Xing Xiumei

  2. Key Laboratory for Bee Genetics and Breeding of Jilin Province, Apiculture Science Institute of Jilin Province, Jilin, People’s Republic of China

    Liu Nannan, Li Xingan, He Jinming & Niu Qingsheng

Authors
  1. Liu Nannan
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Liu Huamiao
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Ju Yan
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Li Xingan
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Li Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Wang Tianjiao
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. He Jinming
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Niu Qingsheng
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Xing Xiumei
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

XX, NQ and LN designed and coordinated the research. LN and HJ conducted laboratory work of wing geometric morphology. LH, LX, JY, LY and WT analyzed the sequencing data. XX, NQ and LN drafted the paper. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Niu Qingsheng or Xing Xiumei.

Ethics declarations

Ethics approval and consent to participate

All animal experiments were performed in accordance with the “Guidelines for Experimental Animals” of the Ministry of Science and Technology (Beijing, China), and approved by the Institutional Animal Care and Use Ethics Committee of Institute of Special Animal and Plant Sciences of CAAS.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Nannan, L., Huamiao, L., Yan, J. et al. Geometric morphology and population genomics provide insights into the adaptive evolution of Apis cerana in Changbai Mountain. BMC Genomics 23, 64 (2022). https://doi.org/10.1186/s12864-022-08298-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-022-08298-x

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us