Research Article
Research Article
A taxonomic re-assessment of the widespread oriental bumblebee Bombus flavescens (Hymenoptera, Apidae)
expand article infoChawatat Thanoosing§|, Michael C. Orr#, Natapot Warrit|, Alfried P. Vogler§, Paul H. Williams
‡ Natural History Museum, London, United Kingdom
§ Imperial College London, London, United Kingdom
| Chulalongkorn University, Bangkok, Thailand
¶ Institute of Zoology, Chinese Academy of Sciences, Beijing, China
# Staatliches Museum für Naturkunde Stuttgart, Stuttgart, Germany
Open Access


Bombus flavescens Smith is one of the most widespread bumblebee species in the Oriental region. Due to colour polymorphisms, this species or species-complex has been a challenge for taxonomy. This study aims to assess the taxonomic status of the flavescens-complex using evidence from COI barcodes and morphology. We then reconstruct its biogeographic history from a phylogenetic analysis of populations across the current range, combining COI with 16S and nuclear PEPCK data. Despite a large range of polymorphisms across its distribution, the results show that B. flavescens is a single species based on algorithmic species delimitation methods, and it is clearly separated from its sister species, B. rotundiceps Friese. We suggest that B. flavescens diverged from its sister lineage in the Himalaya and dispersed into Southeast Asia in the Pleistocene. Conservation of the widespread B. flavescens will need to consider its several unique island populations.


COI, Museum specimens, Polymorphism, Pyrobombus


Bumblebees (genus Bombus Latreille) are well-studied pollinators, especially in temperate regions (Williams et al. 2014; Rasmont et al. 2021). They can also be found in both subtropical and tropical areas in Central and South America, and Asia (Williams 1998), although taxonomic knowledge gaps and the relative rarity of bumblebees in these areas still constrain our understanding of their biogeography and ecology.

The taxonomic studies of bumblebees in Asia, based on morphological evidence only, have been considered highly problematic, due to colour pattern polymorphism within the same species (Huang et al. 2015b; Ding et al. 2019) and cryptic species (Williams et al. 2020, 2022a). Moreover, the same colour pattern is also observed from different species which are co-existed in same locality (Hines and Williams 2012). Therefore, these confusions of taxonomic delimitation of Asian bumblebees require molecular evidence (Williams et al. 2012). Taxonomic status of bumblebees has been revised and various new species have been recently described in Asia, based both morphology and molecular evidence (Williams et al. 2020, 2022a, c; Williams 2022).

Bombus flavescens Smith is a Southern Asian species of the subgenus Pyrobombus Dalla Torre (Williams 1998). While the type locality of B. flavescens is Zhoushan, China (Smith 1852), this taxon is found in many other countries in the Himalaya and Southeast Asia with tropical low montane habitats, including Nepal, India, Bhutan, Myanmar, Thailand, Laos, Vietnam, Malaysia, and the Philippines (Starr 1992; Williams 1998, 2022; Williams et al. 2009; Koch and General 2019).

Most frequently, B. flavescens is a predominantly black pubescence (or hair) bumblebee with an orange tail (the area of hairs covering posterior part of metasoma) and legs in workers and queens. Males often (e.g., in parts of China and in the Himalaya) show a predominantly pale yellow (flavescent) hair pattern with orange legs (Fig. 1), although they sometimes exhibit the widespread dark pattern of the females. However, the colour pattern of hair varies across its wide area of distribution (Fig. 2), which has caused taxonomic problems for more than a century and led to the formal description of local variants (Frison 1934). For example, the name B. alienus Smith in China (Fig. 2E; morphologically synonymised by Williams (2022)) was applied to a specimen with an anterior yellow band on its metasoma (terga 1–2 or T1–2) whereas specimens of from the Philippines (Fig. 2G), with yellow hair covering the side of the mesosoma and metasoma (T1–2) but without the orange tail, have been described as B. baguionensis Cockerell (morphologically synonymised by Williams (1998)). In addition, a form with an entirely body, covered with uniformly orange hairs, has been described as B. rufoflavus Pendlebury (Fig. 2J; provisionally synonymised by Williams (1998), based on morphology). Further names were published for colour patterns resembling the typical B. flavescens: B. mearnsi Ashmead from Mindanao in the Philippines (morphologically synonymised by Pittioni (1949)); B. bakeri Cockerell from Negros (Fig. 2H; morphologically synonymised by Frison (1925)); and B. tahanensis Pendlebury from the Tahan mountains of Malaysia (Fig. 2I; morphologically synonymised by Frison (1934)). Apart from this high level of colour variation, other morphological characters, including shape of labrum, punctures on clypeus, and male genitalia, are relatively similar (Williams 1998; Thanoosing 2022). Yet, at various times through the 20th century different intraspecific names have been applied to describe minor B. flavescens variants (Frison 1925, 1928, 1934; Bischoff 1936; Chiu 1948; Pittioni 1949; Table 1) although they are clearly grouped together because of their morphological similarity (Williams 1998).

Figure 1. 

Holotype male of B. flavescens, deposited in the entomological type collection at the Natural History Museum, London (NHMUK 014025381). Scale bar: 5 mm.

Figure 2. 

Colour variation of B. flavescens workers with their unique specimen identifiers A India (FLA#5) B Nepal (NHMUK 010814958) C Thailand (CT#1024) D Kunming, China (FLA#6) E alienus; Fujian, China (NHMUK 010815023) F Taiwan (FLA#4) G baguionensis; Luzon, the Philippines (CT#932) H bakeri; Negros, the Philippines (ZMA.INS.758600) I tahanensis; Gunung Tahan, Malaysia (NHMUK 010820651) J rufoflavus; Cameron highlands, Malaysia (NHMUK 010814574). Some images are reversed. Scale bar: 5 mm.

Table 1.

List of Bombus flavescens names.

Author Name
Smith (1852) Bombus flavescens
Smith (1854) Bombus alienus
Ashmead (1905) Bombus mearnsi
Friese (1905) Bombus rufocaudatus
Cockerell (1917) Bombus geei
Cockerell (1920) Bombus irisanensis var. baguionensis, Bombus bakeri
Pendlebury (1923) Bombus tahanensis, Bombus rufoflavus
Frison (1925) Bremus mearnsi, Bremus mearnsi var. bakeri, Bremus irisanensis var. baguionensis
Hedicke (1926) Bombus imuganensis
Frison (1928) Bremus (Pratobombus) baguionensis, Bremus (Pratobombus) baguionensis var. imuganensis
Dover (1929) Bremus rufoflavus
Skorikov (1933) Pratobombus flavescens
Frison (1934) Bremus mearnsi var. deflectus, Bremus mearnsi var. ditutus, Bremus mearnsi var. bakeri, Bremus mearnsi var. geei
Bischoff (1936) Bombus (Pratobombus) mearnsi ssp. chekiangensis
Chiu (1948) Bombus mearnsi var. deflectus, Bombus mearnsi var. ditutus, Bombus mearnsi var. bakeri, Bombus mearnsi var. geei, Bombus mearnsi var. subrufus, Bombus mearnsi var. luteus
Pittioni (1949) Bombus flavescens f. dilutior

Colour-pattern variation of bumblebees can be controlled by differential gene expression of, e.g., the Abd-B and nubbin genes affecting pigment generation during pupal development (Rahman et al. 2021). However, the evolutionary forces determining variation in colour pattern among B. flavescens populations remain unclear. Selection for mimicry might take a major role in driving evolutionary differences because the colour pattern of B. flavescens locally matches other bumblebee species in the same areas (Williams 2007). For example, on Luzon lsland, baguionensis shows a similar pattern to B. (Megabombus) irisanensis Cockerell, and in the Malay Peninsula, rufoflavus resembles B. (Megabombus) montivagus Smith with a similar predominantly orange pattern of both pubescence and underlying sclerites, which had been named B. maxwelli Pendlebury (Hines and Williams 2012).

Orange body colour (pubescence and sclerites) and enlarged ocelli in Hymenoptera are often associated with a nocturnal or crepuscular lifestyle, for example, Megalopta bees and Apoica wasps in Central and South America (Roubik 1989; Williams 2007; Warrant 2008). Interestingly, male and female nocturnal carpenter bees, Xylocopa (Nyctomelitta) myops Ritsema, show a pattern closely similar to the orange pattern bumblebees. This carpenter bee is found in Malaysia, Singapore, and Borneo, based on specimens in the Natural History Museum, London (NHMUK) collection and records by Ascher et al. (2022). Similarly, the two orange pattern bumblebees in Malaysia, rufoflavus (B. flavescens) and maxwelli (B. montivagus), might be active nocturnally also, at least in part (Williams 2007). Species of the subgenus Nyctomelitta Cockerell display relatively large ocelli compared to day-flying carpenter bees (Cockerell 1929; Michener 2007). However, the ocelli of the orange pattern bumblebees remain to be assessed in this regard. Observation of a night-flying bumblebee had been claimed in Myanmar (Doria 1886), but this record is a misidentification of the carpenter bee, X. (Nyctomelitta) tranquebarica (Fabricius) (Cameron 1910). Therefore, genuine evidence of nocturnal bumblebees has yet to be recorded in Southeast Asia.

The colour-pattern diversity of B. flavescens must be seen in light of the complex biogeography of its distribution range across the Southeast Asian mainland and the Philippines. The most recent common ancestor (MRCA) of the flavescens species complex has been placed into the Miocene epoch around eight million years ago (Ma), presumably on the mainland (Hines 2008) with two possible dispersal routes to the Philippines (Starr 1989): 1) through the Sundaland route, via Borneo–Sulawesi from the South; or 2) across the Luzon strait via Taiwan from the North. When global temperatures increased after the last glacial maximum (LGM), bumblebees likely became restricted to highlands (e.g., the Cameron Highlands, Malaysia) and diverged in allopatry. At the same time, as the sea level rose and submerged land bridges across Sundaland (Voris 2000; Sathiamurthy and Voris 2006; Woodruff 2010), populations of bumblebees on the Philippines were even further isolated.

