Morphometrical analyses revealed high diversity of the eastern honey bee (Apis cerana) in mountains and islands in China

Abstract

The eastern honey bee (Apis cerana) is of great importance for ecosystems and agriculture in Asia, whose diversity has yet to be explored. Previous studies have found a high diversity of A. cerana in mountainous areas or islands. To better understand the morphological diversity of A. cerana, we focused on the eastern Qinghai-Tibet Plateau (QTP), a region with high mountains, Weizhou Island, Yongxing Island and Hainan Island, much of which lacks systematic investigation of honey bee diversity. Five sampling sites in the plain were also selected as references. A total of 140 colonies of A. cerana were collected from 9 localities in QTP, 3 localities in islands, and 5 localities in plain. We used standard morphometry and geometric morphometrics to assess the level of diversity, in which principal component analysis and discriminant analysis were conducted. The results revealed the A. cerana populations from QTP and islands were differentiated from the plain populations. Moreover, two populations from QTP and one from Yongxing Island showed significant divergence. Our study provided insights into the A. cerana diversity and provides information for the effective conservation of A. cerana genetic resources in China.

Keywords:

• Apis cerana
• diversity
• Qinghai-Tibet Plateau
• islands
• standard morphometry
• geometric morphometrics

Introduction

Pollinators are an integral part of ecosystems because of their essential role in maintaining ecosystem function, wild plant reproduction and crop production (Allen-Wardell et al., 1998; Kremen et al., 2007). It is estimated that 87.5% of angiosperm species need pollinators (Ollerton et al., 2011), while pollination enhances the economic value of world agriculture by an additional US$235–577 billion in 2015 prices annually (Potts et al., 2016). Honey bees are the preferred pollinators for most plants due to their biology and characteristics (Hung et al., 2018; National Research Council, 2007).

The eastern honey bee, Apis cerana, whose natural distribution covers a large variety of geographic regions, ranging from hot to cold temperate climates, from lowlands to high-altitude mountainous areas and from humid coast to arid or semi-arid midlands (Hepburn & Radloff, 2011), is one of the most important honey bee species in Asia. A. cerana has high genetic and phenotypic variations occurring in different spatial scales, which led to the division into several subspecies or ecotypes (Radloff et al., 2010; Ruttner, 1988). Previous research found native honey bees have persistent colony-level behavioral adaptations to the local climate. For example, A. cerana populations in Qinghai-Tibet Plateau (QTP) and tropical Hainan Island show strong nectar collection and strong pollen collection capacity, respectively (China National Commission of Animal Genetic Resources, 2011), reflecting their adaptation to different environments. And in apiaries, colonies with local queens survived obviously longer than colonies with introduced queens (Büchler et al., 2014). Overall, native bees outperformed introduced bees in terms of both ecological and economic value. Therefore, it is particularly important to protect the diversity of A. cerana, whose population has declined in recent decades (Theisen-Jones & Bienefeld, 2016). To protect this, it calls for a better understanding and effective protection of its biodiversity. For this, the crucial prerequisites we need are the classification status and detailed differentiation of A. cerana. The first comprehensive multivariate statistical analysis of this species was published by Ruttner (1988) based on morphological data, which stated that A. cerana has four subspecies: A. c. cerana, A. c. indica, A. c. japonica, A. c. himalaya. The other two comprehensive studies of A. cerana across its full geographical range have been conducted by Smith et al. (2000) using mitochondrial DNA and Radloff et al. (2010) using morphometric characteristics, who identified 4 major groups of haplotypes and 6 morphoclusters of A. cerana, respectively. Compared with the western honey bees (A. mellifera), which were divided into nearly 30 subspecies (Engel, 1999; Ruttner, 1988), A. cerana classification is unclear. And all these 3 studies classified the majority of A. cerana in China as one group.

