Skip to main content
Download PDF

A key gene for the climatic adaptation of Apis cerana populations in China according to selective sweep analysis

BMC Genomics volume 24, Article number: 100 (2023) Cite this article

Abstract

Background

Apis cerana is widely distributed in China and, prior to the introduction of western honeybees, was the only bee species kept in China. During the long-term natural evolutionary process, many unique phenotypic variations have occurred among A. cerana populations in different geographical regions under varied climates. Understanding the molecular genetic basis and the effects of climate change on the adaptive evolution of A. cerana can promote A. cerana conservation in face of climate change and allow for the effective utilization of its genetic resources.

Result

To investigate the genetic basis of phenotypic variations and the impact of climate change on adaptive evolution, A. cerana workers from 100 colonies located at similar geographical latitudes or longitudes were analyzed. Our results revealed an important relationship between climate types and the genetic variation of A. cerana in China, and a greater influence of latitude compared with longitude was observed. Upon selection and morphometry analyses combination for populations under different climate types, we identified a key gene RAPTOR, which was deeply involved in developmental processes and influenced the body size.

Conclusion

The selection of RAPTOR at the genomic level during adaptive evolution could allow A. cerana to actively regulate its metabolism, thereby fine-tuning body sizes in response to harsh conditions caused by climate change, such as food shortages and extreme temperatures, which may partially elucidate the size differences of A. cerana populations. This study provides crucial support for the molecular genetic basis of the expansion and evolution of naturally distributed honeybee populations.

Peer Review reports

Background

Over the past few decades, the numbers of many insects have shown a rapid decline [1], including those of bumble bees [2], caterpillars [3], beetles [4], and honeybees [5]. The sources of insect survival pressure leading to these declines were abundant and included invasive species, habitat loss, pesticide abuse, environmental pollution, and climate change [6]. Of these factors, climate change, or its combination with all other factors, is more likely to be the main cause of large-scale global insect population and diversity declines [7,8,9,10,11,12]. However, evolutionary adaptations can potentially help species cope with survival pressure brought about by climate change [13]. Apis cerana is an important crop-pollinating insect in China. Its long collection cycle and ability to utilize sporadic nectar sources have enabled A. cerana to play a key role in the development of Chinese agriculture. During the long-term natural evolutionary process, many unique phenotypic variations have occurred among A. cerana populations in different regions under different climatic conditions [14, 15]. In general, moving from south to north in geographical distribution, the phenotype of A. cerana gradually shows a trend of weakening migration, slightly larger size, darker body color, and stronger cold resistance [16]. The characteristics of this adaptive radiation distribution under different climatic conditions show exceptional stress resistance, making A. cerana an excellent research subject for the study of climatic adaptation on adaptive evolution [17]. Thus, the study of the climatic adaptation of A. cerana can help effectively conserve A. cerana against climate change challenges. Simultaneously, understanding the genetic basis of these unique phenotypic variations is beneficial for their conservation and provides insight to better harness their genetic resources.

So far, many studies have used different methods to elucidate the basis of the unique phenotypic differences in A. cerana. Based on direct measurements of the phenotypic and morphological characteristics of A. cerana in different regions [18], nine varieties were classified. With the development of molecular biology technology, population genetics studies revealing genetic diversity based on sequencing the individual genomes of different bee species have been widely used [19,20,21]. Once the genome of A. cerana was published in 2015 (GCA_001442555.1) [22], genome research entered a rapid stage (GCA_002290385.1 and GCA_011100585.1) [23, 24]. Studies on the origin and evolution [14, 15, 17, 18, 25] and sequencing of A. cerana with special traits—A. cerana abansis (Aba strain) [26] and the endemic A. cerana with cold resistance [27]—progressed rapidly. These studies have also demonstrated distinct variations in DNA and microsatellites among A. cerana from different regions. A 2020 study identified the role of the leucokinin receptor (Lkr) in influencing the foraging division of labor [18]. However, existing research is far from sufficient to fully elucidate the climatic adaptation during the evolutionary process of A. cerana.

In this study, we collected 100 samples from 10 regions in China (covering the most important climate types) and utilized population genetics methods to identify the key traces in the process of adaptation to different environments at the genome level; we aimed to reveal the genetic basis of the strong environmental adaptability of A. cerana. Simultaneously, we chose collection sites at similar latitudes or longitudes to explore the influence of geographical coordinate factors on the adaptive radiation of A. cerana. Upon genetic structure, genetic differentiation, and population diversity analyses of these 10 populations, as well as through selective sweep analysis and morphometry analysis combination of the populations with high genetic differentiation, we were able to identify the key signaling pathways and genes involved in the process of adaptation to different climatic conditions in China. This study provides a useful reference for the improvement of A. cerana conservation and breeding.

Results

Sampling and sequencing

We collected and resequenced samples from 10 regions of China and included previously published resequencing data from three regions (Table 1, Fig. 1). A total of 429.71 Gb of raw data were obtained from 100 samples. The total data volume after quality control was 418.97 Gb, the average alignment rate with the reference genome was 90.98%, and the average sequencing depth was 19.16 × . The average Q20 was 96.65% and the Q30 averaged 91.59%.

Table 1 Information on sampling sites
Full size table
Fig. 1
figure 1

Sampling sites of Apis cerana. Samples from regions in blue were obtained from previously published data

Full size image

Population genetic structure analysis

The admixture model-based software (Admixture) was used to estimate the population structure. With K = 2, MK, JL, and WC formed an ancestral cluster, while the other populations showed different degrees of mixed lineage. As K increased from 3 to 5, MK, JL, and WC showed distinct lineages from the other populations, with ZZ and MY forming an ancestral cluster (Fig. 2). Simultaneously, we calculated the cross-validation error rate, which was the smallest when K = 4, of genetic structure analysis of 100 samples (Fig. S1). In particular, the SZ and QM groups remained indistinguishable when K increased from 4 to 9 (Fig. 2). The population subdivision pattern classified by principal component analysis (PCA) showed similar results. According to the first and second principal components, MK, JL, and WC were separated clearly (Fig. 3). The results of the group structure analysis and PCA of the 10 regional groups demonstrated that JL, MK, and WC could be separated from the remaining closely related and difficult-to-separate populations. The JL, MK, and WC groups were highly differentiated compared with the other groups. Based on the obtained SNP information, we used neighbor-joining methods to construct a phylogenetic tree (Fig. 4).