This study seeks to clarify genetic patterns within the flavescens-complex, to establish the species status of geographically separated lineages and to infer biogeographic scenarios for bumblebees in Southeast Asia more generally. Species-level entities in bumblebees have been delimited by morphological, molecular (e.g., Williams et al. (2019, 2020)), and chemical approaches, the latter using Cephalic Labial Gland Secretions or CLGS (e.g., Brasero et al. (2021), Ghisbain et al. (2021)). Given the challenges of morphology for taxonomic assessment of the flavescens-complex, molecular data are required to resolve the species status of sub-lineages and their relationships. In addition, this study aims to clarify the question about the nocturnal lifestyle in B. flavescens, using morphological examination of their ocelli.



Museum and institutional specimens of the flavescens-complex (Table 1) were examined for morphological analysis. Collecting information was recorded and made available at the NHMUK Data Portal ( and Additional records were obtained from reliable sources, including: 1) published literature (Williams et al. 2009, 2010); 2) major natural history collections, including the Naturalis Biodiversity Center (NMNL; Bakker F, Creuwels J), the Smithsonian Institution (USNM; Orrell T, Informatics Office), the University of Illinois at Urbana-Champaign (INHS; McElrath T), and the Taiwan Forestry Research Institute (TFRI; Lu S), available on the Global Biodiversity Information Facility (GBIF) ( 2021). Records with unclear locality information and clear geographic outliers were ignored. For example, there is a record of B. flavescens from Japan (e.g., GBIF occurrence 3801818589; McElrath 2022), far outside the established range of B. flavescens, based on the specimen records by Williams (1998, 2022).

COI barcode data of B. flavescens and relatives were obtained from public databases (Suppl. material 1) and newly sequenced from pinned or fresh specimens (Table 2). The samples were selected to represent the geographical distribution and all major colour pattern variants of B. flavescens (Fig. 3). Bombus flavescens belongs to the pratorum-group (Williams 1998). Thus, closely related species in the pratorum-group (Cameron et al. 2007; Williams et al. 2009, 2010), were included in the dataset: B. ardens Smith, B. pratorum (Linnaeus), B. pyrenaeus Pérez, B. modestus Eversmann, B. nursei Friese, B. biroi Vogt, B. wangae Williams et al., and B. rotundiceps Friese. The status of B. nursei as a distinct species was suggested in Williams (2022). Thirty COI sequences (> 600 bp) of B. flavescens and the relatives were downloaded from the Barcode of Life Data System (BOLD; Ratnasingham and Hebert 2007) and the National Center for Biotechnology Information (NCBI) database (GenBank; Sayers et al. 2022) (Suppl. material 1). Nine B. flavescens sequences from China were provided by co-author MO, sequenced from freshly collected specimens. Thirteen further B. flavescens and five B. rotundiceps specimens from museum and institution collections were selected for sequencing (Table 2). The identifier numbers were applied to the studied specimens: 1) “CT#n” for Southeast Asian bumblebee specimens 2) “FLA#n” and “ROT#n” refer to B. flavescens and B. rotundiceps specimens from outside Southeast Asia, respectively.

Figure 3. 

Map showing the distribution of B. flavescens with the pattern of hair colour in the dorsal view. The blue dots represent the records from museum collections, literature, and GBIF database (n = 702). The red dots represent the localities of barcoded specimens.

Table 2.

List of specimens included, species name, ID, deposit place (KKIC = Kasetsart Kamphaeng Sean Insect Collection, PHW = Paul H. Williams Research Collection, CUNHM = Chulalongkorn University Natural History Museum Collection, NMNL = Naturalis Biodiversity Center Collection), sex/caste (w = worker, m = male), locality, collecting date, and COI GenBank accession number.

Species Project ID Collection (specimen ID) Sex/caste Locality Collecting date GenBank accession
B. flavescens CT#662 KKIC w Thailand, Nakorn Pathom? 7/5/2017 OP355718
CT#669 KKIC m Thailand, Loei 13/4/2016 OP355719
CT#926 PHW w Malaysia, Gunung Tahan NA OP355720
CT#1018 NMNL (ZMA.INS.758598) w Philippines, Negros 4/5/1953 OP355722
CT#1023 CUNHM (BSRU-AB-9609) w Thailand, Chiang Mai 15/2/2021 OP355723
CT#1024 CUNHM (BSRU-AB-9610) w Thailand, Chiang Mai 15/2/2021 OP355724
FLA#2 PHW w Philippines, Luzon 27/4/1986 OP355725
FLA#3 PHW w Taiwan, Tai Chung 5/5/1980 OP355726
FLA#4 PHW w Taiwan, Nantou 24/6/1989 OP355727
FLA#5 PHW w India, Uttar Pradesh 6/5/1990 OP355728
FLA#6 PHW w China, Kunming 8/4/2018 OP355729
FLA#7 NMNL (ZMA.INS.773029) w Bhutan, Lungtenphu 10/6/1996 OP355730
FLA#10 PHW w Philippines, Luzon 27/4/1986 OP355731
B. rotundiceps CT#960 CUNHM (BSRU-AB-1268) m Thailand, Chiang Mai 17/6/2019 OP355721
ROT#1 PHW w India, Uttar Pradash 6/5/1990 OP355732
ROT#2 PHW w India, Uttar Pradash 6/5/1990 OP355733
ROT#3 PHW w China, Guangxi 3/6/2016 OP355734
ROT#4 PHW w China, Guangxi 5/6/2016 OP355735

DNA extraction

Pinned and fresh specimens of B. flavescens and B. rotundiceps were selected for DNA extraction. For the recent specimens (< 20 years old), the tissue sources were a right front leg or a whole body for high DNA concentration yield. Each leg sample was ground using a small pestle in a microcentrifuge tube, whereas the whole-body sample was directly put in a microcentrifuge tube without specimen damages. The tissue samples were incubated at 56 °C for 24 hours in the ATL buffer with Proteinase K enzyme. Genomic DNA was extracted using the Qiagen DNeasy Blood and Tissue Kit, following the kit protocol. DNA quality and quantity were assessed using the Nanodrop spectrophotometer (Thermo Scientific ND8000) and the Agilent 2200 TapeStation system (Agilent Technologies, Inc.).

For older specimens (> = 20 years old), genomic DNA was likely to be degraded and all DNA extraction steps were performed inside a laminar flow chamber using dedicated UV sterilised equipment and consumables different from those used with modern bumblebee DNA. The DNA extraction method followed Sproul and Maddison (2017) using the Qiagen QIAamp Micro Kit, but without the use of carrier RNA. For each specimen, the right front leg was cleaned with sterile water, frozen with liquid nitrogen, and then ground with a small pestle. The samples were incubated in ATL buffer with Proteinase K at 55 °C for 24 hours. A blank sample (only the ATL buffer with the enzyme) was included for each extraction batch as a contamination control. Samples were eluted from the column twice with 15 µl AE buffer and incubated at room temperature for 10 minutes. The DNA concentration was measured using the Qubit fluorometer high sensitivity (Invitrogen) and samples were stored at -20 °C.

Primer design, amplification, and sequencing

For the recent specimens, PCR was performed to amplify the full length of COI amplicons (658 bp) using primers LepF1 and LepR1 (Hebert et al. 2004) and an annealing temperature at 50 °C. For older specimens, due to the DNA degradation, three pairs of primers were newly designed for amplifying partial COI contigs (150–300 bp), using Primer3 software (Untergasser et al. 2012). The COI sequence of B. flavescens (GenBank accession number: GU085209) was used as a reference sequence for the first two pairs of primers. The third pair was designed based on the partial mitogenome of B. pratorum (GenBank accession number: KT164684). There are three pairs of new primers for partial COI amplification (contig) in this study (Table 3): 1) FLA1: ParCOI_246_FLA_F1 (Tm° = 66.2 °C) and ParCOI_416_FLA_R1 (Tm° = 63.7 °C), 2) FLA2: ParCOI_349_FLA_F2 (Tm° = 64.2 °C) and ParCOI_662_FLA_R2 (Tm° = 66.6 °C), and 3) PRA1: ParCOI_1816_PRA_F1 (Tm° = 61.7 °C) and ParCOI_1994_PRA_R1 (Tm° = 65.5 °C). These primers were tested for amplification efficiency (Suppl. materials 24). The annealing temperature for the pair FLA1 and FLA2 was set at 65 °C and PRA1 was set at 63 °C.

Table 3.

Newly designed primers were used in this study, including gene, primer name, strand (F = forward, R = reverse) and primer sequence.

Gene Primer name Strand Primer sequence