However, China occupies an extremely vast territory, having various and complex geographic environments and climate zones with a high degree of biodiversity. Previous researchers have studied ecotypes of A. cerana in China and classified A. cerana into nine ecotypes, including Hainan, Yungui Plateau, Tibet, Aba, Changbai mountains, south Yunnan, north China, south China, and central China (China National Commission of Animal Genetic Resources, 2011). Zhu et al. (2022) divided A. cerana in China into 14 morphological clusters, including Hainan, southern Yunnan, Changbai mountains, Taiwan, Bomi (southeastern Tibet), Aba, Batang, Derong, Yajiang, Chuandian, Chuangui, Northwest, Southern and Northern. There also have been several studies with Chinese A. cerana populations suggesting large genetic differences, some of which are at the subspecies level (Chen et al., 2018; Ji et al., 2020; Li et al., 2019; Liu et al., 2016, 2022; Tan et al., 2016; Xu et al., 2013; Yin & Ji, 2013; Zhao et al., 2014) or high morphological variation (Liu et al., 2022; Peng et al., 1989; Radloff et al., 2010; Tan et al., 2003, 2006; Zhou et al., 2016; Zhu et al., 2017), even among populations from small geographic regions. Overall, these researches indicate great diversity of A. cerana in China, which still remains to be explored further.

Advances in understanding the drivers of divergence between honey bee populations have provided clues for further exploration. Notably, Chen et al. (2018) found that geographical isolation such as high mountains and islands is the driving force for populations differentiation. Therefore, the QTP and islands in the South China Sea become interesting candidates. QTP lies in the southwest of China, with an average altitude of 4,000 meters above sea level. Since A. cerana is mainly observed below 3,500 meters (Yu et al., 2019), the eastern and southern parts of QTP could be suitable for the residence of A. cerana. This area is also considered an important biodiversity hotspot (Myers et al., 2000) that is a priority area that requires pollinator protection (Teichroew et al., 2017). The islands, including Hainan, Weizhou and Yongxing, located at the southern end of China in the tropical monsoon climate zone, are rich in flowering plants. Hainan Island, the second largest island in China, has a high degree of bee diversity (Chen et al., 2018; Ji et al., 2020). Weizhou Island and Yongxing Island are small in size (total areas of 27.74 and 2.13 km² respectively). All these 3 islands are located in the South China Sea, with nectar sources to support populations of honey bee. As islands, they are geographically isolated and considered areas worth investigating. We hypothesize that these regions have potentially new A. cerana populations that are differentiated from known populations.

The methods for morphological studies of honey bees mainly include standard morphometry (SM) and geometric morphometrics (GM). As a traditional method, standard morphometry uses 36 morphometric characters and provided good discrimination among honey bee subspecies and populations (Ruttner, 1988). Geometric morphometrics in bees utilizes the wing shape, that is, landmark coordinates which are superimposed by translation, scaling and rotation in order to remove the effect of size (Bookstein, 1991; Rohlf & Marcus, 1993). Compared with the SM method, GM is only based on forewings as a new and labor-saving approach, without the requirements of other body parts. And now, geometric morphometrics were proven to be effective and widely used in the taxonomy of A. mellifera (Francoy et al., 2008; Kandemir et al., 2011; Miguel et al., 2011; Oleksa & Tofilski, 2015; Tofilski, 2008) and A. cerana (Frunze et al., 2022; Ji et al., 2020; Liu et al., 2022; Zhou et al., 2016).

In this study, we collected samples from QTP, Hainan Island, Yongxing Island and Weizhou Island, including regions that were previously thought to lack A. cerana populations. Furthermore, traditionally well-known populations were also sampled for comparison. To understand the variation among the samples, we used both standard morphometry and geometric morphometrics to explore the current situation of A. cerana population differentiation in the studied areas.

Materials and methods

The sampling of Apis cerana bees

A. cerana samples were collected from 17 regions (Figure 1, Table S1). Of these, 12 regions are located in QTP, Hainan Island, Weizhou Island and Yongxing Island, and 5 reference populations are located in areas of relatively low and flat terrain in southern and eastern China. The sampling sites for the references span a large area of ∼600 km from north to south and ∼1,200 km from east to west.

Figure 1. Locations of A. cerana sampling sites.

Worker bees were collected from colonies from wild nests or human-made hives populated with wild-caught swarms. A total of 140 colonies were collected, with five to ten colonies for each sampling site. For each colony, no less than 30 worker bees were collected and preserved in ethanol. Samples were deposited in the Institute of Apicultural Research, Chinese Academy of Agricultural Sciences.

Data acquisition

Ten worker bees per colony were dissected for morphometric analysis according to the methods described by Ruttner (1988), and the dissected body parts were photographed by a LEICA DMS300 digital microscope. For the standard morphometry method, a total of 33 morphometric characters, corresponding to Ruttner (1988) numbers of (1) to (3), (5) to (18), and (21) to (34), and two more characters, number of hamuli (Nh) and cubital index (Ci), were measured (Table S2). For the geometric morphometrics method, 19 landmarks on the forewings described by Meixner et al. (2013) were identified.

