Bumblebees take the high road: climatically integrative biogeography shows that escape from Tibet, not Tibetan uplift, is associated with divergences of present-day Mendacibombus

Many claims that uplift of the Qinghai-Tibetan plateau (QTP) drove the divergences of extant high-elevation biota have recently been challenged. For Mendacibombus bumblebees, high-elevation specialists with distributions centred on the QTP, we examine broader explanations. We extend integrative biogeography to cover multiple contributing factors by using a framework of sequential filters: 1) molecular evidence from four genes is used to estimate phylogenetic relationships, with time calibration from a published estimate; 2) spatial evidence from current distributions is combined with the phylogeny and constrained by a model of short-distance dispersal along mountain corridors to estimate ancestral distributions by both S-DIVA and S-DEC analysis; 3) geological evidence from the literature is used to constrain when high mountain ranges were uplifted to become potential corridors; and 4) climatological evidence from Mendacibombus niche-evolution reconstructions and from palaeoclimate simulations is used to constrain when habitat was suitable in key gaps within corridors. Explanations for Mendacibombus distributions can be identified that require only short-distance dispersal along mountain corridors, commensurate with the limited dispersal ability observed for bumblebees. These explanations depend on the timing of uplift of mountain ranges, regional climate change, and climate-niche evolution. The uplift of the QTP may have contributed to the initial Oligocene divergence of the common ancestor of Mendacibombus from other bumblebees, but for the first two thirds of the history of Mendacibombus, only a single lineage has present-day descendants. Divergence of multiple extant Mendacibombus lineages coincided with the Late Miocene–Pliocene uplift of externally connecting mountains, combined with regional climate cooling. These changes provided greater connectivity of suitable habitat, allowing these bumblebees to disperse out of the western QTP via new high bridges, escaping along the mountain corridors of the Tian Shan and Hindu Kush ranges, reaching eventually far to the west (Iberian Peninsula) and to the north-east (Kamchatka).


Introduction
Over the last decade there have been claims that divergence patterns in many groups of animals and plants have been driven by the high uplift of the Qinghai-Tibetan plateau (QTP), citing similar dates for both within the last 15 Ma (Favre et al. 2015,

Research
Bumblebees take the high road: climatically integrative biogeography shows that escape from Tibet, not Tibetan uplift, is associated with divergences of present-day Mendacibombus Renner 2016). However, growing geological evidence now supports high uplift of the QTP as having occurred much earlier, before 45 Ma (Favre et al. 2015, Renner 2016. We explore patterns of distribution and divergence in one group of specialist high-mountain bumblebees and its coincidence with mountain uplift and climate change in and around the QTP. For these cool-adapted insects, mountain chains, rather than representing barriers to dispersal as they do for some organisms, might act as corridors facilitating dispersal.

Bumblebees and cool distributions
Bumblebees (genus Bombus Latreille) are well suited to studies of distribution (Williams 1985(Williams , 1998 because: 1) they include a manageable number of species (ca 265 world-wide); 2) their distributions are well recorded because they are large, colourful and conspicuous, and active under similar conditions to field biologists; 3) the different species have similar behaviour and ecology as relative floral generalists, with many species co-occurring in the same habitats; and 4) they may be especially poor at long-distance dispersal and establishment (see below), so that their indigenous distributions are not near equilibrium with the global distribution of potentially suitable habitat (as evidenced by the many successful introductions to new regions: Williams 1998, Goulson 2010. Departures from equilibrium can be informative about history. Bumblebees world-wide are associated primarily with habitats in cool climates, especially in mountains, predominantly in the northern hemisphere but also in South America (Williams 1998: his Fig. 15). This preference can be explained by bumblebees: 1) sharing a strongly developed ability for thermoregulation, especially when foraging in cold weather (Heinrich 1979); 2) using nest sites that are thermally insulated by being often underground and by using nest material left from their former occupants, usually small mammals or birds (Williams et al. 2008, their Table III); and 3) being able to hibernate -to pass long predictable annual seasons of adverse climate (usually winters) at relatively low cost in terms of food requirements (cf. the large honey stores required by honeybee colonies).

