Population Structure, Demographic History, and Adaptation of Giant Honeybees in China Revealed by Population Genomic Data

Abstract

There have been many population-based genomic studies on human-managed honeybees (Apis mellifera and Apis cerana), but there has been a notable lack of analysis with regard to wild honeybees, particularly in relation to their evolutionary history. Nevertheless, giant honeybees have been found to occupy distinct habitats and display remarkable characteristics, which are attracting an increased amount of attention. In this study, we de novo sequenced and then assembled the draft genome sequence of the Himalayan giant honeybee, Apis laboriosa. Phylogenetic analysis based on genomic information indicated that A. laboriosa and its tropical sister species Apis dorsata diverged ∼2.61 Ma, which supports the speciation hypothesis that links A. laboriosa to geological changes throughout history. Furthermore, we re-sequenced A. laboriosa and A. dorsata samples from five and six regions, respectively, across their population ranges in China. These analyses highlighted major genetic differences for Tibetan A. laboriosa as well as the Hainan Island A. dorsata. The demographic history of most giant honeybee populations has mirrored glacial cycles. More importantly, contrary to what has occurred among human-managed honeybees, the demographic history of these two wild honeybee species indicates a rapid decline in effective population size in the recent past, reflecting their differences in evolutionary histories. Several genes were found to be subject to selection, which may help giant honeybees to adapt to specific local conditions. In summary, our study sheds light on the evolutionary and adaptational characteristics of two wild giant honeybee species, which was useful for giant honeybee conservation.

Keywords: Apis dorsata; Apis laboriosa; adaptive evolution; demographic history; giant honeybee; population genomics.

Fig. 1.

Geographical distribution and population structure of giant honeybees in China. (A) Geographic distribution of two giant honeybees sampled in this study. A. dorsata were collected from Xishuangbanna, Yunnan province (XS, n = 7); Puer, Yunnan province (PE, n = 3); Qiongzhong, Hainan province (QZ, n = 5); Baoting, Hainan province (BT, n = 4); Chongzuo, Guangxi province (CZ, n = 7); and Nanning, Guangxi province (NN, n = 2). A. laboriosa were collected from Honghe, Yunnan province (HH, n = 9); Baoshan, Yunnan province (BS, n = 6); Diqing, Yunnan province (DQ, n = 4); Linzhi, Tibet (LZ, n = 4); and Rikaze, Tibet (RK, n = 6). Shading represents elevation. (B) Population structure analysis of A. dorsata and A. laboriosa. Samples were clustered into 2–4 groups (K = 2–4), and cross-validation errors are represented by CVe. (C) Scatterplot of principal component analysis (PCA). Principal components 1 and 2 (PC1 and PC2) were calculated for all giant honyebee samples. (D) Maximum likelihood phylogenetic tree using whole-genome SNPs (numbers on branches represent bootstrap value, numbers in the brackets mean time of divergence in million years.

Fig. 2.

Demographic history of giant honeybees in China. PSMC analysis was performed on representative individuals exhibiting inferred variation in N_e over time for (A) A. laboriosa and (B) A. dorsata. For (A), thin lines represent individuals from Yunnan and thick lines represent individuals from Tibet. For (B), thin lines represent individuals from the China Mainland and thick lines represent individuals from the Hainan Island. Horizontal bars reflect three major glaciations: the Naynayxungla glaciation (NG), the penultimate glaciation (PG), and the last glaciation (LG). Red bars represent warm periods. The dynamics of Ne were inferred by SMC++ for (C) A. laboriosa and (D) A. dorsata populations. Generation time (g) = 1 year and transversion mutation rate (µ) = 0.5 × 10⁻⁸ mutations per bp per generation.

Fig. 3.

Identification and annotation of selected genes in giant honeybees. (A) Distribution of log2 (π ratios) and F_ST values calculated using 50 kb sliding windows with 20 kb increments between Eastern Yunnan (including HH populations) and Tibetan (including RK and LZ populations) groups. (B) Distribution of log2 (π ratios) and F_ST values between mainland (including XS, PE, CZ, and NN populations) and island (including BT and QZ populations) groups. For A and B, the red points are genomic regions under selection in Tibetan populations and Hainan Island populations, respectively. (C) Scatterplot of enriched GO terms and KEGG pathways for selected genes of the Tibetan group. (D) Scatterplot of enriched GO terms and KEGG pathways for selected genes of the island group.