1 Introduction

As important pollinators for wild plants and crops, bumblebees provide vital pollination services in natural and agricultural ecosystems (Pywell et al. 2006; Velthuis and van Doorn 2006; Goulson et al. 2008; Julier and Roulston 2009). Bumblebees are widely distributed, and approximately 250 species have been described worldwide (Williams 1998). With diverse landforms and rich vegetation, China harbours the highest richness of bumblebee species in the world (Williams 1998; Williams et al. 2009; Huang and An 2018). To date, 125 bumblebee species have been identified in China, representing 50% of the bumblebee species recorded worldwide (Williams et al. 2017; Huang and An 2018). Although China is the greatest hotspot of bumblebee diversity in the world, the complex taxa of this country are insufficiently described. These taxa include species with diverse intraspecific colour patterns and cryptic species that differ little in colour and morphology. Further detailed species identifications and descriptions are needed to enhance both the conservation of various bumblebee species and our understanding of their roles as pollinators.

A prominent feature of bumblebees is that they are covered with long and brightly coloured hairs that form different colour patterns (Williams 2007). High variation in colour patterns occurs within some species, and the high evolutionary rates of colour patterns even exceed those of mitochondrial genes (Duennes et al. 2012; Hines and Williams 2012; Huang et al. 2015b). In addition, colour patterns often converge among bumblebee species depending on geography and habitat (Williams 2007; Lozier et al. 2013). The exceptional colour-pattern diversity within species and the convergence among species have led to inaccurate descriptions of some bumblebee species, such as the description of the Bombus lucorum complex (Murray et al. 2008; Waters et al. 2011). Studies have revealed that colour patterns cannot be used to identify some bumblebee species (Carolan et al. 2012; Williams et al. 2012). New colour patterns have been discovered from expanded sample collections, and quantitative and systematic studies of colour variations are lacking, especially for those species with multiple colour forms.

The progress in the application of molecular methods has made it much easier and more reliable to identify species and estimate phylogenetic relationships. Nuclear genes and mitochondrial sequences have been used to obtain robust phylogenetic information and have led to revisions in bumblebee systematics (Kawakita et al. 2004; Cameron et al. 2007; Vesterlund et al. 2014). Due to the high rates of sequence changes in mitochondrial cytochrome oxidase subunit I (COI) and constraints on intraspecific divergence, COI barcodes provide an easy and reliable solution for species discrimination and phylogeny reconstruction (Hebert et al. 2003). Studies have demonstrated that COI barcoding is a cost-effective approach to resolve taxonomic uncertainty of bumblebee taxa. It has been used to clarify the taxonomy of the cullumanus group, which comprises three species highly similar in morphology, Bombus semenoviellus, Bombus unicus and Bombus cullumanus (Williams et al. 2013). Furthermore, it has been applied to distinguish the similarly coloured species in the subgenus Mendacibombus (Williams et al. 2016), identify the cryptic species of the subgenus Bombus s. str. (Williams et al. 2012) and recognise species with variable colour patterns, such as Bombus koreanus (Huang et al. 2015b).

Bombus bicoloratus is a medium-sized oriental species with highly variable colour forms and a long tongue. This species was recently reported to belong to the subgenus Megabombus, which is considered very diverse in China (Williams 1998; Williams et al. 2009; An et al. 2014; Huang et al. 2015a). Bombus bicoloratus Smith was first described in 1879 based on specimens from Taiwan that were in the collection of the British Museum. The abdominal terga are entirely ferruginous or fulvous except for a patch of black hairs on the basal middle portion of the first tergite (Frison 1934). Because of the distinct colour morphs of B. bicoloratus from Taiwan, bumblebees from mainland China with yellow bands were previously described as a separate species, Bombus kulingensis, which was grouped into the subgenus Senexibombus with B. bicoloratus (Williams 1998; Cameron et al. 2007). Females of B. bicoloratus from mainland China typically have hairs on the thorax with yellow anterior and posterior bands, sometimes with narrow yellow thoracic bands intermixed with black. Metasomal terga 1–2 and often most of tergum 3 are yellow, and tergum 5 is orange-red (Williams et al. 2009). As a species that is widely distributed over middle-low altitudes in Taiwan, B. bicoloratus has been screened for artificial rearing, and 52% of wild queens can produce small colonies (Ho et al. 2002). Although B. bicoloratus has recently been reported in some regions of China (Williams et al. 2009; An et al. 2014), its distribution patterns, food plants and colour-pattern polymorphism have yet to be investigated in detail.