The PCR products were sequenced in both forward and reverse directions using ABI technology at the NHMUK’s sequencing facility. The sequences were edited using MEGA version 7.0.26 (Kumar et al. 2016:, to trim low-quality bases near the ends of the traces. The COI barcodes have been deposited in GenBank (Table 2).

Phylogenetic analysis and species delimitation

The aligned dataset was analysed to determine the best-fitting nucleotide substitution model using jModelTest software version 2.1.6 (Darriba et al. 2012). Phylogenetic trees were constructed using Bayesian Inference (BI) with MrBayes version 3.2.2 (Ronquist and Huelsenback 2003). Markov Chain Monte Carlo (MCMC) chains were run for 10 million generations and sampled every 1000th generation. The first 10% was discarded as burn-in. The phylogenetic trees were visualised using FigTree version 1.4.3 ( Bumblebee species from closely related subgenera (Cameron et al. 2007) were used as outgroups, including B. (Bombus) terrestris (Linnaeus) and B. (Alpinobombus) polaris Curtis. The outgroup sequences were generated by Thanoosing (2017).

Poisson Tree Process (PTP) analysis is used to identify likely species coalescents. The test establishes the transition from long branches defining between-species diversification to short branches within species (estimated from the number of substitutions on branches) (Zhang et al. 2013; Kapli et al. 2017). The proportion of between- and within-species sampling was in a range where the PTP analysis performs well (Luo et al. 2018). The PTP analysis was run on the bPTP server (; accessed 2021) under default parameters. The analysis was restricted to unique haplotype sequences in order to avoid introducing misleading ‘zero-length’ branches. The input phylogenetic tree was estimated using MrBayes rooted with Bombus terrestris. Results from the PTP analysis were compared to the Generalised Mixed Yule Coalescent (GMYC), which uses the branching rate rather than the number of nucleotide changes to establish the species-to-population transition (Pons et al. 2006). The required ultrametric tree was obtained with BEAST (Bouckaert et al. 2019) using the same dataset as for the PTP analysis. The clock rate was set at 0.0177 using the divergence rate from tenebrionid beetle COI (Papadopoulou et al. 2010). The GMYC analysis was run in the splits package of R (Ezard et al. 2009;

Biogeographic analysis

To reconstruct phylogenetic trees for biogeographic scenarios of B. flavescens, more genetic markers were added to the dataset, including mitochondrial 16S rDNA (16S) and nuclear phosphoenolpyruvate carboxykinase (PEPCK). Due to a lack of fresh specimens, most of 16S and PEPCK sequences for relatives of B. flavescens were retrieved from GenBank, especially from the dataset compiled by Cameron et al. (2007) (Suppl. material 5). 16S and PEPCK were not available from existing sources for three species; B. rotundiceps, B. wangae, and B. nursei. Bombus nursei is rare in inaccessible areas of the western Himalaya (Williams, 1991, 2022), and B. wangae also is a rare bumblebee from Sichuan (Williams et al. 2009). For B. rotundiceps and additional B. flavescens), the 16S and PEPCK were amplified using primers 16SWb (Dowton and Austin 1994) /874–16SIR (Cameron et al. 1992), and FHv4/RHv4 (Cameron et al. 2007) and sequenced in both directions (Table 4).

Table 4.

List of 16S and PEPCK sequences, included species, ID, and the 16S and PEPCK GenBank accession numbers, generated in this study.

Taxa Project ID 16S PEPCK
B. flavescens CT#926 OP354521 -
CT#1023 OP354523 OP382631
FLA#6 OP354524 OP382632
B. rotundiceps CT#960 OP354522 OP382630

An ultrametric tree was constructed based on the combined dataset of COI (unique haplotypes), 16S and PEPCK (Suppl. material 6), using *BEAST (Ogilvie et al. 2017). If not available from the same individual, sequences in these three datasets were concatenated for different representative from the same area of distribution. For PEPCK, the sequence was partitioned for introns and exons. The closely related B. (Pyrobombus) lepidus Skorikov, based on Cameron et al.’s (2007) tree, was used as an outgroup. The nucleotide substitution model for each fragment was selected according to the BIC criterion. The BEAUti software (Bouckaert et al. 2019) was used to generate an input file for BEAST. Trees were estimated under a strict clock and a Calibrated Yule Model as a species tree prior with the birth rate prior drawn from a Gamma distribution. Due to a lack of fossil evidence for the subgenus Pyrobombus, the clock was calibrated from molecular rate estimates based on five genes by Hines (2008), placing the origin to approximately 8.5 +/- 1 Ma. Accordingly, the prior distribution was 10.1–6.86 Ma (5% tails). The MCMC algorithm was run for 500 million generations and sampled every 5000th generation. The trace file was visualilsed using Tracer (Rambaut et al. 2018). A 20% burn-in was discarded using TreeAnnotator (Bouckaert et al. 2019) to construct the maximum clade credibility tree.

S-DIVA analysis (Yu et al. 2010) and BioGEOBEARS with DIVALIKE+J model (Matzke 2013a, 2013b, 2014) in the RASP package (Yu et al. 2020) were used to estimate biogeographic scenarios for the flavescens-complex and their relatives. The biogeographical distributions were divided into 11 principal areas, based on bumblebee global biogeography (Williams 1996), areas of endemism in Southeast Asia and the distributions of flavescens-complex and their relatives (Table 5): 1) Europe (A), 2) Central Asia (B), 3) North and Central China and North Asia (C), 4) Korean peninsula and Japan (D), 5) South China (E), 6) Himalaya (F), 7) Thailand (G), 8) Peninsular Malaysia (H), 9) Taiwan (I), 10) Luzon (J), and 11) Negros (K). These areas were chosen so that the maximum range size for any species was limited to three of these areas, due to acknowledged problems with the methods in inheriting broader distributions (Lamm and Redelings 2009). Then, the dispersal-corridor model was created with possible bumblebee permitted dispersal routes (Fig. 4). According to this model, the possible dispersal areas in this study are AB, ABC, ABF, AC, ACD, ACE, BEF, BF, BFG, CD, CDE, CE, CEF, CEG, CEI, EF, EFG, EFI, EG, EGH, EGI, EI, EIJ, FG, FGH, GH, GHK, HJK, HK, IJ, IJK, and JK. The burnin was set at 20%. Then, the *BEAST trees were used as the input trees.

Figure 4. 

A corridor-dispersal model diagram of the flavescens-complex and their relatives in this study.

Table 5.

Principal areas of endemic distribution of the flavescens-complex and their relatives in this study.

Area Principal ranges included Species recorded
A Europe B. pratorum, B. pyrenaeus
B Central Asia B. biroi
C North China, and North Asia B. modestus
D Korea peninsula and Japan B. ardens, B. modestus
E South China B. flavescens, B. lepidus, B. rotundiceps, B. wangae
F Himalaya B. flavescens, B. lepidus, B. nursei, B. rotundiceps
G Thailand B. flavescens, B. rotundiceps
H Peninsular Malaysia rufoflavus (B. flavescens), tahanensis (B. flavescens)
I Taiwan B. flavescens
J Luzon baguionensis (B. flavescens)
K Negros bakeri (B. flavescens)

Morphological study of ocelli in the flavescens-complex

The idea of the nocturnal lifestyle in the orange pattern bumblebees, rufoflavus (B. flavescens) and maxwelli (B. montivagus) from Peninsular Malaysia, has been introduced by Williams (2007). Although orange hair pattern can be recognised in diurnal bumblebees (e.g., B. humillis Illiger, B. pascuorum (Scopoli), and B. muscorum (Linnaeus)), rufoflavus and maxwelli show the uniformly orange pattern both of pubescence and sclerites uniquely. This orange pattern can be recognised in the female specimens of nocturnal carpenter bees, X. myops, mentioned in Williams (2007), which can be also found in Peninsular Malaysia. For this reason, we selected the female specimens of X. myops and their related species, X. tranquebarica from the NHMUK collection, included in this study to test the idea.

Due to the sexual dimorphism bias and the availability of male specimens, only female or worker specimens were chosen in this study. Female flavescens-complex specimens (n = 107) together with its relative, B. rotundiceps (n = 12), outgroups including the diurnal B. irisanensis (n = 4), the orange pattern B. montivagus (taxon maxwelli) (n = 3), and the nocturnal carpenter bees, X. myops (n = 20) and X. tranquebarica (n = 12), were selected for morphological character measurement, including median ocellus width (MOW), intertegular distance (ID), and head width (HW) (Fig. 5), under a Wild Heerbrugg M4A stereomicroscope, calibrated with ocular and stage micrometres. The specimens of the flavescens-complex were divided into eight groups: 1) Thailand (n = 4), 2) baguionensis (n = 15), 3) Himalaya (n = 16), 4) China (n = 20), 5) Taiwan (n = 6), 6) bakeri (n = 6), 7) rufoflavus (n = 20), and 8) tahanensis (n = 20). The data are available in the Suppl. material 8.

Figure 5. 

The measurement of morphological characters A dorsal aspect: intertegular distance (ID) B frontal aspect: median ocellus width (MOW) C frontal aspect: head width (HW). The specimen shown is B. flavescens, NHMUK 010814574, from the Cameron Highlands, Malaysia.

The morphological values, including MOW, the ratio between MOW and ID (MOW:ID), and the ratio between MOW and HW (MOW:HW), were visualized in box plots using R package ggplot2 (Wichham 2016). To test for differences between species, and within the flavescens-complex, statistical tests were performed in R (R core team 2021). First, a Shapiro-Wilk normality test was used to examine the distribution of data of each of the flavescens-complex taxa and each bee species. Where the data showed a normal distribution, an ANOVA test and Tukey test were then performed. A Kruskal-Wallis test and Dunn’s test were used instead if the data were not normally distributed. The Dunn’s test was conducted using the R package FSA (Ogle et al. 2022).


DNA extraction, amplification, and sequencing