TPS software series (Rohlf, 2008a, 2008b) were adopted to acquire 19 landmarks and size-related characters. Since the 19 landmark coordinates include the vertex coordinates of the wing venation angles, this paper used the landmark coordinates and a custom script in R statistical software (version 4.2.0) (Vienna, Austria) to get the angles. Then a subjective empirical classification of the pigmentation was done, classified into 10 grades ranging from 0 to 9 (from dark to light). Characters of hair were measured with a stereo microscope (ZEISS Stemi 508).

Morphometric analysis

For SM method, colony mean values were calculated for each character. For the GM method, the landmarks were aligned using the Procrustes fit method in R statistical software using the gpagen function in geomorph package (version 4.0.4) and then averaged across colonies. All further analyses were performed using the average values in R statistical software.

Principal component analysis (PCA) and linear discriminant analysis (LDA) of colony means were performed using geomorph and MASS (version 7.3-56) packages and plotted using ggplot2 (version 3.4.0) and factoextra (version 1.0.7) packages. The standard morphometry analysis was performed including all measured characters except (1) due to missing data. Relationships between the 33 morphometric characters and longitude, latitude and altitude were analyzed with the Pearson correlation using the Hmisc package (version 4.7-1). Multiple comparisons were corrected using the Bonferroni method.

Results

Standard morphometry analysis

A principal component analysis was performed to compare the 12 populations from QTP and islands to the 5 reference populations from plains (Figure S1). The island populations occupy the negative axis of PC1 and the positive axis of PC2, with little overlap with other populations. The QTP populations occupy the center and right of the PC1 axis, and the reference populations occupy the center of the PC1 axis and the negative axis of PC2. The island populations, some QTP populations and the reference populations are basically distinguishable from one another, while CYU, another QTP population, overlaps with the reference.

Further examination of the PCs shows PC1 accounted for 33.36% of the total variance while PC2 and PC3 accounted for 9.38% and 7.20%, respectively. PC1 was mainly associated with size-related characters (17), (7), (15), (16), (13), (18), (6), (11), (5), (12), (9) and (10); PC2 was mainly associated with color-related characters (33), (32), (34), cubital index and wing venation angle A4 (21); PC3 was mainly associated with angle-related characters (30) and (27) (Table S3). The characters with high importance were visualized as the quality of representation values in Figure S2. Out of the top five characters, 3 were size-related and 2 were color-related (Figure S2a), and all three pigmentation characters are among the top 10 characters (Figure S2b). The means and standard deviations of the 33 characters for each sampling site are given in Table S4.

A discriminant analysis of the same dataset was also performed. LD1 (37.20% of total variance) were plotted against LD2 (18.47%) in Figure 2. It could be seen that the reference populations (CQ, JX, HS, JH and JXU) were clustered closely together and overlapped with CYU from QTP. Two of the three islands, YX and HN, formed two separate clusters in the lower left corner of the plot, while the third island, WZ, slightly overlaps with DJ from QTP. The QTP populations, which occupied lots of positions in the center and right of LD1 axis and the positive axis of LD2 in the plot, indicate the obvious morphological differences. It was also noticed that the 5 reference populations were packed tightly together and formed much more compact clusters than other populations, indicating less intra-region variation in the reference populations and larger intra-region variation in the QTP and island populations.

Figure 2. Discriminant analysis plot of A. cerana populations based on colony means of 32 characters measured (LD1 vs. LD2), with confidence ellipses at 95%. Each symbol represents one colony. All reference populations (CQ, HS, JH, JX, and JXU) are shown in gray with different shapes for different populations.

Next, we explored the Pearson correlations among the 33 morphometric characters and longitude, latitude and altitude. We observed most hair-related characters and size-related characters are significantly positively correlated with the latitude or altitude, and significantly negatively correlated with the longitude (|r| > 0.25, p < 0.05, Table S5), suggesting with increasing latitude or altitude, the size of wings and body and the length of hairs increase, and with increasing longitude, the body size and wing size decrease. All the color-related characters, including (32), (33) and (34), are significantly positively correlated with the longitude and are significantly negatively correlated with the latitude and altitude, indicating with decreasing longitude or increasing latitude and altitude, the color of tergites become darker.

Geometric morphometrics analysis