In this study, with ample samples of B. bicoloratus collected in China in the last decade, we describe this species’ geographic distribution, summarise the diversity and spatial distribution of its colour patterns and document its food plants in China. The results provide information essential for the conservation of B. bicoloratus and for its potential application in commercial pollination.

2 Materials and methods

2.1 Materials

A systematic survey of bumblebees in China has been ongoing since 2002 (Williams et al. 2017; Huang and An 2018). Bees were collected by sweeping with a nylon hand net. Detailed information, including the name, elevation and location of the collection site, was recorded with a hand-held GPS (Garmin 60CS, China). The collected bumblebees were pinned, labelled, dried and deposited in the Institute of Apicultural Research, Chinese Academy of Agricultural Sciences (IAR-CAAS), Beijing, China. Among these collections, 506 B. bicoloratus–like specimens were preliminarily identified according to the keys described by Williams et al. (2009) and An et al. (2014). The records of B. bicoloratus from Taiwan, which have been confirmed by the authors, were obtained from the website http://gaga.biodiv.tw/new23/cp03_76.htm. Photographs of the food plants visited by bumblebees were also collected, and the food plants were identified from these photographs.

2.2 Genetic analysis

For each of the three castes (queens, workers and males), we selected 1–3 bees from the specimens displaying different colour patterns. A total of 85 bumblebees were selected for DNA barcoding. DNA was extracted from the right foreleg using a Wizard SV 96 kit (Promega Inc.) as described in Huang et al. (2015b). The primers LEP-F1 and LEP-R1 (Hebert et al. 2004) were used to amplify the region of the mitochondrial COI gene. We carried out PCR amplification in a 40-μL volume containing 20 μL of PCR Mix (2×), 0.5 μM of each primer and 2 μL of extracted DNA under the following conditions: 95 °C for 5 min; 35 cycles of 94 °C for 30 s, 53 °C for 30 s and 72 °C for 60 s; and a final extension at 72 °C for 10 min. PCR products were subjected to electrophoresis in a 1.0% agarose gel and visualised under UV light. Positive PCR products were sent for sequencing from both ends (SinoGenoMax, Beijing, China).

The forward and reverse sequence reads were assembled into contiguous sequences and trimmed to a length of 657 bp using BioEdit v7.0.5.2 (Hall 1999). The sequences were deposited into GenBank (Accession numbers: KP671640, KP671641, MH986124 and MH986125). Sequence divergence was calculated using MEGA v7.0.14 (Kumar et al. 2016) based on the Kimura-2 parameter distance model with 1000 bootstrap replications. COI barcoding sequences of related species were downloaded from the NCBI database. We used Clustal W (Thompson et al. 1994) to align these sequences and jModelTest v2.1.7 (Darriba et al. 2012) to select the best nucleotide substitution model according to Akaike’s information criterion (AIC). Species of the subgenera Thoracobombus and Subterraneobombus were used as outgroup taxa. A Bayesian tree was constructed using MrBayes v3.2.5 (Ronquist et al. 2012) with two runs of 100,000,000 generations and a sampling frequency of once every 1000 generations. The phylogenetic tree was displayed and edited using FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).