Genomic DNA was extracted successfully for most specimens of the flavescens-complex and B. rotundiceps. Most of the old samples were degraded and had relatively low DNA concentrations (< 10 ng/µl). The three amplicons obtained with new primer pairs PRA1, FLA1, and FLA2 were assembled into a single contig. However, there was a 25-base pairs gap between the first (PRA1) and second (FLA1) fragments, and the 5’ part of the COI barcode region (88 base pairs) was not amplified at all. COI sequences were generated for 18 collected specimens with top BLAST hits to the subgenus Pyrobombus.

Phylogenetic analysis and species delimitation

The final dataset of the flavescens-complex and relatives comprised 57 COI sequences, including 25 unique haplotypes. The resulting 575 bp alignment was used to construct the BI phylogenic tree under the GTR+Γ model selected by jModelTest (BIC = 4596) (Fig. 6). Bombus rotundiceps, sampled from China, Nepal, India, and Thailand, was the reciprocally monophyletic sister group of the flavescens-complex. The latter included four regional clades: 1) Thailand and Himalaya, 2) China, 3) Malaysia, and 4) the Philippines and Taiwan.

Figure 6. 

Bayesian inference (BI) phylogenetic tree based on a 575 bp fragment of COI in MrBayes. MCMC chains were run for 10 million generations, sampled every 1000th generation. Burn-in fraction is 10%. The posterior probabilities are shown above the branches. The tip of the tree shows the sample label including the sequence length, a taxon name, an identifier code from the database, and its geographic origin. The results of species delimitation methods are shown in the right-hand side columns. The black and grey bars within the same columns indicate the same species.

The unique haplotype COI dataset was used to estimate the input tree for the species-delimitation analysis. The nucleotide substitution model of this dataset was GTR+Γ (BIC = 3877). Species delimitation using PTP returned each of the nine species as a separate entity, with Bayesian support values between 0.14–1.00 (Fig. 6). In contrast, the GMYC lumped the flavescens-complex and B. rotundiceps, as well as the B. modestus and B. wangae species pair (Fig. 6). All analyses recovered the flavescens-group, including B. flavescens s. str., B. baguionensis, B. bakeri, B. rufoflavus, and B. tahanensis, as a single species.

Biogeographic analysis

The phylogenetic analysis was based on three loci (see supplement for missing data in 16S and PEPCK; Suppl. material 6) and 26 terminals, including the outgroup. Tree searches were conducted under partitioning and separate model choice for COI (575 bp), 16S (551 bp), PEPCK introns (483 bp), and PEPCK exons (369 bp). The resulting maximum clade credibility tree showed a deep separation of the B. flavescens/B. rotundiceps clade from their closest relatives, B. pratorum and B. ardens, estimated at 6 +/- 2 Ma (Fig. 7). Within the flavescens-complex the early split in was occupied by the Malaysian group. The extant members of the flavescens-complex diversified during a time interval approximately 1–2 Ma. However, there is some uncertainty about relationships among the four regional groups, as in the COI tree the Thai and Himalayan groups were sister to the China, Malaysia, the Philippines and Taiwan groups.

Figure 7. 

The maximum clade credibility tree of B. flavescens and their closely related species, reconstructed with *BEAST from gene trees for the COI, 16S, and PEPCK genes and calibrated using estimate time from Hines (2008). The MCMC chains were run for 500 million generations, sampled every 5000th generation. Burn-in fraction is 20%. The posterior probabilities are shown under the nodes. Blue bars represent the 95% confidence limits on the estimate dates of divergence. Bombus lepidus is the outgroup.

Although DEC model was the best fit model for the BioGEOBEARS analysis in this study (the highest AICc wt; Suppl. material 7), we selected DIVALIKE+J model instead (the second high AICc wt; Suppl. material 7). The DIVALIKE+J model supports widespread founder-event speciation at nodes (Yu et al. 2020) which fits the distribution changes of bumblebees (Williams et al. 2020, 2022b).

Results from S-DIVA and DIVALIKE+J showed that there was uncertainty in ancestral areas of B. flavescens and relatives, with more than one possible ancestral area suggested in both scenarios (Figs 8, 9). There were only four nodes (III, XI, XIV, and XVI; Fig. 8) for which we identified a single ancestral area in S-DIVA, whereas the DIVALIKE+J result (Fig. 9), only node XI showed a single ancestral area. The results suggested that the ancestral area for the MRCA of the pratorum-group was likely to be in area covered by Europe, Central Asia, and the Himalaya (ABF) in S-DIVA, and Central Asia (B) in DIVALIKE+J with low probability. This group diversified in the late Miocene.

Figure 8. 

Biogeographic scenarios for the flavescens-complex and its close relatives by S-DIVA analysis. The input tree is reconstructed with *BEAST from gene trees for the COI, 16S, and PEPCK genes and calibrated using estimate time from Hines (2008). The internal nodes are represented in Roman numerals. The colours of the nodes represent the possible ancestral areas. The area codes used are given in Table 5. The most likely ancestral areas are represented above the nodes with the probabilities below.

Figure 9. 

Biogeographic scenarios for the flavescens-complex and its close relatives by BioGeoBEARS with DIVALIKE+J analysis. The input tree is reconstructed with *BEAST from gene trees for the COI, 16S, and PEPCK genes and calibrated using estimate time from Hines (2008). The internal nodes are represented in Roman numerals. The colours of the nodes represent the possible ancestral areas. The area codes used are given in Table 5. The most likely ancestral areas are represented above the nodes with the probabilities below.

For the flavescens-complex, the S-DIVA biogeographic scenario illustrated that the MRCA of B. rotundiceps and B. flavescens (Node IX; Fig. 8) might have occurred in the Himalaya (F) around 1.5 Ma during the Pleistocene epoch. Next, the MRCA of B. flavescens in the S-DIVA analysis was inferred to be in the FGH areas, including Himalaya, Thailand, and Peninsular Malaysia (Node X; Fig. 8) at around 1 Ma. However, the DIVALIKE+J biogeographic scenario suggested that the group might have originated in Malaysia (H), but with low probability. The DIVALIKE+J showed higher support for a Peninsular Malaysian origin (Node X; Fig. 9).

The diversification of the B. flavescens populations began after 1 Ma, resulting in three distinct clades: 1) rufoflavus+tahanensis from Malaysia (Node XI, Figs 8, 9); 2) flavescens s. str. from Thailand, Himalaya and China (Node XIII, Figs 8, 9); 3) baguionensis+bakeri+flavescens s. str. from Taiwan (Node XV, Figs 8, 9). According to our corridor-dispersal model diagram (Fig. 4), the S-DIVA showed that the B. flavescens used all possible corridors between E–K, except the corridor H–K, the Sundaland route, whereas the DIVALIKE+J included the corridor H–K as a permitted route.

Morphological study of Bombus flavescens ocelli

The carpenter bees, Xylocopa tranquebarica and X. myops, had a significantly larger median ocellus than the bumblebees (Fig. 10). Bombus irisanensis clearly showed the smallest median ocellus and B. flavescens, B. rotundiceps, and B. montivagus exhibited similar median ocellus sizes (Fig. 10).

Figure 10. 

Box plots showing the median ocellus width (MOW), the ratio between MOW and intertegular distance (MOW:ID), the ratio between MOW and the head width (MOW:HW) of B. irisanensis, B. rotundiceps, B. flavescens, B. montivagus taxon maxwelli, Xylocopa tranquebarica, and X. myops.

For the flavescens-complex (Fig. 11), the Shapiro-Wilk normality test suggested that the MOW was not normally distributed (W = 0.93, p-value < 0.05), whereas the MOW: ID and MOW: HW did not differ significantly from a normal distribution (W = 0.99, p-value = 0.28 and W = 0.99, p-value = 0.45 respectively). There was a significant difference within the MOW of the flavescens-complex (Kruskal-Wallis test, chi-squared = 36.217, p-value < 0.05; Dunn’s test, p-value < 0.05), including the taxa baguionensis-rufoflavus, baguionensis-tahanensis, Himalaya-rufoflavus, and Himalaya-tahanensis. For the MOW: ID, nine pairs were significantly different (ANOVA, F-value = 6.781, p-value < 0.05; Tukey test, p-value < 0.05): bakeri-China, rufoflavus-China, tahanensis-China, bakeri-Himalaya , rufoflavus-Himalaya, tahanensis-Himalaya, Taiwan-bakeri, Taiwan-rufoflavus, and Taiwan-tahanensis. In addition, only two pairs of taxa were significantly different in MOW:HW (ANOVA, F-value = 3.701, p-value < 0.05; Tukey test, p-value < 0.05), rufoflavus-Himalaya, and tahanensis-Himalaya.

Figure 11. 

The box plots of the median ocellus width (MOW), the ratio between MOW and intertegular distance (MOW:ID) and the ratio between MOW and head width (MOW:HW) among Bombus flavescens population: Thailand, China, Himalaya, Taiwan, baguionensis, bakeri, rufoflavus, tahanensis.


Species status of the flavescens-complex

The status of species within the flavescens-complex taxa has been debated for nearly a century. Theodore Frison (1895–1945), an American entomologist wrote in his work in 1934, “This species of bumblebee [B. flavescens] has been a problem to most persons who have attempted to determine…” (Frison 1934). Our results for the COI-barcode tree (Fig. 6) support that the flavescens-complex is a monophyletic group, including the populations in China, the Himalaya, and Southeast Asia. The PTP and GMYC results group the flavescens-complex together as one species. This demonstrates that B. flavescens s. l. includes high intraspecific variation in colour pattern across its distributional range. The species-delimitation analyses confirm the conspecific status of the various geographical and morphological variants named as B. baguionensis, B. bakeri, B. tahanensis and B. rufoflavus by previous authors. Although this study did not include the taxon mearnsi from Mindanao, the taxon mearnsi also is expected to be a synonym of B. flavescens as morphologically suggested by Pittioni (1949), similar to the taxon bakeri from Negros, another population from the southern Philippines islands.