Fig. 2
figure 2

Analysis of population structure

Full size image
Fig. 3
figure 3

Principal component analysis of the 10 populations

Full size image
Fig. 4
figure 4

Phylogenetic tree of Apis cerana

Full size image

Genetic differentiation and genetic diversity

To understand the genetic differentiation between the 10 groups, we calculated pairwise FST (Table 2) and genetic diversity parameters, including Ho, He, and PIC (Table 3). The maximum value of FST was 0.39 (between WC and JL), and the minimum value was less than 0.01 (a value close to 0, rounded to two decimal places is 0.00). These results also showed that JL, MK, and WC were in a state of high genetic differentiation (FST > 0.15), which is consistent with the analysis of population structure and PCA. JL displayed the highest degree of genetic differentiation (average FST = 0.28), followed by WC (average FST = 0.19) and MK (average FST = 0.16). The remaining regions were in a state of low genetic differentiation. The parameter calculation results for genetic diversity were greatest for SZ (PIC = 0.205) and GZ (PIC = 0.206) among the 10 populations. The genetic diversity of JL (PIC = 0.099), MK (PIC = 0.140), and WC (PIC = 0.148) was lower than that of other regions, indicating that bees in these three areas have been subjected to a higher intensity of natural selection.

Table 2 Pairwise FST distance between 10 populations of Apis cerana
Full size table
Table 3 Genetic diversity parameters of Apis cerana
Full size table

Historical effective population sizes

Based on the above population genetic structure analysis results, PCA, and genetic differentiation analysis, JL, MK, and WC have each independently become a subpopulation separated from the remaining populations. Moreover, we noticed that these three populations, in addition to a fourth group comprising the remaining seven populations, each occupy different regional climate types. Jilin is in the northeastern part of China and has a temperate continental monsoon climate; Mangkang is on the Qinghai-Tibet Plateau of China and has a typical plateau climate; Wenchang is on Hainan Island, China, with a unique tropical monsoon island-ridge climate; and the remaining populations are mostly in the subtropical monsoon climate region. The SZ group was selected as the representative subtropical monsoon climate region. To understand the important evolutionary events that occurred during the past adaptation of these four highly differentiated populations of A. cerana, we estimated their historical effective population sizes (Fig. 5). According to the estimated results, the four regions had nearly the same effective colony size during the period 80–200 ka years ago. About 70 ka years ago, the JL group began to show a change in population size that was different from that of the rest of the group. The population size of the three regions of MK, WC, and SZ in the middle and low latitudes showed an upward trend at this time, which seems to indicate that the population in these regions could begin to have continuous differentiation of population structure and mixed blood, resulting in an increase in population size, whereas the JL region did not undergo this process. Overall, the evolution of populations at high latitudes showed different trends than those at low and middle latitudes.

Fig. 5
figure 5

Estimation of historical effective population sizes

Full size image

Morphometric analysis

We measured 10 morphological indicators related to body size on the samples collected from WC, MK, and SZ for morphometric analysis, including the right forewing length (FL), the right forewing width (FB), the sixth sternum length (L6), the sixth sternum width (T6), the third sternum length (S3), the third tergum length and the fourth tergum length (T3 + T4), the femur length (Fe), the tibia length (Ti), the basitarsus length (ML), and the basitarsus width (MT) (Fig. S2). The results of one-way ANOVA on the measured 10 indicators showed that, except for the insignificant difference between the MT of the SZ and WC groups (P > 0.05), all other indicators were significantly different among the three groups (P < 0.01) (Table 4 and Fig. 6). These results confirm that, in terms of body size, MK population is significantly larger than that of SZ, which in turn is significantly larger than that of WC.

Table 4 Summary morphological data of the three geographical populations
Full size table
Fig. 6
figure 6

Boxplot of 10 morphological indicators of Apis cerana among the three geographical populations. A The right forewing length (FL). B The right forewing width (FB). C The sixth sternum length (L6). D The sixth sternum width (T6). E The third sternum length (S3). F The third tergum length and the fourth tergum length (T3 + T4). G The femur length (Fe). H The tibia length (Ti). I The basitarsus length (ML). J The basitarsus width (MT). Note: ** means significant difference at 0.01 level

Full size image

Selective sweep analysis

To further study the adaptive radiation distribution, such as the difference in body size of A. cerana from different regions of China and the genome-level changes that occurred during the evolution of adaptation to the unique climate of each region under natural selection pressure, we estimated pairwise genetic differentiation (FST) and differences in nucleotide diversity (π ratios) from the four different types of populations to identify key selective sweeps. For the 40 samples from JL, MK, WC, and SZ, we performed the linkage disequilibrium (LD) decay analysis. The LD analysis results showed that when r2 was half of the maximum value, the decay rate was MK > JL > WC > SZ (Fig. S3). This also indicates that the bees in MK, JL, and WC were affected by more intense selection pressure than those in SZ; therefore, we selected the SZ region with richer genetic diversity and less selection pressure as the reference population. The JL, MK, and WC populations were scanned and analyzed against the SZ reference population to identify the regions in the A. cerana genome that have been selected under natural selection pressure, thereby revealing the adaptive evolution of A. cerana.

Selective sweep regions were chosen according to the intersection of two indices (FST and π ratios) with a threshold of the top 5% level. The results of the selection signal analysis of the highly differentiated populations under different climate types, with the SZ population (subtropical monsoon oceanic climate) as the reference population, identified 839 candidate regions involving 527 genes (Fig. 7A and Table S3) in the JL population (temperate monsoon climate). In the MK population (plateau climate), 589 candidate regions involving 565 genes were identified (Fig. 7B and Table S4). In the WC population (tropical monsoon island climate), 224 candidate regions involving 311 genes were identified (Fig. 7C and Table S5). In addition, 33 genes were identified in the selection signal analysis results for all three regions (Fig. 7D and Table S6). These 33 genes may play important roles in the adaptation of A. cerana to different climatic conditions in China. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed on these 33 candidate genes. The 33 shared genes were involved in 846 GO and 44 KEGG pathways, with 138 GO and 6 KEGG pathways significantly enriched at a threshold of P < 0.05 (Table S7 and Table S8). GO enrichment analysis revealed that the biological process with the most enriched genes was the single-organism process, the cellular component with the most enriched genes were cell and cell part, and the molecular function with the most enriched genes was binding (Fig. S4). The six signaling pathways were amino sugar and nucleotide sugar metabolism, peroxisome, RNA degradation, AMPK signaling pathway, ubiquitin-mediated proteolysis, and 2-oxocarboxylic acid metabolism. These results indicate that A. cerana responds to natural environmental stress primarily by regulating sugar and protein metabolism. In addition, among the top 20 enriched signaling pathways, we noticed some that may be related to the adaptation of honeybees to regional climates, such as the citrate cycle (TCA cycle) related to energy metabolism in the body, the thermogenesis pathway related to adaptation to changes in ambient temperature, and the Hedgehog signaling pathway involved in regulating gene expression (Fig. S5).

