Morphological and Genetic Variation of the Chinese Honey Bee (Apis cerana cerana Fabricius, 1793) in Wanyuan, Southwest China

TL;DR

This study explores how Chinese honey bees in Wanyuan adapt to different altitudes, showing changes in body traits and maintaining genetic diversity.

Contribution

The study reveals altitudinal morphological and genetic adaptations in a local ecotype of Chinese honey bees.

Findings

Mid-altitude bees showed larger body size, while low-altitude bees had more hindwing hooks.

Genetic diversity was high with 25 haplotypes and strong gene flow between populations.

One site at 900 m showed slight genetic distinctiveness despite overall connectivity.

Abstract

This study investigated how the Chinese honey bee (Apis cerana cerana Fabricius, 1793) adapts to different altitudes in Wanyuan, Sichuan. We collected 656 bees from 15 sites (565–1611 m) and measured 13 body traits and analyzed mitochondrial DNA. We found clear changes in bee shape with altitude: mid-elevation bees were larger in several traits, while low-elevation bees had more hindwing hooks, likely helping them fly in dense vegetation. Key identification traits also varied with elevation. Genetically, the population showed high diversity with 25 haplotypes, and most groups were well connected through gene flow, though one site at 900 m was slightly distinct. Overall, the Wanyuan honey bee displays noticeable physical adjustments to altitude while maintaining strong genetic connectivity across the landscape. These insights can help guide local conservation and sustainable use of this valuable honey bee resource.

Chinese honey bees (Apis cerana cerana Fabricius, 1793) are crucial native pollinators in China, with substantial ecological and economic value. Their morphological traits may vary along altitudinal gradients, particularly in hilly regions such as Wanyuan City, Sichuan Province, which provides typical suitable habitat for a locally thriving ecotype known as the Wanyuan honey bee. To elucidate its adaptive variation across environmental gradients, this study investigated the morphological and genetic diversity of this ecotype along an altitudinal transect in Wanyuan. A total of 656 worker individuals from 15 sampling sites (565–1611 m) were analyzed for 13 morphological traits and mitochondrial DNA (tRNAleu–COII fragment) sequences. Results revealed significant altitudinal clines in morphology: Honey bees from mid-altitude sites exhibited larger body size for several traits, while low-altitude bees possessed a significantly higher number of hindwing hamuli. Key taxonomic indices like the cubital index and proboscis length also varied significantly with altitude. Genetic analysis identified 25 haplotypes with high haplotype diversity and nucleotide diversity, indicating substantial genetic variation. Population differentiation was generally low, with one site (Yinbazhai, 900 m) showing relatively higher distinctiveness. The detected high gene flow suggests frequent genetic exchange among most populations. These findings demonstrate that the Wanyuan honey bee exhibits clear altitudinal adaptation in morphology while maintaining high genetic diversity and connectivity. This study provides a crucial scientific basis for the conservation and sustainable management of this genetic resource by highlighting the importance of its population-specific adaptations and genetic structure.

Figures

Full text

1. Introduction

Apis cerana cerana Fabricius, 1793 is an ecologically and economically important indigenous bee subspecies in China. Compared with widely introduced Western honey bee subspecies (e.g., Apis mellifera ligustica Spinola, 1806), the Chinese honey bee has undergone long-term natural selection, thus evolving strong adaptability to local climates, geographical conditions, and plant communities, facilitating its effective integration into regional ecological networks [1,2]. As a functionally critical pollinator, the Chinese honey bee holds irreplaceable ecological and economic significance [3]: it is widely distributed across China, with hilly regions serving as its primary population strongholds [4]. A notable case is Wanyuan City, which locates in the northeastern part of Sichuan and the hinterland of the Daba Mountains, Southwest China. As the northern margin of the subtropical zone, the city features a mild climate, abundant rainfall and distinct seasons, which are key characteristics of typical optimal habitats for the Chinese honey bee [5]. Additionally, Wanyuan harbors nectar-rich taxa such as Rosaceae and Fabaceae, providing abundant pollen and nectar resources [6]. The local Wanyuan Chinese honey bee (A. c. cerana), a distinct ecotype, exhibits advantageous traits including robust colony viability, low swarming tendency, and high foraging efficiency. These traits not only enhance its local ecological adaptability and apicultural utility [5] but also support population persistence, thereby providing critical support for regional pollination services and the sustainable development of the apicultural industry. Environmental gradient factors, particularly altitude, are key drivers of morphological differentiation within honey bee populations [7,8]. However, systematic investigations into the morphological and genetic traits of the Wanyuan honey bee across different altitudinal gradients remain limited.