There is a discrepancy between the number of species recognised by PTP and GMYC. The GMYC analysis lumped 1) B. flavescens and B. rotundiceps 2) B. modestus and B. wangae, as single species, whereas PTP split them into four species. GMYC requires an ultrametric tree, which distorts the tree and branch lengths (Fujisawa and Barraclough 2013). For this, the closely related sister taxa which establish their own short clade might be identified as the same species. In contrast, the PTP keeps the original tree shape, and should be considered more reliable between the two approaches (Williams 2021), in this instance.

Nevertheless, when we investigate the morphological evidence of these four taxa, their morphological characters are unique (Williams 1998; Williams et al. 2009, 2010). The males of B. wangae and B. modestus are distinct, especially in their genitalia, which the gonostylus of B. modestus is broader than B. wangae. For the B. flavescens and B. rotundiceps, their male genitalia are also different: the gonostylus of B. flavescens is broad with round inner margin, whereas the gonostylus of B. rotundiceps is narrower with straight inner margin. In this case, the PTP results show congruence with the morphological characters whereas GMYC does not. Then, we suggest that the PTP analysis corroborates more consistently morphologically recognisable species for these bumblebees. For further study, if fresh male specimens of B. flavescens are available, chemical evidence, including, CLGS, is recommended as an alternative approach for defining the species recognition, whether it shows agreement with the molecular delimitation in this study.

Fresh specimens of B. flavescens and their relatives are not available from enough samples because of the rarity of these bees. Bumblebees in subgenus Pyrobombus are active particularly early in the year (Williams et al. 2014), so subgenus Pyrobombus sampling can only be conducted in a relatively short period. The accessibility of their habitats in Southeast Asia is also relatively limited and no records have been reported during recent decades, especially for B. flavescens from Gunung Tahan, Malaysia, and from Negros and Mindanao in the Philippines. Moreover, although the flavescens-complex is widely distributed, recorded in at least eleven countries, it is difficult to obtain specimens from some countries, because of limits imposed by local or international regulations. Therefore, old museum specimens are the best available option for clarifying genetic relationships within this group. However, these specimens were collected between 1980 and late 1990. The DNA of those specimens has become degraded through time (Pääbo 1989). It is difficult to extract the DNA from old specimens with standard protocols. However, numerous DNA extraction techniques for old insect specimens have been developed (Gilbert et al. 2007; Thomsen et al. 2009; Ballare et al. 2019). Next generation sequencing would help to obtain DNA sequences from historical specimens (Nazari et al. 2016; Prosser et al. 2016; Sproul and Maddison 2017; Call et al. 2021), but a high DNA quantity is required (Goodwin et al. 2016). According to our results, the DNA concentration of the samples in this study was not high enough for the sequencing requirement. The newly designed primers in this study were successful for the flavescens-complex and B. rotundiceps. The primers can amplify the DNA from specimens up to 70 years old. From our results, fresh specimens or pinned specimens less than 20 years old are recommended for use as DNA sources as the first option. Museum specimens are an alternative way if fresh specimens are not available. However, using museum specimens requires additional molecular processes, for example, extra-sterilisation of equipment and laboratory space, and conducting negative extraction and PCR controls. This is the first time that at least the partial COI barcodes of the flavescens-complex populations from the Philippines from both North and South islands have been retrieved. This information is invaluable because this is the key to resolving cryptic status among this group, but it also provides a workflow for approaching other difficult bee groups for which only old specimens are available.

Orange pattern bumblebees in Peninsular Malaysia

In the Cameron Highlands, Malaysia, several orange pattern bumblebees have been recorded, including the rufoflavus form of B. flavescens (Fig. 2J) and the maxwelli form of B. montivagus. Their fulvous hair and sclerite colour might relate to the nocturnal lifestyle of hymenopterans, similar to the nocturnal carpenter bees in subgenus Nyctomelitta (Williams 2007). In this study, the ocelli of B. flavescens populations were similar in size (Fig. 11) which did not match the hypothesis of enlarged ocelli, whereas the nocturnal carpenter bees have much more distinctively enlarged ocelli (Fig. 10). However, in the context of the variation among B. flavescens populations, the median ocellus of both Malaysian populations, B. flavescens taxon rufoflavus and taxon tahanensis, were significantly larger than the populations at higher latitudes, including from Himalaya and Luzon Island (Fig. 11). The results showed a gradient of ocellar size within the B. flavescens populations from high latitude (smaller ocellus) to low latitude (larger ocellus): MOW, MOWID, and MOWHW differed significantly between low latitude (0°–20°N) and high latitude (>20°N) group (MOW: Mann-Whitney U test, p-value = 0.037; MOWID: ANOVA, F-value = 30.49, p-value < 0.05; MOWHW: ANOVA, F-value = 8.573, p-value < 0.05).

Although there is a significant difference between low-latitude and high-latitude populations of B. flavescens, intraspecific variation of ocellar size in bumblebee populations has been reported before for B. terrestris (Kapustjanskij et al. 2007), without any evidence for nocturnal behavior in this well-known species. This also suggests that ocellar size might not be a good diagnostic character to distinguish species in bumblebees; it may be that it is relatively plastic to local selection regimes. The ocelli are visual organs that detect polarized light in low-light conditions (Wellington 1974; Roubik 1989). Ocellar size shows a negative correlation with light intensity (Kerfoot 1967; Warrant et al. 2006; Kapustjanskij et al. 2007), and there are a high number of nocturnal and crepuscular bees in tropical or low latitude areas (Dorey et al. 2020). The trend seen here towards larger ocellar size in the tropics might be evidence that bumblebees are more active in low-light conditions, for example, in the early morning or late afternoon, to avoid heat stress, as is suggested for B. breviceps Smith and B. haemorrhoidalis Smith (Williams 1991; Thanoosing 2022). Larger ocelli can also be observed in other tropical forest bees, for example, stingless bees, because sunlight is filtered by the tree canopy, thereby reduced in the understory (Streinzer et al. 2016).

Numerous light traps were run at night on the Cameron Highlands recently by researchers (Musthafa et al. 2021) and by tourist services (e.g., Cameron Service:, but no bumblebees have been observed at traps. Consequently, a nocturnal or crepuscular lifestyle of the orange pattern bumblebees in the Cameron Highlands remains unsupported.

Apart from orange B. flavescens, another bumblebee, B. montivagus taxon maxwelli, and a variation of the diurnal hornet Vespa velutina Lepeletier taxon divergens, in the same locality, are also orange in colour (Hines et al. 2012; Perrard et al. 2014). In addition, not only the orange pattern but also the black hair with a red tail pattern of B. flavescens in South China (Fig. 2D) is recognised in the colour variation of V. velutina taxon nigrithorax (Perrard et al. 2014). For this reason, the colour patterns of B. flavescens might potentially be a result of Müllerian mimicry.

Biogeography of Bombus flavescens

Pyrobombus is a bumblebee subgenus in which five species groups are found in the Old and New World (Williams 1998). Bombus flavescens and its relatives are in the ‘pratorum-group’. This group is the sister group to the ‘lepidus-group’ (Cameron et al. 2007). The groups originated in the Palearctic and Oriental regions, in the late Miocene (~11 Ma), and the lepidus-group is entirely restricted to the Oriental region (Williams 1998). In this study, biogeographic analyses show that the origin of the pratorum-group is likely to be in the Oriental region, specifically around the Qinghai-Tibetan Plateau (QTP) or Himalaya. Then, during the late Miocene and early Pliocene, the diversification of this group occurred. This is coincident with the period when the monsoon climate in this area intensified (Zhou et al. 2018), which also accelerated the diversification of bumblebee food plants in the Oriental region (Yu et al. 2015; Ma et al. 2016; Matuszak et al. 2016). Our results show that there were dispersal events between the Oriental (E, F, G, H, I, J, and K; Fig. 8) and the Palearctic region (A, B, C, and D; Fig. 8), for example, the clade of biroi+nursei, pyrenaeus+modestus+wangae, and pratorum+ardens. In addition, there are the lineages of B. nursei and B. wangae that reinvaded the Oriental Region: B. nursei from the west and B. wangae from the east. Bombus nursei prefers subalpine meadow habitats similar to B. biroi (Williams 1991; Williams 2022). This dispersal pattern between the Oriental and Palearctic regions during the Miocene-Pliocene can be observed in other bumblebee subgenera, for example, in subgenera Mendacibombus (Williams et al. 2018) and Melanobombus (Williams et al. 2020).

Only the group rotundiceps+flavescens has been restricted to the Oriental until now. The divergence between B. rotundiceps and B. flavescens was in the Pleistocene (ca 1.5 Ma). The climate in Southeast Asia during the late Pliocene and Pleistocene was cooler than the present day (Morley 2012). This might have given the opportunity for temperate and montane flora to colonise these areas. At the same time, the MRCA of B. flavescens dispersed southward into Southeast Asia or Sundaland, following suitable habitat (A; Fig. 12). Bombus flavescens is likely to have dispersed to and colonised the Philippines islands when the islands or the Pleistocene Aggregate Island Complex (PAIC), including Luzon, Negros-Panay or Visayan, and Mindanao, were connected in the LGM (Brown et al. 2013). This study supports that B. flavescens dispersed to the Philippines by the Taiwan-Luzon route in the same way as B. irisanensis, the other species in the Philippines. The oscillation of sea level and raising of global temperature in the late Pleistocene facilitated the population isolation and higher latitude immigration of B. flavescens (B, C; Fig. 12). First, the populations from the North (Luzon) and the South (Negros and Mindanao) were separated, because the South population (taxon bakeri) is the sister of the North population (taxon baguionensis) and the Taiwan group (Fig. 7). Then, the populations on the North of Philippines and Taiwan were isolated when the sea submerged the bridges, so that the populations in Southeast Asia were captured in montane refuges as biological islands until the present, similar to the populations on the Cameron Highlands and Gunung Tahan mountain in Malaysia (D; Fig. 12).

Figure 12. 

Biogeographic scenarios of Bombus flavescens through time A the most recent ancestor of B. flavescens dispersed through Southeast Asia or Sundaland around 1 Ma (orange) B the Malaysian population was isolated on the highland of the peninsula due to the temperature rising around 1 Ma (blue), at that time, the populations in Taiwan and Philippines were connected to the mainland populations (orange) C the bridge between the mainland and Taiwan was submerged, the Taiwan-Philippines populations were isolated around 0.5 Ma (purple) D at present day, B. flavescens populations distribute in the subtropical of Southeast Asia, China and Himalaya (yellow), on the Highlands of Malay peninsula (blue), and the islands of Taiwan and the Philippines (purple).

Biogeography of bumblebees in Southeast Asia

Many lineages of bumblebees in Southeast Asia might are hypothesised to have originated around the QTP, then dispersing into the region via the Himalaya–Hengduan corridor (Williams 1985; Hines 2008; Williams et al. 2022). The bumblebee fauna at the Northern border and highlands of the Southeast Asian region is similar to the fauna of the Himalaya and Southern China, adjacent neighbours both spatially and genetically (Williams et al. 2009; Williams et al. 2010; Hines and Williams 2012; Streinzer et al. 2019; Williams 2022). The subgenera Orientalibombus Richards and Alpigenobombus Skorikov are mainly restricted to the North of the Southeast Asian region, whereas the subgenera Megabombus Dalla Torre, Melanobombus Dalla Torre, and Pyrobombus are more widely distributed throughout the region (Williams 1998).

The long-faced subgenus Megabombus crown group might have originated around 13 Ma (Hines 2008). Then, two independent lineages of subgenus Megabombus diverged in Southeast Asia approximately between 4.25–1.5 Ma, based on molecular dating (Huang et al. 2015a). The lineages are 1) a trifasciatus-group: B. montivagus, B. albopleuralis Friese, B. burmensis (Skorikov) in the North of Southeast Asia and the Malay Peninsula, and 2) a senex-group: B. irisanensis on Luzon Island, in the Philippines. However, due to a lack of genetic information on B. senex Vollenhoven and B. melanopoda Cockerell, the phylogenetic relationship between mainland and Sumatran subgenus Megabombus remains relatively unresolved (Huang et al. 2015a). Consequently, the biogeography of Sumatran subgenus Megabombus bumblebees remains unclear.