Fig. 7
figure 7

Selective sweep analysis against the reference SZ population. A Results of selective sweep analysis (SZ/JL). B Results of selective sweep analysis (SZ/MK). C Results of selective sweep analysis (SZ/WC). The selected intersect regions in red areas are at the top 5% level by FST and π ratios. D Venn diagram of selected genes

Full size image

Discussion

Our study results reveal an important link between climate type and the adaptive evolution of A. cerana in China. A. cerana with unique phenotypic variations from different climate regions (Jilin, Mangkang, and Wenchang) could have independently formed a subpopulation from the analyses of the population structure, genetic diversity, and genetic differentiation based on whole-genome resequencing data. Furthermore, we confirmed the significant differences in the body size of A. cerana by comparing morphological data. This evidence shows that the JL, MK, and WC populations were highly differentiated compared with other populations, which may be a consequence of different climate types and environmental selection pressures. This is in line with the findings of a study on the adaptive evolution of A. cerana in which isolated groups showed a high degree of differentiation and lower genetic diversity, whereas other groups were less differentiated and had higher genetic diversity [14]. In the present study, we found that the average PIC of the 10 populations was 0.17, and the PIC values of all populations were lower than 0.5, representing low polymorphism and low genetic diversity of A. cerana. Human activities are considered to be an important reason for the decrease in the genetic diversity of A. cerana in recent years. Human activities led to the reduction of the habitat area of A. cerana, which directly led to the reduction of wild resources of A. cerana. 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 A. cerana. In addition, the introduction of A. cerana from other regions also affected the original genetic structure and genetic diversity of local populations. A rich genetic diversity can enhance the adaptability of species to their environment, more stably maintain population balance [28, 29], and improve production performance [30]. Therefore, it is imperative to protect and rationally utilize the genetic resources of A. cerana.

It is also worth noting that the geographic latitude factor had a high degree of influence on the biodiversity of biological groups. This is supported by the estimations of the effective population size. The different change patterns of the effective population size of A. cerana at high and low latitudes also reflect the different impacts of sudden climate change on populations at different latitudes. Our study shows that this discrepancy first originated around 70 to 80 ka years ago. Since the last glacial period of the Quaternary Ice Age [31, 32], different degrees of cooling have occurred throughout China, where the temporal and spatial differences in the cooling rates are remarkable. A significant feature is that cooling is greater in winter than in summer and at high latitudes than at low latitudes. This widespread cooling may have prompted the migration and fusion of high-latitude regions with low-latitude regions. During this process, species at high latitudes could be subjected to stronger natural selection than those at low latitudes.

The main driving forces in the process of species evolution include selection and genetic drift. The joint action of multiple forces may not only make the mutant increase its frequency faster, but also may be that genetic drift slows down the speed of selection, or even reduces the frequency of the mutant. The intensity of genetic drift depends on the effective population size. The larger the effective population size, the weaker the genetic drift; the smaller the effective population size, the higher the probability of genetic drift. According to our estimation of the effective population size, the effective population size of A. cerana in Jilin is lower than that of other populations. This means that A. cerana in Jilin may have experienced serious bottleneck events, and the evolution is more affected by genetic drift. This is consistent with the previously reported situation that human capture has greatly changed the population dynamics of A. cerana in the area around JL in recent decades [22, 25]. Therefore, in order to reduce the false positives caused by the influence of genetic drift and other factors on the SNP scan, we increase the sensitivity of the selection signal by calculating the FST between populations in the sliding window.

In this study, we mainly focus on the relationship between genetic variations and climatic conditions. Hence, we pay attention to the factors like latitude, longitude, and altitude related to climate conditions. It is worth noting that in the process of adaptive evolution, especially in the case of short generation interval of bees, the impact of artificial domestication cannot be ignored. Both the changes in genome structure characteristics caused by artificial selection and natural selection are called selection signals. For these selection signals, it is difficult to distinguish whether their source is artificial selection or natural selection. Therefore, when using the selective sweep analysis to analyze the evolution of species, it is necessary to select appropriate research objects according to the direction of focus. We selected A. cerana as the research object rather than Apis mellifera to study climatic adaptation. Because the domestication history of A. cerana is shorter than that of A. mellifera, and the degree of domestication is weaker than that of A. mellifera. In addition, during the sampling process, the samples we collected were A. cerana living in different geographical environments in the wild or semi-wild state. These samples are relatively less affected by artificial domestication. We believe that such samples are suitable for studying climatic adaptation during the adaptive evolution of A. cerana, and the analysis results based on these samples are representative.

Long-term natural selection often causes changes in the allele frequencies of selected loci and their linked loci. Upon analyzing these traces of natural selection in the genome, we can better understand some of the important events in evolution. We used selective sweep analysis to determine the relationship between climate type and the adaptive evolution of A. cerana. Positive selection will lead to an increase in the frequency of the dominant allele at the locus, which in turn will lead to a decrease in the genetic polymorphism of the selected locus. Therefore, for the studied colonies of A. cerana under the four climatic types, we set the SZ population as the reference, based on the calculated population polymorphism (SZ > WC > MK > JL). An abundance of polymorphisms means that the selection pressure received is less; hence using the SZ population with less selection pressure as the background facilitates the detection of traces of selection in the remaining populations. In addition, we use different selective signal methods (FST and π ratios) for selective signal detection. The adoption of an overlapping detection strategy can further avoid the occurrence of false positives through mutual verification between methods to a certain extent. It is undeniable that the results obtained after these pairwise comparisons are relatively rich, but not all signals are related to adaptation to the climate environment and body size; therefore, we pay more attention to the signals detected in all three selective sweep analyses. These signals will help gain insights into the climatic adaptation mechanisms of A. cerana. According to the results of GO and KEGG pathway analyses, we speculate that the adaptive evolution of A. cerana may be preferentially reflected in its response to the abundance of food resources and changes in temperature. A. cerana actively regulates its metabolism to cope with food shortages in harsh environments and maintain its body temperature. This also indicates that A. cerana energy utilization patterns may differ in different climate types.

