Phylogeny and population genetic analyses reveals cryptic speciation in the Bombus fervidus species complex (Hymenoptera: Apidae)

Bumble bees (Bombus Latrielle) are significant pollinators of flowering plants due to their large body size, abundant setae, and generalist foraging strategies. However, shared setal coloration patterns among closely and distantly related bumble bee species makes identification notoriously difficult. The advent of molecular genetic techniques has increased our understanding of bumble bee evolution and taxonomy, and enables effective conservation policy and management. Individuals belonging to the North American Bombus fervidus species-complex (SC) are homogenous in body structure but exhibit significant body color phenotype variation across their geographic distribution. Given the uncertainty of the genealogical boundaries within the SC, some authors have synonymized all members of the B. fervidus SC within a single taxon, while others propose an alternative two taxa hypothesis. Operating under the phylogenetic species concept, our analysis supports the hypothesis that there are two independent lineages of bumble bees within the B. fervidus SC. With the current evidence, however, it is not possible to assign valid names to either of them, because both lineages include the color phenotypes found in the original species descriptions of B. fervidus and B. californicus. Cryptic speciation does not seem to be the product of Müllerian mimicry between the clades, because diverging coloration patterns are observed when the distribution of the clades overlaps. Furthermore, within each lineage there is evidence for strong population differentiation that is correlated with geographic distribution rather than color phenotype. In our study, we demonstrate the importance of obtaining a broad sample of multiple populations when conducting lower-level phylogenetic analyses. In addition to improving our knowledge of bumble bee diversification patterns, characterizing the evolutionary history of these pollinators provides the foundation needed to guide contemporary conservation assessments and management strategies.