Nonetheless, the link between the Southeast Asian mainland and the Indonesian islands can be explained by the phylogenetic relationship of two subgenus Melanobombus sister species, B. eximius Smith and B. rufipes Lepeletier of the rufipes-group. The MRCA of the rufipes-group lineage diverged from the other subgenus Melanobombus lineages around 16 Ma and was distributed in the Himalaya and QTP (Williams et al. 2020). Speciation within the rufipes-group is likely to have occurred around 1 Ma in the Pleistocene epoch, after the Sundaland land bridge was submerged, isolating part of the rufipes-group on Sumatra and Java (Williams et al. 2020).

Southeast Asian bumblebees in the subgenus Pyrobombus are also distributed both on the mainland and on the islands of the Philippines, for example, B. flavescens. There is no record of subgenus Pyrobombus on the islands of Indonesia and adjacent areas. A member of subgenus Pyrobombus had been recorded on the Andaman Islands, named B. andamanus Gribodo (Gribodo 1882). However, this is a mislabeled specimen of a Nearctic bumblebee, B. bifarius Cresson (Williams 1998).


Genetic information proved crucial for the study of bumblebees in Southeast Asia. This is the first gene-based study to address the taxonomic status of B. flavescens, which is highly variable in hair colour. This study confirms that the populations of B. flavescens on the Asian mainland and on the islands are parts of the same species. Despite this, a trend among B. flavescens populations can be observed towards larger ocelli at lower latitudes. This might only reflect local selection based on foraging in the low-light conditions within dense forest, or in more early morning or twilight activity in warmer environments. Bombus flavescens originated in the Himalaya and dispersed to Southeast Asia during the Pleistocene. Its constituents, various regional colour forms, diversified through an allopatric divergence process. Bombus flavescens is a useful model for studying the biogeography of bumblebees in Southeast Asia, many of which are less known. Nevertheless, more genetic information is required to investigate the conservation of endemic populations of B. flavescens.