We used three methods for species delimitation. The aligned sequences were analysed using Automatic Barcode Gap Discovery (ABGD) software (Puillandre et al. 2012) with a transition/transversion ratio of 2 and simple distance. The Bayesian Poisson tree processes (bPTP) method (Zhang et al. 2013) was employed using the default parameters. We used Beast v1.8.4 (Drummond et al. 2012) to construct the ultrametric and bifurcating Bayesian tree with two runs with chain lengths of 2,000,000,000 generations and a sampling frequency of once per 200,000 generations. The tree was used to delimit species using the generalised mixed Yule coalescent (GMYC) method (Fujisawa and Barraclough 2013).

2.3 Colour-pattern scoring and analysis of climate variables

Following Williams’ (2007) description of bumblebee colour patterns, the dorsum was divided into 24 elements in females and 26 elements in males. Hair colour was scored separately for each element. The colour patterns were coded as pubescence colour, which covers over 50% of the element area. The minority colour was also coded when it formed strongly contrasting fringes or spots.

To explore whether the spatial distributions of the different colour patterns are correlated with variations in climate and elevation, we performed a canonical correspondence analysis (CCA) in Past (Hammer et al. 2001). For this analysis, we included only those factors that are independent of one another and for which high values could reduce food-plant nectar and pollen production (Williams et al. 2015). These factors included mean temperature of the warmest quarter (bio10), irradiance in the warmest quarter (irradiance), ratio of the precipitation in the wettest month to precipitation in the warmest quarter (ratio), water deficit (wdef) and altitude (alt). Raster images of all these factors were imported in ArcGIS v 10.3, and point values of the above factors were extracted for all the localities of B. bicoloratus.

3 Results

3.1 Species delimitation

The COI barcode was sequenced for the 85 B. bicoloratus samples representing different colour patterns. We obtained 83 sequences of the COI gene, which were 657 bp in length, and we identified only 4 unique haplotypes (GenBank accession numbers: KP671640, KP671641, MH986124 and MH986125), with only 4 polymorphic sites. Haplotype MH986124 was widespread and found in 96% of the sequenced samples. The remaining haplotypes were local types that were detected separately only in Sichuan, Guangdong and Hainan Provinces. Moreover, haplotypes KP671640 and KP671641 differed by only 1 to 3 nucleotides from the most common haplotype, MH986124, and both haplotypes were detected in a single specimen. MH986125 was identified in samples from Hainan. The within-group sequence divergence was 0.4%.

Species delimitation using ABGD, bPTP and GMYC analyses identified the same species groups. These results indicated that 451 of B. bicoloratus–like specimens collected between 2009 and 2017 and representing all the identified colour patterns belonged to the species B. bicoloratus (Figure 1).

Figure 1
figure 1

Species delimitation based on COI barcoding sequences. Values at each node are Bayesian posterior probabilities for the groups. The scale bar represents 0.03 substitutions per nucleotide site.

Full size image

3.2 Distribution of B. bicoloratus

Bombus bicoloratus has a broad distribution in southern China, ranging from the subtropical central Qinling Mountains to the tropical southern Hainan Island (Figure 2). The lowest elevation at which the specimens were collected was 64 m, at the southeast coast of Jinhua City, Zhejiang Province, and the highest elevation was 2302 m, at the edge of the southwest Yunnan-Guizhou Plateau in Tianlin county, Guangxi Province. More than 90% of all bumblebees were found at elevations between 286 and 1615 m, within the provinces of Fujian, Zhejiang, Jiangxi, Hainan, Guangxi, Yunnan, Guizhou, Hunan, Hubei, Chongqing, Sichuan, Gansu, Shaanxi and Henan. The mountain regions around Sichuan Basin, especially in Guizhou (35.9%), Sichuan (12.6%), Hubei (10.2%) and Shaanxi (8.0%), had a high abundance of this bumblebee species (Figure 2).

Figure 2
figure 2

The distribution of B. bicoloratus used in this study. Records of specimens are from the Institute of Apicultural Research, Chinese Academy of Agricultural Sciences (IAR-CAAS), Beijing, China, and the website http://gaga.biodiv.tw/new23/cp03_76.htm.

Full size image

3.3 Colour-pattern variation and their spatial distribution