Introduction
Cryptic speciation is the process in which organisms share a nearly identical phenotype but belong to different species [1]. It is a common phenomenon observed across the understudied and numerically dominant insects, and can pose a significant hurdle to effective conservation and management [2]. Biodiversity is rapidly declining on a global scale primarily due to resource extraction activities associated with economic growth and expansion. In fact, it is estimated that the contemporary extinction rate is 1,000 times higher than what has been experienced prior to the global effects of humanity's economic and developmental activities [3]. A major impediment to the effective conservation of biodiversity includes the lack of consensus among scientists and conservation practitioners on the taxonomic resolution appropriate to a conservation or management goal. Without an operational unit that considers the ecology and evolutionary history of a species, efforts to promote species conservation will remain daunting [2].
Bumble bees (Hymenoptera: Apidae, Bombus) are one of the most important native pollinators of North America, contributing to the ecosystem services required by wild and economically important flowering plant species [4,5]. They are dominant pollinators of the northern hemisphere, specifically in alpine and temperate ecosystems [6,7]. Furthermore, wild bumble bee populations have been found to enhance crop productivity through effective pollination [4,8,9]. However, the global decline of wild bumble bee populations due to disease, pesticides, urbanization, and agricultural intensification have prompted state, national, and international efforts to document the diversity and distribution of these iconic bee fauna [10][11][12].
Concurrent efforts to conserve bumble bees are dependent on recognizing operational units, whether they are species, taxonomic, evolutionary, or otherwise [2,11,13]. These units have been useful in unveiling local biotic and abiotic factors that are specific to unique evolutionary lineages of cryptic species [2,14]. Due to the spatial cohabitation of aposematic setal coloration patterns, bumble bees have proven to be difficult to identify to species by both novice and seasoned taxonomists [15][16][17][18][19][20]. The dependence on setal coloration patterns to delineate between closely related species has caused debate on the species status of many of these taxa [15,16,20]. Contemporary phylogenetic investigations using both single and multiple genetic loci, as well as morphology-based taxonomic studies, have resolved some cryptic species complexes across bumble bee subgenera [18,19,21]. It has been demonstrated with a single gene, Cytochrome c oxidase I (COI), that bumble bees exhibiting nearly identical aposematic coloration patterns have been found to be separate species [22][23][24][25]. However, a lack of COI variation between species has also been detected, leading to synonymizations [6,15].
In this study, we examine the evolutionary history of the Bombus fervidus species-group (SC), which contains two species: B. fervidus (Fabricius, 1798) and the nominal sister taxon B. californicus Smith, 1954. These species belong to the globally distributed subgenus Thoracobombus [21,26]. The decline in flowering plants with long corollas due to urbanization and agricultural intensification has been implicated in the decline of European Thoracobombus [27,28]. Additionally, North American B. fervidus and B. californicus have been found to be declining in abundance in both wild and urban environments, relative to historic population abundance estimates [29][30][31][32]. Increased disease detection in wild populations of B. fervidus and another Thoracobombus, B. pensylvanicus, has been hypothesized to be a major contributor to their decline in the wild [10,33].
Bombus fervidus and B. californicus have been recognized to be legitimate species, based on historic and contemporary investigations using taxonomic and comprehensive phylogenetic tools [16,20,21,26]. However, the lack of strong divergence in COI, and exhibition of transitional color patterns following a continuum of variation from mostly black (i.e., B. californicus) to mostly yellow (i.e., B. fervidus) has been suggested to be evidence that they are conspecific [15,34]. Bombus californicus is distributed from the Pacific Coast of North America, east to the Black Hills of South Dakota [16,17,20]. Unlike B. californicus, which is distributed across a broad latitudinal gradient relative to the longitudinal range, B. fervidus has a transcontinental distribution, from the Pacific Coast to the northeastern United States [15,17,20,35]. While both species are sympatric in portions of their range in western North America, Hobbs [36] suggested that B. fervidus and B. californicus differ in nesting habitats in Canada, with B. californicus nesting in wooded areas and the foothills of southern Alberta, and B. fervidus primarily found to be in the prairies [20].
Multiple taxonomic investigations of the two bumble bee species have agreed on one central idea: they are nearly impossible to separate morphologically [15,16,20,21,38]. In regards to distinguishing between B. californicus consanguineus and B. fervidus, W.P. Stephen stated in Bumble Bees of Western North America, "There are no morphological features in either species by which they can be distinguished, and separation is made exclusively on color pattern" [16]. In regards to distinguishing B. californicus (sensu stricto, phenotype 1) from B. fervidus (sensu stricto, phenotype 4), he went on to write, "The species is close morphologically to B. fervidus (Fabr.) and is impossible to separate structurally from that species." [16]. Finally, W.P. Stephen citing Franklin [38], went on to state that "californicus and fervidus may eventually prove to be subspecies of a single species". Twenty-six years later, R. Thorp led the writing of Bumble Bees and Cuckoo Bumble Bees of California (Hymenoptera: Apidae), and expressed a similar sentiment for the lack of variability (outside of setal color) between B. californicus and B. fervidus [20]. Furthermore, a comparison of the male genitalia between B. fervidus and B. californicus found no morphological differences [20]. However, he stated that there were distinct ecological differences between B. californicus and B. fervidus when sympatric, showing no signs of intergradation. At present, there is no biological evidence that B. californicus and B. fervidus have the capacity to breed in the wild, despite historic reports that initially proposed this hypothesis [34]. In a global systematic survey of bumble bees, Cameron et al. [26] inferred a phylogeny based on five genetic loci and found that B. fervidus and B. californicus were separated by substantial branch lengths, suggesting that they might be separate species. However, Williams et al. [15] considered that the lack morphologically diagnostic traits and COI divergence between the two species as evidence that B. fervidus and B. californicus are conspecific.
There are three major hypotheses concerning the species status of B. californicus and B. fervidus. The first hypothesis proposes that B. californicus and B. fervidus are distinct species [16,17,20,21,26]. The second hypothesis proposes that B. californicus and B. fervidus are distinct species, and produce a hybrid subspecies, B. californicus consanguineus [34]. Finally, the third hypothesis proposes that B. californicus and B. fervidus are conspecific [15,16]. Operating under the phylogenetic species concept, our goal in this study is to test all three hypotheses simultaneously. We use data from neutral and adaptive genetic loci to examine their species boundaries. We first infer a phylogeny with three mitochondrial loci: COI, 12s RNA, and 16s RNA with specimens distributed across a broad geographic range, and exhibiting diverse setal phenotypes. Next, we expand our genetic sampling effort of specimens and genotype populations using neutral microsatellite loci to examine potential hybridization and species assignment. We predict that neutral microsatellite loci will have the power to identify introgression between B. fervidus and B. californicus [39,40].