We would like to thank K. Atthasopa and K. Klaithin (ACMU), N. Badruddin (FRIM), F. Bakker and W. van Bohemen (NMNL), N. Chatthanabun and P. Nalinrachatakan (CUNHM), N. Likhitrakarn (MJUMZ), K. Lee (Entopia), D. Notton and J. Monks (NHMUK), N. Pinkaew (KKIC), B. Santos and A. Touret-Alby (MNHN) for their generous help in accessing collections. Thanks are extended to T. Srimaneeyanon for collecting fresh Thai specimens; to I. Barnes and S. Brace (NHMUK) for their suggestion on ancient DNA technique; to M. Musthafa for sharing information on light trapping in the Cameron Highlands; and to P. Nalinrachatakan and S. Yendee for facilitating specimen imaging at CUNHM. Finally, we thank N. Meeprom and J. Satthaphorn for helping with bioinformatic software. This research was supported by a postgraduate scholarship of the Royal Thai Government to CT.


  • Ascher JS, Soh ZWW, Chui SX, Soh EJY, Ho BM, Lee JXQ, Gajanur AR, Ong XR (2022) The bees of Singapore (Hymenoptera: Apoidea: Anthophila): First comprehensive country checklist and conservation assessment for a Southeast Asian bee fauna. The Raffles Bulletin of Zoology 70: 39–64.
  • Ashmead WH (1905) Additions to the recorded Hymenopterous fauna of the Philippine Islands, with descriptions of new species. Proceedings of the United States National Museum 28(1413): 957–971.
  • Ballare KM, Pope NS, Castilla AR, Cusser S, Metz RP, Jha S (2019) Utilizing field collected insects for next generation sequencing: Effects of sampling, storage, and DNA extraction methods. Ecology and Evolution 9(24): 13690–13705.
  • Bischoff H (1936) Schwedisch-chinesische wissenschaftliche Expedition nach den nordwestlichen Provinzen Chinas, unter Leitung von Dr. Sven Hedin und Prof. Sü Ping-chang. Insekten gesammelt vom schwedischen Arzt der Expedition Dr. David Hummel 1927–1930. 56. Hymenoptera. 10. Bombinae. Arkiv för Zoologi 27: 1–27.
  • Bouckaert R, Vaughan TG, Barido-Sottani J, Duchêne S, Fourment M, Gavryushkina A, Heled J, Jones G, Kühnert D, De Maio N, Matschiner M, Mendes FK, Müller NF, Ogilvie HA, du Plessis L, Popinga A, Rambaut A, Rasmussen D, Siveroni I, Suchard MA, Wu C-H, Xie D, Zhang C, Stadler T, Drummond AJ (2019) BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS Computational Biology 15(4): e1006650.
  • Brasero N, Ghisbain G, Lecocq T, Michez D, Valterová I, Biella P, Monfared A, Williams PH, Rasmont P, Martinet B (2021) Resovling the species status of overlooked West-Palaearctic bumblebees. Zoologica Scripta 50(5): 616–632.
  • Brown RM, Siler CD, Oliveros CH, Esselstyn JA, Diesmos AC, Hosner PA, Linkem CW, Barley AJ, Oaks JR, Sanguila MB, Welton LJ, Blackburn DC, Moyle RG, Peterson AT, Alcala AC (2013) Evolutionary Processes of Diversification in a Model Island Archipelago. Annual Review of Ecology, Evolution, and Systematics 44(1): 411–435.
  • Call E, Mayer C, Twort V, Dietz L, Wahlberg N, Espeland M (2021) Museomics: Phylogenomics of the Moth Family Epicopeiidae (Lepidoptera) Using Target Enrichment. Insect Systematics and Diversity 5(2): 1–6.
  • Cameron P (1910) Some Further Notes on Nocturnal Hymenoptera. The Annals of Scottish Natural History 74: 86–87.
  • Cameron SA, Derr JN, Austin AD, Wooley JB, Wharton RA (1992) The application of nucleotide sequence data to phylogeny of the Hymenoptera: A review. Journal of Hymenoptera Research 1(1): 63–79.
  • Chiu SC (1948) Revisional notes on the Formosan bombid-fauna (Hymenoptera). Notes d’entomologie chinoise 12: 57–81.
  • Cockerell TDA (1917) Two new Humble-bees from China. Entomologist 50(655): 265–266.
  • Cockerell TDA (1920) Supplementary notes on the social bees of the Philippine Islands. Philippine Journal of Science 16: 631–632.
  • Darriba D, Taboada GL, Doallo R, Posada D (2012) jModelTest 2: More models, new heuristics and parallel computing. Nature Methods 9(8): e772.
  • Ding G, Zhang S, Huang J, Naeem M, An J (2019) Colour pattern, distribution and food plants of the Asian bumblebee Bombus bicoloratus (Hymenoptera: Apidae). Apidologie 50(3): 340–352.
  • Dorey JB, Fagan-Jeffries EP, Stevens MI, Schwarz MP (2020) Morphometric comparisons and novel observations of diurnal and low-light-foraging bees. Journal of Hymenoptera Research 79: 117–144.
  • Dover C (1929) Wasps and Bees in the Raffles Museum, Singapore. Bulletin of the Raffles Museum 2: 43–70.
  • Dowton M, Austin AD (1994) Molecular phylogeny of the insect order Hymenoptera: Apocritan relationships. Proceedings of the National Academy of Sciences of the United States of America 91(21): 9911–9915.
  • Friese H (1905) Neue oder wenig bekannte Hummeln des russischen Reiches (Hymenoptera). Ezhegodnik Zoologicheskago muzeya 9(1904): 507–523.
  • Frison TH (1925) The bumble bees of the Philippine Islands (Bremidae: Hymenoptera). Philippine Journal of Science 27: 113–121.
  • Frison TH (1928) The bumble bees of the Philippine Islands (Bremidae: Hymenoptera). Philippine Journal of Science 37: 273–281.
  • Frison TH (1934) Records and descriptions of Bremus and Psithyrus from Formosa and the asiatic mainland. Transactions of the Natural History Society of Formosa 24: 150–185.
  • Fujisawa T, Barraclough TG (2013) Delimiting species using single-locus data and the generalized mixed Yule coalescent approach: A revised method and evaluation on simulated data sets. Systematic Biology 62(5): 707–724.
  • Ghisbain G, Martinet B, Wood TJ, Przybyla K, Cejas D, Gérard M, Rasmont P, Monfared A, Valterová I, Michez D (2021) A worthy conservation target? Revising the status of the rarest bumblebee of Europe. Insect Conservation and Diversity 14(5): 661–674.
  • Goodwin S, McPherson J, McCombie W (2016) Coming of age: Ten years of next-generation sequencing technologies. Nature Reviews. Genetics 17(6): 333–351.
  • Gribodo G (1882) Alcune nuove specie e nuove genere di imenotteri aculeati. Annali del Museo Civico di Storia Naturale Giacomo Doria 18: 261–268.
  • Hebert PDN, Penton EH, Burns JM, Janzen DH, Hallwachs W (2004) Ten species in one: DNA barcoding reveals cryptic species in the neotropical skipper butterfly Astraptes fulgerator. Proceedings of the National Academy of Sciences of the United States of America 101(41): 14812–14817.
  • Hedicke H (1926) Beiträge zur Apidenfauna der Philippinen (Hym.). (2. Beitrag zur Kenntnis orientalischer Apiden.). Deutsche Entomologische Zeitschrift 1926: 413–423.
  • Hines HM (2008) Historical Biogeography, Divergence Times, and Diversification Patterns of Bumble Bees (Hymenoptera: Apidae: Bombus). Systematic Biology 57(1): 58–75.
  • Hines HM, Williams PH (2012) Mimetic colour pattern evolution in the highly polymorphic Bombus trifasciatus (Hymenoptera: Apidae) species complex and its comimics. Zoological Journal of the Linnean Society 166(4): 805–826.
  • Huang J, An J, Wu J, Williams PH (2015a) Extreme Food-Plant Specialisation in Megabombus Bumblebees as a Product of Long Tongues Combined with Short Nesting Seasons. PLoS ONE 10(8): e0132358.
  • Huang J, An J, Wu J, Williams PH (2015b) Newly discovered colour-pattern polymorphism of Bombus koreanus females (Hymenoptera: Apidae) demonstrated by DNA barcoding. Apidologie 46(2): 250–261.
  • Kapli P, Lutteropp S, Zhang J, Kobert K, Pavlidis P, Stamatakis A, Flouri T (2017) Multi-rate Poisson tree processes for single-locus species delimitation under maximum likelihood and Markov chain Monte Carlo. Bioinformatics 33(11): 1630–1638.
  • Kapustijanskij A, Streinzer M, Paulus HF, Spaethe J (2007) Bigger is better: Implications of body size for flight ability under different light conditions and the evolution of alloethism in bumblebees. Functional Ecology 21(6): 1130–1136.
  • Kerfoot WB (1967) Correlation between Ocellar Size and the Foraging Activities of Bees (Hymenoptera; Apoidea). American Naturalist 101(917): 65–70.
  • Koch JB, General DEM (2019) A preliminary assessment of bumble bee (Hymenoptera: Apidae) habitat suitability across protected and unprotected areas in the Philippines. Annals of the Entomological Society of America 112(1): 44–49.
  • Kumar S, Stecher G, Tamura K (2016) MEGA7: Molecular evolutionary genetics analysis Version 7.0 for Bigger Datasets. Molecular Biology and Evolution 33(7): 1870–1874.
  • Luo A, Ling C, Ho SYW, Zhu C-D (2018) Comparison of methods for molecular species delimitation across a range of speciation scenarios. Systematic Biology 67(5): 830–846.
  • Ma YP, Xie WJ, Sun WB, Marczewski T (2016) Strong reproductive isolation despite occasional hybridization between a widely distributed and a narrow endemic Rhododendron species. Scientific Reports 6(1): e19146.
  • Matuszak S, Muellner-Riehl AN, Sun H, Favre A (2016) Dispersal routes between biodiversity hotspots in Asia: The case of the mountain genus Tripterospermum (Gentianinae, Gentianaceae) and its close relatives. Journal of Biogeography 43(3): 580–590.
  • Matzke NJ (2013b) Probabilistic historical biogeography: New models for founder-event speciation, imperfect detection, and fossils allow improved accuracy and model-testing. Frontiers of Biogeography 5(4): 242–248.
  • Matzke NJ (2014) Model selection in historical biogeography reveals that founder-event speciation is a crucial process in island clades. Systematic Biology 63(6): 951–970.
  • McElrath T (2022) Data from: Illinois Natural History Survey Insect Collection. Illinois Natural History Survey, Occurrence dataset, GBIF. [Accessed 7 October 2022]
  • Musthafa MM, Abdullah F, Martínez-Falcón AP, de Bruyn M (2021) How mountains and elevations shape the spatial distribution of beetles in Peninsular Malaysia. Scientific Reports 11(1): e5791.
  • Nazari V, Schmidt BC, Prosser S, Hebert PDN (2016) Century-old DNA barcodes reveal phylogenetic placement of the extinct Jamaican Sunset moth, Urania sloanus Cramer (Lepidoptera: Uraniidae). PLoS ONE 11(10): e0164405.
  • Ogilvie HA, Bouckaert RR, Drummond AJ (2017) StarBEAST2 brings faster species tree inference and accurate estimates of substitution rates. Molecular Biology and Evolution 34(8): 2101–2114.
  • Pääbo S (1989) Ancient DNA: Extraction, characterization, molecular cloning, and enzymatic amplification. Proceedings of the National Academy of Sciences of the United States of America 86(6): 1939–1943.
  • Papadopoulou A, Anastasiou I, Vogler AP (2010) Revisiting the Insect Mitochondrial Molecular Clock: The Mid-Aegean Trench Calibration. Molecular Biology and Evolution 27(7): 1659–1672.
  • Pendlebury HM (1923) Four new species of Bombus from the Malay Peninsula. Journal of the Federated Malay States Museums 11: 64–67.
  • Perrard A, Arca M, Rome Q, Muller F, Tan J, Bista S, Nugroho H, Baudoin R, Baylac M, Silvain J, Carpenter JM, Villement C (2014) Geographic Variation of Melanisation Pattern in a Hornet Species: Genertic Differences, Climatic Pressures or Aposematic Constraints? PLoS ONE 9(4): e94162.
  • Pittioni B (1949) Beiträge zur Kenntnis der Bienenfauna SO-Chinas. Die Hummeln und Schmarotzerhummeln der Ausbeute J. Klapperich (1937/38). (Hym., Apoidea, Bombini). Eos 25: 241–284.
  • Prosser SWJ, de Waard JR, Miller SE, Hebert PDN (2016) DNA barcodes from century-old type specimens using next-generation sequencing. Molecular Ecology Resources 16(2): 487–497.
  • R Core Team (2021) R: A Language and environment for statistical computing. R Foundation for Statistical Computing, Vienna.
  • Rahman SR, Terranova T, Tian L, Hines HM (2021) Developmental Transcriptomics Reveals a Gene Network Driving Mimetic Color Variation in a Bumble Bee. Genome Biology and Evolution 13(6): evab080.
  • Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA (2018) Posterior summarisation in Bayesian phylogenetics using Tracer 1.7. Systematic Biology 67(5): 901–904.
  • Rasmont P, Ghisbain G, Terzo M (2021) Bumblebees of Europe and neighbouring regions. N.A.P. Editions, Verrières-le-Buisson, 631 pp.
  • Sathiamurthy E, Voris HK (2006) Maps of holocene sea level transgression and submerged lakes on the Sunda shelf. Tropical Natural History 2: 1–44.
  • Sayers EW, Bolton EE, Brister JR, Canese K, Chan J, Comeau DC, Connor R, Funk K, Kelly C, Kim S, Madej T, Marchler-Bauer A, Lanczycki C, Lathrop S, Lu Z, Thibaud-Nissen F, Murphy T, Phan L, Skripchenko Y, Tse T, Wang J, Williams R, Trawick BW, Pruitt KD, Sherry ST (2022) Database resources of the national center for biotechnology information. Nucleic Acids Research 50(D1): D20–D26.
  • Skorikov AS (1933) Zur Hummelfauna Japans und seiner Nachbarländer. Mushi 6: 53–65.
  • Smith F (1852) Descriptions of some new and apparently undescribed species of hymenopterous insects from north China, collected by Robert Fortune, Esq. The Transactions of the Entomological Society of London 2(2): 33–45.
  • Smith F (1854) Catalogue of Hymenopterous Insects in the Collection of the British Museum. Part II. Apidae. London.
  • Sproul JS, Maddison DR (2017) Sequencing historical specimens:successful preparation of small specimens with low amounts of degraded DNA. Molecular Ecology Resources 17(6): 1183–1201.
  • Starr CK (1992) The bumble bees (Hymenoptera: Apidae) of Taiwan. Bulletin of National Museum of Natural Science 3: 139–157.
  • Streinzer M, Huber W, Spaethe J (2016) Body size limits dim-light foraging activity in stingless bees (Apidae: Meliponini). Journal of Comparative Physiology. A, Neuroethology, Sensory, Neural, and Behavioral Physiology 202(9): 643–655.
  • Streinzer M, Chakravorty J, Neumayer J, Megu K, Narah J, Schmitt T, Bharti H, Spaethe J, Brockmann A (2019) Species composition and elevational distribution of bumble bees (Hymenoptera, Apidae, Bombus Latreille) in the East Himalaya, Arunachal Pradesh, India. ZooKeys 851: 71–89.
  • Thanoosing C (2017) Species status of Bombus (Alpinobombus) kluanensis and its close relatives–an exploration of nuclear markers. MSc Thesis. Imperial College London, London, 45 pp.
  • Thanoosing C (2022) Systematics and ecology of Southeast Asian bumblebees. PHD Thesis. Imperial College London, London, 363 pp.
  • Thomsen PF, Elias S, Gilbert MTP, Haile J, Munch K, Kuzmina S, Froese DG, Sher A, Holdaway RN, Willerslev E (2009) Non-destructive sampling of ancient insect DNA. PLoS ONE 4(4): e5048.
  • Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, Rozen SG (2012) Primer3–New capabilities and interfaces. Nucleic Acids Research 40(15): e115.
  • Warrant EJ (2008) Seeing in the dark: Vision and visual behaviour in nocturnal bees and wasps. The Journal of Experimental Biology 211(11): 1737–1746.
  • Williams PH (1991) The bumble bees of the Kashmir Himalaya (Hymenoptera: Apidae, Bombini). Bulletin of the British Museum (Natural History). Historical Series 60(1): 1–204. [Natural History]
  • Williams PH (1996) Mapping variations in the strength and breadth of biogeographic transition zones using species turnover. Proceedings of the Royal Society B, Biological Sciences 263(1370): 579–588.
  • Williams PH (1998) An annotated checklist of bumble bees with an analysis of patterns of description (Hymenoptera: Apidae, Bombini). Bulletin of The Natural History Museum (Entomology) 67: 79–152.
  • Williams PH (2007) The distribution of bumblebee colour patterns world-wide: Possible significance for thermoregulation, crypsis, and warning mimicry. Biological Journal of the Linnean Society. Linnean Society of London 92(1): 97–118.
  • Williams PH (2021) Not just cryptic, but a barcode bush: PTP re-analysis of global data for the bumblebee subgenus Bombus s. str. supports additional species (Apidae, genus Bombus). Journal of Natural History 55(5–6): 271–282.
  • Williams PH (2022) The Bumblebees of the Himalaya. AbcTaxa, Belgium, 198 pp.
  • Williams PH, Ito M, Matsumura T, Kudo I (2010) The bumblebees of the Nepal Himalaya (Hymenoptera: Apidae). Insecta Matsumurana 66: 115–151.
  • Williams PH, Brown MJ, Carolan JC, An J, Goulson D, Aytekin AM, Best LR, Byvaltsev AM, Cederberg B, Dawson R, Huang J, Ito M, Monfared A, Raina RH, Schmid-Hempel P, Sheffield CS, Šima P, Xie Z (2012) Unveiling cryptic species of the bumblebee subgenus Bombus s. str. worldwide with COI barcodes (Hymenoptera: Apidae). Systematics and Biodiversity 10(1): 21–56.
  • Williams PH, Thorp RW, Richardson LL, Colla SR (2014) An Identification Guide: Bumble bees of North America. Princeton University Press, Princeton, 208 pp.
  • Williams PH, Lobo JM, Meseguer AS (2018) Bumblebees take the high road: Climatically integrative biogeography shows that escape from Tibet, not Tibetan uplift, is associated with divergences of present-day Mendacibombus. Ecography 41(3): 461–477.
  • Williams PH, Altanchimeg D, Byvaltsev A, De Jonghe R, Jaffar S, Japoshvili G, Kahono S, Liang H, Mei M, Monfared A, Nidup T, Raina R, Ren Z, Thanoosing C, Zhao Y, Orr MC (2020) Widespread polytypic species or complexes of local species? Revising bumblebees of the subgenus Melanobombus world-wide (Hymenoptera, Apidae, Bombus). European Journal of Taxonomy 719(1): 1–120.
  • Williams PH, Dorji P, Ren Z, Xie Z, Orr M (2022a) Bumblebees of the hypnorum-complex world-wide including two new near-cryptic species (Hymenoptera:Apidae). European Journal of Taxonomy 847(1): 46–72.
  • Williams PH, Françoso E, Martinet B, Orr MC, Ren Z, Santos Junoir J, Thanoosing C, Vandame R (2022b) When did bumblebees reach South America? Unexpectedly old montane species may be explained by Mexican stopover (Hymenoptera: Apidae). Systematics and Biodiversity 20(1): 1–24.
  • Williams PH, Sung I-H, Lin Y-J, Lu S-S (2022c) Discovering endemic species among the bumblebees of Taiwan (Apidae, genus Bombus). Journal of Natural History 56(5–8): 435–447.
  • Woodruff DS (2010) Biogeography and conservation in Southeast Asia: How 2.7 million years of repeated environmental fluctuations affect today’s patterns and the future of the remaining refugial-phase biodiversity. Biodiversity and Conservation 19(4): 919–941.
  • Yu Y, Harris AJ, He X (2010) S-DIVA (Statistical Dispersal-Vicariance Analysis): A tool for inferring biogeographic histories. Molecular Phylogenetics and Evolution 56(2): 848–850.
  • Yu WB, Liu ML, Wang H, Mill RR, Ree RH, Yang JB, Li DZ (2015) Towards a comprehensive phylogeny of the large temperate genus Pedicularis (Orobanchaceae), with an emphasis on species from the Himalaya-Hengduan Mountains. BMC Plant Biology 15(1): e176.
  • Yu Y, Blair C, He X (2020) RASP 4: Ancestral state reconstruction tool for multiple genes and characters. Molecular Biology and Evolution 37(2): 604–606.
  • Zhou Z-K, Su T, Huang Y-J (2018) Neogene Paleoenvironmental Changes and their Role in Plant Diversity in Yunnan, South‐Western China. In: Hoorn C, Perrigo A, Antonelli A (Eds) Mountains, Climate and Biodiversity. John Wiley & Sons, Hoboken, 449–458.