Colour patterns were scored for the 451 identified B. bicoloratus specimens, which comprised 14 queens, 397 workers and 40 males collected from 102 locations within China. Wide ranges of colour variation were observed in all three castes. Seven different colour patterns were detected among the 14 queens, which were collected from 6 locations, and 12 colour patterns were identified among the 397 workers, collected from 88 sites. The 40 males, collected from 15 locations, exhibited 15 different colour patterns. The workers and queens shared two colour patterns (w3/q3 and w10/q7), and the variation in colour patterns was more pronounced in the females of B. bicoloratus than in the males (Figure 3).

Figure 3
figure 3

The colour patterns of the dorsal pubescence of B. bicoloratus. Twelve colour patterns (w1–w12) were detected in workers, seven (q1–q7) were detected in queens and fifteen (m1–m15) were detected in males.

Full size image

In females, the hairs on the thorax vary from yellow or orange-red anterior and posterior bands to entirely black, and the hairs on the abdomen vary from yellow and orange-red, sometimes mixed with black bands or patches, to completely orange-red. The thorax fringe of some workers is yellow with a sparse admixture of black hairs. Males possess an entirely yellow thoracic dorsum or yellow metanotum with a black mesothorax and a yellow/black prothorax. Metasomal terga 1 and 2 are yellow, whereas tergum 3 is yellow, black or yellow in the middle with a bilateral fringe of black hairs. Metasomal tergum 4 is black, with a posterior fringe of orange-red hairs in some specimens. Terga 5 and 6 are black, orange-red or a mixture of these two colours (Figure 3).

The distribution frequency of the colour patterns was heavily skewed. The most common phenotype among workers was w7, with a yellow patch near the midline of T3 and a black fringe of T4. This phenotype was observed in 50.9% of the worker specimens. Phenotypes w3 and w8 occurred in approximately 10% of the workers, and the remaining colour patterns were present at much lower frequencies (1.3–8.6%) (Figure 4). In most samples (77.3%), metasomal tergum 3 was yellow in the middle with a bilateral fringe of black hairs. For the 10 locations with more than 10 samples, the frequency of colour pattern w7 ranged from 14.2 to 100%. The phenotype w7 predominated in 9 of these 10 locations, in which 1–9 colour patterns were identified.

Figure 4
figure 4

The frequency of the different colour patterns of B. bicoloratus workers deposited in the Institute of Apicultural Research, Chinese Academy of Agricultural Sciences (IAR-CAAS), Beijing, China.

Full size image

Bombus bicoloratus exhibits a striking intraspecific colour-pattern polymorphism across its range, and, in some cases, among samples collected from the same location. The colour pattern of workers collected from Hainan Island (w12) is obviously distinct from the patterns of other specimens from the mainland but identical to the pattern recorded in Taiwan. Worker samples from Wuyi County of Zhejiang Province, Wudang Mountain of Hubei Province, Longli County of Guizhou Province and Yingjing County of Sichuan Province exhibited 5, 5, 7 and 9 different colour patterns, respectively (Figure 5).

Figure 5
figure 5

The distributions of the different colour patterns (w1–w12) detected in B. bicoloratus workers from China. Locality information was obtained from the specimen labels of the Institute of Apicultural Research, Chinese Academy of Agricultural Sciences (IAR-CAAS), Beijing, China, and the website http://gaga.biodiv.tw/new23/cp03_76.htm.

Full size image

The relationships between the colour patterns of workers (w1–w12, from mainland China and Taiwan) and the climatic and topographical variables are shown in Figure 6. The eigenvalues for the first two axes are 0.40 and 0.19. The distinct colour pattern w12, detected from the islands of Hainan and Taiwan, was most positively associated with high values of irradiance, whereas colour patterns w4, w6, w8 and w9 were most positively associated with high values of the topographic factor altitude. Positive associations were also detected between colour patterns w2, w3, w5 and high values of bio10 and between colour pattern w1 and high values of the ratio of precipitation in the wettest month to precipitation in the warmest quarter. The colour patterns w10 and w11 were associated with low values of water deficit. The colour pattern w7 had a wide distribution, and no climatic or topographical variables were highly associated with its spatial distribution (Figure 6).

