Open Access
How to translate text using browser tools
31 December 2022 Diversification Pattern of the Widespread Holarctic Cuckoo Bumble Bee, Bombus flavidus (Hymenoptera: Apidae): The East Side Story
Patrick Lhomme, Sarah D. Williams, Guillaume Ghisbain, Baptiste Martinet, Maxence Gérard, Heather M. Hines
Author Affiliations +

Recent bumble bee declines have made it increasingly important to resolve the status of contentious species for conservation purposes. Some of the taxa found to be threatened are the often rare socially parasitic bumble bees. Among these, the socially parasitic bumble bee, Bombus flavidus Eversmann, has uncertain species status. Although multiple separate species allied with B. flavidus have been suggested, until recently, recognition of two species, a Nearctic Bombus fernaldae (Franklin) and Palearctic B. flavidus, was favored. Limited genetic data, however, suggested that even these could be a single widespread species. We addressed the species status of this lineage using an integrative taxonomic approach, combining cytochrome oxidase I (COI) and nuclear sequencing, wing morphometrics, and secretions used for mate attraction, and explored patterns of color polymorphism that have previously confounded taxonomy in this lineage. Our results support the conspecificity of fernaldae and flavidus; however, we revealed a distinct population within this broader species confined to eastern North America. This makes the distribution of the social parasite B. flavidus the broadest of any bumble bee, broader than the known distribution of any nonparasitic bumble bee species. Color polymorphisms are retained across the range of the species, but may be influenced by local mimicry complexes. Following these results, B. flavidus Eversmann, 1852 is synonymized with Bombus fernaldae (Franklin, 1911) syn. nov. and a subspecific status, Bombus flavidus appalachiensis ssp. nov., is assigned to the lineage ranging from the Appalachians to the eastern boreal regions of the United States and far southeastern Canada.

In the last several decades, bumble bees (Hymenoptera, Apidae, Bombus Latreille) have been observed to be declining in many regions across the globe, with factors driving these declines largely attributed to habitat modification, pathogens, and climate change (Cameron et al. 2011, Goulson et al. 2015, Cameron and Sadd 2020). Bumble bees have a rich taxonomic history (Williams 1998); however, as a result of low morphological variation combined with high color variation, many bumble bee species retain contentious species status (Williams et al. 2020). Recent studies testing species hypothesis have revealed interesting surprises, such as cryptic diversity in the arctic (Potapov et al. 2018, Martinet et al. 2019), speciation hidden by mimicry (Ghisbain et al. 2020a), and the unreliability of color characters (e.g., Carolan et al. 2012, Koch et al. 2018, Ghisbain et al. 2020a). Given recent threats, the need to conserve these bees makes defining species ever more imperative (Ghisbain et al. 2020b).

Among the bumble bee species noted to be declining or threatened are the often rare bumble bees in the subgenus Psithyrus Lepeletier 1832 (Hymenoptera, Apidae) (Colla and Packer 2008, Colla et al. 2012, Suhonen et al. 2015). Psithyrus is a monophyletic clade (28 species worldwide, Williams 1998, Lhomme and Hines 2019) that is unique among bumble bees in consisting exclusively of social parasites (cuckoo bumble bees) of other bumble bees (reviewed in Lhomme and Hines 2019). These social parasites lack a worker caste and have reduced wax glands (Sramkova and Ayasse 2008) and are thus unable to found their own nests. They instead usurp the colonies of other bumble bee species, functionally replacing the host queen and exploiting the worker force (Lhomme and Hines 2018). To do so, cuckoo bumble bees have evolved several chemical and behavioral adaptations to overcome the recognition system of their host, forcing the host workers to feed and care for the parasitic offspring (Kreuter et al. 2012; Lhomme et al. 2012, 2013, 2015; Ayasse and Jarau 2014). They are morphologically distinct in their lack of pollen-collecting corbiculae and display additional defensive adaptations, such as a thickened cuticle and strengthened mandibles, which allow them to fight to enter nests successfully (Fisher and Sampson 1992).

Although Psithyrus are easy to recognize from other bumble bees, their identification at the species level is much more challenging. The variability in size and color patterns in Psithyrus (Reinig 1935, Williams 2008, Lecocq et al. 2011), morphological similarities among the species, and the natural rarity of many species led to many specific and subspecific names in the early years of bumble bee taxonomy that have since been synonymized (Williams 1998). Although 28 species of Psithyrus are currently recognized, more than 350 specific and subspecific names have been proposed (Williams 1998). Nevertheless, taxonomic revisions within the subgenus are scarce (Lecocq et al. 2011) in comparison to other widespread sub-genera (e.g., Williams et al. 2012, 2019, 2020).

Molecular methods have been used to investigate bumble bee species diversity and species delineation (Hines et al. 2006; Cameron et al. 2007; Duennes et al. 2012; Williams et al. 2012, 2015; Brasero et al. 2020). Mitochondrial cytochrome oxidase I (COI) barcoding is a primary tool used to discriminate cryptic species among bumble bees (e.g., Bertsch 2009; Williams et al. 2012, 2019; Ghisbain et al. 2020a); however, reliance solely on the COI barcode can be misleading as the level of divergence required to separate species can greatly differ depending on species group or ecology (Soltani et al. 2017), and introgression, selection, sex-biased histories, numts, and sequence divergence resulting from historical isolation prior to admixture can lead to difficulties in taxonomic interpretations at the species level (Williams et al. 2012, 2019, 2020; Cong et al. 2017; Gueuning et al. 2020). Non-molecular methods such as wing geometric morphometrics have been used to separate bumble bee species (Aytekin et al. 2007, Duennes et al. 2017); however, they have limitations for discriminating between closely related or cryptic species (Lecocq et al 2015a, Gérard et al. 2020). Male cephalic labial gland secretions (CLGS) have also been shown to be useful in separating bumble bee species (Lecocq et al. 2011; Lecocq 2015a,b,c; Martinet et al. 2018, 2019; Brasero et al. 2018, 2020). CLGS play a major role in premating isolation and thus are highly species-specific (reviewed in Valterová et al. 2019); however, it can be difficult to demarcate what level of differences in eco-chemical traits reflect species-level divergence as opposed to subspecific pheromonal dialects (Lecocq et al. 2015b, Brasero et al. 2020). For these reasons, it is best to follow a multi-evidence approach in delimiting species (Dayrat 2005, Schlick-Steiner et al. 2010).

This study aims to shed light on the taxonomic position of two closely related Psithyrus taxa: Bombus flavidus Eversmann, 1852 and Bombus fernaldae (Franklin, 1911). Bombus flavidus was originally described by Eversmann (1852) from specimens of Irkutsk (Russia). The species was first described as a Bombus, even though at that time Bombus and Psithyrus were considered as different genera. The first description of B. flavidus as a cuckoo bumble bee was made by Thompson (1872) (as Apathus lissonurus) from specimens of Lappland (Sweden). Bombus flavidus and Psithyrus lissonurus were considered separate species until Popov (1931) synonymized both taxa. Up to 26 different names (species, subspecies, varieties, and forms; listed in the type material revision section in the results) were attributed to B. flavidus, a problem magnified by the extensive color variability of this species. Bombus flavidus has a broad Palearctic distribution (Fig. 1) spanning from Scandinavia to far eastern Russia (Proshchalykin 2007) with southern limits of the Pyrenees and the Alps in Europe (Rasmont et al. 2015). It is mainly found in sub-alpine/subarctic northern habitats, being primarily distributed across boreal forest/taiga (Løken 1984). Bombus flavidus is reported to be a social parasite of several bumble bee species belonging to the subgenus Pyrobombus (Hymenoptera, Apidae) (Pekkarinen et al. 1981, Rasmont 1988, Pekkarinen and Teräs 1993, Lhomme and Hines 2019), but has been recorded breeding only in one nest of B. (Pyrobombus) jonellus (Kirby) (Brinck and Wingstrand 1951).

Bombus fernaldae was described by Franklin (1911) from specimens of New England, Alaska, Washington (United States), and British Columbia (Canada). Other historically described Nearctic species have also been synonymized with B. fernaldae, including B. tricolor which was determined to be the male of B. fernaldae, and Bombus wheeleri from Oregon and California (synonymized by Frison 1926). Bombus fernaldae has been recognized as the Nearctic sister lineage to B. flavidus given morphological similarities. Its distribution (Fig. 1) is largely restricted to the boreal zone of the Nearctic north, ranging from Alaska across much of Canada, south in the American West through the Rocky, Cascades, and Sierra Nevada mountains, and in the east into the far southern extents of the Appalachian Mountains. Its host is currently unknown as there are no direct observations of this species breeding in host colonies (Williams 2008, Koch et al. 2014, Lhomme and Hines 2019), but main hosts are argued to likely be one or more species of subgenus Pyrobombus, as it is mainly found in geographic association with these species (Wilson et al. 2010) and it is in a Psithyrus subgroup known to be specialized on Pyrobombus hosts (Lhomme and Hines 2019). Bombus fernaldae has been found in nests of B. (Cullumanobombus) rufocinctus Cresson (Hobbs 1965), B. (Bombus) occidentalis Greene (Hobbs 1968), B. (Subterraneobombus) appositus Cresson (Hobbs 1966), and several Pyrobombus species (e.g., Bombus perplexus Cresson; Hobbs 1967, Laverty and Harder 1988). Although this does not necessarily translate to usurpation and control of the nest (they occasionally shelter in nests that they will not usurp), this is suggestive of B. fernaldae being a more generalized social parasite.

Bombus flavidus and B. fernaldae have been considered as distinct species (e.g., Williams 1998) until recent genetic data came to light. Genetic analysis of a single specimen from Scandinavia and a specimen from the western United States using four nuclear markers and one mitochondrial marker showed these two representatives to be nearly identical in sequences, suggestive of the two lineages being conspecific (Cameron et al. 2007). As a result of these findings, New World B. fernaldae has tentatively been synonymized as B. flavidus awaiting further information (Williams et al. 2014). Nevertheless, several authors still recognize B. fernaldae as a valid species (Jacobson et al. 2018, Richardson et al. 2018, Cole et al. 2020). If both taxa are conspecific, B. flavidus would exhibit an unusually large range for bumble bees, including a broad Holarctic distribution spanning across Palearctic, Western, and Eastern Nearctic regions. This is a broader range that has been found in any other bumble bee (Hines 2008). For the purpose of clarity, we hereafter refer to the taxa fernaldae and flavidus as the Nearctic and Palearctic lineages, respectively.

Fig. 1.