Supplementary materials

Supplementary material 1 

List of COI sequences from databases, research and collaboration, including the accession number or ID, the original sequence length and country

Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams

Data type: Sequence ID (word document)

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (16.53 kb)
Supplementary material 2 

Primer testing

Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams

Data type: Text (word document)

Explanation note: Testing newly designed primers in this study.

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (12.71 kb)
Supplementary material 3 

Primer testing with gradient annealing temperatures

Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams

Data type: Experimental (word document)

Explanation note: The samples were B. breviceps (CT#552) and B. flavescens (CT#662). The brightness of gel electrophoresis is presented in symbols: *** = strong, ** = medium, * = weak, NS = unsuccessful.

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (13.79 kb)
Supplementary material 4 

Polymerase chain reaction temperature profile of different pairs of primers in this study

Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams

Data type: Experimental (word document)

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (13.52 kb)
Supplementary material 5 

List of 16S and PEPCK sequences from GenBank with the accession number, the original sequence length, and country

Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams

Data type: Sequence ID (word document)

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (14.23 kb)
Supplementary material 6 

The accession number or ID of the parsing dataset of COI, 16S, and PEPCK sequences for *BEAST analysis

Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams

Data type: Sequence ID (word document)

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (14.96 kb)
Supplementary material 7 

Ancestral area reconstruction model selection, estimated in BioGeoBEARS

Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams

Data type: Model test (word document)

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (13.41 kb)
Supplementary material 8 

List of specimens, including species, collection, specimen ID, origin country and latitude group, and measurement of intertegular distance (ID), head width (HW), median ocelli width (MOW) in millimetres (mm)

Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams

Data type: Morphological (word document)

Explanation note: List of specimens, including species, collection (CUNHM = Chulalongkorn University Natural History Museum, Bangkok, Thailand; KKIC = Kasetsart Kamphaeng Saen Insect Collection, Kasetsart University Kamphaeng Saen Campus, Nakhon Pathom, Thailand; NHMUK = Natural History Museum, London, UK; NMNL = Naturalis Biodiversity Center, Leiden, the Netherlands; PHW = Paul H. Williams research collection, UK), specimen ID, origin country and latitude group (Low = 0°–20°N; High = > 20°N), and measurement of intertegular distance (ID), head width (HW), median ocelli width (MOW) in millimetres (mm) (available at the NHM Data Portal:

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (42.27 kb)
login to comment