Figure 6
figure 6

The canonical correspondence analysis (CCA) biplot graph showing the relationships between colour patterns (w1–w12) of B. bicoloratus workers and the five selected climatic variables. Abbreviations: bio10, mean temperature of the warmest quarter; ratio, ratio of the precipitation in the wettest month to precipitation in the warmest quarter; wdef, water deficit; irradiance, irradiance in the warmest quarter; alt, altitude.

Full size image

3.4 Food plants

Based on the photographs collected (Figure 7), B. bicoloratus was found to forage on flowers of 19 plant species (Helianthus annuus, Buddleja officinalis, Buddleja sp., Hypericum sp., Campylotropis macrocarpa, Astragalus sinicus, Leonurus japonicus, Lagerstroemia indica, Oenothera biennis, Astilbe chinensis, Rubus innominatus, Rubus tephrodes, Amygdalus persica, Rhododendron pulchrum, Nandina domestica, Luffa cylindrica, Jasminum mesnyi, Brassica sp. and Raphanus sativus), belonging to 14 families (Asteraceae, Buddlejaceae, Clusiaceae, Fabaceae, Lamiaceae, Lythraceae, Onagraceae, Saxifragaceae, Rosaceae, Ericaceae, Berberidaceae, Cucurbitaceae, Oleaceae and Brassicaceae).

Figure 7
figure 7

Bombus bicoloratus worker visiting Rubus innominatus (Rosaceae) in Guizhou, China.

Full size image

4 Discussion

Our results showed that in China, the oriental species B. bicoloratus mainly occurs in southern China at low to middle elevations. This species was found to exhibit a wide range of colour variations, with workers having more distinct thoracic dorsum and metasomal terga forms than do males, which exhibited many subtle differences. Bombus bicoloratus foraged on 19 plant species during our collections. In addition, in Taiwan, it has been recorded to visit the flowers of Rhododendron (Ericaceae), Prunus (Rosaceae), Trichodesma (Boraginaceae) and Alpinia (Zingiberaceae) (Sung et al. 2011). These findings indicate that this species is polylectic. Bombus bicoloratus has been reported as a major pollinator of Hosta ventricosa (Liliaceae) in Sichuan (Cao et al. 2011), and efforts have been made to artificially rear this species in Taiwan (Ho et al. 2002). Studies have predicted that large areas of overlapping habitat with imported Bombus terrestris will severely threaten many local bumblebee species, including B. bicoloratus (Naeem et al. 2018). As a native species with potential for practical applications, the conservation of B. bicoloratus is of paramount importance. Moreover, the biology of the B. bicoloratus population from mainland China needs further study, and research on the rearing of this species in captivity is necessary.

Variable levels of intraspecific COI sequence divergence in bumblebees have been identified in previous studies. A high amount of differentiation among distinct haplotypes was observed for Bombus cryptarum (Murray et al. 2008; Carolan et al. 2012), which was supported by nine polymorphisms (Carolan et al. 2012). Seven unique haplotypes with six polymorphic sites were found in B. koreanus (Huang et al. 2015b), whereas B. lucorum demonstrated the lowest levels of intraspecific variation, having only a single morph across different populations (Carolan et al. 2012). Previous work found a level of sequence divergence of 0.4% for B. bicoloratus, which was higher than that identified for Bombus lapponicus (0.1%), Bombus monticola (0.2%), B. lucorum (0.1%) and Bombus magnus (0.1%) and consistent with the level detected in B. cryptarum (0.4%) (Carolan et al. 2012; Gjershaug et al. 2013). As suggested by Carolan et al. (2012), this variation may imply species differences in natural migration rates and, thus, in population structure across their distributions. The pronounced divergence in colour patterns between samples from the islands of Hainan and Taiwan and the other areas might represent very different population genetic properties between them, as observed in other Bombus species, because of geographical isolation (Bombus ignitus: Shao et al. 2004; Bombus muscorum: Darvill et al. 2006; B. terrestris: Schmid-Hempel et al. 2007). The systematic significance of such divergence requires further investigation using appropriate markers, such as microsatellites, or omics.