Global distribution of Bombus flavidus (Palearctic) and Bombus fernaldae (Nearctic). Distributional ranges are approximated based on literature records (Williams et al. 2014; Ascher and Pickering 2020) and thus are not meant to be accurate at a fine scale.


In this article, we seek to resolve the taxonomic relationships between B. fernaldae and B. flavidus through applying an integrative taxonomic approach that combines several data types including semiochemical, genetic, and morphometric evidence. We further explore patterns of color polymorphism that have previously confounded taxonomy in these two taxa. These data refine understanding of the biogeography, evolution of parasitic life-history traits, and conservation status of these clades.

Material and Methods


The natural rarity of the species limited the total number of specimens sampled; however, we were able to examine specimens of the Nearctic (4 females, 149 males) and Palearctic (23 males) lineages spanning the Holarctic, including from across North America (Alaska, Western and Eastern Canada, and United States), Scandinavia, and sites in Central and Far Eastern Russia. Specimens from these three areas were compared on the basis of COI sequences and male color patterns and wing morphometrics, and specimens from North America and Scandinavia were additionally compared for male CLGS. The full collection of examined specimens of both Nearctic and Palearctic lineages and associated traits investigated can be downloaded online as  Supp Table 1 (online only).

Genetic Analysis: COI

We constructed a COI Bayesian phylogeny of 62 specimens of both the Nearctic and Palearctic lineages, including 33 newly sequenced individuals spanning its range, and 29 sequences from North America obtained from GenBank (;  Supp Table 1 [online only]; sequence alignments available at Dryad data repository doi: 10.5061/dryad.cvdncjt3t). To examine the placement of both lineages in relation to their closest known extant relatives, we included exemplar COI sequences of the latter selected based on the comprehensive phylogeny of bumble bees of Cameron et al. (2007): Bombus barbutellus (Kirby) (n = 4), Bombus bohemicus Seidl (n = 4), Bombus norvegicus (Sparre-Schneider) (n = 7), Bombus quadricolor (Lepeletier) (n = 3), Bombus skorikovi (Popov) (n = 1), Bombus sylvestris (Lepeletier) (n = 6), and Bombus vestalis (Geoffroy) (n = 4) along with the slightly more distant Bombus citrinus (Smith) (n = 1) as an outgroup, using both newly sequenced (n = 5) individuals and sequences extracted from GenBank (n = 25) ( Supp Table 1 [online only]; sequence alignments available at Dryad data repository doi: 10.5061/dryad.cvdncjt3t).

For newly sequenced specimens, either thoracic muscle or leg tissues were extracted from specimens and DNA was isolated from the samples following the QIAmp DNA Micro kit or Omega Bio-tek E.Z.N.A. Tissue DNA kit protocols (Omega Bio-tek, Norcross, GA). Polymerase chain reactions (PCR) were then used to amplify the mitochondrial barcoding gene COI for each of the isolated DNA samples. Primers utilized included LepF1 and LepR1 (Herbert et al. 2004) or HCO2198 and LCO1490 primers (Folmer et al. 1994), which amplify the classic barcode fragment used across animals. Each PCR included 0.3 µl each 10 µM primer, 6.9 µl distilled water + DNA (1–2.5 µl), and 7.5 µl AMRESCO Taq polymerase master mix, and the PCR conditions included 2 min at 94°C, an initial cycle of 30 s at 94°C, 40 s at 45°C, 1 min at 72°C, then 35 cycles of 30 s at 94°C, 40 s at 49°C, 1 min at 72°C, and lastly 10 min at 72°C.

PCR products were then purified using ExoSAP (Thermo Fisher Scientific) and sequenced at the Penn State Genomics Core Facility using single-end sequencing. Sequences were edited and aligned using Geneious 2 software, and all single-nucleotide polymorphism (SNP) variants were manually double-checked against chromatograms. COI sequences were edited manually and aligned in Geneious v8.1.9 (Biomatters,, with ends trimmed to remove missing data to yield a 559-bp aligned fragment.

We performed Bayesian inference on COI sequences using the GTR model, the preferred model under the Bayesian information criterion inferred from MEGA X (v10.0.5; Kumar et al. 2018), performed in MrBayes 3.2.6 analyzed in the CIPRES Science Gateway v3.3 servers ( with four independent runs, four chains, 15 million generations and sampling trees every 1,000 generations. Convergence was assessed using likelihood plots (for stationarity) and convergence statistics in MrBayes, and ESS values in Tracer 1.7.1 (Rambaut et al. 2018), which led us to discard the first 25% of the trees as burn-in for all runs to obtain a majority-rule 50% consensus tree.

To examine population-level relationships without forcing bifurcating topologies, a TCS haplotype network (Clement et al. 2002) was generated using the program popART (Leigh and Bryant 2015) using inferred haplotypes from a subset of the COI sequences above that were free of missing data (n = 65 individuals, including 61 B. flavidus s.l. and 4 B. norvegicus). Close relative B. norvegicus was used as an outgroup to root the network. These were trimmed to a length of 601 bp to avoid missing data.

We used a Bayesian implementation of the general mixed Yule-coalescent model (bGMYC) for species delimitation based on the COI phylogenetic tree (Reid and Carstens 2012; see examples of this approach on bumble bees in Williams et al. 2015 and Martinet et al. 2018, 2019). We performed these analyses with the R package ‘bGMYC’ (Reid and Carstens 2012). Because the bGMYC requires an ultrametric tree, we performed a phylogenetic analysis with BEAST 1.8. (Drummond and Rambaut 2007), with Markov chain Monte Carlo (MCMC) = 100 million generations with the first million trees as burn-in, using the maximum clade credibility method and setting the posterior probability limit to 0. A range of probabilities > 0.5 was considered as strong evidence that taxa were conspecific, whereas a range of probabilities < 0.05 suggested that taxa were heterospecific (Reid and Carstens 2012). The range of 0.05–0.5 corresponds to an intermediate zone of uncertainty.

Genetic Analysis: Nuclear Data

Previous nuclear data comparing a specimen each from New Mexico and Sweden suggested that both Nearctic and Palearctic lineages did not differ much in sequence (Cameron et al. 2007; 2 total bp difference between fragments of four genes ArgK [1 bp], PEPCK [0 bp], EF-1a [1 bp], and opsin [0 bp]). To further explore whether different clades inferred from mitochondrial data showed signs of nuclear gene sequence differences, we obtained additional sequence data for the PEPCK gene, a gene previously used to aid species delimitation decisions (e.g., Brasero et al. 2020), including 11 newly sequenced individuals spanning the inferred clades from COI and including the two closest relatives as outgroups, B. skorikovi and B. norvegicus, using sequences for these available on GenBank ( Supp Table 1 [online only]; sequence alignments available at Dryad data repository doi: 10.5061/dryad.cvdncjt3t). We also sequenced taxa from across the clades for the nuclear Internal Transcribed Spacer (ITS) gene, which lies between the large and small subunit rRNA, a gene often used for assessing population-level variation and species delimitation due to its high levels of variation (Hines and Williams 2012). For this we sequenced 10 individuals and the same two outgroup species using the vouchers from Cameron et al. (2007) used to obtain the PEPCK sequences from GenBank ( Supp Table 1 [online only]; sequence alignments available at Dryad data repository doi: 10.5061/dryad.cvdncjt3t). The PEPCK fragment—an ∼847 bp fragment that included both intron and exon sequences and a single indel—used primers FHv4 and RHv4 and respective thermocycler conditions from Cameron et al. (2007). ITS primers ITS1-HHRInt and CAS18SFA (see thermocycler conditions in Hines and Williams 2012) amplified an ∼810 bp fragment with two indels. PCR reagents, purification, and sequencing followed the same procedures as for COI.

Semiochemical Analysis

We examined chemical composition of the most studied reproductive trait involved in the bumble bee premating recognition, the male CLGS (Valterová et al. 2019). CLGS consist of a complex mixture of compounds (mainly aliphatic or isoprenoid), with variable main compounds (Valterová et al. 2019). We extracted CLGS in 400 µl of n-heptane derived from the method described by De Meulemeester et al. (2011). Samples were stored at –40°C prior to the analyses. We extracted CLGS of 23 specimens of the Nearctic lineage from Alaska (n = 4), Yukon (n = 5), Wyoming (n = 1), and Pennsylvania (n = 13) and 11 specimens of the Palearctic lineage from Siberia (n = 8) and Sweden (n = 3;  Supp Table 1 [online only]; data available at Dryad data repository doi: 10.5061/dryad.cvdncjt3t).

The qualitative composition of the CLGS was determined by gas chromatography–mass spectrometry using a Finigan GCQ quadrupole system (GC/MS) with a nonpolar DB-5 ms capillary column (5% phenyl [methyl] polysiloxane stationary phase; column length 30 m; inner diameter 0.25 mm; film thickness 0.25 µm). All samples of CLGS were quantified with a gas chromatograph Shimadzu GC-2010 system (GC-FID) equipped with a nonpolar SLB-5 ms capillary column (5% phenyl [methyl] polysiloxane stationary phase; column length 30 m; inner diameter 0.25 mm; film thickness 0.25 µm) and a flame ionization detector. The composition of CLGS was analyzed according to the protocol used in Martinet et al. (2019). All compounds for which the relative abundance was recorded as less than 0.1% for all specimens were excluded from the analysis (De Meulemeester et al. 2011). The data matrix for each taxon was based on the alignment of each relative proportion of compound between all samples performed with GCAligner 1.0 (Dellicour and Lecocq 2013a,b). To facilitate the alignment of compounds and the identification, before each sample injection, a standard mixture of alkenes (Kovats) from C10 (decane) to C40 (tetracontane) was injected. We calculated Kovats indices with GCKovats 1.0 according to the method described by Dellicour and Lecocq (2013a,b). We performed statistical comparative analyses of the CLGS using R 3.5.1 (R Development Core Team 2020) to detect CLGS differentiations. We transformed [log (x + 1)] and standardized data to reduce the great difference of abundance between highly and lowly concentrated compounds. We used a clustering method computed with the unweighted pair-group method with average linkage (UPGMA) based on correlation distance matrices (relative abundance of each compound; R package ape; Paradis et al. 2004). We assessed the uncertainty in hierarchical cluster analysis using P-values calculated by multiscale bootstrap resampling with 1,000,000 bootstrap replications (significant branch supports > 0.85; R package pvclust, Suzuki and Shimodaira 2011). We also assessed CLGS differentiations between taxa by performing a multiple response permutation procedure (MRPP, R package vegan, Oksanen et al. 2020) based on groups identified by hierarchical cluster analysis. When a significant difference was detected, pairwise multiple comparisons were performed with an adjustment of P-values (Bonferroni correction) to avoid type I errors.

Geometric Morphometrics Analysis

We performed geometric morphometric analysis of wing venation on 103 total male specimens including 94 males of the Nearctic lineage from Oregon (n = 78), Alaska (n = 2), Pennsylvania (n = 12), and Maine (n = 2) and 9 specimens of the Palearctic lineage from Scandinavia (n = 5) and Russia ( Supp Table 1 [online only]; data available at Dryad data repository doi: 10.5061/dryad.cvdncjt3t). As 20 specimens per group is usually recommended to get the full statistical power of geometric morphometric analyses (Dewulf et al. 2014), the present results of geometric morphometrics will have some limitations in interpretation.

We took pictures of the right forewings with the Camera Canon 5D markII (1× magnification), uploaded pictures into tps series software (tpsUTIL 1.74; Rohlf 2013a), and digitized the 18 homologous landmarks ( Supp Fig. 2 [online only]) on the veins and cells of the bumble bee wings (Owen 2012) with 2D Cartesian coordinates (tpsDIG 2.30; Rohlf 2013b). We then superimposed these raw coordinates to align the different landmark configurations using GLS Procrustes superimposition on R version 3.6.1. (Rohlf and Slice 1990, R Development Core Team 2020).

Principal component analysis (PCA) was used to visualize the landmark configuration variability as well as potential clustering and outliers. Linear discriminant analysis (LDA) was then used to assess the potential discriminating power of the wing shape and maximize the interlineage differences. We then assessed the effectiveness of the LDA to identify given taxa using a leave-one-out (LOO) cross-validation procedure. This method assesses the percentage of true positive attributions (i.e., individuals correctly assigned to their own taxon) based on the posterior probabilities of assignment (pp). The posterior probability of assignment represents the probability of an unknown individual being assigned to a given group, compared with all the other groups. The unknown individual is then attributed to the group displaying the highest pp (Huberty and Olejnik 2006, Dehon et al. 2019). Finally, we assessed wing size differences between groups by measuring the centroid size of the right forewing (i.e., square root of the sum of squared distance between all landmarks and the centroid; Gérard et al. 2018). We used type one analysis of variance (ANOVA) to detect differences in wing size between the different taxa. Post hoc pairwise comparisons were made using Tukey's HSD test (R function glht from the R package multcomp; Hothorn et al. 2008).

Color Pattern Analyses

To examine effects of local environment on color pattern, we analyzed a total of 77 males of the Nearctic lineage collected from 16 different localities across Oregon (United States;  Supp Table 1 [online only]). We chose males because they present an especially broad range of color patterns across their range. These specimens clustered into two major regions: the Cascade Mountain Range in the west and the mountain ranges in the Northeast corner of the state. To quantify differences in color pattern across the state, each specimen was observed under a microscope and their individual color patterns were recorded for the head, dorsal thorax, pleuron (lateral thorax), and abdomen. The color patterns were partitioned into 68 total squares for each individual ( Supp Fig. 3 [online only]). Each square was given a score representing the color percentage of the pile (black, yellow, or orange). The percentages of each color were summed, and each individual was given a total body pile percent black, yellow, and orange (data available at Dryad data repository doi: 10.5061/dryad.cvdncjt3t). An unpaired student t-test was performed on percent of yellow color from Western versus Eastern regions (data are normally distributed; Shapiro–Wilk test, P = 0.29). The percentages of yellow coloration for each cluster of localities were mapped to assess the role of geography in coloration. We also performed a regression analysis on the overall body percentage of yellow pile for each individual against elevation to assess whether it was a factor driving variation of color patterns in the Nearctic lineage.

To examine broad variability in color pattern across the clade and how it might contribute to population/species delineation patterns, we analyzed color pattern variation in both the Nearctic and Palearctic lineages using numerous specimens across the range of both taxa ( Supp Table 1 [online only]). Color patterns were recorded using the same method as above. After recording color variation, the extremes and presence of intermediates were noted for each sampled area, and representatives of the range of color variation were plotted on a map of global distribution. Color patterns from Scandinavian specimens were extracted from Løken (1985).

Type Material Examined

We located and examined all type material of B. fernaldae declared by Franklin which includes 17 cotypes (all females), of which 8 were deposited in the United States National Museum (USNM, Smithsonian, Washington, DC), 5 were deposited in University of Massachusetts (UMASS, Amherst, MA; since transferred to the USNM), and 4 were deposited in the University of New Hampshire Insect Collection (UNHC, Durham, NC). One cotype, a specimen from Alaska from USNM, was unofficially declared a lectotype (unpublished, by label only) by H. E. Milliron, 1959.

We also located and examined all type material of B. flavidus declared by Eversmann consisting of two male specimens deposited at the Zoological Institute of Russian Academy of Sciences (ZISP, Saint-Petersburg, Russia). One specimen was later declared a lectotype and the other one was declared a paralectotype (unpublished, by label only) by Yu. A. Pesenko, 1988.

Species Delineation Framework

Following the unified theoretical concept of de Queiroz (2007) that considers species as independently evolving lineages, we assess here support for species using an approach integrating multiple evidence of divergence based on genetic, chemical, and morphometric traits. In this context, in considering species status we assess whether lineages are 1) genetically distinct for COI (bGMYC, P < 0.05), and nuclear genes 2) monophyletic, 3) significantly differentiated in its semiochemical traits (MRPP, P < 0.05; bootstrap value > 0.85) with different main compounds, and 4) differentiated in wing shape (LOO cross-validation procedure). These patterns will be considered together to make decisions about species status.