Mendacibombus as high mountain specialists
The bumblebees of the subgenus Mendacibombus are a group with 12 known species that are recognised as the sistergroup to all of the other extant bumblebees (Williams 1985, Kawakita et al. 2004, Cameron et al. 2007). Mendacibombus species have been described as being closely associated with high mountains (Skorikov 1923(Skorikov , 1931, where they occur but are uncommon in flower-rich subalpine and alpine grassland (meadow) habitats (Rasmont 1988, Williams 1991, Rasmont and Flagothier 1996, Williams et al. 2009).
For much of the history of Mendacibombus species, and despite minor variations, montane meadow can be expected to have been suitable habitat. A recent revision of Mendacibombus species world-wide includes many new samples, especially from Asia (Williams et al. 2016: Fig. 1a, with mountain ranges labelled in Fig. 1b). For terrestrial organisms, individuals within any one species usually share relatively similar climate requirements. This results in a characteristic shift in the elevations a species inhabits when the species occurs at different latitudes, as it tracks similar climates (Cogbill and White 1991). As expected, Mendacibombus bumblebees show this pattern within species (Fig. 2). They also show a similar pattern among species, with most species at any particular latitude overlapping broadly in elevation (Fig. 2). Therefore they probably shared similar climate niches for much of their shared history. Thus for this group, mountains would not have been barriers to dispersal in the same way that they are for other, more lowland bumblebees (e.g. B. impatiens Cresson in North America: Williams et al. 2014;B. pascuorum (Scopoli) in Europe: Widmer and Schmid-Hempel 1999).

Range growth: dispersal and establishment
For bumblebees, dispersal across short distances (here referred to as short-distance dispersal, SDD) within a suitable habitat (and at least close to an existing population) should be distinguished from the process of long-distance dispersal (LDD) that requires a jump over a substantial barrier of unsuitable habitat. The distinction is important because dispersal requires not only moving but also the founding and establishment of a new isolated population. Bumblebee species' distributions have often been documented to change, often spreading year after year along corridors of suitable habitat (Macfarlane and Gurr 1995, Smissen van der and Rasmont 1999, Macdonald 2001, Schmid-Hempel et al. 2013. Some authors have believed that bumblebees were less able to undergo LDD with establishment (Panfilov 1957, Ito and Sakagami 1980, Ito 1987, Pekkarinen and Teräs 1993, Macfarlane and Gurr 1995, Estoup et al. 1996, Lecocq et al. 2016), but others disagree (see below).
For bumblebees it is the dispersal ability of queens in spring that is relevant to range growth by the establishment of new populations. Bumblebees reproduce in immobile colonies, which in most cases persist for only a few months and are founded by solitary queens in spring (Sakagami 1976). Queens also leave colonies later in the summer to hibernate, but at this time they spend less time outside the nest than in spring (plots of monthly activity in Prys-Jones andCorbet 1987, Williams et al. 2014) and they are then much heavier with food reserves (Alford 1975). Despite the ability of spring queen bumblebees to fly long distances, sometimes even across sea (Haeseler 1974, Mikkola 1984, it is unknown but perhaps unlikely that they would arrive fit enough to found colonies (Williams et al. 2015b). Furthermore, because sex determination in bumblebees depends on heterozygosity (Crook and Crozier 1995), mating between siblings from isolated founder colonies would result in the rearing of many diploid males (Duchateau et al. 1994). This is likely to be exacerbated because most bumblebee species are monandrous (Estoup et al. 1995, Schmid-Hempel andSchmid-Hempel 2000). Males do not forage for the colony and therefore a high proportion of males in early broods would inhibit colony development and represent a severe cost to colony reproductive success. Therefore unless many queens in a fit physiological state arrive in an isolated new area simultaneously, the result of LDD is likely to be starvation and . The three largest current disjunctions in distribution ('gaps') are indicated by grey lines and by the letters a (Aladagh mountains), b (Balkan-Dinaric mountains), and g (Stanovoy-Okhotsk mountains), which follow presumed earlier paths of least resistance for dispersal at higher elevations (cf. minimum straight-line distances between the nearest records: a, 1520 km; b, 1450 km; g, 2300 km). (b) Map with black labels showing the principal mountain ranges named in Table 1 or discussed in the text; grey lines show the areas of endemism recognised in Table 1 (grey large capitals for area abbreviations) with dotted lines following presumed earlier paths of least resistance for dispersal at higher elevations (in addition to where areas of endemism are adjoining); grey small labels show selected deserts. Map projected onto a sphere in ArcGIS using World_Shaded_Relief basemap © 2014 ESRI.
inbreeding depression, making long-term establishment of an isolated new population less likely.
Assessing the potential for LDD by bumblebees empirically, attempts to introduce the relatively invasive B. terrestris (Linnaeus) to islands have demonstrated that as many as 100-150 queens in one nesting season have been needed in order to guarantee establishment (Macfarlane and Gurr 1995). In an example that appears at first sight contrary to this expected difficulty in establishment, Schmid-Hempel et al. (2007) estimated that the population of B. terrestris in Tasmania is descended from perhaps as few as two successful queens, probably transported there on ships. These queens were unusually pathogen-free and apparently pre-selected from 'above-average families' (Schmid-Hempel et al. 2007). Nonetheless, the early generations in Tasmania still produced a high frequency of diploid males (Buttermore et al. 1998). In another well-known case, B. hypnorum (Linnaeus) has recently become established in Britain (Goulson and Williams 2001). But rather than flying unaided across the English Channel from the European mainland (Crowther et al. 2014), it may be that this species, which has a stronger preference for areas inhabited by people than many other north European bumblebees (Løken 1973), could have arrived in Britain by assisted transport, for example by groups of queens that were hibernating in flower pots being transported by people from garden-plant nurseries in Europe (B. terrestris has been found hibernating in the soil in flower pots: F. Adamson in McCluskey 2012). It may also be significant that B. hypnorum is very unusual among bumblebees for being polyandrous (Estoup et al. 1995, Schmid-Hempel andSchmid-Hempel 2000). When establishing new small isolated populations, polyandry should reduce inbreeding depression, placing B. hypnorum at a particular advantage for introduction.
Some recent genetic studies have agreed that LDD over wide barriers appears to be rare for bumblebees (Darvill et al. 2006, Goulson et al. 2011. Other recent genetic studies have claimed support for frequent LDD for bumblebees, but actually they have dealt either with dispersal within regions with widespread broadly suitable habitat requiring only SDD steps (Lepais et al. 2010, Cameron et al. 2011, or with dispersal to islands with pre-existing conspecific populations that should remove the problems of inbreeding depression (Kraus et al. 2009). Therefore in general, it may be reasonable to expect that bumblebee species rarely establish in areas that lack existing populations if they arrive by unaided LDD across wide barriers of severely unsuitable habitat. Thus SDD through contiguous suitable habitat is expected to be the more important biogeographic process of range growth for bumblebees.

Assessing changing distributions and divergence
We seek a simple explanation for the current distribution and divergence of all present-day Mendacibombus species. We have no fossils of Mendacibombus from which to study speciation and extinction of all lineages and the group is too small for estimating these from phylogeny (Morlon et al. 2011). We can integrate complementary approaches that represent some of the many factors associated in time with the Figure 2. Plot of elevation (y-axis) against latitude (x-axis) for the 4416 specimen records of the 12 species of the bumblebee subgenus Mendacibombus world-wide updated from Williams et al. (2016). White letters show locations of the three gaps (a b g) from Fig. 1a; large light grey triangles show the ranges of latitude and the maximum elevations of mountain ranges in their higher sections along the grey curves in Fig. 1a; large dark grey triangles show the maximum elevations of mountain ranges in their lower sections (bottlenecks) along the grey curves in Fig. 1a. A key to the symbols for species is shown to the left. distribution and divergence of present-day species, including the uplift of the QTP and of the associated mountains. This requires that we assess whether corridors of habitat suitable for species' SDD are sufficient to explain their current distribution across Asia and Europe. It includes considering whether climate conditions are likely to have been suitable for SDD in the appropriate past periods within the three current large 'gaps' (a b g) in their distribution in Fig. 1a. To integrate the major governing factors, we use a framework with sequential filters -to constrain progressively the explanations to those most likely: 1) evidence from genes is used to estimate phylogeny and dates for divergence events between extant species; 2) evidence from current distributions is used with the phylogeny to estimate ancestral distributions; 3) evidence from the geological literature is used to estimate when mountain ranges showed high uplift; and 4) evidence from Mendacibombus niche-evolution reconstructions and from palaeoclimate simulations is used to estimate whether suitable habitat was present in these corridors within the appropriate temporal windows.

Estimating and dating phylogeny
A time-calibrated phylogeny for species within the subgenus Mendacibombus was estimated with a Bayesian analysis of four genes using BEAST (ver. 1.8.2, accessed 2015: < www.beast. bio.ed.ac.uk >). This linked trees (Heled and Drummond 2010) for the mitochondrial genes COI and 16S and for the nuclear genes PEPCK (phosphoenolpyruvate carboxykinase) and opsin (long-wavelength rhodopsin copy 1). Methods for the extraction, amplification, sequencing, and alignment for the gene sequences are described in Williams et al. (2016). This study examined 4416 specimens and extracted DNA from 171 specimens. We reduced the dataset to a single example sequence for each gene for each of the species, using the species recognised from an integrative analysis of COI-coalescents and morphology (Williams et al. 2016). Details of primers are in Supplementary material Appendix 1 Table A1 and sequence accession numbers are in Supplementary material Appendix 1 Table A2. As outgroups, we include several close relatives of Mendacibombus, using a selection of species from other bumblebee subgenera (Williams et al. 2008) for which relationships are strongly supported (Cameron et al. 2007 (Alpigenobombus) nobilis Friese. Sequences were aligned using the MUSCLE procedure and amino acid translation for coding regions was checked within MEGA (ver. 6.06, accessed 2014: < www.megasoftware.net >). Different nucleotidesubstitution models were used for different genes: the best model for our COI-barcode and 16S data according to the Bayesian information criterion (BIC) obtained from MEGA is the general time-reversible model with a gamma-frequency distribution of changes among sites (GTR  G); the best model for the PEPCK exon (boundaries and alignment follow Cameron et al. 2007) is HKY (Hasegawa-Kishino-Yano); for the opsin exon is HKY  G; and for the PEPCK and opsin introns is HKY. For B. makarjini, B. marussinus and B. himalayanus, we could not acquire sequences of 16S, opsin, or PEPCK. The clock prior was set to the uncorrelated lognormal relaxed clock, the tree speciation prior was a Yule process, and the chain length for the Markov-chain Monte Carlo (MCMC) algorithm was set to two billion generations, with sampling of the resulting trees every 200 000 generations. The sample of trees was examined using Tracer (ver. 1.6.0, accessed 2013: Drummond and Rambaut 2007), which showed that stationarity had been achieved within the first 1% of the MCMC generations. The large number of MCMC generations was needed with these data to increase effective sample sizes. A Bayesian maximum-clade-credibility tree was obtained from the post burn-in tree sample after rejecting the first 1% of sampled trees using TreeAnnotator (ver. 1.8.2, as for BEAST).
The phylogeny had to be calibrated with a date from a molecular study because there are no known fossils of Mendacibombus species. Hines (2008) estimated the age of divergence between the subgenus Mendacibombus and all of the remaining extant bumblebees to be at 40-25 Ma, although this estimate was then fixed at 34 Ma for most of her subsequent analyses (consistent with earlier suggestions by Skorikov 1923, Williams 1985. We used 34 Ma for the mean and initial age estimate for the most recent common ancestor of Mendacibombus plus all other bumblebees (tMRCA), assuming a normal distribution with a standard deviation of 4 Myr, which places the 5% tails of the tMRCA distribution at 41-27 Ma (late Eocene to mid Oligocene). We truncated the distribution at 50 Ma and 20 Ma to prevent extreme values.

Ancestral distributions
To assess the roles of SDD, vicariance, and extinction in explaining current distributions from the phylogeny, we used Bayesian methods implemented in the RASP package (ver. 3.0, accessed 2015: Yu et al. 2014) with a corridor model. Statistical divergence-vicariance analysis (S-DIVA: Yu et al. 2010) assumes that all distribution changes happen only at nodes (associated with speciation). This contrasts with a more gradualist model of distribution changes happening along branches in parametric methods (Ree and Sanmartin 2009), such as the analogous statistical dispersal-extinctioncladogenesis analysis (S-DEC, also known as Bayes-Lagrange analysis: Yu et al. 2014). Both procedures take account of the uncertainty in the estimate of the species' phylogeny by using a sample of 10 000 trees from the phylogenetic analysis of the combined gene data (after a 1% burn-in). S-DIVA was set to allow 'extinctions' (area-unit extirpations).
Area units were chosen by agglomerating the 4416 specimen records for Mendacibombus species (Fig. 1a) Fig. 1b). The number of areas recognised within a species' range was kept low because the biogeographic analyses we employ do not have mechanisms by which two daughter species can both inherit ranges consisting of multiple areas (Lamm andRedelings 2009, Ree andSanmartin 2009). The range of B. margreiteri was coded without including the narrow disjunct southern presence in the Burhan Budai mountains ( Fig. 1a, b) in order to keep the coded ranges small (this is interpreted as a relatively recent dispersal, see Discussion). Ancestral ranges are constrained to be a maximum of three contiguous areas in order: 1) to exclude LDD jumps between non-contiguous areas (Introduction); 2) to not greatly exceed the size coded for current ranges of species; 3) to offset the bias of analyses towards combining all areas from daughter species into widespread ancestral ranges; and 4) to prevent reconstructed ranges from being highly disjunct (Ronquist 1996, Ree and Smith 2008, Lamm and Redelings 2009, Ree and Sanmartin 2009).
To seek an explanation for current distributions using only SDD, a simple model of potential mountain corridors (assuming suitable climates in appropriate periods) was defined for this analysis (Fig. 3, Fig. 1b). This summarises a set of SDD-based dispersal events (excluding LDD 'jumps') permitted between neighbouring areas (in either direction) assuming that suitable habitat is present. The dispersal model is implemented within both analyses by constraining the set of permitted narrow contiguous distributions for each lineage (GH GHP GP GPT GPC HP HPT HPC PT PTM PTC PC PCA TM CA, from Fig. 3).

Ancestral niches and palaeoclimate suitability
To model the potential climate niche of Mendacibombus species, seven climate variables important in governing the distribution of vegetation were included initially as potential predictors of Mendacibombus distribution: annual precipitation; annual variation in precipitation; maximum monthly precipitation; maximum monthly temperature; mean annual temperature; minimum monthly precipitation; and minimum monthly temperature (these variables also reflect the effects of elevation). We would prefer to exclude any measurements that include winter and nocturnal conditions, times when these bumblebees are insulated underground, but corresponding variables are unavailable for palaeoclimates. Based on the seven initial variables, we also estimated two climate indices: aridity and continentality, calculated using Map showing mountain-range names in Fig. 1b. Authors of species' names: avinoviellus (Skorikov), convexus Wang, defector Skorikov, handlirschianus Vogt, himalayanus (Skorikov), makarjini Skorikov, margreiteri Vogt, marussinus Skorikov, mendax Gerstaecker, superbus (Tkalcu), turkestanicus Skorikov, waltoni Cockerell. Figure 3. Diagram representing a corridor model, encompassing a set of the short-distance dispersal events permitted (in either direction) between the areas defined in Table 1, based on their geographical proximity and the likely disposition of suitable habitat within mountain ranges and favourable climates in the past (see the text). 467 the formulae provided by Valencia- Barrera et al. (2002). Maps for these nine variables were compiled both for the present and for past epochs. Climate data for current conditions (1950-2000) were extracted from the WorldClim database (< www.worldclim.org >; Hijmans et al. 2005) and managed at a resolution of 2.5  2.5 minute grid cells (ca 5  5 km). Hadley Centre coupled ocean-atmosphere general circulation models (HadCM3L) were used for the Miocene and Pliocene. These models incorporate the effect of changes in atmospheric CO 2 concentration (Beerling and Royer 2011). The simulations were generated at a 2.5°  3.75° grid resolution, corresponding to: 1) a 400 ppm CO 2 Late Miocene simulation representing the early warm period of the Miocene and the Mid Miocene Climate Optimum around 15 Ma; 2) a 280 ppm CO 2 Miocene simulation representing the cold and dry conditions prevailing around 8 Ma; and 3) a 560 ppm CO 2 Pliocene simulation of the conditions around the Mid Pliocene warming event (3.6-2.6 Ma) (Beerling et al. 2009and Bradshaw et al. 2012. These time periods were chosen based on the results of the dated phylogeny to cover the principal periods of Mendacibombus dispersal. Contemporary site records of Mendacibombus bumblebees (n  493 sites) were grouped into 378 cells on a 2.5  2.5 minute grid. This occurrence information and the nine climate variables were used to discern which variables have a high probability of explaining the current distribution of Mendacibombus. In order to minimize the effects of historical events and dispersal limitations (Acevedo et al. 2012), we first delimited an appropriate compromise extent for the study area (geographical background), choosing the river basins of level 2 (subregions) that have presence observations and exclude areas disjunct by more than 1000 km from extant observations. Watersheds from the WaterBase project (< www. waterbase.org >) were used to estimate the relevant climate variables for the accessible regions. We removed highly correlated climate variables within the geographical background area, eliminating sequentially the variables with a variance inflation factor (VIF) lower than 10, following Neter et al. (1990). We selected the variables with the highest capacity for discriminating between the environmental conditions in the cells with recorded presences from the conditions prevailing in the geographical background area by applying a generalised linear model (GLZ, with a binomial distribution, logit link) according to the use-availability procedure of the resource-selection functions (Johnson 1980): the Mendacibombus occurrences (use) were related to a random sample of 10 000 cells from the geographical background (availability). A type 3 likelihood-ratio (type 3LR) test was applied in order to estimate the independent predictive capacity of each variable, controlling for the effects of the other variables. After a Bonferroni correction for multiple comparisons, the statistically significant variables were selected as those with a higher probability of being relevant to explain the occurrence of these species. All analyses were carried out in STATISTICA (StatSoft, ver. 12.0) and ModestR (García-Roselló et al. 2013; < www.ipez.es/modestr/ >).
Selected climate variables were used to estimate climate suitability for both extant species and their ancestors. We calculate the scale-invariant Mahalanobis distance of current and simulated palaeoclimate conditions in each grid cell from the centroid of the climate conditions existing in the cells with extant Mendacibombus records. The upper quartile of the Mahalanobis distances calculated from the current records was used as the threshold to delimit favourable areas and this threshold was used for all of the scenarios. Mapping these Mahalanobis distances allows us to recognise present and past climatically suitable areas according to the current conditions occupied by Mendacibombus (i.e. assuming niche conservatism).
In addition, we estimated the Mahalanobis distance of the palaeoclimate in each grid cell from the reconstructed ancestral climate niches for selected nodes in the phylogenetic tree. Ancestral climate niches were reconstructed by comparing the fit of each one of the previously selected climate variables (mean values) to different models of character evolution (Brownian, Ornstein-Uhlenbeck, white noise, and early burst), using the Akaike information criterion with a correction for small sample size (AICc) with the function Fit-Continuous in 'Geiger' (Pennell et al. 2014). We ran this 1) assuming no measurement error and 2) using an estimate of the error. We inferred rates of evolution (d) and associated ancestral values for each climate trait over the Mendacibombus chronogram in 'Phytools' (Revell 2012). From these results, we made three extrapolations for the Palaearctic background area: 1) the climate niche estimated for the most recent common ancestor (MRCA) of Mendacibombus was used to calculate the Mahalanobis distance from each grid cell in the Early Miocene palaeoclimate simulation for 15 Ma; 2) the climate niche for the MRCA of the makarjini-mendax group was used to calculate the Mahalanobis distance from each grid cell in the Late Miocene palaeoclimate simulation for 8 Ma; and 3) the climate niche for the MRCA of the margreiterimendax group was used to calculate the Mahalanobis distance from each grid cell in the Pliocene palaeoclimate simulation for 3.6-2.6 Ma. By using palaeoclimates and climate niches for the ancestors at the times in which they diverged, we assess the extent and connectivity of areas suitable for them without assuming niche conservatism.

Estimating and dating phylogeny
The best estimate of the dated phylogenetic relationships for all 12 species of the subgenus Mendacibombus shows that divergence of extant Mendacibombus species is likely to have begun in the Late Miocene epoch (ca 9.4 Ma, 95% confidence 13.6-6.1 Ma, Fig. 4).

Ancestral distributions
The most likely solution for distributions of ancestral Mendacibombus lineages from the S-DIVA analysis is shown in Fig. 5a and Supplementary material Appendix 1 Table A3. Results from the S-DEC analysis are similar ( Fig. 5b and Supplementary material Appendix 1 Table A4), but 'cost' slightly more in terms of the number of dispersal/vicariance/ extinction steps (15) than S-DIVA (13), making the S-DIVA solution (Fig. 5a) the simpler explanation. The S-DIVA solution includes eight principal events, which are mapped in Fig. 6. Therefore all of these events could be resolved using SDD as represented in Fig. 3, if suitable habitat were available at the appropriate times in all of the corridors. Three of the transitions in Supplementary material Appendix 1 Table A3-A4 require more than one SDD step, but these could have taken place over periods of as long as 3-5 Myr.

Ancestral niches and palaeoclimate suitability
In assessing variables for the potential climate niche of Mendacibombus species, VIF analysis eliminated four variables.
From the remaining five variables, the presence/background analysis selected aridity, annual variation in precipitation, maximum monthly temperature, and minimum monthly precipitation, as significant at the level of p  0.001 (Bonferroni correction for multiple comparisons, p  0.05/5  0.01). Reconstructions of ancestral states for these variables using 'Phytools' are shown in Supplementary material Appendix 1 Table A5 and Fig. A1.
Annual variation in precipitation and maximum monthly temperature showed a significant fit to a Brownian model, while aridity and minimum monthly precipitation were better explained by a white noise model (Supplementary material Appendix 1 Table A6). Estimation of potentially suitable current climates for all species of Mendacibombus shows the existence of favourable climate conditions scattered across the distribution of the group (Fig. 7a), including broad favourable areas in the northeast of Asia (gap g from Fig. 1a), but much less favourable areas in the gaps a and b   Table 1 and Fig. 1b; letter combinations in grey below terminals show species' current distributions; letter combinations above the nodes show the most likely optimised reconstructions of ancestral distributions, with the probabilities of these solutions shown below the nodes (alternative reconstructions with lower probabilities are shown in parenthesis). Numbers in circles at nodes identify the events discussed in the text. Grey bars represent mountain ranges uplifted to high elevation (see the text for references).  Table A3) projected onto present-day mountain ranges. (a) Vicariance and divergence between the most recent common ancestors (MRCAs) of (solid) Mendacibombus and (dashed) its sister group (outgroup) of all other extant bumblebees; (b) event 1, vicariance and divergence between the MRCAs of the superbus-convexus group and the makarjini-mendax group; (c) event 2, vicariance and divergence between the MRCAs of B.superbus and the waltoni-convexus group; (d) event 3, dispersal and divergence between the MRCAs of B.makarjini and the himalayanus-mendax group; (e) event 4, dispersal and divergence between the MRCAs of the himalayanus-avinoviellus group and the margreiteri-mendax group; (f ) event 5, dispersal and divergence between the MRCAs of the himalayanus-defector group and the marussinus-avinoviellus group; (g) event 8, dispersal, vicariance and divergence between the MRCAs of B.margreiteri and of the handlirschianus-mendax group; (h) dispersal of ancestral B.margreiteri to Kamchatka and Qinghai. Coloured ellipses encompass (overestimate) the approximate reconstructed ranges that ancestral lineages came to occupy after numbered events (Supplementary material Appendix 1 Table A3; colours do not correspond to the species in Fig. 1a). Doubleheaded arrows show dispersal/vicariance and divergence events and single-headed arrows show dispersal events. Relief map with hill shading, projected onto a sphere and viewed as though from above 46°N, 64°E using ArcGIS and World_Shaded_Relief basemap © 2014 ESRI. (Fig. 1a). During the warm Middle Miocene Climate Optimum, the projection based on the reconstruction of ancestral climate niche for the MRCA of Mendacibombus shows that favourable conditions in the south of Asia and Europe were widely scattered with multiple gaps (Fig. 7c). The extrapolation using the climate niche derived from current observations shows that only the QTP and the mountains central Asia and of the far northeast of Asia would be favourable at this time (Fig. 7b). However, during the cold and dry conditions of the Late Miocene, the reconstruction of the ancestral climate niche for the makarjini-mendax group shows a continuous band of favourable conditions along the southern Palearctic region extending as far west as Turkey, bridging gap a (Fig. 7e, cf. Fig. 1a). The extrapolation using the current climate niche again shows that only the QTP and the mountains of central Asia and of the far northeast of Asia would be favourable (Fig. 7d). At the time of the Pliocene ancestor of the margreiteri-mendax group (Fig. 4), the reconstructed ancestral climate niche again show a gap in suitability across gap a (Fig. 7g, cf. Fig. 1a). In general, extrapolations based on estimated ancestral climate niches appear to have a higher explanatory capacity than current climate niches in accounting for the biogeographic history of Mendacibombus, supporting the importance of climate-niche evolution for this group.

Discussion
Previous reconstructions for bumblebees have suggested that their ancestors diverged from other social bees in the mountains of Asia and dispersed throughout much of Asia and Europe; then via Beringia from Asia to North America; and further via the Panamanian isthmus from North to South America (Williams 1985, Kawakita et al. 2004, Hines 2008. Our results are consistent with this interpretation, although we seek to improve the resolution of the reconstructions and to compare them eventually with results for other high-elevation organisms (Favre et al. 2015, Renner 2016. Uncertainties in the dated phylogeny (Fig. 4) and in the distribution of Mendacibombus species are discussed in Williams et al. (2016). The nodes with lowest support in Fig. 4 are near the terminal branches of the tree and consequently have little effect on our reconstructions of ancestral distributions in Fig. 5. Uncertainties in the dating of mountain uplifts and in the changing climate are discussed below.

Integrating sequential filters for explaining distributions and gaps
For our first two filters, 1) genetic and 2) distribution evidence shows that the current distribution of Mendacibombus bumblebees could be explained by SDD through contiguous mountain ranges if suitable habitat were present (Fig. 5). The question is whether the habitat suitable for SDD was present at appropriate times or whether unsuitable habitat would have required LDD, as in some regions at present (especially for gaps a b). The window for dispersal through the gaps a and b (Fig. 1a) has its older bound set by the origin of the margreiteri-mendax group at event 4 (ca 5.3 Ma, Early Pliocene), because there is no evidence of Mendacibombus distributions in Europe from earlier periods (Fig. 5). Dispersal through gap g has its older bound set most likely by the divergence of the ancestor of B. margreiteri at event 8 (ca 3.1 Ma, Late Pliocene). This is because populations of B. margreiteri from Kamchatka and Mongolia are very similar even in the fast-evolving COI gene and yet divergent from the handlirschianus-mendax group (Williams et al. 2016). Therefore dispersal of B. margreiteri across gap g is unlikely to have pre-dated event 8 (Fig. 5). For dispersal across gap a, the younger bound is likely to have been set by the divergence between the ancestors of B. margreiteri and of the handlirschianus-mendax group also at event 8, because this is when these lineages diverged on either side of gap a. For gap b, the younger bound is likely to have been set by the divergence between the ancestors of B. handlirschianus and of B. mendax in the Early Pleistocene, because this is when these lineages diverged on either side of gap b. For gap g, there is no known recent constraint other than climate. Integration of the temporally broader evidence for our last two filters (3: geology, and 4: changes in palaeoclimate and in climate niches) is discussed below.

Distributions of early bumblebees
The divergence between the MRCA of Mendacibombus and the MRCA of all other bumblebees (Fig. 4) is likely to have coincided (Williams 1985, Hines 2008) with a rapid and pronounced cooling in global climate at the beginning of the Oligocene at ca 34 Ma (Zachos et al. 2008). An ability to thrive in cool climates may have facilitated the SDD of the earliest common ancestor of Mendacibombus at high elevations southwards into the early QTP region, probably initially into the older higher region of the Tanggula Shan (Yin et al. 2008, Yin 2010 (Fig. 6a; mountains labelled in Fig.  1b). The earliest-diverging extant species of the sister group of all other bumblebees (Cameron et al. 2007, their Fig. 1) still occur just to the north in the central Asian and Mongolian highlands (B. confusus and B. soroeensis : Skorikov 1931, Tkalcu 1974IAR collection). Consequently the two oldest sister-lineages of bumblebees (Cameron et al. 2007) have vicariant distributions in the mountains on either side (north and south) of the relatively low-lying arid Dzungar -Tarim (Taklamakan) -Qaidam basins (Fig. 1b, 6a). These basins now form a barrier of severe aridity (deserts) across central Asia (Yin et al. 2008) that is continuous with the Gobi Desert to the east. This arid zone developed especially with the Early Oligocene cooling of the climate and the westward retreat of the Tethys Sea (Dupont-Nivet et al. 2007) and remains the sharpest boundary between the large Palaearctic (northern) and Oriental (southern) bumblebee faunas in Asia today (Williams 1996). Available reconstructions of Late Eocene and Oligocene mountain ranges in these regions show stronger high-elevation connections with more potential for habitats suitable for SDD by bumblebees between the northern and southern mountains in the east (via the Great Khingan mountains: Fig. 1b) than in the west, where western highelevation connections were not yet available as the Tethys retreated (Scotese 2003, Markwick 2006, Blakey 2014.

The early (QTP) phase for Mendacibombus
The first phase of the history of the subgenus Mendacibombus is characterised by a long period (ca 34-8 Ma) of what appears to be relative stability in distribution (Fig. 4), centred on the QTP (Fig. 5). Part of the explanation for apparent stability might be that global temperatures became warmer than at present (Zachos et al. 2001(Zachos et al. , 2008, so that the early Mendacibombus may have become even more restricted to the higher mountains than now. The mountains of the Tian Shan, Altai, and some sections of the Alpide ranges might still have been relatively low and still unsuited as corridors by which Mendacibombus could have dispersed outwards from southern Asia by SDD. For the first ca 24 Myr of this phase (ca 34-9.4 Ma), only one lineage of Mendacibombus with present-day descendants is known (possibly a 'phylogenetic fuse': Cooper and Fortey 1998; or alternatively resulting from high rates of extinction: Condamine andHines 2015, Sanmartin andMeseguer 2016). The ancestral distributions reconstructed with S-DIVA and S-DEC (Fig. 5) place this lineage in the areas GHP (Table 1, Fig. 1b). However, the GP (and associated) mountains are older, with the more peripheral mountains to the north, south, and east (including the outer Himalaya and Hengduan, area H) uplifted later and progressively towards the east (Searle 1995, Ducea et al. 2003, Harris 2006, Wang et al. 2008, Yin 2010, so that GP may be a more realistic distribution for the earliest common ancestor of all Mendacibombus. By ca 9.4 Ma, global temperatures had declined again to near where they are now (Zachos et al. 2001(Zachos et al. , 2008. Although the central QTP was already arid before the Oligocene (DeCelles et al. 2007), our reconstruction of the QTP climate suggests that there was some habitat suitable for contemporary Mendacibombus (Fig. 7e). However, the distribution of the MRCA of Mendacibombus became divided by vicariance (event 1, Supplementary material Appendix 1 Table A3-A4), separating the superbus-convexus group with an ancestral distribution in the east QTP (area G, centred on the Tanggula: Fig. 5, 6b), from the makarjini-mendax group with an ancestral distribution to the west (area P) in the Pamir. The ancestor of B. superbus in central QTP (area G) then diverged following vicariance or SDD (event 2, Supplementary material Appendix 1 Table A3-A4) from the waltoni-convexus group to the east, which dispersed to the outer Himalaya and Hengduan mountains (area H: Fig. 5, 6c) as their uplift continued later.

The middle (dispersal) phase for Mendacibombus
In contrast, the second phase of Mendacibombus history (ca 8-3 Ma) is a period in which these bees dispersed much more broadly through other mountain ranges across Asia. This phase coincided with a long-term trend towards cooling in global temperatures (Zachos et al. 2008), which is likely to have expanded the area with suitable habitats, allowing Mendacibombus to show SDD to lower elevations. Simultaneously, many of the other Asian mountain ranges were being uplifted towards high elevations, bringing them within the climate zone suitable for Mendacibombus and increasing the connectivity between suitable habitats in previously isolated mountains: 1) uplift continued in the peripheral mountain ranges of the QTP (Wang et al. 2008, Yin 2010; 2) uplift continued of key high links from the QTP to (and within) the Alpine (or Alpide) mountains, and forming a continuous mountain belt from Afghanistan westwards into southern and central Europe (Schmid et al. 2004, Yin 2010, joining together some highland areas that had previously been separated by sea, e.g. in the Armenian highlands (Scotese 2003, Markwick 2006; and 3) further north, while uplift of the Tian Shan and Altai ranges may have begun or been reactivated in the Early Miocene (Yin 2010), high uplift accelerated in the Mid to Late Miocene (Yin et al. 2008, Yin 2010, Miao et al. 2012. In this middle phase there were three main divergence events within Mendacibombus (events 3-5: Supplementary material Appendix 1 Table A3-A4, Fig. 5) centred within the Pamir but followed by SDD into the uplifting contiguous mountain systems: 1) northwards into the Tian Shan (event 3: Fig. 6d); 2) westwards via the Hindu Kush and Alborz (event 4, Fig. 6e); and 3) south-eastwards into the Himalaya as well as north-eastwards into the Tian Shan (event 5: Fig. 6f ). With climate cooling, the ancestor of the margreiteri-mendax group may have been able to disperse by SDD westwards at slightly lower elevations of the Alpine or Alpide mountains: through the Hindu Kush, via the Aladagh across gap a (Fig. 1), and into the Alborz ranges (after event 4: Fig. 6e). Our reconstructions of ancestral climate niches and palaeoclimate simulations (Fig. 7b-g) support the idea that there could have been greater connectivity of potentially suitable habitat available for SDD by Mendacibombus across gap a in the cold Late Miocene (Fig. 5, 7e; although only if reconstructed ancestral climate niches are used, cf. Fig. 7d). Subsequently, gap a has become too warm and too arid at its elevation for Mendacibombus (Fig. 1a, 2, 7g). The postulated temporary existence of a cooler, less arid corridor across gap a is consistent with the small isolated and apparently relict pockets with bumblebees now present in some river valleys in the Aladagh mountains (in Khorasan). These bumblebees include species (e.g. B. fragrans (Pallas)) that are known today otherwise only as rare species from distant (by  900 km) regions of Turkey and central Asia to either side (Williams et al. 2011). LDD jumps over such long distances and across deserts to these small pockets in gap a seems unlikely for bumblebees (Introduction). There is also evidence for two later vicariance events between the Pamir and the Himalaya (events 6-7: Supplementary material Appendix 1 Table  A3-A4, Fig. 5).

The late (fragmentation) phase for Mendacibombus
The third phase of Mendacibombus history (ca 3 Ma to the present) is characterised by the break-up of populations within the margeiteri-mendax group. This phase coincided with a continuing global trend towards even cooler and increasingly fluctuating climates, including the Pleistocene glaciations (Zachos et al. 2001). The cooling trend is likely to have provided broader areas with suitable habitat at lower elevations, with more opportunities for species to show SDD to reach suitable areas that were previously isolated. In combination with continuing mountain uplift, cooling could have contributed towards increasingly broad Mendacibombus distributions in Fig. 6b-h.
Early Pleistocene cold periods could have provided corridors of suitable habitat at lower elevations through the European and Middle Eastern Alpide mountains (Fig. 1a,  gap b). There are no climate simulations available for this cold period from which to assess suitability quantitatively. However, mapping suitability for the much more recent Last Glacial Maximum by considering the climate conditions of current observations and following the same protocol (Supplementary material Appendix 1 Fig. A3) suggests that such Alpide corridors could have been likely in an earlier cold period. Warmer Pleistocene periods may then have contributed to the loss of Mendacibombus populations from the mountains of western Turkey, the Dinaric Alps, and much of southern France. These regions are currently only slightly too low at their latitudes (too warm) to have large areas of suitable habitat (Fig. 2, gap b).
For B. margreiteri, the modern Transbaikal and Kamchatka populations may once have been connected, extending through the older Stanovoy -Okhotsk mountains of northeastern Russia (Fig. 6h). The possibility of a single broad population of B. margreiteri throughout this region is supported by the lack of divergence in fast-evolving COI sequences between the present-day Mongolian and Kamchatkan samples (Williams et al. 2016). It may have taken only minor climate cooling to render habitat suitable for these bees to cross gap g by SDD (Fig. 1a, 2, 7a). However, subsequently this population has become fragmented into widely disjunct populations (by up to ca 2300 km: Fig. 1a, gap g), probably as a result of the more extreme climate fluctuations with warm (or very cold) periods in the Pleistocene. Around the time of the Last Glacial Maximum, much of this region is likely to have had similar steppe-tundra vegetation (Ray and Adams 2001). There were likely to have been periodic local increases in the Kamchatka ice cap (Ehlers and Gibbard 2007), with ice coverage reaching a maximum around 26.5-19.5 ka (Clark et al. 2009). The persistence of B. margreiteri as an isolated relict in south-eastern Kamchatka through cold periods may have been aided by locally ameliorating geothermal effects (Williams et al. 2016) during the last 50 ka (Braitseva et al. 1995). When sufficient fresh material of B. margreiteri is available, a more sensitive phylogeographic analysis (cf. Wang et al. 2011, Lecocq et al. 2013, Dellicour et al. 2016 should provide greater resolution of the history of this species. According to the more recent reconstructions (Wang et al. 2011: their Fig. 4;Lehmkuhl and Owen 2005), glaciations probably had a more limited effect on the drier QTP, even at the extreme elevation ( 5000 m) of the Tanggula Shan (Colgan et al. 2006). In contrast, after the end of the last glaciation, there is evidence of a wetter period with more extensive woodland between the Mongolian mountains and the Qilian mountains (Yu et al. 1998a, b). This may have been one of perhaps several periods that could have provided a bridge of suitable habitat in the intervening hills to allow B. margreiteri to show SDD southwards directly from the Altai or Tian Shan and into the Burhan Budai mountains (Fig. 6h). An eastern-bridge interpretation is supported by similar narrowly disjunct distributions in Qinghai that are known for several other bumblebee species, all of which are more widespread in Mongolia (Panfilov 1957, Peters and Panfilov 1968, Williams et al. 2015a Table 3) or in Neimenggu (Inner Mongolia: An et al. 2014). These species are more likely to have dispersed southwards into the Qinghai mountains at the same time as B. margreiteri, because in each case their sister species or closest relatives (Cameron et al. 2007) occur in the Palaearctic region (Williams 1998), to the north. In recent decades, the QTP has experienced a distinct warming trend, associated in its central and southwestern regions with aridification and further expansion of the Qiangtang desert (Zhang et al. 2013). This might be further reducing bumblebee distributions within the QTP (Williams et al. 2015a).

Conclusions
Our results show that the divergence of the present-day Mendacibombus species, rather than being associated in time with the early uplift of the QTP, is associated with the later Miocene and Pliocene uplift of mountain ranges connecting with the QTP and with dispersal out of the QTP. Biogeographic analyses show that the current distribution of Mendacibombus can be explained in part by species' SDD for range growth, although they identify gaps in the distribution with habitat that currently would be unsuitable and would require LDD. Nonetheless, we find corroboration for a more general role for SDD by Mendacibombus from evidence of temporary corridors of habitat in these gaps that would have been suitable in the appropriate periods in the past, resulting from changes in the climate and from changes in climate niches. Mendacibombus species appear subsequently to have been extirpated from these gaps, reducing their range sizes. Our results are consistent with those from previous studies of bumblebees, which have suggested in broad terms that Asian mountains are likely to have played an important role in early bumblebee evolution (Skorikov 1931, Williams 1985, Hines 2008, Williams et al. 2011. Formal diversification analysis applied across all bumblebees associated the only increase in net diversification rate among all bumblebees with climate cooling during the Pliocene (Condamine and Hines 2015), which is paralleled by our finding of increased divergence in this period. The strongest increase in extinction rate was associated with Pleistocene climate oscillations, which is also paralleled by our finding of the timing of extirpations in the gaps b and g. We seek to provide greater spatial and temporal resolution for this one group of bumblebees, but detailed comparisons of shared patterns with other more species-rich mountain-bumblebee groups and with other groups of organisms would be informative as more data become available. Changing geological and climate conditions on Asian mountains have been claimed to impact the evolution 475 of many other groups of mountain biota, although care is needed in assessing the evidence (Favre et al. 2015, Renner 2016. We would expect that the patterns shown by Mendacibombus are most likely to be shared by other cool-adapted groups from similar alpine and subalpine QTP habitats.