Based on our findings, we speculate that the selection of RAPTOR could help A. cerana adapt to the complex climatic environment in China. A total of 33 candidate genes related to climatic adaptation were screened using selective elimination analysis. Upon annotation results combination with existing research reports, we identified RAPTOR as a candidate gene related to body size in climate adaptation. RAPTOR is the intersection of the AMPK and mTOR signaling pathways. The evidence to date suggests that the mTOR signaling pathway is involved in the development of A. mellifera [33,34,35,36]. Studies have shown that the mTOR signaling pathway plays a key role in the development of the queen bee [37,38,39,40,41], and the knockdown of the TOR encoding gene can affect the fate of the queen [37]. Under nutrient deprivation, AMPK, an important kinase that regulates energy homeostasis, directly phosphorylates RAPTOR to inhibit the mTORC1, thereby inhibiting cell growth [42]. The mTORC1 participates in glucose metabolism indirectly in response to insulin secreted by glucose entering the blood or is directly activated by amino acids [43, 44]. The mTOR pathway is also involved in regulating autophagy. When insufficient nutrient supply causes insufficient mTORC1 activity, mTORC1 initiates the inhibition of autophagy and activates lysosomes to degrade relatively uncritical proteins and organelles to provide material and energy to maintain the basic survival needs of cells [45, 46]. Bees absorb nutrients (glucose and amino acids) by feeding on plant pollen and nectar [23]. Native plants are affected by local climatic conditions, exhibiting shorter flowering periods and lower abundance in colder regions than in warmer regions [26, 47]. As a result, the degree of difficulty in finding food and absorbing nutrients for A. cerana in different regions is inconsistent. The selection of RAPTOR at the genome level could be in response to changes in climate and food sources. Early studies also reported the functions of the RAPTOR and TOR pathways in Drosophila [43, 46, 48, 49]. Decreased TOR activity inhibits ecdysone release and leads to prolonged development time and increased body size, whereas activating TOR can reverse stunting caused by nutrient restriction [50]. In addition, RAPTOR was upregulated in a study of cold resistance in Drosophila [51]. A recent study on thermal and oxygen flight sensitivity in Drosophila showed that downregulation of TOR activity produces smaller flies with smaller wing epidermal cells, and flies with small cells can maintain superior performance in metabolically demanding activities, such as flying under hypoxic conditions [52]. Whether A. cerana has adopted the same strategy as other insects to combat different climate types requires further investigation.

In summary, our findings suggest that climate type, geographic location, and food resources are all related to the adaptive radiation distribution and unique geographic phenotypes of A. cerana. It is highly probable that food and nutrition deficits caused by diverse climates are the predominant driving forces for unique geographic phenotypes. Further research is necessary to verify the molecular mechanism between RAPTOR and the body size.

Conclusions

Our study shows that the genetic structure and genetic differentiation of A. cerana are related to the distribution of climate and environment in China and are strongly affected by latitude. A key gene, RAPTOR, plays an important role in this process. The selection of RAPTOR at the genomic level helps A. cerana respond to nutritional problems caused by harsh environments by modulating TOR activity to alter body size, which partly explains the difference in body size among A. cerana distributed in China. These results help us understand the genetic basis of the climatic adaptation of A. cerana and provide a reference for the protection and utilization of germplasm resources and future genetic improvement.

Methods

Sample collection

A total of 100 honeybee samples were studied in this experiment, 70 of which were collected from Tibet, Hubei, Anhui, Jiangsu, Guangdong, and Hainan provinces in China. The collection sites for each sample were centered in Suzhou, Jiangsu, and were located at similar latitudes and different longitudes or similar longitudes and different latitudes. In addition, 30 samples from sites with similar longitude and different latitude as Suzhou were included for comprehensive data analysis [14, 18]. Ten colonies of 3–5 bees were collected from each sampling point. The collected bee samples included wild and semi-wild local bees. Samples were placed in 75% alcohol and then stored at -20 °C for future use.

Whole-genome sequencing, quality control, and clean reads mapping

Total genomic DNA was extracted from samples, and at least 3 µg genomic DNA was used to construct paired-end libraries with an insert size of 500 bp using a Paired-End DNA Sample Prep kit (Illumina Inc., San Diego, CA, USA). These libraries were sequenced using the HiSeq X10 NGS platform (Illumina Inc., San Diego, CA, USA). Raw reads were processed to obtain high-quality clean reads according to two stringent filtering standards: 1) removing reads containing > 50% of low-quality bases (Q < 20) or > 10% unidentified nucleotides (Ns); and 2) removing reads aligned to the barcode adapter. The Burrows-Wheeler Aligner (BWA) was used to align the clean reads against the reference genome (GCA_011100585.1) with the settings ‘mem 4 -k 32 -M’ [53]. Duplicates were marked using Picard 2.18.7 (http://broadinstitute.github.io/picard/).

Variant identification and annotation

Variant identification was performed using the Genome Analysis Toolkit (GATK) [54]. SNPs were filtered using GATK Variant Filtration with proper standards (-Window 4, -filter "QD < 2.0 || FS > 60.0 || MQ < 40.0", -G_filter "GQ < 20"). SNPs were filtered according to two stringent filtering standards: 1) missing ratio < 20%; and 2) minor allele frequency (MAF) > 5%. The ANNOVAR software [55] was used to annotate SNPs.

Principal component analysis (PCA), population structure analysis, and phylogenetic analysis

The admixture model-based software, Admixture V1.3.0 [56], was used to estimate the population structure. The tested K was set from 1 to 9, and the optimal K was determined based on the lowest cross-validation error. The population subdivision pattern was preliminarily classified using PCA in the GCTA software [57]. We constructed a phylogenetic tree using the neighbor-joining method with TreeBeST software [58]. The bootstrapped confidence interval was based on 1000 replicates.

Genetic diversity and FST statistics

Ho and He for each group or population were calculated using Plink [59]. Ne and PIC were calculated using the Perl script. The pairwise FST matrix was calculated using Genepop software [60].