Over the past two decades, with advances in Apis (honey bee) research methodologies, scholars have gradually developed an integrated research strategy based on the combination of “improved morphological identification methods (e.g., geometric morphometrics) and molecular techniques (e.g., mitochondrial DNA fragment analysis)” (e.g., [9,10,11]). Morphological traits of the honey bee serve as key indicators for elucidating its physiological characteristics, behavioral adaptations, and subspecies delineation [12,13]; they further provide critical morphological support for interpreting the evolutionary history of honey bee populations and implementing genetic resource conservation [14,15].

Meanwhile, mitochondrial DNA (mtDNA) possesses biologically distinct characteristics, including a small genome size, simple structural organization, high copy number per cell [16,17], and strict maternal inheritance [18,19]. Its evolutionary trajectory is primarily driven by base substitutions, with insertions/deletions (indels) being rare, and it exhibits a relatively rapid evolutionary rate. These attributes collectively position it as an ideal molecular marker for population-level phylogenetic studies of the honey bee [20,21], one widely applied in phylogenetic reconstruction, inference of Chinese honey bee evolutionary histories, and characterization of Chinese honey bee population genetic structure [22,23]. Notably, research leveraging the aforementioned approaches has indicated that even within a single Chinese honey bee ecotype, significant morphological variation and genetic differentiation may arise among geographically discrete populations [24,25,26], a finding that lays a robust theoretical and methodological foundation for in-depth studies of region-specific Chinese honey bee populations.

Building on this, the present study focuses on the Wanyuan honey bee, conducting quantitative measurements and subsequent statistical analyses of morphological traits relevant to Chinese honey bee genetic resource identification, including proboscis length, forewing length, and other key traits [27], and performing molecular analyses of its mitochondrial tRNAleu–COII fragment region. Specifically, this study aims to elucidate the population-specific morphological phenotypes of the Wanyuan honey bee across different altitudinal gradients. In doing so, it seeks to provide a theoretical basis and empirical data support for the scientific conservation, rational development, and sustainable utilization of this Wanyuan honey bee genetic resource, thereby enriching the body of research on the overall genetic diversity of Chinese honey bees.

2. Materials and Methods

2.1. Sample Collection

A total of 656 individuals of the Wanyuan honey bee were collected from 15 sampling sites across Wanyuan City, Sichuan Province (Figure 1, Table S1). These sites included Heishan, Taiping Town, Jiuyuan Town, Fengtong Town, Guandu Town, Baisha Town, Zhuyu Town, and Zixi Town. At each site, 42–45 colonies were randomly sampled with 15–20 worker bees per colony. These honey bee individuals were preserved in 75% ethanol and stored at −20 °C until further use.

2.2. Experimental Treatment and Sample Processing

2.2.1. Morphological Measurements

Following the standard method [12,28], we compared 13 morphological characteristics of the Wanyuan honey bee workers to evaluate this subspecies’ morphological parameters. These characteristics specifically include proboscis length, length of forewing, width of forewing, cubital index, length of femur, length of tibia, length of metabasitarsus, width of metabasitarsus, length of tergum 3, length of tergum 4, length of sternum 6, width of sternum 6, and number of hindwing hamuli. The worker bees were dissected and measured under a stereomicroscope (Carl Zeiss, Stemi 305, Oberkochen, Germany).

2.2.2. DNA Extraction

