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 (https://doi.org/10.5519/qd7f4uuw and https://doi.org/10.5519/isxh6saw). 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) (GBIF.org 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.

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.

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 2–4). The annealing temperature for the pair FLA1 and FLA2 was set at 65 °C and PRA1 was set at 63 °C.

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: https://www.megasoftware.net/), to trim low-quality bases near the ends of the traces. The COI barcodes have been deposited in GenBank (Table 2).

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 1000^th generation. The first 10% was discarded as burn-in. The phylogenetic trees were visualised using FigTree version 1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/). 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 (https://species.h-its.org; 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; https://r-forge.r-project.org/projects/splits/).

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).

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 5000^th 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.

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).

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.

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 1000^th 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.

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 5000^th 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.

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.

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.

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: https://cameronservice.blogspot.com/), 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.

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).

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).