Historical effective population sizes

SMC +  + 1.13.1 [61] was used to estimate the changes in Ne over the past one million years (Ma). The mutation rate was set to 5.27 × 10−9 per base pair per generation, following a divergence estimate of 7 Ma between A. mellifera and A. cerana [62]. The generation time was assumed to be one year. The polarization error was set to 0.5.

Morphometric analysis

A total of 130 worker bees (30 from MK, 50 from WC, and 50 from SZ) were used for morphometric measurements. The bees were dissected to make samples of each tissue, observed and photographed under a stereo microscope, and measured by a measurement system (M-Shot Image Analysis System V1.1.4). The measure of each sample in parallel was repeated 3 times and the mean value was taken.

Linkage disequilibrium (LD)

To evaluate the LD pattern, we estimated the squared allele frequency correlation (r2) using Haploview 4.2 [63]. LD decay graphs were plotted using the R script.

Selective sweep analysis

For the 40 samples from JL, MK, WC, and SZ, we performed selective sweep analysis, with SZ as the reference population. We estimated pairwise genetic differentiation (FST) [64] and differences in nucleotide diversity (π ln-ratio) [65] from the four different populations to identify key selective regions using the PopGenome software [66] with the sliding window approach. We set the window size to 100 kb and the step size to 10 kb. Selective sweep regions were selected according to the interception of two indices (FST and π ratios), with a threshold of the top 5% level. All related graphs were drawn using R scripts. Candidate genes within sweep regions were extracted for GO and KEGG enrichment analysis.

Availability of data and materials

Sequencing data of 70 samples in this study have been submitted to the NCBI Short Read Archive (SRA) under the BioProject accession number PRJNA876655. The morphometry data and other data are included in this paper.