We selected 3 individuals from each of the 15 sampling sites, resulting in a total of 45 individuals for DNA extraction. Head and thorax tissues were sampled from these selected individuals, and DNA extraction was performed using a DNA extraction kit (Tiangen Biotech, Beijing, China) in strict compliance with the manufacturer’s instructions.

2.2.3. PCR Amplification

Amplification of the mitochondrial tRNAleu–COII region of the Wanyuan honey bee was performed using specific primers (E2: 5′-GGCAGAATAAGTGCATTG-3′, H2: 5′-CAATATCATTGATGACC-3′; [29]). The tested nucleotide sequences of mtDNA correspond to the 5′ end of the tRNAleu gene (71 bp) and the 3′ end of the COII gene (239 bp) fragments, i.e., a total of 310 bp [29]. A 20 μL PCR reaction mixture was used, which contained approximately 50 ng of genomic DNA, 2.0 μL of 10× buffer, 1 μL of 200 μmol/L dNTPs, 1 μL of 10 pmol forward primer, 1 μL of 10 pmol reverse primer, 0.25 μL of Taq DNA polymerase, and double-distilled water (ddH_2_O) to bring the total volume to 20 μL. The PCR thermal profile was set as follows: initial denaturation at 96 °C for 3 min, followed by 35 cycles of denaturation at 94 °C for 45 s, annealing at 50 °C for 45 s, and extension at 70 °C for 120 s, and a final extension at 72 °C for 5 min.

2.3. Statistical Analyses

PCR products were subjected to 1.5% agarose gel electrophoresis for quality detection, after which they were sent to BGI Biotech Co., Ltd. (Beijing, China) for sequencing. Obtained sequences were aligned using Clustal X software v2.1, with subsequent manual checking of corresponding sequencing electropherograms to identify sequencing errors and trim primers. Sequences from individuals with interindividual nucleotide variations (e.g., substitutions) in the mitochondrial tRNAleu–COII fragment were re-sequenced to validate sequence accuracy. DNASP 5.0 software [30] was employed to calculate the number of nucleotide differences, nucleotide diversity, and haplotype diversity for each sampling site. Meanwhile, Mega 6.0 software was used to determine genetic differentiation coefficients between sampling sites, and these coefficients were further used to compute gene flow parameters.

SPSS software v30.0 was employed to assess morphological differences among the 15 sampling sites. One-way Analysis of Variance (One-way ANOVA) was conducted to test the significance of trait variation across the 15 populations, and morphological characteristics of each population were described using mean values and standard deviations (mean ± SD). Post hoc multiple comparisons were performed using the Least Significant Difference (LSD) method for data with homogeneous variances, and Tamhane’s T2 test for data with heterogeneous variances, with different letters assigned to values to denote significant differences at p < 0.05.

GraphPad Prism (Version 10; GraphPad Software, LLC, La Jolla, CA, USA) software was used to analyze morphological data aggregated from the 15 sampling sites into three altitude gradients (low: <900 m, medium: 900–1165 m, high: >1165 m). Prior to overall significance testing via One-way ANOVA, a variance homogeneity test was conducted: standard One-way ANOVA was applied for data with homogeneous variances (p > 0.05), while the Brown-Forsythe test and Welch’s ANOVA were used for corrected analysis of data with heterogeneous variances (p ≤ 0.05). If the overall ANOVA result indicated p < 0.05 (significant differences between at least two altitude groups), further pairwise comparisons were carried out; the analysis was concluded directly if p ≥ 0.05 (no significant differences across all altitude groups).

3. Results

3.1. Morphological Analysis

By quantifying 13 morphological traits of 656 worker bees from 15 sampling sites, we identified significant inter-regional differences in these morphological traits (Table S2). To simplify complex data and intuitively visualize group differences, we constructed a parallel coordinate plot to illustrate the morphological differences across distinct sampling sites (Figure 2).