This article and the nomenclatural act(s) it contains have been registered in Zoobank (, the official register of the International Commission on Zoological Nomenclature. The LSID (Life Science Identifier) number of the publicationis:urn:lsid:zoobank. org:pub:96311861-3558-4B27-AEB4-6541970CAB34.


Mitochondrial Phylogeny and Haplotype Network

The Bayesian inference of mitochondrial COI showed the most distinct divergence between two lineages: an Eastern Nearctic lineage distributed across the Appalachians northward into the southeastern boreal zone of Canada and the northeastern United States, and a lineage comprised the remaining western Nearctic and all Palearctic specimens (Fig. 2A and B). The haplotype network suggests the same two larger haplotype groups (Fig. 2A), with nine unique SNPs separating the two broader groups (1.7% COI divergence) and only a few SNPs separating individuals within these groups. These two lineages were found to meet in northern Michigan, United States, where both mitochondrial haplotypes from both lineages were found in the western Upper Peninsula, but were otherwise not found to be admixed.

The haplotype network supports 28 nucleotide changes separating our ingroup from the clearly distinct sibling lineage B. norvegicus (in black; Fig. 2A; 4.7% divergence). In both phylogeny and network, Palearctic + Western Nearctic lineage is more closely related to the outgroup, whereas the Eastern Nearctic lineage is a group derived from the Palearctic + Western Nearctic lineage.

Fig. 2.

Relationships of Bombus flavidus lineages inferred with cytochrome oxidase I (COI). (A) Haplotype network for the sequenced individuals, highlighting the connections between the haplotypes shown on the map (right). (B) Bayesian tree of B. flavidus and related species. Clade support values are Bayesian posterior probabilities.The right zoomed-in part of the figure highlights population structure. (C)Visualization of the major haplotypes mapped at their respective geographical locations.


The bGMYC analysis (Fig. 3) supports homospecificity within each of the Palearctic + Western Nearctic and Eastern Nearctic lineages. In respect to each other however, they do not pass the heterospecificity threshold of P < 0.05, but instead fall in an intermediate zone of uncertainty (P = 0.05–0.5). The analysis suggests a more conservative delimitation of eight prospective ingroup species in our complete data set: B. barbutellus, B. bohemicus, B. norvegicus, B. quadricolor, B. skorikovi, B. sylvestris, B. vestalis, and a single species including both Nearctic and Palearctic lineages of our ingroup (P < 0.05).

Fig. 3.

Bayesian general mixed Yule-coalescent model (bGMYC) results based on COI phylogenetic tree. The colored matrix corresponds to the pairwise probabilities of conspecificity (color scale below the figure).


Nuclear Data

The nuclear gene PEPCK contained no SNP variation among all our specimens sampled across the Holarctic (except for a single specimen from Pennsylvania), confirming previous results that there is no nuclear variation to distinguish these lineages in this gene. For the ITS marker, variation was found in only a single SNP, with Oregon specimens and a W. Siberian specimen having one SNP that matched the outgroups (Fig. 4).

Chemical Analyses

For both taxa, 112 compounds were identified in the pheromonal bouquet. Although the same CLGS compounds were found across samples including main compounds (hexadecen-1-ol, ethyl tetradecanoate, ethyl tetradecenoate), these differed in relative amount. Our chemical analyses mirrored the analysis based on COI: the Eastern Nearctic samples clustered separately from the Western Nearctic + Palearctic specimens; however, the extent of these differences was not statistically significant (MRPP, A = 0.09, P = 0.11; Fig. 5).

Wing Morphometric Analyses

The Nearctic lineage could not be separated from the Palearctic lineage based on wing shape. In the PCA, the three different groups (Palearctic, Western Nearctic, Eastern Nearctic) have slight to strong overlaps (Fig. 6A), where as the LDA shows fairly good discrimination of the groups (Fig. 6B). To further explore the discriminating power of wing morphometrics, the leave-one-out cross-validation procedure was applied to show a high hit-ratio of specimens from the Western Nearctic (HR = 92.5%), whereas specimens from Eastern America and Palearctic displayed very low hit-ratios (HR = 50% and 33.3%, respectively;  Supp Table 4 [online only]). The most common misidentification was the attribution of Eastern Nearctic and Palearctic specimens to the Western Nearctic lineage (4 and 5 individuals, respectively). The low discriminating power for two of the lineages as well as the high proportion of misattribution to Western Nearctic taxa could be the consequence of the unbalanced groups (i.e., many more specimens from Western Nearctic than from the two other areas).

Fig. 4.

Median-joining network of haplotypes based on nuclear genes ITS and PEPCK. Circle sizes are proportional to frequencies of haplotypes. Colors of haplotypes refer to geographic areas (Fig. 3). Black bars on lines represent the number of mutation(s) between two close haplotypes. Black circles are outgroups (B. skorikovi and B. norvegicus).


Fig. 5.

Dendrogram based on cephalic labial gland secretions of Bombus flavidus with represented Eastern Nearctic (purple) and Palearctic + western Nearctic (green) lineages. This cluster was obtained by hierarchical clustering using an unweighted pair-group method with arithmetic mean (UPGMA) based on correlation distance matrices calculated from the relative abundance of the compounds present in cephalic labial gland secretions (CLGS). The values near nodes represent multiscale bootstrap resampling values. Multiple response permutation procedure (MRPP) indicates level of CLGS differentiation between taxa.


One-way ANOVA on centroid size revealed significant differences between the Eastern Nearctic, Western Nearctic, and Palearctic specimens concerning wing size (P < 0.001). Although no differences were found between Palearctic specimens and the other two groups, specimens from the Eastern Nearctic group were significantly smaller in wing size, a proxy for body size (Gérard et al. 2015), than specimens from the Western Nearctic (Tukey's HSD test; P < 0.001; Fig. 7).

Color Data

On a global scale, males from both the Nearctic and Palearctic lineages are highly color polymorphic. Although all specimens are yellow in the first thoracic and fourth abdominal segments, specimens exhibited remarkable variation in the degree of yellow as opposed to black coloration in the remainder of tergal and pleural segments (Fig. 8). This variation appeared to be continuous, and the full range from mostly yellow to mostly black in the color variable segments often occurs within a given geographic region. The amount of yellow fluctuates across the Holarctic and is not sorted by mitochondrial haplotype in northern Michigan where both COI mitochondrial lineages are found in sympatry, which has mostly black and mostly yellow individuals.

Fig. 6.

(A) Ordination of individuals of Bombus flavidus along the first two principal components (PCA1 and PCA2) of the principal components analysis of the landmark configuration variability. (B) Ordination of individuals of Bombus flavidus along the two first axes (LDA1 and LDA2) of the linear discriminant analysis of the landmark configuration variability.


Fig. 7.

Wing size differences between Palearctic, Eastern Nearctic, and Western Nearctic specimens of Bombus flavidus assessed by comparing centroid size. Box plots show the median and 25–75% percentiles. Whiskers show all data excluding outliers. Outliers (dots) are values being more than 1.5 times box length from upper and lower edge of respective box. The different letters indicate significant differences in centroid size between populations (type one ANOVA,Tukey's HSD test, P < 0.05).


Although the global analysis revealed high variability even within regions, our examination of regional patterns in Oregon revealed local partitioning of the degree of yellow-black coloration (Fig. 9A). Specimens collected from more western locations have a lower total average percent of yellow pile than specimens from more eastern locations (unpaired two-tailed t-test, df = 53.691, t = 6.0509, P < 0.0001; Fig. 9B).

Type Material Revision

Bombus fernaldae (Franklin, 1911)

  • Apathus insularis Cresson, 1863: 113

  • Psithyrus insularis Ashmead, 1902: 130

  • Psithyrus fernaldae Franklin, 1911: 164

  • Psithyrus tricolor Franklin, 1911: 167

  • Psithyrus wheeleri Bequaert and Plath, 1925: 265

  • Lectotype (USNM). United States, Alaska (Sitka), lectotype ♀ by designation of Milliron, 1959 (by label only). This lectotype has not been published and is accepted here.

  • Cotypes (16 Specimens: 12 in USNM and 4 in UNHC). United States, Alaska (Sitka, Nushagak), Washington (Mt. Rainier), New York (Ithaca), New Hampshire (Webster, Conway, Crawford, Durham); Canada, British Columbia (Kaslo, Metlakatla). Cotypes, all ♀, by original designation of Franklin (1911).

  • Bombus flavidus Eversmann, 1852

  • Bombus autumnalis Zetterstedt, 1838: 474

  • Bombus flavidus Eversmann, 1852: 131

  • Apathus lissonurus Thomson, 1872: 49

  • Psithyrus quadricolor var. lissonurus Pérez, 1884: 266

  • Psithyrus quadricolor var. lutescens Pérez, 1890: 154

  • Psithyrus quadricolor st. lissonurus Dalla Torre, 1896: 570

  • Psithyrus lissonurus race alpium Richards, 1928: 357

  • Psithyrus lissonurus var. alpium Richards, 1928: 358

  • Psithyrus lissonurus var. atricolor Richards, 1928: 356

  • Psithyrus lissonurus var. lutescens Richards, 1928: 357

  • Psithyrus flavidus alpium Popov 1931: 198

  • Psithyrus flavidus frisoni Popov, 1931: 199

  • Psithyrus flavidus var. leucochromus Popov, 1931: 166

  • Psithyrus flavidus var. leucochrous Popov, 1931: 198

  • Psithyrus flavidus var. lissonurus Popov, 1931: 199

  • Psithyrus flavidus var. maculinotus Popov, 1931: 199

  • Psithyrus flavidus alpium f. frey-gessneri Pittioni, 1939: 112

  • Psithyrus flavidus flavidus f. thomsoni Pittioni, 1939: 111

  • Psithyrus flavidus alpium f. analirufescens Pittioni, 1942: 206

  • Psithyrus flavidus alpium f. paradoxus Pittioni, 1942: 207

  • Psithyrus flavidus alpium f. quasiquadricolor Pittioni, 1942: 207

  • Psithyrus flavidus alpium f. rufiolutescens Pittioni, 1942: 207

  • Psithyrus flavidus alpium f. thomsoniformis Pittioni, 1942: 207

  • Psithyrus flavidus flavidus f. intermedius Pittioni, 1942: 206

  • Psithyrus flavidus flavidus f. latofasciatus Pittioni, 1942: 206

  • Psithyrus flavidus flavidus f. superlissonurus Pittioni, 1942: 207

  • Two subspecies of B. flavidus are currently recognized, B. flavidus lutescens (Pérez, 1890) restricted in the Pyrenees and B. flavidus alpium (Richards 1928) restricted in the Alps. Our study did not include Pyrenees and Alps specimens to assess the taxonomic relationships of the latter two taxa.

  • Lectotype (ZISP). Russia (Irkutsk), lectotype ♂ by designation of Pesenko, 1988 (by label only). This lectotype has not been published and is accepted here.

  • Paralectotype (ZISP). Russia (Irkutsk), paralectotype ♂ by designation of Pesenko, 1988 (by label only). This paralectotype has not been published and is accepted here.

  • Photographs of all type specimens of B. fernaldae and B. flavidus and associated labels can be downloaded online as  Supp File 5 (online only).

  • Fig. 8.

    Global distribution of color patterns of Bombus flavidus.The extremes and intermediates of color observed are shown for each region. Color patterns with asterisks (Scandinavia) were extracted from Løken (1985).


    Fig. 9.

    Geographic distribution of B. flavidus color patterns in Oregon. (A) Average percent of individual yellow pile; (B) boxplots comparing total body pile percent yellow between individuals from Western and Eastern localities. The asterisks indicate a significant different between boxplots (unpaired two-tailed t-test, ***P < 0.0001).



    The data presented here are largely congruent in supporting two lineages within a single broadly distributed bumble bee species: an Eastern Nearctic lineage represented by specimens collected from the Eastern Nearctic boreal and Appalachian region, and a Palearctic + Western Nearctic lineage covering the large expanse of the Palearctic and the boreal and mountainous regions of the Western Nearctic. These patterns are supported by CLGS and COI data and to some degree by our morphometric analysis. However, CLGS and COI species delimitation analyses, as well as the inability to clearly distinguish taxa with wing morphometrics and nuclear data examined thus far, suggest that these lineages are not distinct enough to represent separate species. We thus regard these species as conspecific and parts of a single widespread species named Bombus flavidus Eversmann (the oldest available name for the species). We provide a new subspecific name for the Eastern Nearctic population, mostly distributed across the Appalachian Mountains and differing from all other studied subspecies in the studied traits, Bombus flavidus appalachiensis Lhomme and Hines ssp. nov. The choice to declare a subspecies is motivated by the distinction of the lineage. It also distinguishes this lineage for conservation and will facilitate future research on host use and speciation.

    These data point toward a west–east separation of historic Nearctic B. flavidus populations. Several molecular investigations of North American boreal species, from birds (Drovetski 2003, Weir and Schluter 2004) and fishes (Bernatchez and Wilson 1998, Taylor and McPhail 1999) to mammals (Arbogast and Kenagy 2001, Stone et al. 2002), have such west–east genetic structure due to Pleistocene glaciation. This reproductive isolation separating boreal eastern and western zones is thought to have occurred 1.8–0.8 Myr ago during the Early to Mid-Pleistocene (Weir and Schluter 2004), a timing concordant with inferred COI divergence between B. f. flavidus and B. f. appalachiensis (0.85–1 mya inferred using a roughly 1.5–2% divergence per million years; Brower 1994, Quek et al. 2004). Both populations probably reconnected as they expanded their ranges after glaciation, although still keeping some degree of differentiation. The area where both subspecies seem to co-occur, surrounding the Great Lakes, is also the zone where eastern and western boreal regions finally reconnected after glaciation (Barendregt and Duk-Rodkin 2011).

    The COI data also suggest potential repeated movements between the Nearctic and Palearctic, as Western Nearctic specimens show the closest affinity to the Palearctic outgroups (all sibling species are Old World), with the eastern Nearctic and Palearctic lineages derived from the western Nearctic population. The inferred patterns most parsimoniously suggest that B. flavidus originally separated from Palearctic or Oriental sister lineages (B. norvegicus/B. skorikovi) through Palearctic/Nearctic vicariance and that the current Palearctic populations of B. flavidus could be derived from later migration of Nearctic populations back to the Palearctic. This would explain why this single species has such a broad Holarctic range without speciation. Such an event would most likely have occurred during a recent ice age considering relative divergence compared to B. f. appalachiensis.

    Both taxa could not be separated based on wing shape, suggesting either that wing morphology is not a useful tool for differentiating lineages at the species level in this particular group or that these lineages are conspecific. A strong differentiation in wing morphology is not generally expected at the subspecies level due to stabilizing selection on wing shape (Dockx 2007). Nevertheless, we found that B. f. appalachiensis tend to have smaller wings. Most likely differences in wing size is a result of smaller overall body size in these bees, as we saw similar size patterns in genitalia measurements on these bees (Williams 2018, data not presented). More data would be needed to determine whether B. f. appalachiensis tends to be smaller bodied or whether this could be a product of the populations sampled. Body size in bumble bees has been correlated with temperature (Ramírez-Delgado et al. 2016; Scriven et al. 2016) and habitat conditions (Gérard et al. 2018), and Psithyrus species are prone to considerable size plasticity (Plath 1922, Medler 1962). Bombus flavidus appalachiensis is more rare than western Nearctic B. f. flavidus and, given near complete faunal turnover between western and eastern bumble bee species, likely has different hosts, both factors that could affect body size (Plowright and Jay 1977).

    Regarding color pattern variability, B. flavidus specimens display an unusually high degree of sustained local polymorphism in relative amounts of yellow and black pile in their body dorsum across their range. This coloration did not sort distinctly by population, lending further support to the unreliability of using color alone for species delimitation in bumble bees (Carolan et al. 2012; Williams et al. 2019, 2020; Ghisbain et al. 2020a). The broader sylvestris species-group of Psithyrus, in which B. flavidus is contained, is especially prone to high variability in color patterns and other morphological features, leading to historical confusion in species delimitation (Løken and Framstad 1983, Løken 1984). We did observe trends at a local level in Oregon. Oregon includes a transition zone between Müllerian mimicry complexes in bumble bees (Ezray et al. 2019, Tian et al. 2019), where typically paler (usually more red) Rocky Mountain bumble bees occur in eastern parts of the state, and darker forms occur in a mimicry complex along the coast. Bombus flavidus may be under selection to favor darker forms near the coast where its color pattern matches the local mimetic color form. Although the eastern form does not exhibit the red abdominal color typically found in the Rocky Mountain species, this color transition matches well to the same color transition observed for B. fervidus/californicus species complex (Koch et al. 2018). Therefore, the B. fervidus (Fabricius) complex and B. flavidus may converge onto Rocky Mountain patterns by exhibiting paler color forms. A geographic study in Scandinavia found similar localized frequency shifts with a higher incidence of the black form on the western coastal region and more yellow form inland and to the east in B. flavidus (Løken 1984), perhaps suggestive of an adaptive role of climate in coloration shifts as well (Williams 2007). Several species of bumble bees tend to exhibit darker forms in northern oceanic areas with cold and wet springs (Løken 1973, P. Rasmont, personal communication).

    The conspecificity of Nearctic and Palearctic lineages is quite fascinating given that the distribution of this species now spans across the Western United States, Canada, and the Old World, giving B. flavidus, a social parasite, a distribution broader than is achieved by any other bumble bee species. Only two other bumble bee species come close: B. jonellus, ranging from Iceland to Western Canada, and B. cryptarum Fabricius, ranging from Ireland to Western Canada. Although their higher trophic level and potential host specialization might lead one to infer that social parasites should be more locally distributed than host species, B. flavidus may have achieved such an exceptional range because of its socially parasite lifestyle. Cuckoo bumble bees may be less constrained by ecological factors, as they leave nest initiation and most foraging to their hosts and have a shorter seasonality than host species. Their distribution is likely most constrained by the distribution of their host (Suhonen et al. 2016). Although the hosts of B. flavidus have not been clearly identified, it is considered most likely to be a generalist species, which may have promoted its broad range (Suhonen et al. 2016). Closely related species B. sylvestris and B. norvegicus are both host specialists on Pyrobombus species, and because there is a general trend for species belonging to the same Psithyrus subgroup to have closely related hosts (Lhomme and Hines 2019), it is expected that hosts of B. flavidus should be among the Pyrobombus. Pyrobombus is most abundant in the Nearctic, but also has Palearctic relatives, which would have promoted transcontinental movement. Only one nest parasitized by B. flavidus has been described so far, belonging to B. jonellus (Brinck and Wingstrand 1951); however, although the geographical range of both species match well in the Palearctic zone (Pekkarinen et al. 1981, Løken 1984), B. jonellus is limited to the far Northwest Nearctic (Koch et al. 2014). In the Northern Palearctic, the range of B. flavidus is very close to that of Bombus (Bombus) sporadicus Nylander, Bombus (Pyrobombus) cingulatus Wahlberg, B. (Pyrobombus) lapponicus (Fabricius), B. (Pyrobombus) monticola Smith, and B. (Alpinobombus) balteatus Dahlbom, leading them to be proposed as potential host species (Sparre-Schneider 1906, Friese 1923, Pittioni 1942, Løken 1973, Pekkarinen et al. 1981, Rasmont 1988). These species are either absent in the Nearctic or generally only observed at very high altitude or latitude, outside the range of B. flavidus. Bombus flavidus could have broadened its host range in the New World, where Pyrobombus hosts were more abundant and diverse, which could also explain why B. flavidus achieves more abundance in the Nearctic (Lhomme and Hines 2019). Further research on host specificity differences among the newly defined subspecies, examining how host species shift across these regions of endemicity, are needed.

    This study has clarified the contentious species status of this socially parasitic bumble bee. Although B. flavidus flavidus is abundant in certain parts of its range, B. f. appalachiensis is not abundant and is among those taxa that have been documented in the eastern Nearctic to be in decline (Colla and Packer 2008, Colla et al. 2012). It may thus be a population of concern. This study also revealed a cryptic pattern in diversity: that the most distinct divisions in this clade are not between Nearctic and Palearctic, but rather distinguish a population in the eastern Nearctic. Understanding speciation and historical biogeography in this lineage would be improved by examining admixture in the Great Lakes region at a genomic scale.

    Impact on Nomenclature

    Taken together, these data support the conspecificity between both taxa, however with a distinct lineage confined to the eastern Nearctic, that can be considered as a separate subspecies from B. f. flavidus. As B. fernaldae was originally described from both western and eastern Nearctic regions, and the labeled lectotype is from Alaska, thus clearly within the range of B. f. flavidus, this name is best not applied to the new eastern subspecies. Following these results, B. flavidus Eversmann, 1852 is synonymized with B. fernaldae (Franklin, 1911) syn. nov. and a subspecific status, B. f. appalachiensis ssp. nov., is assigned to the distinct lineage ranging from the Appalachians to the eastern boreal regions of the United States and far southeastern Canada.

    A holotype specimen of Bombus flavidus appalachiensis Lhomme and Hines (male; Labels: “USA: Maine, Steuban, 11-vii-2016, 44.44116, -69.9047”, “flavME116, DNA voucher, Hines, Flavidus Project”, “116”, PSUC_FEM 000071169) is deposited at the Frost Entomological Museum at Pennsylvania State University (University Park, PA) (PSUC). The genitalia have been detached and imaged. Aside from its location, this subspecies can be distinguished from other Bombus flavidus by its COI barcode sequence (MW711776 for the holotype; MW711754 reflects the sequence of the paratypes). A photograph of the type specimen, the genitalia of the specimen, and its respective labels are provided in  Supp File 5 [online only]. In addition, five paratype specimens confirmed with COI barcodes as B. flavidus appalachiensis, are also deposited, including: (1) Female, USA: North Carolina, Swain Co., DNA voucher NCflav15115 [photos in Supp File 5], (2) 2 Males from USA: North Carolina, Smoky Mts., DNA vouchers NCflav534 and NCflav539, (3) Male, USA: Michigan, Tuttle Marsh, DNA voucher JSfern1, and (4) Male, USA: Michigan, Stickley, DNA voucher JSfern3; the first three are deposited at PSUC, the later two are housed at the University of Manitoba Wallis Roughley Museum of Entomology (JBWM). Additional information on the paratypes are provided in  Supp Table 1 [online only].

    Supplementary Data

    Supplementary data are available at Insect Systematics and Diversity online.


    We thank Eric Rayfield, Jason Gibbs, Kalyn Bickerman-Martens, Pierre Rasmont, Sydney Cameron, Maxim Proshchalykin, and Arkady Lelej for providing specimens. Thanks to Li Tian for help with imaging, Briana Ezray for advice, and Irena Valterová for assistance with chemical data. G.G. is supported by a grant ‘Aspirant’ from the Fonds National de la Recherche Scientifique (F.R.S.-FNRS, Brussels, Belgium). B.M. is supported by a F.R.S.-FNRS grant ‘Chargé de Recherches’. M.G. is supported by the Swedish Research Council grant number 2018-06238. Funds for this research were provided by the National Science Foundation grant NSF CAREER DEB-1453473 to H.M.H., and a Schreyer Honors College grant to S.W. The authors declare that they have no conflicts of interest. Original data files are available at Dryad data repository, doi: 10.5061/dryad.cvdncjt3t.

    Author Contributions

    H.M.H.: Initiated, Conceptualized, and Funded the project; P.L., S.W., G.G., and H.M.H.: Generated and processed molecular data; G.G., B.M., and H.M.H.: Conducted phylogenetic analyses; P.L. and S.W.: Generated morphometric data; P.L. and M.G.: Conducted morphometric analyses; P.L., B.M., and H.M.H.: Generated chemical data; B.M.: Conducted chemical analyses; P.L., S.W., and H.M.H.: Generated color pattern data; S.W. and H.M.H.: Conducted color pattern analyses; All authors interpreted the data; P.L. and S.W.: Led the writing of the manuscript, with contributions by G.G., B.M., and H.M.H.; All authors contributed to figures, critically edited drafts, and approved the final manuscript.

    References Cited


    Arbogast, B. S., and G. J. Kenagy . 2001. Comparative phylogeography as an integrative approach in biogeography. J. Biogeogr. 28: 819–825. Google Scholar


    Ascher, J. S., and J. Pickering . 2020. Discover life bee species guide and world checklist (Hymenoptera: Apoidea: Anthophila). ( (accessed November 2020). Google Scholar


    Ashmead, W. H. 1902. Papers from the Harriman Alaska Expedition XXVIII. Hymenoptera. Proc. Wash. Acad. Sci. 4: 117–274. Google Scholar


    Ayasse, M., and S. Jarau . 2014. Chemical ecology of bumble bees. Annu. Rev. Entomol. 59: 299–319. Google Scholar


    Aytekin, A. M., M. Terzo, P. Rasmont, and N. Cagatay . 2007. Landmark based geometric analysis of wing shape in Sibiricobombus Vogt (Hymenoptera: Apidae: Bombus Latreille). Ann. Soc. Entomol. Fr. 43: 95–102. Google Scholar


    Barendregt, R. W., and A. Duk-Rodkin . 2011. Chronology and extent of late Cenozoic ice sheets in North America: a magnetostratigraphic assessment, pp. 419–426. In J. Ehlers, P. L. Gibbard, and P.D. Hughes (eds.), Quaternary glaciations: extent and chronology; a closer look. Elsevier, Amsterdam, The Netherlands. Google Scholar


    Bequaert, J., and O. E. Plath . 1925. Description of a new Psithyrus, with an account of Psithyrus laboriosus, and notes on bumblebees. Bull. Mus. Comp. Zool. Harv. 67: 265–288. Google Scholar


    Bernatchez, L., and C. C. Wilson . 1998. Comparative phylogeography of nearctic and Palearctic fishes. Mol. Ecol. 7: 431–452. Google Scholar


    Bertsch, A. 2009. Barcoding cryptic bumblebee taxa: B. lucorum, B. crytarum and B. magnus, a case study. Beitr. Entomol. 59: 287–310. Google Scholar


    Brasero, N., B. Martinet, T. Lecocq, P. Lhomme, P. Biella, I. Valterová, K. Urbanová, M. Cornalba, H. M. Hines, and P. Rasmont . 2018. The cephalic labial gland secretions of two socially parasitic bumblebees Bombus hyperboreus (Alpinobombus) and Bombus inexspectatus (Thoracobombus) question their inquiline strategy. Insect Sci. 25: 75–86. Google Scholar


    Brasero, N., B. Martinet, D. Michez, T. Lecocq, I. Valterová, and P. Rasmont . 2020. Taxonomic revision of the sylvarum-group of bumblebees using an integrative approach. Syst. Biodivers. 18: 12–28. Google Scholar


    Brinck, P., and K. G. Wingstrand . 1951. The mountain fauna of the Virihaure area in Swedish Lapland. Act. Uni. Lund. 45: 1–69. Google Scholar


    Brower, A. V. Z. 1994. Rapid morphological radiation and convergence among races of the butterfly Heliconius erato inferred from patterns of mitochondrial DNA evolution. Proc. Natl. Acad. Sci. USA 91: 6491–6495. Google Scholar


    Cameron, S. A., and B. M. Sadd . 2020. Global trends in bumble bee health. Annu. Rev. Entomol. 65: 209–232. Google Scholar


    Cameron, S. A., H. M. Hines, and P. H. Williams . 2007. A comprehensive phylogeny of the bumble bees (Bombus). Biol. J. Linn. Soc. 91: 161–188. Google Scholar


    Cameron, S. A., J. D. Lozier, J. P. Strange, J. B. Koch, N. Cordes, L. F. Solter, and T. L. Griswold . 2011. Patterns of widespread decline in North American bumble bees. Proc. Natl. Acad. Sci. USA 108: 662–667. Google Scholar


    Carolan, J. C., T. E. Murray, Ú. Fitzpatrick, J. Crossley, H. Schmidt, B. Cederberg, L. McNally, R. J. Paxton, P. H. Williams, and M. J. Brown . 2012. Colour patterns do not diagnose species: quantitative evaluation of a DNA barcoded cryptic bumblebee complex. PLoS One 7: e29251. Google Scholar


    Clement, M., Q. Snell, P. Walker, D. Posada, and K. A. Crandall . 2002. TCS: estimating gene genealogies. In First IEEE International Workshop on High Performance Computational Biology (HiCOMB), 15 April 2002, Fort Lauderdale, FL. Google Scholar


    Cole, J. S., R. B. Siegel, H. L. Loffland, E. A. Elsey, M. B. Tingley, and M. Johnson . 2020. Plant selection by Bumble bees (Hymenoptera: Apidae) in montane riparian habitat of California. Environ. Entomol. 49: 220–229. Google Scholar


    Colla, S. R., and L. Packer . 2008. Evidence for decline in eastern North American bumblebees (Hymenoptera: Apidae), with special focus on Bombus affinis Cresson. Biodiv. Conserv. 17: 1379–1391. Google Scholar


    Colla, S. R., F. Gadallah, L. Richardson, D. Wagner, and L. Gall . 2012. Assessing declines of North American bumble bees (Bombus spp.) using museum specimens. Biodiv. Conserv. 21: 3585–3595. Google Scholar


    Cong, Q., J. Shen, D. Borek, R. K. Robbins, P. A. Opler, Z. Otwinowski, and N. V. Grishin . 2017. When COI barcodes deceive: complete genomes reveal introgression in hairstreaks. Proc. R. Soc. B 284: 20161735. Google Scholar


    Cresson, E. T. 1863. List of the North American species of Bombus and Apathus. Proc. Entomol. Soc. Phila. 2: 83–116. Google Scholar


    Dalla Torre, C. G. 1896. Catalogus Hymenopterorum, hucusque descriptorum, systematicus et synonymicus. Apidae (Anthophila). Lipsiae, Sumptibus Guilelmi Engelmann. Google Scholar


    Dayrat, B. 2005. Towards integrative taxonomy. Biol. J. Linn. Soc. 85: 407–415. Google Scholar


    De Meulemeester, T., P. Gerbaux, M. Boulvin, A. Coppée, and P. Rasmont . 2011. A simplified protocol for bumble bee species identification by cephalic secretion analysis. Insectes Soc. 58: 227–236. Google Scholar


    de Queiroz, K. 2007. Species concepts and species delimitation. Syst. Biol. 56: 879–886. Google Scholar


    Dehon, M., M. S. Engel, M. Gérard, A. M. Aytekin, G. Ghisbain, P. H. Williams, P. Rasmont, and D. Michez . 2019. Morphometric analysis of fossil bumble bees (Hymenoptera, Apidae, Bombini) reveals their taxonomic affinities. Zookeys 891: 71–118. Google Scholar


    Dellicour, S., and T. Lecocq . 2013a. GCALIGNER 1.0 and GCKOVATS 1.0: manual of a software suite to compute a multiple sample comparison data matrix from eco-chemical datasets obtained by gas chromatography. University of Mons, Mons, Belgium. Google Scholar


    Dellicour, S., and T. Lecocq . 2013b. GCALIGNER 1.0: an alignment program to compute a multiple sample comparison data matrix from large eco-chemical datasets obtained by GC. J. Sep. Sci. 36: 3206–3209. Google Scholar


    Dewulf, A., T. De Meulemeester, M. Dehon, M. S. Engel and D. Michez . 2014. A new interpretation of the bee fossil Melitta willardi Cockerell (Hymenoptera, Melittidae) based on geometric morphometrics of the wing. ZooKeys 389: 35–48. Google Scholar


    Dockx, C. 2007. Directional and stabilizing selection on wing size and shape in migrant and resident monarch butterflies, Danaus plexippus (L.), inCuba. Biol. J. Linn. Soc. 92: 605–616. Google Scholar


    Drovetski, S. V. 2003. Plio-Pleistocene climatic oscillations, Holarctic biogeography and speciation in an avian subfamily. J. Biogeogr. 30: 1173–1181. Google Scholar


    Drummond, A. J., and A. Rambaut . 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 7: 214. Google Scholar


    Duennes, M. A., J. D. Lozier, H. M. Hines, and S. A. Cameron . 2012. Geographical patterns of genetic divergence in the widespread Mesoamerican bumble bee Bombus ephippiatus (Hymenoptera: Apidae). Mol. Phylogenet. Evol. 64: 219–231. Google Scholar


    Duennes, M. A., C. Petranek, E. P. D. de Bonilla, J. Mérida-Rivas, O. Martinez-López, P. Sagot, R. Vandame, and S. A. Cameron . 2017. Population genetics and geometric morphometrics of the Bombus ephippiatus species complex with implications for its use as a commercial pollinator. Conserv. Genet. 18: 553–572. Google Scholar


    Eversmann, E. 1852. Fauna Hymenopterologica Volgo-Uralensis (Continuatio). Byull. Moskovsk. Obshch. Isp. Prir. 25: 1–137. Google Scholar


    Ezray, B. D., D. C. Wham, C. Hill, and H. M. Hines . 2019. Unsupervised machine learning reveals mimicry complexes in bumble bees occur along a perceptual continuum. Proc. R. Soc. B 286: 20191501. Google Scholar


    Fisher, R. M., and B. J. Sampson . 1992. Morphological specializations of the bumblebee social parasite Psithyrus ashtoni (Cresson) (Hymenoptera: Apidae). Can. Entomol. 124: 69–77. Google Scholar


    Folmer, O., M. Black, W. Hoeh, R. Lutz, and R. Vrijenhoek . 1994. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol. Mar. Biol. Biotechnol. 3: 294–299. Google Scholar


    Franklin, H. J. 1911. New North American Bombidae. Trans. Am. Ent. Soc. 37: 157–168. Google Scholar


    Friese, H. 1923. Hymenoptera. Apidae, pp. 3–9. In A. W. Brøggers boktrykkeria/s (eds.), Norwegian expedition to Novaya Zemlya 1921. Report of the scientific results, No. 14. Videnskapsselskapet i Kristiania, Oslo, Norway. Google Scholar


    Frison, T. H. 1926. Descriptions and records of North American Bremidae, together with notes on the synonymy of certain species (Hymenoptera). Trans. Am. Ent. Soc. 52: 129–145. Google Scholar


    Gérard, M., D. Michez, D. Fournier, K. Maebe, G. Smagghe, J. C. Biesmeijer, T. De Meulemeester . 2015. Discrimination of haploid and diploid males of Bombus terrestris (Hymenoptera; Apidae) based on wing shape. Apidologie 46: 644–653. Google Scholar


    Gérard, M., D. Michez, V. Debat, L. Fullgrabe, I. Meeus, N. Piot, O. Sculfort, M. Vastrade, G. Smagghe, and M. Vanderplanck . 2018. Stressful conditions reveal decrease in size, modification of shape but relatively stable asymmetry in bumblebee wings. Sci. Rep. 8: 15169. Google Scholar


    Gérard, M., B. Martinet, K. Maebe, G. Smagghe, L. Marshall, N. J. Vereecken, S. Vray, P. Rasmont, and D. Michez . 2020. Shift in size of bumblebee queens over the last century. Glob. Chang. Biol. 26: 1185–1195. Google Scholar


    Ghisbain, G., D. J. Lozier, S. R. Rahman, B. D. Ezray, L. Tian, J. Strange, J. M. Ulmer, S. D. Heraghty, J. P. Strange, P. Rasmont , et al. 2020a. Substantial genetic divergence and lack of recent gene flow support cryptic speciation in a color polymorphic bumble bee (Bombus bifarius) species complex. Syst. Entomol. 45: 635–652. Google Scholar


    Ghisbain, G., D. Michez, L. Marshall, P. Rasmont, and S. Dellicour . 2020b. Wildlife conservation strategies should incorporate both taxon identity and geographical context – further evidence with bumblebees. Divers. Distrib. 26: 1741–1751. Google Scholar


    Goulson, D., E. Nicholls, C. Botías, and E. L. Rotheray . 2015. Bee declines driven by combined stress from parasites, pesticides, and lack of flowers. Science 347: 1255957. Google Scholar


    Gueuning, M., J. E. Frey, and C. Praz . 2020. Ultraconserved yet informative for species delimitation: Ultraconserved elements resolve long-standing systematic enigma in Central European bees. Mol. Ecol. 29: 4203–4220. Google Scholar


    Herbert, P. D. N., E. H. Penton, J. M. Burns, D. H. Janzen, and H. Hallwachs . 2004. Ten species in one: DNA Barcoding reveals cryptic species in the Neotropical skipper butterfly Astraptes fulgerator. Proc. Natl. Acad. Sci. USA 101: 14812–14817. Google Scholar


    Hines, H. M. 2008. Historical biogeography, divergence times, and diversification patterns of bumble bees (Hymenoptera: Apidae: Bombus). Syst. Biol. 57: 58–75. Google Scholar


    Hines, H. M., and P. H. Williams . 2012. Mimetic colour pattern evolution in the highly polymorphic Bombus trifasciatus (Hymenoptera: Apidae) species complex and its comimics. Zool. J. Linn. Soc. 166: 805–826. Google Scholar


    Hines, H. M., S. A. Cameron, and P. H. Williams . 2006. Molecular phylogeny of the bumble bee subgenus Pyrobombus (Hymenoptera: Apidae: Bombus) with insights into gene utility for lower-level analysis. Invertebr. Syst. 20: 289–303. Google Scholar


    Hobbs, G. A. 1965. Ecology of species of Bombus Latr. (Hymenoptera: Apidae) in southern Alberta: III. Subgenus Cullunanobombus Vogt. Can. Entomol. 97: 1293–1302. Google Scholar


    Hobbs, G. A. 1966. Ecology of species of Bombus Lam. (Hymenoptera: Apidae) in southern Alberta: V. Subgenus Subterramobombus Vogt. Can. Entomol. 98: 288–294. Google Scholar


    Hobbs, G. A. 1967. Ecology of species of Bombus (Hymenoptera: Apidae) in southern Alberta: VI. Subgenus Pyrobombus. Can. Entomol. 99:1271–1292. Google Scholar


    Hobbs, G. A. 1968. Ecology of species of Bombus (Hymenoptera: Apidae) in southern Alberta: VI1. Subgenus Bombus. Can. Entomol. 100: 156–164. Google Scholar


    Hothorn, T., F. Bretz, and P. Westfall . 2008. Simultaneous inference in general parametric models. Biom. Z. 50: 346–363. Google Scholar


    Huberty, C. J., and S. Olejnik . 2006. Applied MANOVA and discriminant analysis, 2nd edn. John Wiley and Sons Inc., Hoboken, NJ. Google Scholar


    Jacobson, M. M., E. M. Tucker, M. E. Mathiasson, and S. M. Rehan . 2018. Decline of bumble bees in northeastern North America, with special focus on Bombus terricola. Biol. Conserv. 217: 437–445. Google Scholar


    Koch, J., J. P. Strange, and P. Williams . 2014. Bumble bees of the Western United States. United States Forest Service, Washington, DC. Google Scholar


    Koch, J. B., J. Rodriguez, J. P. Pitts, and J. P. Strange . 2018. Phylogeny and population genetic analyses reveals cryptic speciation in the Bombus fervidus species complex (Hymenoptera: Apidae). PLoS One 13: e0207080. Google Scholar


    Kreuter, K. R., E. Bunk, A. Lückemeyer, R. Twele, W. Francke, and M. Ayasse . 2012. How the social parasitic bumblebee Bombus bohemicus sneaks into power of reproduction. Behav. Ecol. Sociobiol. 66: 475–486. Google Scholar


    Kumar, S., G. Stecher, M. Li, C. Knyaz, and K. Tamura . 2018. MEGAX: molecular evolutionary genetics analysis across computing platforms. Mol. Biol. Evol. 35: 1547–1549. Google Scholar


    Laverty, T. M., and L. D. Harder . 1988. The bumble bees of eastern Canada. Can. Entomol. 120: 965–987. Google Scholar


    Lecocq, T., P. Lhomme, D. Michez, S. Dellicour, I. Valterová, and P. Rasmont . 2011. Molecular and chemical characters to evaluate species status of two cuckoo bumblebees: Bombus barbutellus and Bombus maxillosus (Hymenoptera, Apidae, Bombini). Syst. Entomol. 36: 453–469. Google Scholar


    Lecocq, T., S. Dellicour, D. Michez, M. Dehon, A. Dewulf, T. De Meulemeester, N. Brasero, I. Valterová, J. Y. Rasplus, and P. Rasmont . 2015a. Methods for species delimitation in bumblebees (Hymenoptera, Apidae, Bombus): towards an integrative approach. Zool. Scr. 44: 281–297. Google Scholar


    Lecocq, T., A. Coppée, T. Mathy, P. Lhomme, M. C. Cammaerts-Tricot, K. Urbanová, I. Valterová, and P. Rasmont . 2015b. Subspecific differentiation in male reproductive traits and virgin queen preferences, in Bombus terrestris. Apidologie 46: 595–605. Google Scholar


    Lecocq, T., N. Brasero, T. De Meulemeester, D. Michez, S. Dellicour, P. Lhomme, R. De Jonghe, I. Valterová, K. Urbanová, and P. Rasmont . 2015c. An integrative taxonomic approach to assess the status of Corsican bumblebees: implications for conservation. Anim. Conserv. 18: 236–248. Google Scholar


    Leigh, J., and D. Bryant . 2015. POPART: full-feature software for haplotype network construction. Methods Ecol. Evol. 6: 1110–1116. Google Scholar


    Lhomme, P., and H. M. Hines . 2018. Reproductive dominance strategies in insect social parasites. J. Chem. Ecol. 44: 838–850. Google Scholar


    Lhomme, P., and H. M. Hines . 2019. Ecology and evolution of cuckoo bumble bees. Ann. Entomol. Soc. Am. 112: 122–140. Google Scholar


    Lhomme, P., M. Ayasse, I. Valterová, T. Lecocq, and P. Rasmont . 2012. Born in an alien nest: how do social parasite male offspring escape from host aggression? PLoS One 7: e43053. Google Scholar


    Lhomme, P., A. Sramkova, K. Kreuter, T. Lecocq, P. Rasmont, and M. Ayasse . 2013. A method for year–round rearing of cuckoo bumblebees (Psithyrus). Ann. Soc. Entomol. Fr. 49: 117–125. Google Scholar


    Lhomme, P., M. Ayasse, I. Valterová, T. Lecocq, and P. Rasmont . 2015. A scent shield to survive: identification of the repellent compounds secreted by the male offspring of the cuckoo bumblebee Bombus vestalis. Entomol. Exp. App. 157: 263–270. Google Scholar


    Løken, A. 1973. Studies on Scandinavian bumblebees (Hymenoptera, Apidae). Norsk Entomol. Tidsk. 20: 1–218. Google Scholar


    Løken, A. 1984. Scandinavian species of the genus Psithyrus (Hymenoptera: Apidae). Entomol. Scand. 23: 1–45. Google Scholar


    Løken, A. 1985. Norske insekttabeller. 9. Humler. Norsk Entomologisk Forening, Oslo, Norway. Google Scholar


    Løken, A., and E. B. Framstad . 1983. Contribution to the taxonomy of Psithyrus (Fernaldaespsithyrus) (Hymenoptera: Apidae). Acta Entomol. Fenn. 42: 46–50. Google Scholar


    Martinet, B., T. Lecocq, N. Brasero, P. Biella, K. Urbanová, I. Valterová, M. Cornalba, J. O. Gjershaug, D. Michez, and P. Rasmont . 2018. Following the cold: geographic differentiation between interglacial refugia and speciation in Arcto-Alpine species complex Bombus monticola (Hymenoptera: Apidae). Syst. Entomol. 43: 200–217. Google Scholar


    Martinet, B., T. Lecocq, N. Brasero, M. Gérard, K. Urbanová, I. Valterová, J. O. Gjershaug, D. Michez, and P. Rasmont . 2019. Integrative taxonomy of an arctic bumblebee species complex highlights a new cryptic species (Apidae: Bombus). Zool. J. Linn. Soc. 187: 599–621. Google Scholar


    Medler, J. T. 1962. Morphometric studies on bumble bees. Ann. Entomol. Soc. Am. 55: 212–218. Google Scholar


    Oksanen, J., F. G. Blanchet, R. Kindt, P. Legendre, P. R. Minchin, R. B. O'Hara, G. L. Simpson, P. Solymos, M. H. H. Stevens, and H. Wagner . 2020. Vegan: community ecology package. R Package Version 2.5–6. ( ). Google Scholar


    Owen, R. E. 2012. Applications of morphometrics to the Hymenoptera, particularly bumblebees (Bombus, Apidae). In C. Wahl (ed.), Morphometrics. IntechOpen, London. Google Scholar


    Paradis, E., J. Claude, and K. Strimmer . 2004. APE: analyses of phylogenetics and evolution in R language. Bioinformatics 20: 289–290. Google Scholar


    Pekkarinen, A., and I. Teräs . 1993. Zoogeography of Bombus and Psithyrus in northwestern Europe (Hymenoptera, Apidae). Ann. Zool. Fenn. 30: 177–208. Google Scholar


    Pekkarinen, A., I. Teräs, J. Viramo, and J. Paatela . 1981. Distribution of bumblebees (Hymenoptera, Apidae: Bombus and Psithyrus) in eastern Fennoscandia. Not. Entomol. 61: 71–89. Google Scholar


    Pérez, J. 1884. Contribution à la faune des Apiaires de France, 2ème partie: Parasites. Actes Soc. Linn. Bordeaux 37: 205–380. Google Scholar


    Pérez, J. 1890. Catalogue des Mellifères du Sud-Ouest. Actes Soc. Linn. Bordeaux 44: 133–200. Google Scholar


    Pittioni, B. 1939. Die hummeln und Schmarotzerhummeln der Balkan-Halbinsel. Mitt. Königl. Naturwissensch. Inst. Sofia 12: 49–115. Google Scholar


    Pittioni, B. 1942. Die boreoalpinen Hummeln und Schmarotzerhummeln (Hymen., Apidae, Bombinae). Mitt. Königl. Naturwissensch. Inst. Sofia 15: 155–218. Google Scholar


    Plath, O. E. 1922. Notes on Psithyrus with records of two new American hosts. Biol. Bull. Mar. Biol. Lab. Woods Hole. 43: 23–44. Google Scholar


    Plowright, R. C., and S. C. Jay . 1977. On the size determination of bumble bee castes (Hymenoptera: Apidae). Can. J. Zool. 55: 1133–1138. Google Scholar


    Popov, V. B. 1931. Zur Kenntnis der paläarktischen Schmarotzerhummeln (Psithyrus Lep.). Eos. Madr. 7: 131–209. Google Scholar


    Potapov, G. S., A. V. Kondakov, V. M. Spitsyn, B. Y. Filippov, Y. S. Kolosova, N. A. Zubrii, and I. N. Bolotov . 2018. An integrative taxonomic approach confirms the valid status of Bombus glacialis, an endemic bumblebee species of the High Arctic. Polar. Biol. 41: 629–642. Google Scholar


    Proshchalykin, M. Y. 2007. Families Colletidae, Andrenidae, Melittidae, Megachilidae, Apidae, pp. 897–908. In A. S. Leley (ed.), Key to the insects of the Russian Far East 4 (5). Dal'nauka, Vladivostok, Russia. Google Scholar


    Quek, S. P., S. J. Davies, T. Itino, and N. E. Pierce . 2004. Codiversification in an ant-plant mutualism: stem texture and the evolution of host use in Crematogaster (Formicidae: Myrmicinae) inhabitants of Macaranga (Euphorbiaceae). Evolution 58: 554–570. Google Scholar


    R Development Core Team. 2020. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. ( ). Google Scholar


    Rambaut, A., A. J. Drummond, D. Xie, G. Baele, and M. A. Suchard . 2018. Posterior summarisation in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67: 901–904. Google Scholar


    Ramírez-Delgado, V. H., S. Sanabria-Urbán, M. A. Serrano-Meneses, and R. Cueva Del Castillo . 2016. The converse to Bergmann's rule in bumblebees, a phylogenetic approach. Ecol. Evol. 6: 6160–6169. Google Scholar


    Rasmont, P. 1988. Monographie écologique et zoogéographique des Bourdons de France et de Belgique (Hymenoptera, Apidae, Bombinae). Ph.D. thesis dissertation, Faculté des Sciences agronomiques de l'Etat, Gembloux, Belgium. Google Scholar


    Rasmont, P., M. Franzén, T. Lecocq, A. Harpke, S. P. M. Roberts, K. Biesmeijer, L. Castro, B. Cederberg, L. Dvoř ák, Ú. Fitzpatrick , et al. 2015. Climatic risk and distribution Atlas of European bumblebees. BioRisk 10: 1–246. Google Scholar


    Reid, N. M., and B. C. Carstens . 2012. Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed Yule-coalescent model. BMC Evol. Biol. 12: 196. Google Scholar


    Reinig, W. F. 1935. On the variation of Bombus lapidarius L. and its cuckoo, Psithyrus rupestris Fabr., with notes on mimetic similarity. J. Genet. 30: 321–356. Google Scholar


    Richards, O. W. 1928. A revision of the European bees allied to Psithyrus quadricolor Lepeletier (Hymenoptera, Bombidae). Trans. Entomol. Soc. Lond. 76: 345–365. Google Scholar


    Richardson, L. L., K. P. McFarland, S. Zahendra, and S. Hardy . 2018. Bumble bee (Bombus) distribution and diversity in Vermont, USA: a century of change. J. Insect Conserv. 23: 45–62. Google Scholar


    Rohlf, F. J. 2013a. tpsUTIL, version 1.56. Department of Ecology and Evolution, State University of New York, Stony Brook, NY. Google Scholar


    Rohlf, F. J. 2013b. tpsDig, version 2.17. Department of Ecology and Evolution. State University of New York, Stony Brook, NY. Google Scholar


    Rohlf, F. J., and D. Slice . 1990. Extensions of the Procrustes method for the optimal superimposition of landmarks. Syst. Zool. 39: 40–59. Google Scholar


    Schlick-Steiner, B. C., F. M. Steiner, B. Seifert, C. Stauffer, E. Christian, and R. H. Crozier . 2010. Integrative taxonomy: a multisource approach to exploring biodiversity. Ann. Rev. Entomol. 55: 421–438. Google Scholar


    Scriven, J. J., P. R. Whitehorn, D. Goulson, and M. C. Tinsley . 2016. Bergmann's body size rule operates in facultatively endothermic insects: evidence from a complex of cryptic bumblebee species. PLoS One 11: e0163307. Google Scholar


    Soltani, G. G., D., Bénon, N. Alvarez, and J. C. Praz . 2017. When different contact zones tell different stories: putative ring species in the Megachile concinna species complex (Hymenoptera: Megachilidae). Biol. J. Linn. Soc. 121: 815–832. Google Scholar


    Sparre-Schneider, J. 1906. Sydherö. Et lidet bidrag til kundskaben on den arktiske skjaergaards malakologiske og entomologiske fauna. Tromsö Mus. Arsh. 27: 170–205. Google Scholar


    Sramkova, A., and M. Ayasse . 2008. Psithyrus females do possess wax glands. Insectes Soc. 55: 404–406. Google Scholar


    Stone, K. D., W. Flynn, and J. A. Cook . 2002. Post-glacial colonization of northwestern North America by the forest associated American marten (Martes americana, Mammalia Carnivora: Mustelidae). Mol. Ecol. 11: 2049–2063. Google Scholar


    Suhonen, J., J. Rannikko, and J. Sorvari . 2015. The rarity of host species affects the co-extinction risk in socially parasitic bumblebee Bombus (Psithyrus) species. Ann. Zool. Fenn. 52: 236–242. Google Scholar


    Suhonen, J., J. Rannikko, and J. Sorvari . 2016. Species richness of cuckoo bumblebees is determined by the geographical range area of the host bumblebee. Insect. Cons. Div. 9: 529–535. Google Scholar


    Suzuki, R., and H. Shimodaira . 2011. Pvclust: hierarchical clustering with P-values via multiscale bootstrap resampling. Contributed package. Version 1–1.10. Vienna, Austria: R Foundation for Statistical Computing. ( ). Google Scholar


    Taylor, E. B., and J. D. McPhail . 1999. Evolutionary history of an adaptive radiation in species pairs of threespine sticklebacks (Gasterosteus aculeatus): insights from mitochondrial DNA. Biol. J. Linn. Soc. 66: 271–291. Google Scholar


    Thomson, C. G. 1872. Hymenoptera Scandinaviae. 2. Berling, Lund, Sweden. Google Scholar


    Tian, L., S. R. Rahman, B. D. Ezray, L. Franzini, J. P. Strange, P. Lhomme, and H. M. Hines . 2019. A homeotic shift late in development drives mimetic color variation in a bumble bee. Proc. Natl. Acad. Sci. USA 116: 11857–11865. Google Scholar


    Valterová, I., B. Martinet, D. Michez, P. Rasmont, and N. Brasero . 2019. Sexual attraction: a review of bumblebee male pheromones. Z. Naturforsch. C. J. Biosci. 74: 233–250. Google Scholar


    Weir, J. T., and D. Schluter . 2004. Ice sheets promote speciation in boreal birds. Proc. R. Soc. B 271: 1881–1887. Google Scholar


    Williams, P. H. 1998. An annotated checklist of bumble bees with an analysis of patterns of description (Hymenoptera: Apidae, Bombini). Bull. Br. Mus. Nat. Hist. Entomol. 67: 79–152. Google Scholar


    Williams, P. H. 2007. The distribution of bumblebee colour patterns worldwide: possible significance for thermoregulation, crypsis, and warning mimicry. Biol. J. Linn. Soc. 92: 97–118. Google Scholar


    Williams, P. H. 2008. Do the parasitic Psithyrus resemble their host bumblebees in colour pattern? Apidologie 39: 637–649. Google Scholar


    Williams, P. H., M. J. F. Brown, J. C. Carolan, J. An, D. Goulson, A. M. Aytekin, L. R. Best, A. M. Byvaltsev, B. Cederberg, R. Dawson , et al. 2012. Unveiling cryptic species of the bumblebee subgenus Bombus s. str. worldwide with COI barcodes (Hymenoptera: Apidae). Syst. Biodivers. 10: 21–56. Google Scholar


    Williams, P. H., R. W. Thorp, L. L. Richardson, and S. R. Colla . 2014. Bumble bees of North America: an identification guide. Princeton University Press, Princeton, NJ. Google Scholar


    Williams, P. H., A. M. Byvaltsev, B. Cederberg, M. V. Berezin, F. Ødegaard, C. Rasmussen, L. L. Richardson, H. Jiaxing, C. S. Sheffield, and S. T. Williams . 2015. Genes suggest ancestral colour polymorphisms are shared across morphologically cryptic species in arctic bumblebees. PLoS One 10: e0144544. Google Scholar


    Williams, P. H., M. V. Berezin, S. G. Cannings, B. Cederberg, F. Ødegaard, C. Rasmussen, L. L. Richardson, J. Rykken, C. S. Sheffield, C. Thanoosing , et al. 2019. The arctic and alpine bumblebees of the subgenus Alpinobombus revised from integrative assessment of species' gene coalescents and morphology (Hymenoptera, Apidae, Bombus). Zootaxa 4625: 1–68. Google Scholar


    Williams, P. H., D. Altanchimeg, A. Byvaltsev, R. De Jonghe, S. Jaffar, G. Japoshvili, S. Kahono, H. Liang, M. Mei, A. Monfared , et al. 2020. Widespread polytypic species or complexes of local species? Revising bumblebees of the subgenus Melanobombus world-wide (Hymenoptera, Apidae, Bombus). Eur. J. Taxon. 719: 1–120. Google Scholar


    Williams, S. D. 2018. Species delimitation and the evolution of highly variable traits in the broad holarctic socially parasitic bumble bee, Bombus flavidus. B. S. thesis dissertation, The Pennsylvania State University, State College, PA. Google Scholar


    Wilson, J. S., L. E. Wilson, L. D. Loftis, and T. Griswold . 2010. The montane bee fauna of north central Washington, USA, with floral associations. West. North Am. Nat. 70: 198–207. Google Scholar


    Zetterstedt, J. W. 1838. Bombus Fabr. columns 473–474. In J. W. Zetterstedt (ed.), Insecta lapponica descripta. Voss, Lipsiae (Leipzig), Germany. Google Scholar
    © The Author(s) 2021. Published by Oxford University Press on behalf of Entomological Society of America.
    Patrick Lhomme, Sarah D. Williams, Guillaume Ghisbain, Baptiste Martinet, Maxence Gérard, and Heather M. Hines "Diversification Pattern of the Widespread Holarctic Cuckoo Bumble Bee, Bombus flavidus (Hymenoptera: Apidae): The East Side Story," Insect Systematics and Diversity 5(2), 1-15, (31 December 2022).
    Received: 14 December 2020; Accepted: 14 February 2021; Published: 31 December 2022
    Integrative taxonomy
    social parasite
    Back to Top