We performed principal component analysis on the average wing morphology of the 17 populations, and obtained values for 34 principal components (the number of landmarks times the number of dimensions minus the degrees of freedom [4 for 2D], which is 19*2 − 4 = 34). In Figure S3, where principal component 1 (27.86%) is plotted against principal component 2 (12.69%), the QTP populations occupy the left and middle of the scatter plot, the reference populations occupy the middle, and the island populations occupy the right. And part of QTP populations, island populations and reference populations overlapped with each other.

Figure 3 showed a scatter plot of discriminant analysis of 34 principal component values in PCA of GM, in which LD1 (27.86%) is plotted against LD2 (16.85%). And the same as Figure S3, the QTP populations occupy the left and middle of the scatter plot, the reference populations occupy the middle, and the island populations occupy the right while the three kinds of populations have large overlaps with each other. It was obvious that YX formed a separate cluster on the right of the plot, and CYA and BY formed two separate clusters in the top left corner.

Figure 3. Discriminant analysis plot of A. cerana populations based on wing shape (LD1 vs. LD2), with confidence ellipses at 95%. Each symbol represents one colony. All reference populations (CQ, HS, JH, JX, JXU) are shown in gray with different shapes for different populations.

Discussion

Our study gained a new understanding of A. cerana populations in China. We made special efforts to collect samples from the QTP, which we hypothesize to have a high diversity according to previous studies (Chen et al., 2018; Ji et al., 2020; Liu et al., 2022; Tan et al., 2016; Zhu et al., 2017, 2022). However, these studies either only include two to four commonly known sampling sites in the QTP, or sampling sites are skewed toward a certain region such as Sichuan Province, one of six provinces where the QTP is located, leaving many regions in QTP unexplored. In order to overcome the limitation of the lack of samples on bee diversity in the QTP, the sampling regions were further expanded in this study. Notably, we were able to collect samples from BY, a region that was previously reported to be free of A. cerana (Zhu et al., 2017). In addition, island samples, which are also isolated populations, are collected. We collected island samples in Hainan Island (HN), Weizhou Island (WZ) and Yongxing Island (YX). Among them, there have been a lot of studies on A. cerana in Hainan Island, which was considered a distinct ecotype (Wang et al., 2021; Zhao et al., 2014, 2017; Zhou et al., 2018), and very few studies in WZ and YX. To our knowledge, this is the first study to report the existence of A. cerana in BY, and the first morphometrical characterization of populations of CYA, CYU, DJ, BY, WZ and YX.

This study found the high morphometrical variation of A. cerana from QTP and islands, which can be distinguished by both standard morphometry analysis and geometric morphometrics analysis. Forty-five colonies from 5 plain regions in China were selected as reference populations, including CQ from an adjacent region of QTP and JXU from an adjacent region of islands. The reference populations are chosen in such a way that they covered a large area of south and east China, where the terrain is relatively flat and the altitudes are relatively low, without major physical barriers. Both SM and GM analyses showed large variation in QTP and island populations. The reference populations formed much more compact clusters, while QTP or island populations formed clusters with relatively less or no overlaps with each other, indicating the observed morphological differences within the reference populations were significantly smaller than those within QTP populations or island populations (Figures 2 and 3), even though the distribution range of the reference populations is about the same as that of QTP populations and much larger than that of island populations. In addition, the large overlaps among QTP populations reflect our dense sampling sites along the edge of QTP, because the smaller the distance between sample sites is, the less distinct the morphoclusters are (Radloff & Hepburn, 1998).

Our analyses identified three new morphological clusters that were previously unknown. Previous studies divided A. cerana in China into 9 ecotypes and 14 morphological clusters respectively (China National Commission of Animal Genetic Resources, 2011; Zhu et al., 2022). Four of these populations are also used in this study, including HN (Hainan ecotype, Hainan cluster), DR (Derong cluster), MEK (Aba ecotype, Aba cluster), and BM (Bomi cluster), which are considered to have a high differentiation and could also serve as valuable references. In Figure 2, the discriminant analysis indicated YX formed a more distinctive cluster even than HN. In Figure 3, the three most obvious clusters are formed by CYA, BY and YX. All these indicated that the morphology of honey bees in CYA, BY and YX had higher morphometric differentiation than others and could probably be new morphological clusters.