Among the 15 sites, we found that the Wanyuan honey bee from Sampling Site 11 (altitude: 1167 m, Table S1) exhibited the largest values in 5 out of the 13 morphological indices compared to those from other sites. These indices include proboscis length (mean: 5.23 mm), width of forewing (mean: 2.98 mm), cubital index (mean: 4.75), length of tergum 4 (mean: 1.93 mm), length of tergum 3 (mean: 1.97 mm). Four additional indices showed the second-largest values, with no significant difference from the largest values: length of hind leg femur (mean: 2.46 mm), length of hind leg metabasitarsus (mean: 2.00 mm), length of sternum 6 (mean: 2.37 mm), and width of sternum 6 (mean: 2.83 mm). However, the number of hindwing hamuli of the Wanyuan honey bee at this sampling site was the smallest among all sites. In contrast, the Wanyuan honey bee from Sampling Site 15 (altitude: 1611 m, Table S1) had the smallest values in 5 indices, including width of forewing (mean: 2.84 mm), length of hind leg tibia (mean: 2.61 mm), width of hind leg metabasitarsus (mean: 1.02 mm), length of tergum 3 (mean: 1.84 mm), and width of sternum 6 (mean: 2.63 mm).

To further investigate the association between the morphological structure of the Wanyuan honey bee and altitude, we stratified the sampling area into three altitude groups (low, medium, and high) and conducted a quantitative comparative analysis of the Wanyuan honey bee’s morphological indices. The results showed that morphological traits exhibited differential variation patterns across the three altitude groups (Figure 3). Proboscis length and cubital index, key taxonomic indices for honey bees, with differences commonly used to differentiate distinct taxa, exhibited no continuous variation along the altitude gradient (Figure 3). Specifically, proboscis length exhibited a significant difference solely between the low- and high-altitude groups, with no significant differences between the low-medium and medium-high groups. This variation pattern may be associated with differences in the corolla tube lengths of nectar-pollen plants across altitudes, indirectly reflecting the potential for co-evolution between the Wanyuan honey bee and plants. For the cubital index, a highly significant difference was detected only between the low- and high-altitude groups, with no significant differences between adjacent groups.

Body size-related indices, including length of forewing, width of forewing, length of hind leg femur, length of hind leg tibia, length of tergum 3, length of tergum 4, length of sternum 6, and width of sternum 6, showed significant differences among groups, but the variation patterns were index-dependent. For example, indices such as length of forewing and width of forewing exhibited highly significant differences between the low-altitude group and both the medium- and high-altitude groups; length of tergum 3 showed a significant difference only between the medium- and high-altitude groups; and length of tergum 4, length of sternum 6, and width of sternum 6 exhibited significant differences solely between the low- and high-altitude groups. Notably, width of hind leg metabasitarsus, directly related to motor function, did not follow the aforementioned discrete variation pattern; instead, it showed gradually significant differences with increasing altitude, reflecting continuous adaptive adjustments to different altitude environments. In addition, analysis of function-specific traits indicated that the number of hindwing hamuli (linked to flight adhesion function) was significantly higher in the low-altitude group than in the medium- and high-altitude groups. This phenomenon may be related to the more complex vegetation structure and greater flight adhesion requirements in low-altitude areas. In summary, most morphological traits of the Wanyuan honey bee exhibited significant differences across different altitude groups, and the variation patterns were distinctly differentiated based on trait function. These patterns not only demonstrate the relative stability of taxonomic indices but also reflect the specific responses of body size-related traits and targeted adaptations of functional traits, collectively presenting a multi-faceted adaptive strategy of the Wanyuan honey bee to the altitude environment.

3.2. Haplotype Diversity

A total of 25 haplotypes were identified, designated as Hap-1–Hap-25. Among the 45 samples of the Wanyuan Chinese honey bee, Hap-3 was detected in 10 samples, accounting for 22% of the total samples; Hap-12 and Hap-14 were identified in 3 and 4 samples, respectively, representing 6.7% and 8.9% of the total samples. Hap-1, Hap-5, Hap-7, Hap-9, and Hap-15 were each detected in 2 samples, while the remaining 18 haplotypes exhibited relatively low frequencies. Notably, each individual Wanyuan Chinese honey bee across all groups carried only one haplotype.

