Apis cerana is known as an excellent pollinator of various economic crops with a wide foraging range (Wang et al., 2021). The native habitat of A. cerana, the Republic of Korea, is a peninsula that is separated by the ocean from other countries in the south, west, and east, having only a border with North Korea in the north region that restricts any random introduction of honey bees outside. The Republic of Korea is located on the east coast of the Eurasian Continent and adjacent to the Western Pacific. It shows climate characteristics with an annual mean temperature from 10 to 16°C that provided a successful environment for the vegetation of honey plants from April to November (https://www.weather.go.kr/eng/biz/climate_01.jsp).
The subspecies determination of honey bees in the Republic of Korea had a complicated path. In this location, Lew (1992) reported that the A. c. indica honey bee was kept in the first part 17th century, but their subspecies traits were not described. Later in the 20th century, the population of A. cerana species in Korea was not be identified to the subspecies level, because the general scientific recommendation was not accepted to estimate the A. cerana subspecies neither by morphometric or genetic methods (Lee and Choi, 1986;Kwon and Huh, 1992, Park et al., 2015). The majority of A. cerana subspecies have no clear taxonomic status by morphometric characters yet (Sarah et al., 2005;Radloff et al., 2010;Ilyasov et al., 2019, 2021; Proshchalykin et al., 2020). However, the differences between A. cerana honey bees from plains and hills regions were found more than in the north and south regions which were studied by morphological and biochemical characteristics (Lee, 1993). Moreover, the A. cerana honey bee in Korea (Seoul) had a distinct morphocluster out from A. cerana in Japan by the 12 morphometric characters (Radloff et al., 2005;2010).
Unfortunately, in the first part of the 20th century in Korea the number and distribution of A. cerana colonies was decreased, because a number of A. mellifera colonies were introduced to increase the beekeeping resources (Ai et al., 2012).
Additionally, in Korea, the sacbrood virus (SBV) disease was the reason for more than 90% of eastern honey bees’ mortality in 2008 (Choi et al., 2010;Vung et al., 2017). However, the A. cerana honey bee resistant to SBV was selected from the remaining honey bees and bred in Korea (Vung et al., 2020). Researchers utilized the hygienic behavior in honey bees to confer colony-level resistance against brood disease for artificial selection of A. cerana colonies (R, H) and checked the brood survival rate after colonies were SBV-inoculated (Vung et al., 2020). Also, some Korean beekeepers manage the A. cerana honey bees which survived by natural selection after extensive eastern honey bees’ mortality.
These honey bees were identified for the first time by not morphometry, but by the genetic method. Ilyasov et al. (2018) used modern genomic technologies to discover and prove the taxonomic uniqueness of Korean honey bee subspecies. So, novel subspecies A. c. koreana resistant to SBV virus in Korea (Ilyasov et al., 2018), has recently been identified and it was uploaded to the DDBJ/Genbank database (AP018431). Scientists performed a comparative analysis of the complete mtDNA sequence of the Korean, Japanese, Chinese, and Taiwan A. cerana population and performed phylogenetic analysis (Ilyasov et al., 2018). So, report above confirmed the differences on genetic level of the A. c. koreana breed (Vung et al., 2020) out from other A. cerana. Although, the differences on the colony level between A. cerana honey bees was not studied in Korea, it was found between A. mellifera honey bees after breeding in closed population for 4 years by Szabo and Lefkovitch (1988).
Similar achievements were studied by morphometric methods in our research. To differentiate the A. c. koreana, we investigated the 22 classic morphological characteristics, the royal jelly secretion, the weight of workers, queen, and drone bees of A. c. koreana bred in the Republic of Korea. To define the result of selection, we additionally applied the geometric morphometric analysis. Specifically, the aim of this research is to provide phenotypic information that can compare A. c. koreana with other A. cerana subspecies from open resources and determine breeding effect by morphometric traits. This study provides the understanding of the A. cerana honey bees’ morphometric response to the geographic isolation and selection. The characteristic of honey bees is important to identify, protect, and preserve of honey bee stock surviving in the new hazards arising in climate change cases.
Materials and Methods
Honey bees
A. c. koreana honey bees (R, H, RH) against sacbrood virus disease were bred by artificial selection on the breeding station using a method of pure breeding from 2015 year (Vung et al., 2020). The honey bees (L) of natural selection were keep in the breeding station without any treatments. The honey bee breeding lines have been controlled through the production of the queens from a selected colony and mate using the artificial insemination and later natural mating method by isolating station on the island.
The names of honey bee colonies and lines L, R, H, RH are at random. A total of 15 A. c. koreana colonies were sampled at the same time in July 2020 from the breeding station in National Institute of Agricultural Science (Wanju-gun) in the south-west of the Republic of Korea.
Classic morphometric method
Thirty-five worker bees were collected from each assigned colony and killed immediately with hot water (higher than 80°C) to extend their proboscises and were then preserved in 75% ethanol until use. 24-32 worker bees from each line were dissected and quantified for morphological analysis using the methods described by Ruttner (1988) and Carreck et al. (2013). A total of 22 characteristics associated with production performance were quantified from 107 sampled honey bees using a Leica MZ16 A stereomicroscope and the TCapture computer program.
For classic morphometric method, the twenty-two characters are the Length of body, Length of abdomen, Width of abdomen, Length of tergite 2, Length of tergite 3, Length of tergite 4, Length of head, Width of head, Length of proboscis, Length of tongue, Length of antenna, Tarsal index, Length of femur, Length of tibia, Length of metatarsus, Width of metatarsus, Length of forewing, Width of forewing, Length of hind wing, Width of hind wing, Ci Goetze, and the Hamuli of hind wing. Means of colony sample, standard deviation, and Coefficient of Variance (CV) were computed for each character from the samples. Eight morphological characters were common to the databases of A. cerana from Korea, China, India, Sri Lanka, Philippines, and Japan. This eight characters were the Cubital index, Length of proboscis, Length of femur, Length and Width of metatarsus, Length and Width of forewing, and the Length of tibia. Three morphological characters were common to the database of A. cerana honey bee studied earlier in Korea by Lee and Choi (1986). This three characters were the Tarsal index, Length of forewing, and the Cubital index. Multivariate statistical analysis of the classic morphometric data included Principal Components Analysis (PCA) to identify possible splitting of A. c. koreana. The hierarchical cluster analysis was done to identify the colonies within A. c. koreana and with data of A. cerana honey bees. The linear discriminant analysis carried out to determine the characters which had high linear discriminant function as result of highest contribution in determination of colonies, and defined the percentage of correct classification of colonies in A. c. koreana honey bees.
Geometric morphometric method
For the geometric morphometric method, the 50 right forewing images for every A. c. koreana R, H, L honey bee lines were photographed and stored in the computer. A total of 19 landmarks on the forewings were identified (Fig. 1). Tps files were prepared using tpsUtil 32, and landmarks were digitized into images using tpsDig232 (Ferdous and Armbruster, 2012). Later, landmarks were superimposed using a generalized leastsquare algorithm in MorphoJ (Klinberger, 2011;Berezin, 2019). This landmark-based morphometric method removes all non shape variation that can be attributed to differences in the location, orientation, and scale of the specimens (Koca and Kandemir, 2013). A multivariate analysis of variance (MANOVA) and pairwise tests were carried out on the landmark data by using MorphoJ in order to compare A. c. koreana honey bees. Superimposed x, y coordinate data were than used as the data set for multivariate statistical analysis of A. c. koreana honey bees. To assess the variation between the studied honey bees, multivariate statistical analysis (Canonical Variate Analysis, PCA) were carried out on geometrics morphometric data collected from the same set of samples taken for classic morphometric analysis.
Weights of Apis cerana koreana workers, queens, and drones and royal jelly secretion by worker bees
To study the weight, the live worker and drone bees were collected from A. c. koreana colonies in totals of 470 worker bees, 78 drones from queenright colonies, and 151 drones from laying worker colonies. The virgin queens were reared from A. c. koreana colonies accordingly to recommendations for A. cerana honey bees using queen cells on the frame (Vung et al., 2017). The queen pupae in cells (11 days from grafting larvae) were put in the queen cage individually and were placed in an incubator at 35°. The next day the weight of the live virgin queens were measured by the scales. Fertile queens were replaced from queenright honey bees’ colonies and measured immediately.
To study the royal jelly secretion, the standard methods for A. mellifera L. royal jelly research were used (Hu et al., 2017) with modification for A. c. koreana honey bees by Vung et al. (2017) from June to September in 2020. Frames with grafted larvae were removed from rearing colonies after 48 hours and taken to the bee breeding laboratory for royal jelly collection and measurements on the electric scale balance.
Statistical analysis
Data were analyzed by using the MS Excel with the XLSTAT (https://www.xlstat.com/en/), AtteStat application (http://offext.ru/ library/science/development/87.aspx), program R, SAS Enterprise Guide 7.1 program. The one-way ANOVA and Tukey post-hoc test (in case for three and more groups) or t-test (in case for two comparing groups) provided the testing of the significant differences between the multiply means of the morphometric characters of honey bees. In all stages, the significant levels were set at ɑ = 0.05.
Results
Classic morphometric method
The group means, standard deviation, and CV of colony means for a total of 22 morphological characteristics were presented in Table 1. CV were calculated and shown the inconsistent variations among the different morphological characteristics of A. c. koreana honey bees (Table 1). The CV of Cubital index the Length of abdomen was the highest among all morphological characteristics, and the average CV was 13.165 and 10.370% respectively. By contrast, the CVs of other morphological characteristics did not exceed 9% and the Length of forewing had the lowest CV (1.719%). It should be noted that low CV demonstrates the extent of variability of morphometric characters in relation to the mean of the population, we suppose that morphometric characters with the low CV is the candidate to assess the breeding treatments as the most accurate characters.
To check breeding clusters, a linear discriminant analysis was used. Lastly, four groups were entered for the following 22 characters with definition of their discriminatory powers (Table 2). The results classified 100% of data for both A. c. koreana R (n = 25) and A. c. koreana H (n = 24), and 96.88% A. c. koreana RH (n = 32) with one data misclassifying into artificially selected A. c. koreana H group of honey bees; 92.31% of the naturally selected A. c. koreana L honey bees group (n = 26) but two data misclassified into A. c. koreana H group.
To test the equality of the groups means for the characters used in the discrimination function, Wilks lambda approximated by the F-statistic was determined. A significant difference between the means of the four groups was detected (λ = 0.008, F66.00 = 14.80, P < 0.0001). Next, the eight characters with scores more than 300 were extracted which could allow a comprehensive assessment of A. c. koreana: the Length of tergite 4, Length of antenna, Length of metatarsus, Length of hind wing, Length of tergite 3, Length of forewing, Width of head, and the Width of forewing (in decreasing order of the linear discriminant function) (bolded characters in Table 2).
Additionally, the hierarchical cluster analysis by these eight characters with high score of linear discriminate function (Fig. 2A) indicated that these can slightly discriminate A. c. koreana honey bees on branch of artificial (A. c. koreana R, RH, H) and natural (A. c. koreana L) selection. Selected A. c. koreana honey bees confirmed significant differences (one-way ANOVA, Tukey post-hoc test, λ=0.05) only in three morphological characteristics with high score of linear discriminant function. These characters were the Length and Width of forewing, and the Length of hind wing (Table S1). Although the breeding results of A. c. koreana honey bees were identified by these three traits, we could not differentiate the lines A. c. koreana by general 22 classic morphometric characters.
Besides, we compared the morphometric characters of A. cerana honey bees from different countries to define the peculiarities of A. c. koreana honey bees in Korea. A. c. koreana is distinct from other A. cerana (Fig. 2B) by the hierarchical cluster analysis based on our research and literature data of eight total morphometric characters (Table S2). A. c. koreana have been included in branch MI and did not match to MII and MV morphoclusters, according to the classification of A. cerana honey bees (Radloff et al., 2010). So, we found the proximity A. c. koreana, A. cerana in China, and A. c. indica honey bees in one brunch by the eight total classic characters which, perhaps, relay the origin of the A. cerana honey bees. In this study, means of the Length of proboscis, Length of forewing A. c. koreana were higher than A. c. indica (Rajkumari et al., 2020), and A. cerana Sri Lanka (Szabo, 1990), but lower than A. cerana in China (Zhu et al., 2017). Moreover, the mean Cubital index of artificially selected A. c. koreana (RH, R) honey bees was lower than those of other A. cerana honey bees in the research. Hence, A. c. koreana honey bees differed from A. cerana subspecies by honey bee morphology due to geographic isolation and adaptation to breeding environment in the Korean peninsula.
Additionally, to test a result of selection by the differences of the group A. c. koreana with A. cerana studied by Lee and Choi (1986) in the same place in Korea, we studied and found that there were increase in the Tarsal index, Length of forewing, and the stability in the Cubital index at 56.537 and 53.530%, 8.625 and 8.274 mm, 3.30 and 3.31 respectively.
Geometric morphometric method
Shape coordinates were superimposed to successfully eliminate the size effect, which was apparent from Procrustes analysis. The deformed wireframe was drawn on the shape among three lines of honey bees (A. c. koreana R, H, L) to interpret shape changes that support the Relative Warp (RW) analysis. This illustrated deformation in shape from the control group A. c. koreana L that corresponds to selected positions in the ordination. Canonical Variate Analysis (CVA) extracted Mahalanobis and Procrustes distances among groups found to be highly significant (p < 0.0001) (Tables 3, 4). Canonical variate 1 and canonical variate 2 explained 50.710% and 36.014% of the total variance, respectively, while canonical variate 3 explained only 13.276% of the total variance. Classification results of CVA indicated that all the specimens of each group were allotted to their respective groups. So, wing vein geometric morphometrics were able to clearly identify differentiation between samples of A. c. koreana honey bees R, H, L by CVA (Fig. 3).
Body weights of Apis cerana koreana workers, queens, and drones and royal jelly secretion by worker bees
The differences in Body weight of honey bees can be caused by many factors such as temperature and species of plants. The mean weight of worker bees A. c. koreana was lower than A. cerana in China but higher than A. cerana in Indonesia, showing 69.29 ± 15.32 mg in Korea, 73.95 ± 0.55 mg in China, and 46.060 ± 7.015 mg in Indonesia, respectively (Yue et al., 2017;Widiawati et al., 2020).
The means weight of drones from the queenright colony were higher than that of the drones from laying worker colony A. c. koreana, drones A. cerana in China (Phokasem et al., 2021), providing 118.17 ± 12.58, 91.56 ± 11.69, and 103.56 ± 12.51 mg, respectively. The means weight of virgin and fertile queen A. c. koreana were 150.42 ± 15.02 and 185.33 ± 20.64 mg, respectively. The weight of fertile queen A. c. koreana in this research was higher than that reported for fertile queen bee A. c. koreana (169.00 ± 0.16 mg) at the same location in Korea (Vung et al., 2018) (Table 5).
Despite that the A. cerana honey bee is not used as royal jelly producers, the amount of secreted royal jelly by worker bees influences the honey bee colony development such as the size of the honey bees (Szabo and Lefkovitch, 1988). The weight of body of the artificially and naturally selected A. c. koreana honey bees did not have significant differences (t-test, P > 0.01), but the royal jelly secretion of the artificially selected honey bees was higher (t-test, P < 0.01) than naturally selected honey bees, providing 132.96 ± 20.70 (n = 24) and 113.00 ± 25.35 mg (n = 22), respectively.
Discussion
The results presented here show that geometric morphometric yielded marginally better discrimination of the artificial and naturally selected Apis cerana honey bees than classic morphometry, which also was marked by Tofilski to differentiate Apis mellifera honey bee subspecies (2008). Better discrimination on geometric morphometric traits compared to classic morphometry is not surprising, as the former used 19 variables as opposed to the 8 variables from 22 classic morphometry. Last 8 characters were selected by discriminant analysis for the classification model only in the variables that contribute to the discrimination. This solved the problem, when a large number of variables are used, although some of them become redundant because of the correlations between them (Tofilski, 2008). However, from eight extracted characters confirmed three with significant differences (one-way ANOVA, Tukey post-hoc test, α = 0.05) were the Length of hind wing, Length of forewing, and the Width of forewing. It can indicate the low identification of artificial and naturally bred A. cerana koreana honey bees by the classic morphometric method which was survived bred after mortality in 2008. Also, the significantly intensive royal jelly secretion were found at artificially bred against naturally survived and bred A. c. koreana honey bees in Korea, indicating the high potential development of artificially bred honey bees.
The aim of this research is to provide phenotypic information that can compare A. c. koreana with other A. cerana subspecies and determine breeding results by morphometric traits. Increase trends were noted for the A. c. koreana honey bee based on results of current research and on literature of the 1980s concerning A. cerana sizes in Korea measured by the Tarsal index, Length of forewing, and the Cubital index. The cluster analysis revealed the proximity of A. c. koreana, A. cerana in China, and A. c. indica honey bees by the eight total classic characters which, perhaps, relay the origin of the A. cerana honey bees. If the A. cerana worker bees in China having higher values of the Length of proboscis, and Length of forewing than A. c. koreana honey bees, but last one had higher these two morphometric characters, than A. c. indica honey bees. The differences between artificially and naturally selected and bred A. c. koreana honey bees were found by 19 landmarks of the geometric morphometric method, but not by 22 characters of the classic morphometric method. However, the artificially selected honey bees secreted significantly more royal jelly than naturally selected A. c. koreana honey bees, which positively influence the development of the honey bee colony through enhanced nutrition of the larva. The results of this study are helpful to elucidate the phenotypic evolution of A. cerana honey bees and provide the information for identification, protection, and preservation of honey bee stocks in Korea by analyzing the morphometric characteristic of newly defined subspecies A. c. koreana honey bees to. Also, this study expands our understanding of how honey bees will respond to natural and artificial selection in response to the new climate condition that we are faced with.