Our results further support the hypothesis that physical barriers are one of the main driving force of population differentiation (Chen et al., 2018), as indicated by the high variation of isolated populations located in islands or mountainous areas. These barriers may cause population differentiation through the reduction of gene flow. Genetic studies to quantify gene flow among the populations are needed to gain a further understanding of the differentiation of A. cerana populations.

In this study, we used characters of length, pigmentation, and angels. In the PCA of SM, characters of pigmentation showed high weight inPC2 and moderate weight in PC1 (Figure S2, Table S3). Previous studies have used pigmentation to distinguish between A. cerana and other bee species (Hadisoesilo et al., 1995; Lo et al., 2010; Oldroyd et al., 2006). However, the role of color in the classification of A. cerana subspecies remains unclear. Ruttner et al. (1978) and Ruttner (1988) described the measurement standard for measuring pigmentation, which was traditionally developed for A. mellifera. Rinderer et al. (1989) found neither A. cerana nor A. koschevnikovi followed the exact patterns indicated for A. mellifera. Wang et al. (2006) and Wang (2007) made new criteria for scutellum, tergites and labrum for A. cerana, but we found it is not applicable to many of our samples in this study. Meanwhile, many studies did not include pigmentation characters. Another drawback of pigmentation characters is the subjectivity of measurement, especially when a new sample does not fit readily into any of the existing patterns. A recent study by Qiu et al. (2023) promoted the use of wing morphological geometrics because wing vein angles were less associated with environmental factors. Here, our results support this approach in that the main pattern revealed by SM and GM is in accordance.

We observed correlations between morphometrical characters and environmental factors. Correlations of hair length, wing size, body size and pigmentation with latitude or altitude are observed, which is consistent with previous observations in A. cerana (Hepburn et al., 2001; Tan et al., 2003, 2006; Verma et al., 1994; Zhu et al., 2017, 2022). This could be related to the capability of nectar collection in low temperature environments, with greater resistance to the cold, foraging range (Greenleaf et al., 2007). It is worth noting that correlations were also observed between longitude and some morphometrical characters, however, this could be the result of the strong correlation between longitude and altitude (r = −0.81, p < 0.001) in our sampling sites. The terrain of China is high in the west and low in the east (You & Yang, 2013), resulting in the QTP samples collected at high altitudes with low longitude, and island or reference samples at low altitudes with high longitude.

At present, it is urgent to protect the native colonies of A. cerana in China. In the past decades, the number and distribution of A. cerana in China have declined by more than 75% (Yang, 2009), possibly due to parasites, pesticides, lack of flowers and the introduction of A. mellifera. One more thing to note is that the local A. cerana resources in different regions may also be at risk. Many colonies or queens of A. cerana from other regions were introduced by beekeepers for economic benefits, which may lead to the loss of local populations and therefore the loss of desirable traits adapted to the local environment. Previous studies highlighted the importance of promoting the use of endemic honey bee populations in apiculture (De la Rúa et al., 2013; Harpur et al., 2013). Here, we also call for the promotion of the use of local A. cerana populations. Studies found there are many small populations of A. cerana in QTP, whose effective population size is generally lower than 500 colonies (Yu et al., 2019; Zhu et al., 2017), making these populations unique and precious genetic resources but with a high risk of extinction. Honey bee populations in Weizhou Island and Yongxing Island are also small, due to the very small size of the islands. Our study showed the high level of morphological diversity of A. cerana from QTP, Hainan Island, Weizhou Island and Yongxing Island. Therefore, these particular A. cerana populations need to be protected as valuable genetic resources.

Conclusions

Through multivariate statistical analyses of both standard morphometry and geometric morphometrics, this study characterizes A. cerana populations to provide insights into the differentiation of A. cerana in China. We sampled A. cerana in a more comprehensive region in China, including areas where previously thought to be free of A. cerana, or where A. cerana is not previously morphologically analyzed, and large morphometrical variations were observed in the studied areas, especially CYA, BY and YX. The results of this study provide information that could elucidate the differentiation of A. cerana in different environments and could be used for the effective protection of A. cerana diversity in China.

Acknowledgments

The authors thank Yuyang Shen, Lijiao Gao and Gulinuer Tulaxi for their laboratory assistance, and beekeepers and local helpers for their help in collecting the samples.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was supported by the National Natural Science Foundation of China (Grant No. 31902219), the earmarked fund for China Agriculture Research System (Grant No. CARS-44-KXJ1), and the Agricultural Science and Technology Innovation Program (CAAS-ASTIP-2023-IAR).