Taxa examined
We included a total of 320 specimens associated with the B. fervidus SC, including the nominal B. californicus. We made an effort to include a diversity of setal color phenotypes associated with the B. fervidus SC (Fig 1) [15]. Exemplars of B. weisi (Thoracobombus) and B. insularis (Psithyrus) were selected as outgroup taxa based on recent Bombus phylogenies [21,26]. Ingroup taxa, exclusive to females were sampled throughout a major portion of their range in North America. We recorded setal color pattern data and locality information associated with queen and worker castes (S1 Table). We categorized specimens into four broad phenotype groups (Fig 1). These phenotype groups are based on previous taxonomic assessments of the B. fervidus SC [16,17,20]. Assignment of setal color patterns to specimens follow the schematic diagram presented in [15] and [17]. In addition to phenotypes, we assigned specimens to either B. fervidus and the nominal B. californicus following the species diagnoses from [20], [16], and [17]. In brief, B. californicus is much more variable than B. fervidus, and has been assigned three predominant setal phenotypes. Fig 1 presents the phenotypes (P) as P1, P2, and P3. Bombus fervidus is not considered to be variable in setal bandings, and is presented as P4 in Fig 1. Forty-nine field sites did not require specific permissions to survey bumble bees as the surveys were conducted in public spaces that had no specific geographic protections (Table 1). Furthermore, collection permits are not necessary as members of the B. fervidus SC are not protected under any state or federal laws. Four sites did require specific permits for bumble bee surveys as they took place in U.S. National Parks. The sites and corresponding permits are as follows: Pop Table 1). Permit details may be retrieved from National Park Service Research and Reporting System (https://irma.nps.gov/rprs/). Collection data associated with specimens used for this study have been digitized and deposited in the United States National Pollinating Insect Collection at Utah State University in Logan, Utah, U.S.A (S1 Table).

DNA extraction, amplification, and gene sequencing
We extracted genomic DNA from the mid-leg of a specimen using a modified Chelex 10 protocol following Strange et al [41]. DNA extracted in this manner was primarily used for microsatellite genotyping (i.e., Fragment Analysis), and was not especially successful when used in PCR aimed at amplifying gene fragments >500 base pairs. In this case, we also extracted genomic DNA using the Roche High Pure Template Preparation Kit (Roche Diagnostics GmbH, Germany) to obtain high quality genomic DNA suitable for downstream amplicon sequencing.

Phylogenetic analysis
The mitochondrial genes were examined separately and combined into a single partitioned dataset (1758 nucleotides) to infer a phylogeny with a Bayesian likelihood-based approach. Models of molecular evolution for each mitochondrial locus and codon position (COI) were first investigated with PartitionFinder v1.0.1 [46]. We implemented the model HKY+Gamma for 12S and 16S, HKY+I for COI first codon position, F81 for COI second codon position and HKY for COI third codon position. The Bayesian single-gene and concatenated phylogenies were estimated with MrBayes v3.2.1 [47] using two independent runs with three heated chains  [45]. Because we were potentially working with two different species in our study, we elected to use a universal bin set when scoring alleles for all specimens. This approach ensured that alleles were being consistently called with the appropriate microsatellites motifs with no a priori assumptions of species identity. Our method did not yield any ambiguous allele calls nor did we observe any "bin creep" [53], suggesting that the genotypes discovered in this study were suitable for downstream analyses.

Population genetic analysis
A Bayesian clustering method implemented in Structure v2.3.4 [54] was used to assign individuals to populations a priori. This method ensured that we did not base species identifications on the setal color phenotype the specimen displayed (Fig 1). We predicted that specimens that were grouped together based on microsatellite genotypes composed distinct genetic clusters separate from specimens in other predicted groups. The Structure algorithm in this way has been found to be useful in identifying distinct genetic clusters in other studies of bumble bees with cryptic phenotypes and evolutionary histories [18]. We used the admixture model in Structure, which assumes that individuals comprise K unknown genetic clusters, to which an individual can be fractionally assigned. This allowed us to group specimens based on their genotype without prior delineation to a population or species. In this case, the inferred population represents a genetic cluster and would illuminate any contemporary admixture of genes. The alternative to the admixture model would be to set the modelling scheme to "no admixture" which would assume that populations are discrete, where genotypes were assigned to a genetic cluster in full (i.e., no fractional assignment). As we are testing whether B. fervidus and B. californicus are conspecific with gene flow among populations, incorporating admixture into the modelling framework would allow for fractional assignment to K population(s). Furthermore, the admixture model would allow us to detect if any hybridization at the microsatellite loci between the two species was evident in areas where the two color phenotypes are sympatric. We set the admixture model to run with 20,000 burnin steps and 100,000 samples, with 10 iterations for each K, where K ranged from 1 to 10. Testing a wide range of K ensured that we did not bias the assignment of genotypes to only one or two species.
To determine the optimal K (i.e., populations/species or genetic lineages), the distributions of the probability of the data (ln P(D)) and ΔK, as described by Earl and von Holdt [55] and Evanno and vonHoldt [55], were visualized with the web-based software program Structure Harvester [55]. To account for multimodality associated with individual Structure simulations, we averaged each individual's admixture proportions over the 10 replicates for the best K using Clumpp v1.1.2 [56]. Finally, in addition to Structure analyses, we combined the 13 microsatellite loci into a principal components analysis to determine if significant clustering of similar genotypes could be inferred.
After determining the appropriate species assignments and number of K genetic clusters, the probability of null alleles was estimated with the software program MicroChecker [57]. We then estimated pairwise linkage disequilibrium (LD) and deviations from Hardy-Weinberg equilibrium (HWE) across populations and loci with the web-based software program Genepop v 4.0.10 using default parameters [58]. Based on the genetic clusters inferred by Structure, we performed an analysis of molecular variance (AMOVA) to test for differences in genetic structure with Arlequin v3.5 [59]. We then tested for a correlation between pairwise estimates of fixation based on allele frequencies with geographic distance (Isolation by Distance) within the genetic clusters inferred from the Structure analysis with GeneAlEx v6.5 [60].

Phylogenetic analysis
Our inferred phylogeny based on the concatenated gene sequences recovered two distinct monophyletic groups with strong support (Bayesian Posterior Probability, BPP = 1.0) (Fig 2A). Our data recover a paraphyletic B. californicus sensu lato and a polyphyletic B. californicus sensu stricto and B. fervidus, but support the hypothesis that there are two phylogenetically distinct species-clade b and clade c-due to fairly long branch lengths separating them. Single gene investigations revealed similar topologies to the full evidenced set but with lower support for clades b and c, specifically, BPP COI = 0.89, BPP 12s = 0.84, and BPP 16s = 0.86. All three genes contributed to the inferred Bayesian phylogeny and were retained in all analyses. Examination of sequence divergence between clades b and c, revealed the COI gene to have 861 identical sites (95.7%) with an average sequence divergence of 1.67% between clades; 16s revealed 473 identical sites (97.1%) with an average sequence divergence of 1.66%; and 12s revealed 348 identical sites (94.8%) with an average sequence divergence of 5.04%. GenBank accession numbers for the three mitochondrial gene fragments of the 64 specimens are found in S2 Table. Species descriptions of B. californicus by Smith (1859) and B. fervidus by Fabricius (1798) did not capture the phenotype (setal color) variability associated with lineages inferred in our well supported phylogeny. While setal color variability has been documented in both species, taxonomic keys and diagnoses by Thorp et al. [20], Stephen [16], Mitchell [35], Koch et al. [17], and others do not account for the shared setal color polymorphisms uncovered in this study. Coloration patterns from the holotypes of both species have been recovered in the two clades, which impedes us from assigning taxonomic names to them. Clade c includes the least color variability, which has traditionally been assigned to B. fervidus. This clade contains individuals from phenotype 1 from the Coastal/South Sierra California, forming the subclade i ( Table 1, sites 39 and 40). Within the Intermountain West + Pacific Northwest h clade, individuals that exhibited no signs of admixed black setae on the dorsal regions of terga two and three of the metasoma were detected, which is typically attributed to B. californicus consanguineus (Fig 2A) [20].
Within the respective b and c clades, we found a degree of support for geographic structuring across lineages (Fig 2A). Specifically, within clade b, we found strong support (BPP = 1.0) for a Rocky Mountain clade g as sister to the populations distributed in the Intermountain West + Black Hills and the Pacific West clade d. An exception was a South Dakota specimen (CusterSD, DD13197) that was found within the Intermountain West clade d, but it was preceded by a node with poor support (BPP = 0.71). Within clade c, we found strong support for the Coastal/South Sierra California clade i as sister to a lineage that comprises specimens distributed from North Sierra California + Intermountain West to Eastern USA clade e. Within clade e we found low support (BPP = 0.61) for the sister relationship between the North Sierra California populations and the populations that comprises the Intermountain West to Eastern USA.

Population genetic analysis
Microsatellite genotypes corroborate the existence of two monophyletic groups inferred from the multi-gene phylogeny within the B. fervidus SC (Fig 2B). Structure analysis of the available genotypes revealed two major genetic clusters within the B. fervidus + B. californicus clade a (Fig 2A). The estimate of the optimal cluster is based on a Structure Harvester analysis that found the highest log likelihood of the inferred models of K to occur at K = 2 ( Table 2; Mean LnP(K|2) = -14577.2). Significantly less explanatory power was gained by additional clusters (ΔK = 954.68) [55] (Table 2). Finally, at six localities in our study, we found sympatric populations of clades b and c as evidenced by distinct microsatellite genotypes (Fig 3A), and the inferred phylogeny (Fig 2A).
In total, 93.8% of the 373 genotyped individuals were matched to the species identifications performed by the authors (B. fervidus or B. californicus) based on the classic setal color phenotypes found in taxonomic keys and field guides (Fig 1) [16,17,20]. Of the 209 specimens exhibiting the B. californicus phenotypes (phenotypes 1, 2, and 3) (Fig 1), 5.26% were assigned to genetic cluster 2 (Fig 2B). While the genotypes of 10 specimens were assigned to genetic cluster 2, they exhibited phenotype 1 (B. californicus sensu stricto). Eight of the specimens were collected in Pinnacles National Park and two of the specimens were collected in Yosemite National Park (S1 Table). An additional individual assigned to genetic cluster 2 exhibited phenotype 3 (B. californicus sensu lato), and was collected in Owyhee County, Idaho. Of the 172 specimens exhibiting the B. fervidus phenotype (phenotype 4), 7.3% were assigned to genetic cluster 1 (Fig 2B). While the genotypes of 12 specimens were assigned to genetic cluster 1, they exhibited phenotype 4 (B. fervidus sensu stricto). Eight specimens were collected in the Toiyabe Range in Lander County, Nevada, one specimen was collected in the Bitterroot Valley in Missoula County, Montana, two specimens were collected in Logan Canyon in Cache County, Table 2. Table of four probabilities of model fit implemented with the Evanno method associated with different values of K (i.e., genetic clusters) based on 13 microsatellites implemented in Structure Harvester. Bold text represents the indices that suggests the value of K that best predicts the microsatellite genotypes assigned in the Structure analysis. Utah, and one specimen was collected in Guardsman Pass in Wasatch County, Utah (S1 Table). Principal components analysis estimated 202 principal components for the 13 genetic loci used in our study. Principal component 1 explained 4% of the variance in the genotype data and principal component 2 explaining 6% of the variance in the genotype data ( Fig 3B). While the number of principal components is large, visual inspection of principal components 1 plotted against principal components 2 revealed two distinct clusters associated with the genotype assignments inferred from the Structure analyses ( Fig 2B). Furthermore, AMOVA results found that 14.66% of the genetic variation was partitioned among the two major genetic clusters, 14.10% among individuals within populations, and 71.24% among individuals within sites (Table 3). Overall F ST among populations is 0.15 (P < 0.001) and F IS is 0.17 (P < 0.001). Microsatellite genotype data is available at https://doi.org/10.6084/m9.figshare.6972518.v1.   To determine HWE and LD associated across populations within each genetic cluster (i.e., clade), we first separated out individuals based on genetic cluster assignment supported by Structure analysis. After partitioning the specimens by genetic clusters, we used Micro-checker to determine if any loci by population combinations exhibited evidence of null alleles or stuttering. From our analyses of population within the genetic cluster 1 (clade b), we elected to remove BTMS0044 as it was found to be in LD with BTERN02. Finally, BL15 and B124 did not amplify in several specimens in genetic cluster 2, and were not used in any further analyses. After the removal of problematic loci, we retained the following eight loci for further analyses with specimens assigned to genetic cluster 1: BT10, B96, BTERN02, B124, BL15, BT28, BTMS0086, BTMS0066, and the following eight loci for specimens assigned to genetic cluster 2: B126, BT10, B96, BTERN02, BTERN01, BTMS0044, BT28, BTMS0066.

Source of Variation
Across genetic cluster 2 (clade c) we detected a strong effect of geographic distance on patterns of allelic fixation (Mantel Tests, r = 0.39, P = 0.03), with estimates of pair-wise linearized F ST ranging from 0 to 0.26 (Fig 4A). We also detected a strong effect of geographic distance on patterns of allelic fixation within genetic cluster 1 (clade b) (Mantel Tests, r = 0.56, P = 0.01), with estimates of pairwise linearized F ST ranging from 0 to 0.53 (Fig 4B).

Discussion
Globally, there are more than 260 species of described Bombus [26]. Bumble bees are typically regarded as well studied relative to other Hymenoptera given that they represent the only extant genus in the tribe Bombini (Apidae), particularly in North America. In our study, we uncovered two well-supported lineages made up of populations that exhibit shared setal color polymorphisms across clades b and c in the B. fervidus SC (Fig 2A). Cluster assignment of 13 microsatellite loci corroborates the results of the inferred phylogeny, specifically, that two distinct genetic lineages are present in areas where the species are broadly sympatric (Figs 2B  and 3).

Phylogeny and population genetic structure
The recent and rapid diversification within the B. fervidus SC was likely driven by climate change and glacial oscillations associated with the late Pleistocene [18,61,62]. Simple pairwise examination of the average levels of divergence across COI between clade b and c is 1.67%. The observed level of divergence is below the 2% level that is often considered reflective of what delimits a species [15]. This suggests that divergence from a common ancestor likely occurred less than~1 million years ago based on estimates of mitochondrial divergence with respect to time [63]. However, in addition to COI, we considered that distinct 16s and 12s haplotypes are characteristic of individuals associated with both clades (Fig 2). Because there are no morphological differences other than setal coloration [16,17,20], and color patterns found in the holotypes of B. californicus and B. fervidus are represented in the two clades recovered, there is no evidence to assign species names to either clade. Future studies including molecular and/or morphological data from molecular specimens and the holotypes will be crucial to establish the species boundaries within this SC.
Despite shared setal color polymorphisms in the B. fervidus SC, we reject the hypothesis that the complex is composed of a single species. However, due to the lack of evidence other than setal coloration, we cannot reject the hypothesis that B. fervidus and B. californicus are conspecific, because the type specimen of both could be included in a single lineage. Microsatellites are powerful molecular tools that have the capacity to uncover introgression between cryptic species [39,40,64]. Our microsatellite data found sympatric populations within the B. fervidus SC to be reproductively isolated, with no evidence of introgression (Figs 2B and 3). Therefore, we reject the hypothesis that B. californicus and B. fervidus produce the hybrid B. californicus consanguineus (phenotype 3).
Williams et al. [15] and Stephen [16] state that current taxonomic tools are not useful for differentiating some closely related species, including B. californicus and B. fervidus. We agree with both Williams et al. [15] and Stephen [16] that B. californicus and B. fervidus cannot be identified to species using setal color in certain parts of the geographic distribution based on the results generated in this study (Figs 2 and 3). For example, based on the data in this study, we found that specimens, which would be identified as B. californicus in southern California (sites 39 and 40) (Fig 1), could be assigned to clade c which are made up of populations exhibiting the "B. fervidus" phenotype (Figs 1, 2A and 3A). Furthermore, populations distributed in the Toiyabe Range in Nevada, and the Bear River Mountain Range in Utah identified as B. fervidus based on the absence of black setae on the dorsum of the metasoma [20] could be assigned to clade b which are made up of populations exhibiting the "B. californicus" phenotype (Fig 2A). Given the results of our study, setal color patterns appear to be of limited taxonomic use. Despite the degree of crypsis associated within the B. fervidus SC, we assigned 89% of the B. californicus specimens used in our study to clade b (genetic cluster 1) with microsatellite genotypes (Fig 2A). With B. fervidus, we assigned 93% of the specimens to clade c (genetic cluster 2) with microsatellite genotypes based on recognized phenotypes of the SC (Figs 2 and 3). If the type specimens were samples, and the clades could be assigned to correct names, species assignment to either B. californicus or B. fervidus based on current taxonomic tools could be possible in some areas of North America. Future research on the B. fervidus SC should evaluate taxonomic characters like wing venation as it has been useful method for distinguishing between cryptic bumble bee species [65].

Cryptic speciation and mimicry
Cryptic speciation is found in a diversity of bumble bee clades, as well as other invertebrate and vertebrate taxa. For example, species in the B. lucorum SC (B. lucorum, B. magnus, and B. cryptarum) are indistinguishable from each other using taxonomic methods of identification (i.e., setal color patterns), and can only be diagnosed to species using molecular techniques such as barcoding [22]. Müllerian mimicry is a well-documented phenomenon where sympatric species share a common aposematic phenotype to warn predators of their noxious chemical composition [66]. Like the B. lucorum SC, the B. trifasciatus SC in east Asia is another example where species identification based on setal color pattern fails to differentiate between species. Genetic divergence among the B. trifasciatus SC is hypothesized to be generated by Himalaya orogeny with Müllerian mimicry the likely factor shaping cryptic speciation among unrelated bumble bees [66]. Other cryptic species complexes among the bumble bees include the B. patagiatus and B. hypocrita SCs of Asia [67], and the B. ephippiatus SC of Mesoamerica [18,65]. Examples of cryptic speciation facilitated by Müllerian mimicry is observed in butterflies (Heliconius spp.) [68], spider wasps (Pompilidae) [69], velvet ants (Mutillidae) [70], and frogs (Dendrobates spp.) [71].
Our phylogenetic and population genetic analytical framework discovered two distinct lineages exist in the B. fervidus SC, and that they can occur in the same habitat space. The results of our study suggest that when both species are sympatric they appear to be phenotypically divergent (Figs 2, 3A), which would indicate they are not mimicking each other. However, while they can be some sympatric, some authors have suggested that B. fervidus and B. californicus inhabit different habitat niches [36]. Future research could examine how climate, mimicry, and floral niche might contribute to their ability to coexist in some portions of their range, but not in others [27,28,65,72]. For example, Pleistocene climate variation has been hypothesized to not only drive genetic divergence in B. huntii, but also differences in bioclimatic niche, and potentially in setal color variation [72]. In addition to B. huntii, there is evidence that historic climate variation has shaped patterns of genetic divergence and habitat partitioning across closely related bumble bee species, and is hypothesized to have also cascaded down to changes in setal color patterns [6,18]. In the B. fervidus SC, shared setal coloration patterns between the two clades is potentially a result of Müllerian mimicry where the model is not in the SC. For example, at MPG Ranch in Missoula, Montana (site 11), populations belonging to genetic cluster 1 (clade b) and 2 (clade c) are detected, and exhibit divergent phenotypes (Figs 1 and 3A). Other sites where B. fervidus SC species are sympatric (both genetic clusters 1 and 2 are detected), exhibit divergent phenotypes, and show no evidence for introgression include the North Sierra Nevada Mountains (sites 2, 3, 4), the Bear River Range (site 25), the Toiyabe Range (site 13), southeastern Oregon (site 14), and the Black Hills (site 23) (Fig 3A).
Among bumble bees, Müllerian mimicry is a common phenomenon, and has been documented across a diversity of communities around the globe [18,19,65,66]. For example, bumble bees in eastern North America share similar yellow and black setal coloration patterns that can make it difficult to correctly identify some individuals to species [15]. Outside of bumble bees, participation in a Müllerian mimicry ring with species of a completely different, or closely related taxonomic group is also common. For example, there is strong phylogenetic evidence that spider wasps (Pompilidae) and velvet ants (Mutillidae) have exhibit similar phenotypes when they sympatric [69].

Conservation implications and future work
Bumble bees are well regarded for their value in agricultural ecosystems as they are efficient pollinators of a diversity of crops [4,5,73]. However, there is global concern for bumble bee decline due to economic activities associated with human growth and expansion, namely the shuffling of Hymenopteran disease due to movement of bee colonies to meet pollination demands, as well as increased urbanization and agricultural intensification [10,[74][75][76]. B. fervidus in particular has been associated with decline at regional scales [30,32], and has been found to be highly susceptible to a suite of pathogens [77]. Despite its co-distribution with B. fervidus throughout western North America, B. californicus does not appear to be associated with high levels of pathogen incidence [74,77].
Cryptic speciation in bumble bees is well documented [22,65,67]. Application of a phylogenetic and population genetic analytical framework has revealed that using setal color patterns as a way to diagnosis species might not be useful in discriminating between closely related species [22,67]. The inability to discriminate between species due to cryptic speciation has significant implications to both biodiversity conservation and agriculture. For example, Carolan [22] discovered that B. hypocrita and B. patagiatus exhibit a similar phenotype, but are phylogenetically distinct. The widespread Russian B. patagiatus are reared by commercial greenhouse growers for pollination of food crops. Because they can be indistinguishable from B. hypocrita, which are found in Japan, there is potential for misidentification and ultimately, the unintended movement of B. patagiatus and B. hypocrita between continental Asia and Japan. The movement of non-native species or populations has the potential to displace native bumble bee species or populations [78,79], cause a reproductive disturbance with native species [80], and potentially facilitate the spread of disease [76,81,82].
A prevailing hypothesis associated with bumble bee decline includes the introduction of novel pathogens or pathogen strains [10,11]. Given the differences in pathogen prevalence between B. fervidus and B. californicus, we suggest that researchers treat the two species in the B. fervidus SC differently in the context of conservation, ecology, and evolution. Our results show that the two lineages are phylogenetically distinct (Fig 2A), with no evidence for introgression when sympatric (Figs 2B and 3). Given the pronounced genetic differences in the species, treating them as separate will allow for a more robust assessment of their conservation needs and disease profiles.
Despite the inability to identify the individuals to species based on current taxonomy, there is potential for alternative, non-destructive ways to ensure correct species identification [83,84]. Specifically, we found that microsatellite genotypes have the capacity to differentiate species, even when they are sympatric (Figs 2 and 3). While we propose that a synoptic collections of the bumble bee community be created when conducting ecological research, we have found that taking a tarsal clipping from the mid-leg for DNA extraction and subsequent genotyping is possible, which avoids sacrificing the whole individual, allowing it to continue with its contribution to the nest economy [83]. Knowledge about the evolutionary processes associated with the formation of a species is required in conservation biology [2,14,85]. In this study, we demonstrate that populations that compose B. fervidus SC lineages are cryptic, yet form well supported clades. To reduce the complex to a single species based on the inability to identify them to species using morphological traits will likely obscure the host-pathogen dynamics associated with the species, and ultimately hinder effective action on their conservation and management.