JHR 96: 507-541 (2023) Ae eS,” JOURNAL OR. She tentonteatiertt doi: 10.3897/jhr.96. 104715 RESEARCH ARTICLE ) I Tymenopter a % https://jhr.pensoft.net The Inarasional Society of Hymenopeeriss, RESEARCH A taxonomic re-assessment of the widespread oriental bumblebee Bombus flavescens (Hymenoptera, Apidae) Chawatat Thanoosing'*?, Michael C. Orr*°, Natapot Warrit’, Alfried P. Vogler'*, Paul H. Williams! | Department of Life Sciences, Natural History Museum, Cromwell Road, London, SW7 5BD, UK 2 Department of Life Sciences, Silwood Park Campus, Imperial College London, Buckhurst Road, Ascot, SL5 7PY, UK3 Center of Excellence in Entomology and Department of Biology, Faculty of Science, Chulalongkorn University, Bangkok 10330, Thailand 4 Institute of Zoology, Chinese Academy of Sciences, 1 Beichen West Road, Chaoyang District, Beijing, 100101, China § Entomologie, Staatliches Museum ftir Naturkunde Stuttgart, Stuttgart, Germany Corresponding author: Chawatat Thanoosing (chawatat.thanoosing 16@imperial.ac.uk; t.chawatat@gmail.com) Academic editor: Michael Ohl | Received 11 April 2023 | Accepted 8 June 2023 | Published 22 June 2023 Attps://zoobank. org/2ACO9EB7- 1B86-48 18-A803-8B 1 170FE8373 Citation: Thanoosing C, Orr MC, Warrit N, Vogler AP, Williams PH (2023) A taxonomic re-assessment of the widespread oriental bumblebee Bombus flavescens (Hymenoptera, Apidae). Journal of Hymenoptera Research 96: 507— 541. https://doi.org/10.3897/jhr.96.104715 Abstract 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 morphol- ogy. 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 polymor- phisms 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. Keywords COI, Museum specimens, Polymorphism, Pyrobombus Copyright Chawatat Thanoosing et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 508 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) Introduction Bumblebees (genus Bombus Latreille) are well-studied pollinators, especially in tem- perate 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 mor- phology 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; Wil- liams 1998, 2022; Williams et al. 2009; Koch and General 2019). Most frequently, B. flavescens is a predominantly black pubescence (or hair) bum- blebee 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 descrip- tion 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. bagui- onensis 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 typi- cal B. flavescens: B. mearnsi Ashmead from Mindanao in the Philippines (morphologi- cally synonymised by Pittioni (1949)); B. bakeri Cockerell from Negros (Fig. 2H; mor- phologically synonymised by Frison (1925)); and B. tahanensis Pendlebury from the Tahan mountains of Malaysia (Fig. 21; morphologically synonymised by Frison (1934)). Apart from this high level of colour variation, other morphological characters, includ- ing shape of labrum, punctures on clypeus, and male genitalia, are relatively similar The oriental bumblebee Bombus flavescens 509 hs NHMUK 014025381 flavescens omm FLA# 92: Figure |. Holotype male of B. flavescens, deposited in the entomological type collection at the Natural History Museum, London (NHMUK 014025381). Scale bar: 5 mm. D E | oF 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 re- versed. Scale bar: 5 mm. (Williams 1998; Thanoosing 2022). Yet, at various times through the 20" century dif- ferent intraspecific names have been applied to describe minor B. flavescens variants (Fri- son 1925, 1928, 1934; Bischoff 1936; Chiu 1948; Pittioni 1949; Table 1) although they are clearly grouped together because of their morphological similarity (Williams 1998). 510 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) Table |. 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 vat. baguionensis, Bombus bakeri Pendlebury (1923) | Bombus tahanensis, Bombus rufoflavus Frison (1925) Bremus mearnsi, Bremus mearnsi var. bakeri, Bremus irisanensis vat. baguionensis Hedicke (1926) Bombus imuganensis Frison (1928) Bremus (Pratobombus) baguionensis, Bremus (Pratobombus) baguionensis vat. imuganensis Dover (1929) Bremus rufoflavus Skorikov (1933) Pratobombus flavescens Frison (1934) Bremus mearnsi vat. deflectus, Bremus mearnsi vat. ditutus, Bremus mearnsi vat. bakeri, Bremus mearnsi vat. geei Bischoff (1936) Bombus (Pratobombus) mearnsi ssp. chekiangensis Chiu (1948) Bombus mearnsi vat. deflectus, Bombus mearnsi vat. ditutus, Bombus mearnsi vat. bakeri, Bombus mearnsi vat. 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 ex- pression of, e.g., the Abd-B and nubbin genes affecting pigment generation during pu- pal 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 Island, baguionensis shows a similar pat- tern 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. max- welli 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 The oriental bumblebee Bombus flavescens pti 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 bioge- ography 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 Cam- eron 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 estab- lish 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 Se- cretions or CLGS (e.g., Brasero et al. (2021), Ghisbain et al. (2021)). Given the chal- lenges 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. fla- vescens, using morphological examination of their ocelli. Methods Sampling Museum and institutional specimens of the flavescens-complex (Table 1) were exam- ined for morphological analysis. Collecting information was recorded and made avail- able 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, in- cluding: 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 Re- search Institute (TFRI; Lu S), available on the Global Biodiversity Information Fa- cility (GBIF) (GBIEorg 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 estab- lished 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). 5112, Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) 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 Univer- sity Natural History Museum Collection, NMNL = Naturalis Biodiversity Center Collection), sex/caste (w = worker, m = male), locality, collecting date, and COI GenBank accession number. Species ProjectID Collection (specimen ID) Sex/caste Locality Collecting GenBank date accession B. flavescens CT#662 KKIC w Thailand, Nakorn = 7/5/2017 OP355718 Pathom? CT#669 KKIC m Thailand, Loei 13/4/2016 OP355719 CT#926 PHW w Malaysia, NA OP355720 Gunung Tahan CT#1018 NMNL (ZMA.INS.758598) w Philippines, Negros 4/5/1953. OP355722 CT#1023 CUNHM (BSRU-AB-9609) w Thailand, 15/2/2021 .OP355725 Chiang Mai CT#1024 CUNHM (BSRU-AB-9610) Ww Thailand, 15/2/2021 OP355724 Chiang Mai FLA#2 PHW w Philippines, Luzon 27/4/1986 OP355725 FLA#3 PHW Ww 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) Ww Bhutan, Lungtenphu 10/6/1996 OP355730 FLA#10 PHW Ww Philippines, Luzon 27/4/1986 OP355731 B. rotundiceps CT#960 CUNHM (BSRU-AB-1268) m Thailand, 17/6/2019 OP355721 Chiang Mai ROT#1 PHW w India, Uttar Pradash 6/5/1990 OP355732 ROT#2 PHW Ww India, Uttar Pradash 6/5/1990 OP355733 ROT#3 PHW Ww China, Guangxi 3/6/2016 OP355734 ROT#4 PHW w China, Guangxi 5/6/2016 OP355735 The samples were selected to represent the geographical distribution and all major col- our pattern variants of B. flavescens (Fig. 3). Bombus flavescens belongs to the pratorum- group (Williams 1998). Thus, closely related species in the pratorum-group (Cam- eron 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. The oriental bumblebee Bombus flavescens 513 oO ae) =) —_ — © — Ss 70°E 80°E 90°E 100°E 110°E 120°E 130°E Longitude 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. 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 sam- ples 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, fol- lowing 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 Maddi- son (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 514 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) with the enzyme) was included for each extraction batch as a contamination control. Samples were eluted from the column twice with 15 yl 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 ampli- cons (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 (Gen- Bank 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. pra- torum (GenBank accession number: K 1164684). There are three pairs of new primers for partial COI amplification (contig) in this study (Table 3): 1) FLA1: ParCOIl_246_ FLA“ Io Cim? “= -66,2°°C) and, ParCOl “416 FLA Rh (im? =°63,7.°C), 2)-FLA2: ParCOI_349 FLA F2 (Tm® = 64.2 °C) and ParCOI_662_FLA_R2 (Tm° = 66.6 °C), and 3) PRAI: ParCOI_1816_PRA_F1 (Im° = 61.7 °C) and ParCOI_1994_PRA_ RI (TIm® = 65.5 °C). These primers were tested for amplification efficiency (Suppl. materi- als 2-4). The annealing temperature for the pair FLA1 and FLA2 was set at 65 °C and PRA 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). 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 (Ron- Table 3. Newly designed primers were used in this study, including gene, primer name, strand (F = for- ward, R = reverse) and primer sequence. Gene Primer name Strand Primer sequence Partial COI ParCOI_246_FLA_F1 lf 5'— CCT GACATAGCPITEGCCACGA -3' ParCOI_416_FLA_R1 5’—- TGCAATATCAACTGAAGGTGATG -3’ ParCOI_349_FLA_F2 5'- CAGGATGAACTGTT FACCETCCT =3° ParCOI_662_FLA_R2 5’ TGGATCACCT CCT CCTATTGGA.=3' ParCOI_1816_PRA_F1 5’- TTCGTATAGAATTAAGTCATCCTGGT -3’ ParCOI_1994_PRA_R1 5’- TCGTGGAAAAGCTATATCAGGTGAT -3’ A mA TA The oriental bumblebee Bombus flavescens S15 quist and Huelsenback 2003). Markov Chain Monte Carlo (MCMC) chains were run for 10 million generations and sampled every 1000" generation. The first 10% was discarded as burn-in. The phylogenetic trees were visualised using FigTree ver- sion 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. (Bom- bus) 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 diversi- fication 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; ac- cessed 2021) under default parameters. The analysis was restricted to unique haplotype sequences in order to avoid introducing misleading ‘zero-length’ branches. ‘The in- put phylogenetic tree was estimated using MrBayes rooted with Bombus terrestris. Re- sults 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 ultra- metric 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/). 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 re- trieved 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 516 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) 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 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" 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+] 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 bio- geographical 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) Cen- tral 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 cho- sen 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 pos- sible bumblebee permitted dispersal routes (Fig. 4). According to this model, the pos- sible dispersal areas in this study are AB, ABC, ABE, AC, ACD, ACE, BEF, BE BFG, CD, CDE, CE, CEE CEG, CEI, EK EFG, EFI, EG, EGH, EGI, EI, EJ, FG, FGH, GH, GHK, HJK, HK, JJ, IJK, and JK. The burnin was set at 20%. Then, the *BEAST trees were used as the input trees. Morphological study of ocelli in the flavescens-complex The idea of the nocturnal lifestyle in the orange pattern bumblebees, rufoflavus (B. fla- vescens) and maxwelli (B. montivagus) from Peninsular Malaysia, has been introduced by The oriental bumblebee Bombus flavescens Say 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 ke 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) C North Asia/ Central Asia Figure 4. A corridor-dispersal model diagram of the flavescens-complex and their relatives in this study. Williams (2007). Although orange hair pattern can be recognised in diurnal bumblebees (e.g., B. humillis liger, B. pascuorum (Scopoli), and B. muscorum (Linnaeus)), rufoflavus and maxwelli show the uniformly orange pattern both of pubescence and sclerites unique- ly. 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 Ma- laysia. For this reason, we selected the female specimens of X. myops and their related spe- cies, 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 fe- male 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 518 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) 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. morphological character measurement, including median ocellus width (MOW), inter- tegular distance (ID), and head width (HW) (Fig. 5), under a Wild Heerbrugg M4A ster- eomicroscope, calibrated with ocular and stage micrometres. The specimens of the flaves- cens-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. 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 dis- tributed. The Dunn’s test was conducted using the R package FSA (Ogle et al. 2022). Results 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/ul). 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 oriental bumblebee Bombus flavescens 519 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+I° 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. 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+I (BIC = 3877). Species delimitation using PTP returned each of the nine spe- cies 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. ta- hanensis, 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 maxi- mum 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. [he 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. Although DEC model was the best fit model for the BioGEOBEARS analysis in this study (the highest AlICc wt; Suppl. material 7), we selected DIVALIKE+J model instead (the second high AlCc wt; Suppl. material 7). The DIVALIKE+J model sup- ports 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+] showed that there was uncertainty in ances- tral 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+] 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+] with low probability. This group diversified in the late Miocene. 520 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) 0.58 0.51 0.3 0.47 0.01 1 687 polaris |4856|CANADA-Yukon 658 terrestris | TF1|UK 664 nursei|6879H11|INDIA-Kashmir 604 biroi|6877G05|KYRGYZSTAN 604 biroi|6877G06| KYRGYZSTAN 602 wangae|3263D07 | CHINA-Sichuan 1 601 wangae|3262D07| CHINA-Sichuan 601 wangae|3262D08|CHINA-Sichuan 596 wangae | 3263D08|CHINA-Sichuan 607 modestus | ASBEE139-08| MONGOLIA 657 modestus | MF361435 | SOUTH KOREA 657 modestus | ASBEEO45-08 | MONGOLIA 0.93 667 modestus | PCHELO25-09|RUSSIA-Kirovskaya Oblast 593 modestus | ASBEE140-08| RUSSIA-Ussuriland 601 modestus | ASBEE043-08 | CHINA-Shanxi 605 modestus | ASBEE137-08 | CHINA-Shanxi 653 ardens|MF361549|SOUTH KOREA 653 ardens|MF361530|SOUTH KOREA 653 ardens| MF361482|SOUTH KOREA 653 ardens|MF361483|SOUTH KOREA 686 pratorum| BMNH970335-BEEEE120-15| UK | 617 pratorum|ACUFI1211-13| FINLAND 659 pratorum|POLLE1300-19 | FRANCE 658 pratorum|FBHAP827-09| GERMANY 658 pyrenaeus | FBHAP828-09|GERMANY 658 pyrenaeus | KJ837876 | FRANCE 674 rotundiceps | CTROT4 | CHINA-Guangxi 674 rotundiceps | CTROT3 | CHINA-Guangxi 671 rotundiceps |CT960| THAILAND 667 rotundiceps | 182-09 | NEPAL 527 rotundiceps|CTROT2-assem|INDIA 527 rotundiceps|CTROT1-assem| INDIA — 667 rotundiceps|181-09| NEPAL — 527 flavescens|CTFLA5-assem|INDIA 527 flavescens | CTFLA7-assem|BHUTAN 610 flavescens|CT662| THAILAND - 527 flavescens|CT669-assem| THAILAND 634 flavescens|CT1024| THAILAND 629 flavescens|CT1023 | THAILAND 658 flavescens | BOFTW091-08 | THAILAND 658 flavescens | BOFTH434-09 | THAILAND 0.57 667 flavescens | XL2-2018-6-25 | CHINA-Yunnan 667 flavescens | Q27-XL4-2018-6-25 | CHINA-Yunnan 0.63 672 flavescens |FLA6|CHINA-Yunnan 667 flavescens | 25-XL9-2018-6-24| CHINA-Yunnan 667 flavescens | 7-XL5-2018-6-24| CHINA-Yunnan 667 flavescens | 19-XL3-2018-6-25 | CHINA-Yunnan 667 flavescens | XL2-2018-6-24 | CHINA-Yunnan 0.521 667 flavescens | Q28-XL4-2018-6-24 | CHINA-Yunnan 667 flavescens | 36-XL1-2018-6-25 | CHINA-Sichuan 667 flavescens | 17-XL8-2018-6-24|CHINA-Yunnan 0.98 658 rufoflavus | 1550E08 | MALAYSIA | 658 tahanensis | 1550E07 | MALAYSIA 0.76 | 9-65 527 tahanensis|CT926-assem|MALAYSIA | 0.93 527 flavescens|CTFLA4-assem | TAIWAN 0.98 | 0.23 } 0.32 527 flavescens| CTFLA3-assem| TAIWAN 545 bakeri|CT1018-assem|PHILIPPINES-Negros 0.91 527 baguionensis | CTFLA2-assem | PHILIPPINES-Luzon 527 baguionensis | CTFLA10-assem | PHILIPPINES-Luzon f+) iy “ay. a eS ee a ao, 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" 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. 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, includ- ing Himalaya, Thailand, and Peninsular Malaysia (Node X; Fig. 8) at around 1 Ma. However, the DIVALIKE+J biogeographic scenario suggested that the group might The oriental bumblebee Bombus flavescens D2Al baguionensis flavescens (Taiwan) bakeri d96 flavescens (Thailand) flavescens (Himalaya) G.84 flavescens (China) — rr ie fa rufoflavus 1.00 1.00 L tahanensis Tad rotundiceps ie pratorum | O57 | ardens modestus 1.00 wangae a 0.99 pyrenaeus biroi 1.00 nursei lepidus Millions of years ago 11 10 9 8 7 6 5 4 3 2 1 0 Figure 7. The maximum clade credibility tree of B. flavescens and their closely related species, recon- structed 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" 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. have originated in Malaysia (H), but with low probability. The DIVALIKE+] 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). Accord- ing 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 permit- ted 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). 522 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) (J) baguionensis @ A- Evrope @ 6 - Thailand ® B-CAsia @ H - Malaysia 1.00 Xvi © c-NAsia+NChina @® t-Taiwan UK aad (1) flavescens (Taiwan) @ D - Korea + Japan @ J- Luzon 0.43 XV @ £-Schina @ k- Negros (K) bakeri @) F - Himalaya WJ oop xt (G) flavescens (Thailand) FG 1.00 XIV aes (F) flavescens (Himalaya) cH “tS ag XIII 029/ * (E) flavescens (China) IX (H) rufoflavus a > (H) tahanensis FO vu (EFG) rotundiceps (A) pratorum ‘ (D) ardens (CD) modestus CDE 0.50 vl ABE ACD (E) wangae (A) pyrenaeus BEF gm (B) biroi (F) nursei (EF) /epidus 8 7 6 5 4 3 2 1 0 Millions of years ago 7 ‘Miocene 7 Pliocene _ | Pleistocene ] 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. 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-val- ue = 6.781, p-value < 0.05; Tukey test, p-value < 0.05): bakeri-China, rufoflavus-Chi- na, tahanensis-China, bakeri-Himalaya , rufoflavus-Himalaya, tahanensis-Himalaya, Taiwan-bakeri, Taiwan-rufoflavus, and Taiwan-tahanensis. In addition, only two pairs The oriental bumblebee Bombus flavescens 523 @ 4- Europe @ 6-Thaitana Pie (J) baguionensis 0.50 @ 8-casia @ b- Malaysia kal (I) flavescens (Taiwan) © c-NAsia+NChina @ 1- Taiwan oas a XV @ D - Korea + Japan @ J- Luzon (K) bakeri 7 4 | 0.30 = : as @ «Negros (E) flavescens (China) - Himalaya = @ Xill , 049 (F) flavescens (Himalaya) L@xiv 0.60 a5 (G) flavescens (Thailand) P (H) rufoflavus H 0.04 Xl at (H) tahanensis 0.23 vill (EFG) rotundiceps (A) pratorum A amy "7 (D) ardens (CD) modestus 0.53 vl B Ad . (E) wangae (A) pyrenaeus 7 Ty (B) biroi 0.03 (F) nursei (EF) /epidus 8 7 6 5 4 3 2 1 O Millions of years ago Miocene Pliocene | Pleistocene 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. E 0.16 0.12 = = a = = = = : = 0.10 = 6 0.12 s = = S o (eo) 4 (e) 5 i = 0.08 Cc G 0.08 or xo) © => 0.06 0.04 \ oy . “ e ot e Figure 10. Box plots showing the median ocellus width (MOW), the ratio between MOW and inter- tegular 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. 524 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) = 0.080 £ s Q = ce) = < 0.075 > = = rn) O O =) = = = 8 2 ° 0.070 fay] —_— oO 5 = = ro] 0.065 ® = 0.060 & soo go © S32 of S o 19 W 2o & gv Ss eo sv wv Cy 0 22 gh SP ay cs PAN oS 3 WP CK aS Ks om A Pag! +S FOF Paw ae SES FF Pag? S Na es ae we ae) os & ad SS Na of eS eo Sad Pe ON «6 Oi vege {0 ROS & us vo K RO & ‘“e ee By \ » e mo 2 & Figure I |. 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. of taxa were significantly different in MOW:HW (ANOVA, F-value = 3.701, p-val- ue < 0.05; Tukey test, p-value < 0.05), rufoflavus-Himalaya, and tahanensis-Himalaya. Discussion 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 distri- butional 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 in- clude 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 The oriental bumblebee Bombus flavescens 525 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 orig- inal 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. ro- tundiceps, 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 morpho- logical characters whereas GMYC does not. Then, we suggest that the PTP analysis cor- roborates 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-com- plex 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 clarify- ing 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 (Paabo 1989). It is difficult to extract the DNA from old specimens with standard protocols. However, numerous DNA extraction techniques for old in- sect 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 526 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) 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 Nyc- tomelitta (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 tahan- ensis, 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 popu- lations 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 bumble- bees; 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 noctur- nal 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 after- noon, 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 re- searchers (Musthafa et al. 2021) and by tourist services (e.g., Cameron Service: https:// cameronservice.blogspot.com/), but no bumblebees have been observed at traps. Con- The oriental bumblebee Bombus flavescens j27 sequently, 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 max- welli, 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. ve- lutina taxon nigrithorax (Perrard et al. 2014). For this reason, the colour patterns of B. flavescens might potentially be a result of Miillerian 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 ‘/epidus-group’ (Cameron et al. 2007). The groups originated in the Palearctic and Oriental regions, in the late Miocene (~11 Ma), and the /epidus-group is entirely restricted to the Oriental region (Williams 1998). In this study, biogeographic analyses show that the origin of the pra- torum-group is likely to be in the Oriental region, specifically around the Qinghai-Ti- betan 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, FE, 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. \n 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 re- gions 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 528 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) A 15Ma Bima C 0.5Ma D present Figure |2. 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 popu- lation 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 popula- tions 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). Taiwan-Luzon route in the same way as B. irisanensis, the other species in the Philip- pines. 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 sub- merged 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). Biogeography of bumblebees in Southeast Asia Many lineages of bumblebees in Southeast Asia might are hypothesised to have origi- nated around the QTP, then dispersing into the region via the Himalaya-Hengduan corridor (Williams 1985; Hines 2008; Williams et al. 2022). The bumblebee fauna The oriental bumblebee Bombus flavescens 529 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 di- verged in Southeast Asia approximately between 4.25—1.5 Ma, based on molecular dat- ing (Huang et al. 2015a). The lineages are 1) a trifasciatus-group: B. montivagus, B. albo- pleuralis 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. melan- opoda Cockerell, the phylogenetic relationship between mainland and Sumatran subge- nus 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 Melanobom- bus 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). Conclusion 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. flave- scens 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, 530 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) 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 conserva- tion of endemic populations of B. flavescens. Acknowledgements We would like to thank K. Atthasopa and K. Klaithin (ACMU), N. Badruddin (FRIM), E Bakker and W. van Bohemen (NMNL), N. Chatthanabun and P. Nalinra- chatakan (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. Srimaneey- anon 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. References 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. https://doi.org/10.26107/RBZ-2022-0004 Ashmead WH (1905) Additions to the recorded Hymenopterous fauna of the Philippine Is- lands, with descriptions of new species. Proceedings of the United States National Museum 28(1413): 957-971. https://doi.org/10.5479/si.00963801.28-1413.957 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. https://doi.org/10.1002/ece3.5756 Bischoff H (1936) Schwedisch-chinesische wissenschaftliche Expedition nach den nordwestli- chen Provinzen Chinas, unter Leitung von Dr. Sven Hedin und Prof. Sii Ping-chang. In- sekten gesammelt vom schwedischen Arzt der Expedition Dr. David Hummel 1927-1930. 56. Hymenoptera. 10. Bombinae. Arkiv for Zoologi 27: 1-27. Bouckaert R, Vaughan TG, Barido-Sottani J, Duchéne S, Fourment M, Gavryushkina A, Heled J, Jones G, Kiithnert D, De Maio N, Matschiner M, Mendes FK, Miiller NE, 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 plat- form for Bayesian evolutionary analysis. PLoS Computational Biology 15(4): e1006650. https://doi.org/10.1371/journal.pcbi. 1006650 The oriental bumblebee Bombus flavescens 531 Brasero N, Ghisbain G, Lecocq T, Michez D, Valterova I, Biella PR. Monfared A, Williams PH, Rasmont P, Martinet B (2021) Resovling the species status of overlooked West-Palaearctic bumblebees. Zoologica Scripta 50(5): 616-632. https://doi.org/10.1111/zsc. 12486 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 Archipel- ago. Annual Review of Ecology, Evolution, and Systematics 44(1): 411-435. https://doi. org/10.1146/annurev-ecolsys-110411-160323 Call E, Mayer C, Twort V, Dietz L, Wahlberg N, Espeland M (2021) Museomics: Phylog- enomics of the Moth Family Epicopeiidae (Lepidoptera) Using Target Enrichment. Insect Systematics and Diversity 5(2): 1-6. https://doi.org/10.1093/isd/ixaa021 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 nucleo- tide sequence data to phylogeny of the Hymenoptera: A review. Journal of Hymenoptera Research 1(1): 63-79. Cameron SA, Hines HM, Williams PH (2007) A comprehensive phylogeny of the bumble bees (Bombus). Biological Journal of the Linnean Society. Linnean Society of London 91(1): 161-188. https://doi.org/10.1111/j.1095-8312.2007.00784.x 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. Phil- ippine Journal of Science 16: 631-632. Cockerell TDA (1929) XXXII.—Descriptions and records of bees.—CXIX. Annals & Magazine of Natural History 4(21): 296-304. https://doi.org/10.1080/00222932908673058 Darriba D, Taboada GL, Doallo R, Posada D (2012) jModelTest 2: More models, new heuristics and parallel computing. Nature Methods 9(8): e772. https://doi.org/10.1038/nmeth.2109 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. https://doi.org/10.1007/s13592-019-00648-1 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 Re- search 79: 117-144. https://doi.org/10.3897/jhr.79.57308 Doria G (1886) A nocturnal Hymenoptera of the genus Bombus. Nature 33(852): e392. https:// doi.org/10.1038/033392a0 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: Apocri- tan relationships. Proceedings of the National Academy of Sciences of the United States of America 91(21): 9911-9915. https://doi.org/10.1073/pnas.91.21.9911 Ezard T, Fujisawa T, Barraclough TG (2009) splits: SPecies’ Llmits by Threshold Statistics [online]. http://R-Forge.R-project.org/projects/splits/ 532 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) 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 simu- lated data sets. Systematic Biology 62(5): 707-724. https://doi.org/10.1093/sysbio/ syt033 GBlEorg (2021) Data from: GBIF Occurrence, GBIE https://doi.org/10.15468/dl.fpw9zm [Accessed 31 January 2021] Ghisbain G, Martinet B, Wood TJ, Przybyla K, Cejas D, Gérard M, Rasmont P, Monfared A, Valterova I, Michez D (2021) A worthy conservation target? Revising the status of the rar- est bumblebee of Europe. Insect Conservation and Diversity 14(5): 661-674. https://doi. org/10.1111/icad.12500 Gilbert MTP, Moore W, Melchior L, Worobey M (2007) DNA extraction from dry museum beetles without conferring external morphological damage. PLoS ONE 2(3): e272. https:// doi.org/10.1371/journal.pone.0000272 Goodwin S, McPherson J, McCombie W (2016) Coming of age: Ten years of next-gener- ation sequencing technologies. Nature Reviews. Genetics 17(6): 333-351. https://doi. org/10.1038/nrg.2016.49 Gribodo G (1882) Alcune nuove specie e nuove genere di imenotteri aculeati. Annali del Mu- seo 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. https://doi.org/10.1073/pnas.0406166101 Hedicke H (1926) Beitrage zur Apidenfauna der Philippinen (Hym.). (2. Beitrag zur Kenntnis orientalischer Apiden.). Deutsche Entomologische Zeitschrift 1926: 413-423. https://doi. org/10.1002/mmnd.48019260506 Hines HM (2008) Historical Biogeography, Divergence Times, and Diversification Patterns of Bumble Bees (Hymenoptera: Apidae: Bombus). Systematic Biology 57(1): 58-75. https:// doi.org/10.1080/10635 150801898912 Hines HM, Williams PH (2012) Mimetic colour pattern evolution in the highly polymor- phic Bombus trifasciatus (Hymenoptera: Apidae) species complex and its comimics. Zoo- logical Journal of the Linnean Society 166(4): 805-826. https://doi.org/10.1111/j.1096- 3642.2012.00861.x 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. https://doi.org/10.1371/journal.pone.0132358 The oriental bumblebee Bombus flavescens 533 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. https://doi.org/10.1007/s13592-014-0319-9 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. https://doi.org/10.1093/ bioinformatics/btx025 Kapustijanskij A, Streinzer M, Paulus HE 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. https://doi.org/10.1111/j.1365-2435.2007.01329.x Kerfoot WB (1967) Correlation between Ocellar Size and the Foraging Activities of Bees (Hyme- noptera; Apoidea). American Naturalist 101(917): 65-70. https://doi.org/10.1086/282470 Koch JB, General DEM (2019) A preliminary assessment of bumble bee (Hymenoptera: Api- dae) habitat suitability across protected and unprotected areas in the Philippines. Annals of the Entomological Society of America 112(1): 44-49. https://doi.org/10.1093/aesa/say046 Kumar S, Stecher G, Tamura K (2016) MEGA7: Molecular evolutionary genetics analysis Ver- sion 7.0 for Bigger Datasets. Molecular Biology and Evolution 33(7): 1870-1874. https:// doi.org/10.1093/molbev/msw054 Lamm KS, Redelings BD (2009) Reconstructing ancestral ranges in historical biogeography: Properties and prospects. Journal of Systematics and Evolution 47(5): 369-382. https:// doi.org/10.1111/j.1759-6831.2009.00042.x 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. https://doi.org/10.1093/sysbio/syy011 Ma YP, Xie WJ, Sun WB, Marczewski T (2016) Strong reproductive isolation despite occa- sional hybridization between a widely distributed and a narrow endemic Rhododendron species. Scientific Reports 6(1): 19146. https://doi.org/10.1038/srep19146 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, Gen- tianaceae) and its close relatives. Journal of Biogeography 43(3): 580-590. https://doi. org/10.1111/jbi.12617 Matzke NJ (2013a) BioGeoBEARS: bioGeography with Bayesian (and likelihood) evolutionary analysis in R scripts. R package, version 0.2.1. [published July 27, 2013] http://CRAN.R- project.org/package=BioGeoBEARS Matzke NJ (2013b) Probabilistic historical biogeography: New models for founder-event spe- ciation, imperfect detection, and fossils allow improved accuracy and model-testing. Fron- tiers of Biogeography 5(4): 242-248. https://doi.org/10.21425/F55419694 Matzke NJ (2014) Model selection in historical biogeography reveals that founder-event spe- ciation is a crucial process in island clades. Systematic Biology 63(6): 951-970. https://doi. org/10.1093/sysbio/syu056 McElrath T (2022) Data from: Illinois Natural History Survey Insect Collection. Illinois Natu- ral History Survey, Occurrence dataset, GBIF. https://doi.org/10.15468/eol0pe [Accessed 7 October 2022] 534 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) Michener CD (2007) The bees of the world (2"4 edn.). John Hopkins University Press, Balti- more, 953 pp. https://doi.org/10.56021/9780801885730 Morley CK (2012) Late Cretaceous—Early Palaeogene tectonic development of SE Asia. Earth- Science Reviews 115(1—2): 37-75. https://doi.org/10.1016/j.earscirev.2012.08.002 Musthafa MM, Abdullah FE Martinez-Falcén AP, de Bruyn M (2021) How mountains and el- evations shape the spatial distribution of beetles in Peninsular Malaysia. Scientific Reports 11(1): e5791. https://doi.org/10.1038/s41598-021-84965-5 Nazari V, Schmidt BC, Prosser S$, Hebert PDN (2016) Century-old DNA barcodes reveal phylo- genetic placement of the extinct Jamaican Sunset moth, Urania sloanus Cramer (Lepidoptera: Uraniidae). PLoS ONE 11(10): 0164405. https://doi.org/10.1371/journal. pone.0 164405 Ogilvie HA, Bouckaert RR, Drummond AJ (2017) StarBEAST2 brings faster species tree infer- ence and accurate estimates of substitution rates. Molecular Biology and Evolution 34(8): 2101-2114. https://doi.org/10.1093/molbev/msx126 Ogle D, Doll J, Wheeler P, Dinno A (2022) FSA: Simple Fisheries Stock Assessment Methods. [online] https://cran.r-project.org/web/packages/FSA/index.html Paabo 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. https://doi.org/10.1073/pnas.86.6.1939 Papadopoulou A, Anastasiou I, Vogler AP (2010) Revisiting the Insect Mitochondrial Molecu- lar Clock: The Mid-Aegean Trench Calibration. Molecular Biology and Evolution 27(7): 1659-1672. https://doi.org/10.1093/molbev/msq05 1 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 E, Tan J, Bista S, Nugroho H, Baudoin R, Baylac M, Sil- vain 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. https://doi.org/10.1371/journal.pone.0094162 Pittioni B (1949) Beitraége zur Kenntnis der Bienenfauna SO-Chinas. Die Hummeln und Sch- marotzerhummeln 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. https://doi.org/10.1111/1755-0998.12474 R Core Team (2021) R: A Language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. https://www.R-project.org/ 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. https://doi.org/10.1093/gbe/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. https://doi. org/10.1093/sysbio/syy032 Rasmont P, Ghisbain G, Terzo M (2021) Bumblebees of Europe and neighbouring regions. N.A.P. Editions, Verriéres-le-Buisson, 631 pp. The oriental bumblebee Bombus flavescens 535 Ratnasingham S, Hebert PDN (2007) BARCODING, BOLD: The Barcode of Life Data System (www.barcodinglife.org). Molecular Ecology Notes 7(3): 355-364. https://doi. org/10.1111/j.1471-8286.2007.01678.x Ronquist FE, Huelsenbeck JP (2003) MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19(12): 1572-1574. https://doi.org/10.1093/bioinformatics/ bte180 Roubik DW (1989) Ecology and natural history of tropical bees. Cambridge University Press, Cambridge, 514 pp. https://doi.org/10.1017/CBO9780511574641 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. https://doi.org/10.1093/nar/gkab1112 Skorikov AS (1933) Zur Hummelfauna Japans und seiner Nachbarlander. 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 En- tomological Society of London 2(2): 33-45. https://doi.org/10.1111/j.1365-2311.1852. tb02208.x 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. https://doi.org/10.1111/1755-0998.12660 Starr CK (1989) Bombus folsomiand the origin of Philippine bumble bees (Hymenoptera: Apidae). Systematic Entomology 14(4): 411-415. https://doi-org/10.1111/j.1365-3113.1989. tb00294.x 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 sting- less bees (Apidae: Meliponini). Journal of Comparative Physiology. A, Neuroethology, Sensory, Neural, and Behavioral Physiology 202(9): 643-655. https://doi.org/10.1007/ s00359-016-1118-8 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. https://doi.org/10.3897/zookeys.85 1.32956 Thanoosing C (2017) Species status of Bombus (Alpinobombus) kluanensis and its close rela- tives—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. 536 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) Thomsen PF, Elias S, Gilbert MTP, Haile J, Munch K, Kuzmina S, Froese DG, Sher A, Holda- way RN, Willerslev E (2009) Non-destructive sampling of ancient insect DNA. PLoS ONE 4(4): e5048. https://doi.org/ 10.137 1/journal.pone.0005048 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. https:// doi.org/10.1093/nar/gks596 Voris HK (2000) Maps of Pleistocene sea levels in Southeast Asia: Shorelines, river systems and time durations. Journal of Biogeography 27(5): 1153-1167. https://doi.org/10.1046/ j.1365-2699.2000.00489.x 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. https://doi.org/10.1242/ jeb.015396 Warrant EJ, Kelber A, Wallén R, Wcislo WT (2006) Ocellar optics in nocturnal and diur- nal bees and wasps. Arthropod Structure & Development 35(4): 293-305. https://doi. org/10.1016/j.asd.2006.08.012 Wellington WG (1974) Bumblebee Ocelli and Navigation at Dusk. Science 183(4124): 550- 551. https://doi.org/10.1126/science.183.4124.550 Wickham H (2016) ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag, New York. https://doi.org/10.1007/978-3-319-24277-4 Williams PH (1985) A preliminary cladistic investigation of relationships among the bum- ble bees (Hymenoptera, Apidae). Systematic Entomology 10(2): 239-255. https://doi. org/10.1111/}.1365-3113.1985.tb00529.x Williams PH (1991) The bumble bees of the Kashmir Himalaya (Hymenoptera: Apidae, Bom- bini). 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 transi- tion zones using species turnover. Proceedings of the Royal Society B, Biological Sciences 263(1370): 579-588. https://doi.org/10.1098/rspb.1996.0087 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 sig- nificance for thermoregulation, crypsis, and warning mimicry. Biological Journal of the Linnean Society. Linnean Society of London 92(1): 97-118. https://doi.org/10.1111/ j-1095-8312.2007.00878.x 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). Jour- nal of Natural History 55(5—6): 271-282. https://doi.org/10.1080/00222933.2021.1900444 Williams PH (2022) The Bumblebees of the Himalaya. AbcTaxa, Belgium, 198 pp. Williams PH, Tang Y, Yao J, Cameron S (2009) The bumblebees of Sichuan (Hymenoptera: Apidae, Bombini). Systematics and Biodiversity 7(2): 101-190. https://doi.org/10.1017/ S$1477200008002843 The oriental bumblebee Bombus flavescens 537 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, Sima PB, Xie Z (2012) Unveiling cryptic species of the bumblebee subgenus Bombus s. str. worldwide with COI barcodes (Hymenoptera: Apidae). Systematics and Bio- diversity 10(1): 21-56. https://doi.org/10.1080/14772000.2012.664574 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. https://doi. org/10.1111/ecog.03074 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. https://doi.org/10.5852/ejt.2020.719.1107 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. https://doi.org/10.5852/ejt.2022.847.1981 Williams PH, Francoso E, Martinet B, Orr MC, Ren Z, Santos Junoir J, Thanoosing C, Van- dame 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. https://doi.org/10.1080/14772000.2022.2092229 Williams PH, Sung I-H, Lin Y-J, Lu S-S (2022c) Discovering endemic species among the bum- blebees of Taiwan (Apidae, genus Bombus). Journal of Natural History 56(5—8): 435-447. https://doi.org/10.1080/00222933.2022.2052991 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 remain- ing refugial-phase biodiversity. Biodiversity and Conservation 19(4): 919-941. https://doi. org/10.1007/s10531-010-9783-3 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. https://doi.org/10.1016/j-ympev.2010.04.011 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. https://doi.org/10.1186/s12870-015-0547-9 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. https://doi.org/10.1093/ molbev/msz257 538 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) Zhang J, Kapli P, Pavlidis PR, Stamatakis A (2013) A general species delimitation method with applications to phylogenetic placements. Bioinformatics 29(22): 2869-2876. https://doi. org/10.1093/bioinformatics/btt499 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 material | List of COI sequences from databases, research and collaboration, including the accession number or ID, the original sequence length and country Authors: Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams Data type: Sequence ID (word document) Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODDbL) 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. Link: https://doi.org/10.3897/jhr.96.104715.suppl1 Supplementary material 2 Primer testing Authors: 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. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODDbL) 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. Link: https://doi.org/10.3897/jhr.96.104715.suppl2 The oriental bumblebee Bombus flavescens D2? Supplementary material 3 Primer testing with gradient annealing temperatures Authors: 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, ** = me- dium, * = weak, NS = unsuccessful. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODDbL) 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. Link: https://doi.org/10.3897/jhr.96.104715.suppl3 Supplementary material 4 Polymerase chain reaction temperature profile of different pairs of primers in this study Authors: Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams Data type: Experimental (word document) Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODDbL) 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. Link: https://doi.org/10.3897/jhr.96.104715.suppl4 540 Chawatat Thanoosing et al. / Journal of Hymenoptera Research 96: 507-541 (2023) Supplementary material 5 List of 16S and PEPCK sequences from GenBank with the accession number, the original sequence length, and country Authors: Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams Data type: Sequence ID (word document) Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). 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. Link: https://doi.org/10.3897/jhr.96.104715.suppl5 Supplementary material 6 The accession number or ID of the parsing dataset of COI, 16S, and PEPCK se- quences for *BEAST analysis Authors: Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams Data type: Sequence ID (word document) Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODDbL) 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. Link: https://doi.org/10.3897/jhr.96.104715.suppl6 Supplementary material 7 Ancestral area reconstruction model selection, estimated in BioGeoBEARS Authors: Chawatat Thanoosing, Michael C. Orr, Natapot Warrit, Alfried P. Vogler, Paul H. Williams Data type: Model test (word document) Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODDbL) 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. Link: https://doi.org/10.3897/jhr.96.1047 15.suppl7 The oriental bumblebee Bombus flavescens 541 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) Authors: 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 = Chu- lalongkorn University Natural History Museum, Bangkok, ‘Thailand; KKIC = Ka- setsart 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 inter- tegular distance (ID), head width (HW), median ocelli width (MOW) in millime- tres (mm) (available at the NHM Data Portal: https://doi.org/10.5519/wknn9sd2). Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODDbL) 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. Link: https://doi.org/10.3897/jhr.96.1047 15.suppl8