References

  1. Hallmann CA, Sorg M, Jongejans E, Siepel H, Hofland N, Schwan H, Stenmans W, Muller A, Sumser H, Horren T, Goulson D, de Kroon H. More than 75 percent decline over 27 years in total flying insect biomass in protected areas. PLoS ONE. 2017;12(10):e0185809. https://doi.org/10.1371/journal.pone.0185809.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Soroye P, Newbold T, Kerr J. Climate change contributes to widespread declines among bumble bees across continents. Science. 2020;367(6478):685–8. https://doi.org/10.1126/science.aax8591.

    Article  CAS  PubMed  Google Scholar 

  3. Salcido DM, Forister ML, Lopez HG, Dyer LA. Loss of dominant caterpillar genera in a protected tropical forest. Sci Rep. 2020;10(1):422. https://doi.org/10.1038/s41598-019-57226-9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Harris JE, Rodenhouse NL, Holmes RT. Decline in beetle abundance and diversity in an intact temperate forest linked to climate warming. Biol Conserv. 2019;240:108219. https://doi.org/10.1016/j.biocon.2019.108219.

    Article  Google Scholar 

  5. van der Zee R, Pisa L, Andonov S, Brodschneider R, Charriere JD, Chlebo R, Coffey MF, Crailsheim K, Dahle B, Gajda A, Gray A, Drazic MM, Higes M, Kauko L, Kence A, Kence M, Kezic N, Kiprijanovska H, Kralj J, Kristiansen P, Hernandez RM, Mutinelli F, Nguyen BK, Otten C, Ozkirim A, Pernal SF, Peterson M, Ramsay G, Santrac V, Soroker V, Topolska G, Uzunov A, Vejsnaes F, Wei S, Wilkins S. Managed honeybee colony losses in Canada, China, Europe, Israel and Turkey, for the winters of 2008–9 and 2009–10. J Apic Res. 2012;51(1):91–114. https://doi.org/10.3896/Ibra.1.51.1.12.

    Article  Google Scholar 

  6. Scheffers BR, De Meester L, Bridge TCL, Hoffmann AA, Pandolfi JM, Corlett RT, Butchart SHM, Pearce-Kelly P, Kovacs KM, Dudgeon D, Pacifici M, Rondinini C, Foden WB, Martin TG, Mora C, Bickford D, Watson JEM. The broad footprint of climate change from genes to biomes to people. Science. 2016;354(6313):aaf7671. https://doi.org/10.1126/science.aaf7671.

    Article  CAS  PubMed  Google Scholar 

  7. Halsch CA, Shapiro AM, Fordyce JA, Nice CC, Thorne JH, Waetjen DP, Forister ML. Insects and recent climate change. Proc Natl Acad Sci U S A. 2021;118(2):e2002543117. https://doi.org/10.1073/pnas.2002543117.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Warren R, Price J, Graham E, Forstenhaeusler N, VanDerWal J. The projected effect on insects, vertebrates, and plants of limiting global warming to 1.5 degrees C rather than 2 degrees C. Science. 2018;360(6390):791–5. https://doi.org/10.1126/science.aar3646.

    Article  CAS  PubMed  Google Scholar 

  9. Roman-Palacios C, Wiens JJ. Recent responses to climate change reveal the drivers of species extinction and survival. Proc Natl Acad Sci U S A. 2020;117(8):4211–7. https://doi.org/10.1073/pnas.1913007117.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Urban MC. Accelerating extinction risk from climate change. Science. 2015;348(6234):571–3. https://doi.org/10.1126/science.aaa4984.

    Article  CAS  PubMed  Google Scholar 

  11. Garcia RA, Cabeza M, Rahbek C, Araujo MB. Multiple Dimensions of Climate Change and Their Implications for Biodiversity. Science. 2014;344(6183):1247579. https://doi.org/10.1126/science.1247579.

    Article  CAS  PubMed  Google Scholar 

  12. Boggs CL. The fingerprints of global climate change on insect populations. Curr Opin Insect Sci. 2016;17:69–73. https://doi.org/10.1016/j.cois.2016.07.004.

    Article  PubMed  Google Scholar 

  13. Hoffmann AA, Sgro CM. Climate change and evolutionary adaptation. Nature. 2011;470(7335):479–85. https://doi.org/10.1038/nature09670.

    Article  CAS  PubMed  Google Scholar 

  14. Chen C, Wang HH, Liu ZG, Chen X, Tang J, Meng FM, Shi W. Population Genomics Provide Insights into the Evolution and Adaptation of the Eastern HoneyBee (Apis cerana). Mol Biol Evol. 2018;35(9):2260–71. https://doi.org/10.1093/molbev/msy130.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Shi P, Zhou J, Song H, Wu Y, Lan L, Tang X, Ma Z, Vossbrinck CR, Vossbrinck B, Zhou Z, Xu J. Genomic analysis of Asian honeybee populations in China reveals evolutionary relationships and adaptation to abiotic stress. Ecol Evol. 2020;10(23):13427–38. https://doi.org/10.1002/ece3.6946.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Ilyasov RA, Park J, Takahashi J, Kwon HW. Phylogenetic Uniqueness of Honeybee from the Korean Peninsula Inferred from The Mitochondrial, Nuclear, and Morphological Data. J Apic Sci. 2018;62(2):189–214. https://doi.org/10.2478/jas-2018-0018.

    Article  CAS  Google Scholar 

  17. Ilyasov RA, Youn HG, Lee ML, Kim KW, Proshchalykin MY, Lelej AS, Takahashi J, Kwon HW. Phylogenetic Relationships of Russian Far-East Apis Cerana with Other North Asian Populations. J Apic Sci. 2019;63(2):289–314. https://doi.org/10.2478/Jas-2019-0024.

    Article  Google Scholar 

  18. Ji Y, Li X, Ji T, Tang J, Qiu L, Hu J, Dong J, Luo S, Liu S, Frandsen PB, Zhou X, Parey SH, Li L, Niu Q, Zhou X. 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 

  19. Wallberg A, Schoning C, Webster MT, Hasselmann M. Two extended haplotype blocks are associated with adaptation to high altitude habitats in East African honeybees. PLoS Genet. 2017;13(5):e1006792. https://doi.org/10.1371/journal.pgen.1006792.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Wallberg A, Pirk CW, Allsopp MH, Webster MT. Identification of Multiple Loci Associated with Social Parasitism in Honeybees. PLoS Genet. 2016;12(6):e1006097. https://doi.org/10.1371/journal.pgen.1006097.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Chen C, Liu ZG, Pan Q, Chen X, Wang HH, Guo HK, Liu SD, Lu HF, Tian SL, Li RQ, Shi W. Genomic Analyses Reveal Demographic History and Temperate Adaptation of the Newly Discovered HoneyBee 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 

  22. Park D, Jung JW, Choi BS, Jayakodi M, Lee J, Lim J, Yu Y, Choi YS, Lee ML, Park Y, Choi IY, Yang TJ, Edwards OR, Nah G, Kwon HW. Uncovering the novel characteristics of Asian honeybee, Apis cerana, by whole genome sequencing. BMC Genomics. 2015;16:1. https://doi.org/10.1186/1471-2164-16-1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Diao QY, Sun LX, Zheng HJ, Zeng ZJ, Wang SY, Xu SF, Zheng HQ, Chen YP, Shi YY, Wang YZ, Meng F, Sang QL, Cao LF, Liu F, Zhu YQ, Li WF, Li ZG, Dai CJ, Yang MJ, Chen SL, Chen RS, Zhang SW, Evans JD, Huang Q, Liu J, Hu FL, Su SK, Wu J. Genomic and transcriptomic analysis of the Asian honeybee Apis cerana provides novel insights into honeybee biology. Sci Rep. 2018;8:1–14. https://doi.org/10.1038/s41598-017-17338-6.

    Article  CAS  Google Scholar 

  24. Wang ZL, Zhu YQ, Yang Q, Yan WY, Zheng HJ, Zeng ZJ. A Chromosome-Scale Assembly of the Asian Honeybee Apis cerana Genome. Front Genet. 2020;11:279. https://doi.org/10.3389/fgene.2020.00279.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Liu NN, Liu HM, Ju Y, Li XA, Li Y, Wang TJ, He JM, Niu QS, Xing XM. Geometric morphology and population genomics provide insights into the adaptive evolution of Apis cerana in Changbai Mountain. BMC Genomics. 2022;23:64. https://doi.org/10.1186/s12864-022-08298-x.

    Article  CAS  Google Scholar 

  26. Lan L, Shi P, Song H, Tang X, Zhou J, Yang J, Yang M, Xu J. De Novo Genome Assembly of Chinese Plateau Honeybee Unravels Intraspecies Genetic Diversity in the Eastern Honeybee, Apis cerana. Insects. 2021;12(10):891. https://doi.org/10.3390/insects12100891.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Xu K, Niu Q, Zhao H, Du Y, Jiang Y. Transcriptomic analysis to uncover genes affecting cold resistance in the Chinese honeybee (Apis cerana cerana). PLoS ONE. 2017;12(6):e0179922. https://doi.org/10.1371/journal.pone.0179922.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Oldroyd BP, Fewell JH. Genetic diversity promotes homeostasis in insect colonies. Trends Ecol Evol. 2007;22(8):408–13. https://doi.org/10.1016/j.tree.2007.06.001.

    Article  PubMed  Google Scholar 

  29. Jones JC, Myerscough MR, Graham S, Oldroyd BP. Honeybee nest thermoregulation: diversity promotes stability. Science. 2004;305(5682):402–4. https://doi.org/10.1126/science.1096340.

    Article  CAS  PubMed  Google Scholar 

  30. Mattila HR, Seeley TD. Genetic diversity in honeybee colonies enhances productivity and fitness. Science. 2007;317(5836):362–4. https://doi.org/10.1126/science.1143046.

    Article  CAS  PubMed  Google Scholar 

  31. Batchelor CL, Margold M, Krapp M, Murton D, Dalton AS, Gibbard PL, Stokes CR, Murton JB, Manica A. The configuration of Northern Hemisphere ice sheets through the Quaternary. Nat Commun. 2019;10:3713. https://doi.org/10.1038/s41467-019-11601-2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Clark PU, Dyke AS, Shakun JD, Carlson AE, Clark J, Wohlfarth B, Mitrovica JX, Hostetler SW, McCabe AM. The Last Glacial Maximum. Science. 2009;325(5941):710–4. https://doi.org/10.1126/science.1172873.

    Article  CAS  PubMed  Google Scholar 

  33. Ronai I, Vergoz V, Oldroyd BP. The Mechanistic, Genetic, and Evolutionary Basis of Worker Sterility in the Social Hymenoptera. Adv Study Behav. 2016;48:251–317. https://doi.org/10.1016/bs.asb.2016.03.002.

    Article  Google Scholar 

  34. Cardoen D, Wenseleers T, Ernst UR, Danneels EL, Laget D, De Graaf DC, Schoofs L, Verleyen P. Genome-wide analysis of alternative reproductive phenotypes in honeybee workers. Mol Ecol. 2011;20(19):4070–84. https://doi.org/10.1111/j.1365-294X.2011.05254.x.

    Article  CAS  PubMed  Google Scholar 

  35. Chen H, Wu GA, Zhou H, Dai XY, Steeghs NWF, Dong XL, Zheng L, Zhai YF. Hormonal Regulation of Reproductive Diapause That Occurs in the Year-Round Mass Rearing of Bombus terrestris Queens. J Proteome Res. 2021;20(5):2240–50. https://doi.org/10.1021/acs.jproteome.0c00776.

    Article  CAS  PubMed  Google Scholar 

  36. Corby-Harris V, Snyder L, Meador C. Fat body lipolysis connects poor nutrition to hypopharyngeal gland degradation in Apis mellifera. J Insect Physiol. 2019;116:1–9. https://doi.org/10.1016/j.jinsphys.2019.04.001.

    Article  CAS  PubMed  Google Scholar 

  37. Patel A, Fondrk MK, Kaftanoglu O, Emore C, Hunt G, Frederick K, Amdam GV. The Making of a Queen: TOR Pathway Is a Key Player in Diphenic Caste Development. PLoS ONE. 2007;2(6):e509. https://doi.org/10.1371/journal.pone.0000509.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Mutti NS, Dolezal AG, Wolschin F, Mutti JS, Gill KS, Amdam GV. IRS and TOR nutrient-signaling pathways act via juvenile hormone to influence honeybee caste fate. J Exp Biol. 2011;214(23):3977–84. https://doi.org/10.1242/jeb.061499.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Wheeler DE, Buck NA, Evans JD. Expression of insulin/insulin-like signalling and TOR pathway genes in honeybee caste determination. Insect Mol Biol. 2014;23(1):113–21. https://doi.org/10.1111/imb.12065.

    Article  CAS  PubMed  Google Scholar 

  40. Corona M, Libbrecht R, Wheeler DE. Molecular mechanisms of phenotypic plasticity in social insects. Curr Opin Insect Sci. 2016;13:55–60. https://doi.org/10.1016/j.cois.2015.12.003.

    Article  PubMed  Google Scholar 

  41. Chen X, Hu Y, Zheng HQ, Cao LF, Niu DF, Yu DL, Sun YQ, Hu SN, Hu FL. Transcriptome comparison between honeybee queen- and worker-destined larvae. Insect Biochem Mol Biol. 2012;42(9):665–73. https://doi.org/10.1016/j.ibmb.2012.05.004.

    Article  CAS  PubMed  Google Scholar 

  42. Mihaylova MM, Shaw RJ. The AMPK signalling pathway coordinates cell growth, autophagy and metabolism. Nat Cell Biol. 2011;13(9):1016–23. https://doi.org/10.1038/ncb2329.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Geminard C, Rulifson EJ, Leopold P. Remote Control of Insulin Secretion by Fat Cells in Drosophila. Cell Metab. 2009;10(3):199–207. https://doi.org/10.1016/j.cmet.2009.08.002.

    Article  CAS  PubMed  Google Scholar 

  44. Duran A, Amanchy R, Linares JF, Joshi J, Abu-Baker S, Porollo A, Hansen M, Moscat J, Diaz-Meco MT. p62 is a key regulator of nutrient sensing in the mTORC1 pathway. Mol Cell. 2011;44(1):134–46. https://doi.org/10.1016/j.molcel.2011.06.038.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Wei Y, Lilly MA. The TORC1 inhibitors Nprl2 and Nprl3 mediate an adaptive response to amino-acid starvation in Drosophila. Cell Death Differ. 2014;21(9):1460–8. https://doi.org/10.1038/cdd.2014.63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Takats S, Varga A, Pircs K, Juhasz G. Loss of Drosophila Vps16A enhances autophagosome formation through reduced Tor activity. Autophagy. 2015;11(8):1209–15. https://doi.org/10.1080/15548627.2015.1059559.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Cleland EE, Chuine I, Menzel A, Mooney HA, Schwartz MD. Shifting plant phenology in response to global change. Trends Ecol Evol. 2007;22(7):357–65. https://doi.org/10.1016/j.tree.2007.04.003.

    Article  PubMed  Google Scholar 

  48. Guertin DA, Guntur KVP, Bell GW, Thoreen CC, Sabatini DM. Functional Genomics identifies TOR-regulated genes that control growth and division. Curr Biol. 2006;16(10):958–70. https://doi.org/10.1016/j.cub.2006.03.084.

    Article  CAS  PubMed  Google Scholar 

  49. Lee G, Chung J. Discrete functions of rictor and raptor in cell growth regulation in Drosophila. Biochem Bioph Res Co. 2007;357(4):1154–9. https://doi.org/10.1016/j.bbrc.2007.04.086.

    Article  CAS  Google Scholar 

  50. Layalle S, Arquier N, Leopold P. The TOR pathway couples nutrition and developmental timing in Drosophila. Dev Cell. 2008;15(4):568–77. https://doi.org/10.1016/j.devcel.2008.08.003.

    Article  CAS  PubMed  Google Scholar 

  51. Telonis-Scott M, Hallas R, McKechnie SW, Wee CW, Hoffmann AA. Selection for cold resistance alters gene transcript levels in Drosophila melanogaster. J Insect Physiol. 2009;55(6):549–55. https://doi.org/10.1016/j.jinsphys.2009.01.010.

    Article  CAS  PubMed  Google Scholar 

  52. Szlachcic E, Czarnoleski M. Thermal and Oxygen Flight Sensitivity in Ageing Drosophila melanogaster Flies: Links to Rapamycin-Induced Cell Size Changes. Biology. 2021;10(9):861. https://doi.org/10.3390/biology10090861.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. 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 

  54. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, DePristo MA. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20(9):1297–303. https://doi.org/10.1101/gr.107524.110.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. 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 

  56. 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 

  57. 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 

  58. Vilella AJ, Severin J, Ureta-Vidal A, Heng L, Durbin R, Birney E. EnsemblCompara GeneTrees: Complete, duplication-aware phylogenetic trees in vertebrates. Genome Res. 2009;19(2):327–35. https://doi.org/10.1101/gr.073585.107.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, Sham PC. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81(3):559–75. https://doi.org/10.1086/519795.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Diniz JAF, Soares TN, Lima JS, Dobrovolski R, Landeiro VL, Telles MPD, Rangel TF, Bini LM. Mantel test in population genetics. Genet Mol Biol. 2013;36(4):475–85. https://doi.org/10.1590/S1415-47572013000400002.

    Article  Google Scholar 

  61. 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 

  62. Wallberg A, Han F, Wellhagen G, Dahle B, Kawata M, Haddad N, Simoes ZLP, Allsopp MH, Kandemir I, De la Rua P, Pirk CW, Webster MT. A worldwide survey of genome sequence variation provides insight into the evolutionary history of the honeybee Apis mellifera. Nat Genet. 2014;46(10):1081–8. https://doi.org/10.1038/ng.3077.

    Article  CAS  PubMed  Google Scholar 

  63. Barrett JC, Fry B, Maller J, Daly MJ. Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics. 2005;21(2):263–5. https://doi.org/10.1093/bioinformatics/bth457.

    Article  CAS  PubMed  Google Scholar 

  64. Hudson RR, Slatkin M, Maddison WP. Estimation of levels of gene flow from DNA sequence data. Genetics. 1992;132(2):583–9. https://doi.org/10.1093/genetics/132.2.583.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Lin T, Zhu GT, Zhang JH, Xu XY, Yu QH, Zheng Z, Zhang ZH, Lun YY, Li S, Wang XX, Huang ZJ, Li JM, Zhang CZ, Wang TT, Zhang YY, Wang AX, Zhang YC, Lin K, Li CY, Xiong GS, Xue YB, Mazzucato A, Causse M, Fei ZJ, Giovannoni JJ, Chetelat RT, Zamir D, Stadler T, Li JF, Ye ZB, Du YC, Huang SW. Genomic analyses provide insights into the history of tomato breeding. Nat Genet. 2014;46(11):1220–6. https://doi.org/10.1038/ng.3117.

    Article  CAS  PubMed  Google Scholar 

  66. Pfeifer B, Wittelsburger U, Ramos-Onsins SE, Lercher MJ. PopGenome: An Efficient Swiss Army Knife for Population Genomic Analyses in R. Mol Biol Evol. 2014;31(7):1929–36. https://doi.org/10.1093/molbev/msu136.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We are deeply grateful to all donors who participated in this program. We would also like to thank Editage (www.editage.cn) for their language editing services.