In this study, we identified 17 colour patterns among female B. bicoloratus specimens. Although diverse phenotypes were noted, most of the samples had black hairs on the T3 segment and a yellow patch near the midline, and all samples had red tails, which might be the most representative morphological characteristic of B. bicoloratus. Extensive colour-pattern polymorphism has been reported for many bumblebee species (Williams 2007), such as Bombus trifasciatus (Hine and Williams 2012) and B. koreanus (Huang et al. 2015b). However, detailed studies of colour-pattern distributions and research to determine whether some patterns predominate among the abundant colour varieties are needed to better understand the morphological characteristics of different bumblebee species. In the present study, the spatial distributions of all the colour patterns of worker samples except w7 are correlated with climate variables (Figure 6). Climate variables have previously been shown to be associated with the bumblebee species distribution within the Tibetan region (Williams et al. 2015), but our results demonstrate that climatic variations also are correlated with colour formation in bumblebees. The darker colour of w12, present on islands that are close to the equator, is associated with strong irradiance. However, the associations between colour patterns and climatic variables need to be examined further.

COI barcoding is an effective method for studying closely related bumblebee taxa and recognising species with variable colour patterns (Williams et al. 2012; Williams et al. 2013; Huang et al. 2015b). Consistent with previous studies, our results support the view that colour patterns cannot reliably be used to identify bumblebee species (Carolan et al. 2012), and there is a lack of covariation between COI and colour-pattern variations (Hines and Williams 2012; Huang et al. 2015b). The sympatric and allopatric samples exhibit discordant colour patterns yet have almost identical sequences. The presence of different colour morphs in the same microhabitat suggests that colour polymorphism cannot be strictly induced by the environment and that even if colour is plastic, not all individuals respond equally to the environmental cues. Therefore, it is likely that there is a genotypic or a genotypic and environmental interaction that determines the colour variation. Classical genetic research on Bombus melanopygus and Bombus rufocinctus showed that either a single locus of two alleles or more complex polygenic inheritance could govern colouration (Owen and Plowright 1988). Recently, highly conserved Hox genes and novel post-Hox specialisation were suggested as potentially regulating adaptive colour-pattern variations in bumblebees (Rapti et al. 2014). In addition, analysis of an RNA sequencing data set identified the Xanthine dehydrogenase/oxidase-like gene as a strong candidate for the genetic basis of abdominal pigment variations in Bombus bifarius (Pimsler et al. 2017). The underlying genetic regulation of the formation of colour patterns remains poorly understood despite extensive descriptions of colour-pattern polymorphism. Novel insights could be obtained from the expansion of genomic data sets of whole-genome sequencing, population genomics, gene expression and functional analyses (Lozier et al. 2016).

Resemblance of colour patterns among polymorphic bumblebee species has been surveyed, and Müllerian mimicry has been suggested to explain regional colour-pattern groups (Williams 2007; Hines and Williams 2012). Sexual automimicry has been proposed to explain the selective advantage that male bumblebees, which lack the defence mechanism of females, gain by resembling their noxious female counterparts. Furthermore, thermoregulatory efficiency might explain the colour divergence in males and females (Stiles 1979). Our results showed that the colour patterns of workers were associated with climatic variables. In addition, the difference among colour patterns was found to be more pronounced in females of B. bicoloratus than in males. However, with detailed descriptions of colour polymorphism only in B. bicoloratus and the limited number of male samples, our study is unable to test the above predictions. Additional comprehensive studies of highly polymorphic bumblebee species and of the similarities in the pubescence characteristics of males and females are needed to determine the roles of these proposed selective pressures in colour polymorphism.