The Wanyuan honey bee exhibited an average haplotype diversity of 0.9112 and an average nucleotide diversity of 0.0328. For nucleotide diversity specifically, values were highest in individuals from Sampling Site 4 (altitude: 859 m, Table S1) and lowest in those from Sampling Site 3 (altitude: 750 m, Table S1) (Table 1).

3.3. Population Differentiation and Genetic Structure

The overall genetic differentiation coefficient (FST) among all sampling points of the Wanyuan Chinese honey bee ranged from 0.005 to 0.071. Specifically, pairwise FST values between Sampling Site 5 (Yinbazhai Village, Baisha Town; altitude: 900 m) and other sampling sites were comparatively elevated, spanning 0.062 to 0.071. In contrast, the lowest FST values (0.005–0.008) were observed between specific site pairs (e.g., Sampling Site 13 vs. 14, 13 vs. 2, 14 vs. 2), with the majority of pairwise FST values falling within the range of 0.015 to 0.07. The gene flow parameter (Nm) derived from the genetic differentiation coefficients between pairs of sampling points ranged from 3.251 to 50.513 (Table 2).

It should be noted that Nm values are calculated based on the island model (FST = 1/(4Nm + 1)), which assumes evolutionary equilibrium and no human-mediated dispersal. However, these assumptions may not hold for semi-managed Wanyuan honey bee (A. cerana) populations [31]. So Nm should be interpreted as an approximate indicator of genetic connectivity rather than direct evidence of natural dispersal.

4. Discussion

Through in-depth analysis of the morphological traits and genetic diversity of the Wanyuan Chinese honey bee in Wanyuan City, Sichuan Province, China, this study has characterized the morphological features and genetic differentiation of this subspecies across altitudinal gradients. Results indicate that worker bees from different altitudes exhibit significant morphological variation, a pattern potentially linked to ecological adaptation to distinct altitudinal environments or inter-population gene flow dynamics (Figure 3); this differentiation is likely facilitated by frequent genetic exchange between populations (Table 2). While it is generally postulated that honey bee body size increases with rising altitude [32], our experimental results deviate from this conventional pattern. Specifically, worker bees from the mid-high altitude sampling site (Site 11) exhibit the largest body size (Table S2), whereas those from low-altitude regions have a higher number of hindwing hamuli (Figure 3). This observed discrepancy could be attributed either to the unique ecological and geographical characteristics of the Wanyuan region or to the frequent gene flow among Wanyuan Chinese honey bee populations, an inference supported by the small FST values and large Nm values (Table 2) that reflect high genetic connectivity between populations. Such a morphological variation pattern is also consistent with ecological adaptation to distinct altitudinal environments [7,33]. Furthermore, the cubital index and proboscis length, key metrics for honey bee population differentiation and habitat adaptation [13,34], display significant differences across altitudinal populations, thereby reflecting the function-specific adaptive responses of Wanyuan Chinese honey bee morphological traits to altitudinal gradients.

After clarifying the morphological and genetic differentiation characteristics of the Wanyuan Chinese honey bee, further analysis reveals that the high genetic connectivity exhibited by the population is not driven by a single factor. Notably, the high genetic connectivity inferred from low FST and high Nm values (which underpin the frequent genetic exchange mentioned above) may not solely result from natural dispersal of honey bees [35]. As a locally important apicultural resource, the Wanyuan Chinese honey bee is semi-managed by local beekeepers, who may conduct colony exchange, queen introduction, or transhumance between adjacent regions [28,36]. These anthropogenic activities could promote gene flow among populations, supplementing natural dispersal. Therefore, the observed strong genetic connectivity is likely a combined result of natural dispersal and human-mediated colony movement.