Funding

This research was funded by the National Natural Science Foundation of China (32272935, 31902220, ZL), the Natural Science Foundation of Jilin Province (20180101022JC, ZW, and QN), the Science and Technology Development Project of Jilin Province (20200402042NC, ZW, and QN), the Earmarked Fund for Modern Agro-industry Technology Research System (CARS-44, TJ, QN, FG, and ZW), and the Lvyangjinfeng Program of Yangzhou (YZLYJFJH2021YXBS155, ZL).

Author information

Authors and Affiliations

  1. College of Animal Science and Technology, Yangzhou University, Yangzhou, 225009, China

    Yi Zhang, Minqi Cai, Kang Wang, Zheguang Lin & Ting Ji

  2. Sericultural Research Institute, Anhui Academy of Agricultural Science, Hefei, 230061, China

    Hao Xu

  3. Apiculture Science Institute of Jilin Province, Jilin, 132108, China

    Zhi Wang & Qingsheng Niu

  4. Jinzhong Agriculture and Rural Affairs Bureau, Jinzhong, 030601, China

    Haoliang Jie

  5. Mudanjiang Branch of Heilongjiang Academy of Agricultural Sciences, Mudanjiang, 157043, China

    Fuchao Gao

  6. College of Animal Science (College of Bee Science), Fujian Agriculture and Forestry University, Fuzhou, 350002, China

    Dafu Chen & Rui Guo

