Diversification Pattern of the Widespread Holarctic Cuckoo Bumble Bee, Bombus flavidus (Hymenoptera: Apidae): The East Side Story

Abstract 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, 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) that is unique among bumble bees in consisting exclusively of social parasites (cuckoo bumble bees) of Research 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 . 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. 2012Lhomme et al. , 2013Lhomme et al. , 2015Ayasse 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 subgenera (e.g., Williams et al. 2012Williams et al. , 2019Williams et al. , 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. 2012Williams et al. , 2015Brasero 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. 2012Williams et al. , 2019Ghisbain 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. , 2020Cong 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. 2018Martinet et al. , 2019Brasero et al. 2018Brasero et al. , 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 . It is mainly found in subalpine/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), 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, 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 . 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 andHarder 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 ). 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.
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.

Sampling
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).
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 Biotek 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, http://www.geneious.com), 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 (https://www.phylo.org) 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 Yulecoalescent 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. 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 . CLGS consist of a complex mixture of compounds (mainly aliphatic or isoprenoid), with variable main compounds . 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 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. 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 andOlejnik 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.

Nomenclature
This article and the nomenclatural act(s) it contains have been registered in Zoobank (www.zoobank.org), the official register of the International Commission on Zoological Nomenclature. The LSID (Life Science Identifier) number of the publication is: 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.
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).

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), whereas 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).
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.
Photographs of all type specimens of B. fernaldae and B. flavidus and associated labels can be downloaded online as Supp File 5 (online only).

Discussion
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 andSchluter 2004) and fishes (Bernatchez andWilson 1998, Taylor andMcPhail 1999) to mammals (Arbogast andKenagy 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 Williams et al. 2019Williams et al. , 2020Ghisbain 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 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 , 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 , 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 (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 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)

Supplementary Data
Supplementary data are available at Insect Systematics and Diversity online.

Acknowledgments
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