It is further important to point out that the high genetic connectivity at the overall level does not obscure the uneven nature of gene flow (Table 2). Specifically, the Nm value between Site 13 and Site 14 reaches an extremely high 50.513, indicating intense genetic exchange between these two populations. In contrast, most other population pairs have significantly lower Nm values (ranging from 3.251 to around 34), reflecting relatively limited genetic connectivity. More importantly, gene flow is generally restricted among the three predefined altitude groups (low: <900 m, medium: 900–1165 m, high: >1165 m) in this study. This pattern is consistent with the significant morphological differences observed among altitude groups: gene flow typically reduces phenotypic divergence by homogenizing genetic variation, but the limited genetic exchange between altitude groups allows strong diversifying selection driven by altitudinal environmental gradients to maintain or even amplify morphological differences, such as more hindwing hamuli in low-altitude populations and larger body size in mid-altitude populations.

Mitochondrial DNA tRNAleu-COII fragment analysis identified 25 distinct haplotypes in the Wanyuan Chinese honey bee, with an average haplotype diversity of 0.9112, evidence of substantial genetic diversity within this population. This high genetic diversity not only reflects the subspecies’ long-term adaptation to heterogeneous natural environments but also acts as a critical buffer for population health and long-term persistence [37,38]. Further analyses detected moderate genetic differentiation among geographically distinct Wanyuan Chinese honey bee populations. Of particular note is that the Yinbazhai sampling site (900 m a.s.l.) displays relatively greater genetic differentiation from other populations (FST: 0.062–0.071), in contrast to FST values of ≤0.04 observed between all other pairs of sampling sites. This pattern may be attributed to the unique local habitat characteristics of Yinbazhai: the site is located in a semi-isolated hilly area with fragmented vegetation cover, which may have slightly restricted gene flow with surrounding populations. Additionally, subtle differences in local apicultural practices (e.g., relatively limited exchange of breeding colonies among local beekeepers) could have further reinforced this weak genetic differentiation, leading to its distinct genetic signature compared to other populations. These findings highlight the need to account for geographical and genetic structural variation when developing conservation and management strategies for Wanyuan Chinese honey bee resources, thereby justifying the implementation of targeted protective measures [39].

Based on the core findings of population genetic structure, gene flow characteristics, and genetic diversity mentioned above, the conservation and sustainable utilization of the Wanyuan Chinese honey bee hold significant ecological and economic importance. While moderate genetic differentiation was detected, the relatively high gene flow values indicate that these populations maintain close genetic connectivity, a pattern that is crucial for sustaining the genetic stability and adaptive potential of the overall population [40,41]. Notably, as a key ecological resource in a national-level ecological demonstration city, the conservation and utilization of Wanyuan Chinese honey bees are pivotal for maintaining local ecological balance and fostering regional economic development. Thus, enhanced monitoring and targeted research efforts are warranted to safeguard the security and stability of this subspecies’ germplasm resources.

Although this study has initially clarified the adaptive characteristics and genetic patterns of the Wanyuan Chinese honey bee, several research limitations require objective clarification. As mitochondrial DNA is maternally inherited, this study cannot fully capture the genetic dynamics of the nuclear genome, especially regions linked to morphological traits, meaning the specific impact of heterogeneous gene flow on the whole genome and particularly on regions regulating key morphological traits such as proboscis length, cubital index and hindwing hamuli number remains unclear. Given the aforementioned genetic limitations of mitochondrial DNA and the constraints on available data types, we did not test the correlation between genetic and morphological distances, nor verify the isolating effect of altitude or geographical distance on population differentiation. Additionally, the lack of quantitative analysis on the consistency between phenotypic variation and genetic structure has resulted in separate rather than integrated discussions of morphology and genetics, which limits the depth of evolutionary interpretation in this study. To address these shortcomings, future research could leverage whole-genome sequencing to explore the genetic basis of morphological adaptation and clarify the associations and intrinsic mechanisms among genetic divergence, phenotypic variation, and environmental factors, as well as the impacts of anthropogenic activities on gene flow dynamics [42,43,44], thereby advancing our understanding of the adaptive evolutionary rules of the Wanyuan honey bee and providing scientific support for the targeted conservation and sustainable utilization of this species.