Authors
  1. Yi Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Hao Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Zhi Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Haoliang Jie
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Fuchao Gao
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Minqi Cai
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Kang Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Dafu Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Rui Guo
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Zheguang Lin
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Qingsheng Niu
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Ting Ji
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conception and design of the study: TJ, QN, FG, ZL, and YZ; Collection of the samples: TJ, QN, FG, DC, ZW, and HX; Acquisition of the morphometry data: YZ, HX, HJ, KW, and MC; Analysis of the sequencing data: YZ, ZW, RG, and HX; Drafting of the paper: YZ, ZL, and TJ; All authors contributed to and approved the manuscript.

Corresponding authors

Correspondence to Qingsheng Niu or Ting Ji.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

All authors declare 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

Additional file 1:

Figure S1. The cross-validation error rate of genetic structure analysis of 100 samples. Figure S2. The diagram of 10 morphological indicators of A.cerana. (A) The right forewing length (FL) and width (FB). (B) The sixth sternum length (L6) and width (T6). (C) The third sternumlength (S3). (D) The third tergum length (T3). (E) The fourth tergum length (T4). (F) The femur length (Fe), the tibia length (Ti), the basitarsus length (ML), and the basitarsus width (MT). Figure S3. Analysis of linkage disequilibrium. Figure S4. GO classification of candidate genes. Figure S5. The top 20 enriched KEGG pathways.

Additional file 2:

Table S1. Sequencing data production and quality statistics of 100 samples. Table S2. Summary of sequencing data depth, coverage, and mapping rate of 10 sampling sites. Table S3. 527 selected genes of A. cerana in Jilin. Table S4. 565 selected genes of A. cerana in Mangkang. Table S5. 311 selected genes of A. cerana in Wenchang. Table S6. 33 common selected genes of the three geographical populations. Table S7. GO enrichment of the 33 common selected genes of A. cerana (P < 0.05). Table S8. KEGG enrichment of the 33 common selected genes of A. cerana (P <0.05).

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

Zhang, Y., Xu, H., Wang, Z. et al. A key gene for the climatic adaptation of Apis cerana populations in China according to selective sweep analysis. BMC Genomics 24, 100 (2023). https://doi.org/10.1186/s12864-023-09167-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-023-09167-x

Keywords

  • Apis cerana
  • Adaptive radiation distribution
  • Climate change
  • Population genomics
  • Selective sweep analysis

BMC Genomics

ISSN: 1471-2164

Contact us