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

Paul H. Williams, Jorge M. Lobo, Andrea Sanchez Meseguer

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Paul H. Williams, Jorge M. Lobo, Andrea Sanchez Meseguer. Bumblebees take the high road: climatically integrative biogeography shows that escape from Tibet, not Tibetan uplift, is associated with divergences of present-day Mendacibombus. Ecography, 2018, 41 (3), pp.461-477. 10.1111/ecog.03074. hal-02623734

HAL Id: hal-02623734

https://hal.inrae.fr/hal-02623734

Submitted on 28 Jun 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
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*

Paul H. Williams, Jorge M. Lobo and Andrea S. Meseguer

P. H. Williams (http://orcid.org/0000-0002-6996-5682) (paw@nhm.ac.uk), Dept of Life Sciences, Natural History Museum, London, UK. - J. M. Lobo, Dept of Biogeography and Global Change, Museo Nacional de Ciencias Naturales, Madrid, Spain. - A. S. Meseguer, INRA-UMR 1062 CBGP (INRA, IRD, CIRAD, Montpellier SupAgro), Montferrier-sur-Lez, France, and CNRS, UMR 5554 Inst. des Sciences de l'Evolution (ISEM, Univ. de Montpellier), Montpellier, France.

Ecography

41: 461–477, 2018

doi: 10.1111/ecog.03074

Subject Editor: Corrie Moreau

Editor-in-Chief: Miguel Araújo

Accepted 1 March 2017

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, comensurate 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,

© 2017 Natural History Museum, London, UK. Ecography © 2017 Nordic Society Oikos

<figure><img src="image_3.png" /></figure>
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, 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 sister-group 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, 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 and Corbet 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 and
<figure><img src="image_4.png" /><figcaption>Figure 1. (a) Distributions of all 12 extant *Mendacibombus* species (separated as differently coloured spots: key shown to the left) from the 4416 specimen records updated from Williams et al. (2016) with additional Mongolian records (A. Dorjoo in litt.). The three largest current disjunctions in distribution ('gaps') are indicated by grey lines and by the letters α (Aladagh mountains), β (Balkan-Dinaric mountains), and γ (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: α, 1520 km; β, 1450 km; γ, 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.</figcaption></figure>

Schmid-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
<figure><img src="image_5.png" /><figcaption>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 (α β γ) 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.</figcaption></figure>

inbreeding depression, making long-term establishment of an isolated new population less likely.

Assessing the potential for LDD by bumblebees empiri-
cally, 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 pro-
duced 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 stron-
ger 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 sig-
nificant that *B. hypnorum* is very unusual among bumblebees
for being polyandrous (Estoup et al. 1995, Schmid-Hempel

and Schmid-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 dis-
persal 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 reason-
able 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 rep-
resent some of the many factors associated in time with the
distribution and divergence of present-day species, includ-
ing the uplift of the QTP and of the associated mountains.
This requires that we assess whether corridors of habitat suit-
able 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 cur-
rent large 'gaps' (α β γ) in their distribution in Fig. 1a. To
integrate the major governing factors, we use a framework
with sequential filters – to constrain progressively the expla-
nations 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) evi-
dence 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 suit-
able habitat was present in these corridors within the appro-
priate temporal windows.

**Methods**

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): *B. (Bombias) confusus* Schenck, *B. (Kallobombus) soroeensis* (Fabricius), *B. (Subterraneobombus) difficillimus* Skorikov, *B. (Megabombus) gerstaeckeri* Morawitz, *B. (Alpinobombus) balteatus* Dahlbom, *B. (Pyrobombus) bimaculatus* Cresson, and *B. (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 nucleotide-substitution 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 + Γ); 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 + Γ; 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 uncorre-
lated 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 gen-
erations, 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 com-
mon 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-extinction-
cladogenesis 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 speci-
men records for *Mendacibombus* species (Fig. 1a) updated
<table><caption>Table 1. Principal areas of endemic distribution identified for their combinations of bumblebee species of the subgenus *Mendacibombus*.</caption><thead><tr><th>Area</th><th>Principal mountain ranges included</th><th>Species recorded</th></tr></thead><tbody><tr><td>G</td><td>Tanggula - Hohxil Shan</td><td>superbus waltoni</td></tr><tr><td>H</td><td>Himalaya - Hengduan - Qilian Shan</td><td>avinoviellus convexus himalayanus waltoni</td></tr><tr><td>P</td><td>Pamir - Hindu Kush</td><td>makarjini marussinus turkestanicus</td></tr><tr><td>T</td><td>Alai - Tian Shan</td><td>defector makarjini margreiteri turkestanicus</td></tr><tr><td>M</td><td>Altai - Khangai (C Mongolia) - Kamchatka</td><td>margreiteri</td></tr><tr><td>C</td><td>Alborz - Armenia - Caucasus</td><td>handlirschianus</td></tr><tr><td>A</td><td>Alps - Pyrenees</td><td>mendax</td></tr></tbody></table>

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.

from Williams et al. (2016) into a few large areas (Table 1, 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 and Redelings 2009, Ree and Sanmartin 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

<figure><img src="image_6.png" /><figcaption>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).</figcaption></figure>
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₂ concentration (Beerling and Royer 2011). The simulations were generated at a 2.5° × 3.75° grid resolution, corresponding to: 1) a 400 ppm CO₂ 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₂ Miocene simulation representing the cold and dry conditions prevailing around 8 Ma; and 3) a 560 ppm CO₂ Pliocene simulation of the conditions around the Mid Pliocene warming event (3.6–2.6 Ma) (Beerling et al. 2009, 2011, 2012 and 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 (δ) 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.

## Results

### 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
<figure><img src="image_7.png" /><figcaption>Figure 4. Dated estimate of phylogeny for all extant species of the subgenus *Mendacibombus* by linked Bayesian (BEAST) analysis of mitochondrial (COI, 16S) and nuclear (PEPCK, opsin) genes from two billion generations with a 1% burn-in (Williams et al. 2016). Values above the nodes are Bayesian posterior probabilities for groups. Values below the nodes are estimates of the dates of divergence in Ma (millions of years before the present) based on the *Mendacibombus* root-date estimate of ca 34 Ma (shown to the left). The grey bars represent the 95% confidence limits for date estimates.</figcaption></figure>

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 γ from Fig. 1a), but much less favourable areas in the gaps α and β
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 5. Most likely biogeographic reconstruction for species of the subgenus *Mendacibombus* from (a) S-DIVA analysis and (b) S-DEC analysis using the short-distance dispersal (no jumps) model in Fig. 3 and the estimate of phylogeny in Fig. 4 with other bumblebees as the outgroup. Letters represent the areas listed in 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).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 6. Maps of the dispersal and speciation history for species of the subgenus *Mendacibombus* from the most likely reconstructed ancestral distributions in the S-DIVA results in Fig. 5a (Supplementary material Appendix 1 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). Double-headed 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.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 7. Reconstruction of climate suitability of the Palaearctic since the Mid Miocene based on the current observations and estimated ancestral climate niche of *Mendacibombus* through time. (a) Suitability of the current climate for all current *Mendacibombus* species shown as a grey-scale (scale bar right) on a 2.5 × 2.5 minute grid for Mahalanobis distances from the centroid of all species' current records to current climates. Maps of suitability of palaeoclimates shown as a grey-scale on a 2.5° × 3.75° grid for Mahalanobis distances from the centroid of all *Mendacibombus* species' current records to: (b) the Mid Miocene (warm) palaeoclimate simulation; (c) Late Miocene (cold) palaeoclimate simulation; (d) Pliocene (warm) palaeoclimate simulation; and for (e) the ancestral climate niche for the most recent common ancestor (MRCA) of all *Mendacibombus* species projected to the Mid Miocene (warm) palaeoclimate simulation; (g) the ancestral climate niche for the MRCA of the *makarjini-mendax* group projected to the Late Miocene (cold) palaeoclimate simulation; and (h) the ancestral climate niche for the MRCA of the *margreiteri-mendax* group projected to the (warm) Pliocene palaeoclimate simulation. Only Mahalanobis distances below the upper quartile for cells with occurrence data are shown. All lowland areas (elevation < 500 m) and areas disjunct (by > 1000 km) from extant observations are excluded. Image created in IDRISI (ver. 17.02, Clark Labs).</figcaption></figure>

(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 $\alpha$ (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 α (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 α β). The window for dispersal through the gaps α and β (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 γ 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 γ is unlikely

to have pre-dated event 8 (Fig. 5). For dispersal across gap α, 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 α. For gap β, 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 β. For gap γ, 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 1974; IAR 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 high-elevation 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, 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 extinc-
tion: Condamine and Hines 2015, Sanmartin and Meseguer
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 *Men-
dacibombus*. By ca 9.4 Ma, global temperatures had declined
again to near where they are now (Zachos et al. 2001, 2008).
Although the central QTP was already arid before the Oligo-
cene (DeCelles et al. 2007), our reconstruction of the QTP
climate suggests that there was some habitat suitable for
contemporary *Mendacibombus* (Fig. 7e). However, the dis-
tribution 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. Simulta-
neously, 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 continu-
ous 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: Supplemen-
tary 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 α (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 α in the cold Late Miocene (Fig. 5, 7e; although only if
reconstructed ancestral climate niches are used, cf. Fig. 7d).
Subsequently, gap α 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
α 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 α 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 β). 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 β).

For *B. margreiteri*, the modern Transbaikal and Kamchatka populations may once have been connected, extending through the older Stanovoy – Okhotsk mountains of north-eastern 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 γ 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 γ), 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: their 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 β and γ. 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
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 *Mendaci-
bombus* are most likely to be shared by other cool-adapted
groups from similar alpine and subalpine QTP habitats.

*Acknowledgements* – We thank all who have collected specimens and maintained collections for study. We thank the Biodiversity Inst. of Ontario for DNA sequencing. We thank L. Bailey, A. Purvis, S. Renner, A. Favre, H. Schneider and Y. for discussion and comments.

References

Acevedo, P. et al. 2012. Delimiting the geographical background in species distribution modelling. - J. Biogeogr. 39: 1383-1390.

Alford, D. V. 1975. Bumblebees. - Davis-Poynter.

An, J.-D. et al. 2014. The bumblebees of north China (Apidae, *Bombus* Latreille). - Zootaxa 3830: 1-89.

Beerling, D. J. and Royer, D. L. 2011. Convergent CO2 history. - Natl Geosci. 4: 418-420.

Beerling, D. J. et al. 2009. Methane and the CH4 related greenhouse effect over the last 400 million years. - Am. J. Sci. 309: 97-113.

Beerling, D. J. et al. 2011. Enhanced chemistry-climate feedbacks in past greenhouse worlds. - Proc. Natl Acad. Sci. USA 108: 9770-9775.

Beerling, D. J. et al. 2012. Ecosystem CO2 starvation and terrestrial silicate weathering: mechanisms and global-scale quantification during the late Miocene. - J. Ecol. 100: 31-41.

Blakey, R. 2014. Paleogeography library/rectangular global maps/ Oligocene (35 Ma). - <http://cpgeosystems.com/paleomaps.html>.

Bradshaw, C. D. et al. 2012. The relative roles of CO2 and palaeogeography in determining late Miocene climate: results from a terrestrial model-data comparison. - Clim. Past 46: 171-175.

Braitseva, O. A. et al. 1995. Ages of calderas, large explosive craters and active volcanoes in the Kuril-Kamchatka region, Russia. - Bull. Volcanol. 57: 383-402.

Buttermore, R. E. et al. 1998. Assessment of the genetic base of Tasmanian bumble bees (*Bombus terrestris*) for development as pollination agents. - J. Apicult. Res. 37: 23-25.

Cameron, S. A. et al. 2007. A comprehensive phylogeny of the bumble bees (*Bombus*). - Biol. J. Linn. Soc. 91: 161-188.

Cameron, S. A. et al. 2011. Patterns of widespread decline in North American bumble bees. - Proc. Natl Acad. Sci. USA 108: 662-667.

Clark, P. U. et al. 2009. The last glacial maximum. - Science 325: 710-714.

Cogbill, C. V. and White, P. S. 1991. The latitude-elevation relationship for spruce-fir forest and treeline along the Appalachian mountain chain. - Vegetatio 94: 153-175.

Colgan, P. M. et al. 2006. Cosmogenic radionuclide evidence for the limited extent of last glacial maximum glaciers in the Tanggula Shan of the central Tibetan Plateau. - Quat. Res. 65: 336-339.

Condamine, F. L. and Hines, H. M. 2015. Historical species losses in bumblebee evolution. - Biol. Lett. 11: e20141049.

Cooper, A. and Fortey, R. 1998. Evolutionary explosions and the phylogenetic fuse. - Trends Ecol. Evol. 13: 151-155.

Crook, J. M. and Crozier, R. H. 1995. Sex determination and population biology in the Hymenoptera. - Trends Ecol. Evol. 10: 281-286.

Crowther, L. P. et al. 2014. Habitat and forage associations of a naturally colonising insect pollinator, the tree bumblebee *Bombus hypnorum*. - PLoS One 9: e107568.

Darvill, B. et al. 2006. Population structure and inbreeding in a rare and declining bumblebee, *Bombus muscorum* (Hymenoptera: Apidae). - Mol. Ecol. 15: 601-611.

DeCelles, P. G. et al. 2007. Late Cretaceous to middle Tertiary basin evolution in the central Tibetan Plateau: changing environments in response to tectonic partitioning, aridification, and regional elevation gain. - Geol. Soc. Am. Bull. 119: 654-680.

Dellicour, S. et al. 2016. Ecological niche modelling and coalescent simulations to explore the recent geographical range history of five widespread bumblebee species in Europe. - J. Biogeogr. doi: 10.1111/jbi.12748

Drummond, A. J. and Rambaut, A. 2007. BEAST: Bayesian evolutionary analysis by sampling trees. - BioMed Central Evol. Biol. 7: 214.

Ducea, M. N. et al. 2003. Building the Pamirs: the view from the underside. - Geology 31: 849-852.

Duchateau, M. J. et al. 1994. Diploid males in the bumble bee *Bombus terrestris*. - Entomol. Exp. Appl. 71: 263-269.

Dupont-Nivet, G. et al. 2007. Tibetan plateau aridification linked to global cooling at the Eocene–Oligocene transition. – Nature 445: 635–638.

Ehlers, J. and Gibbard, P. L. 2007. The extent and chronology of Cenozoic global glaciation. - Quat. Int. 164-165: 6-20.

Estoup, A. et al. 1995. Monandry and polyandry in bumble bees (Hymenoptera, Bombinae) as evidenced by highly variable microsatellites. - Mol. Ecol. 4: 89-93.

Estoup, A. et al. 1996. Genetic differentiation of continental and island populations of *Bombus terrestris* (Hymenoptera: Apidae) in Europe. - Mol. Ecol. 5: 19-31.

Favre, A. et al. 2015. The role of the uplift of the Qinghai-Tibetan Plateau for the evolution of Tibetan biotas. - Biol. Rev. 90: 236-253.

García-Roselló, E. et al. 2013. ModestR: a software tool for managing and analyzing species distribution map databases. - Ecography 36: 1202-1207.

Goulson, D. 2010. Impacts of non-native bumblebees in western Europe and North America. - Appl. Entomol. Zool. 45: 7-12.

Goulson, D. and Williams, P. 2001. *Bombus hypnorum* (Hymenoptera: Apidae), a new British bumblebee? - Br. J. Entomol. Nat. Hist. 14: 129-131.

Goulson, D. et al. 2011. Population structure, dispersal and colonization history of the garden bumblebee *Bombus hortorum* in the Western Isles of Scotland. - Conserv. Genet. 12: 867-879.

Haeseler, V. 1974. Aculeate Hymenopteren über Nord- und Ostsee nach Untersuchungen auf Feuerschiffen. - Insect Syst. Evol. 5: 123-136.

Harris, N. 2006. The elevational history of the Tibetan Plateau and its implications for the Asian monsoon. - Palaeogeogr. Palaeoclimatol. Palaeoecol. 241: 4-15.

Heinrich, B. 1979. Bumblebee economics. - Harvard Univ. Press.

Heled, J. and Drummond, A. J. 2010. Bayesian inference of species trees from multilocus data. - Mol. Biol. Evol. 27: 570-580.

Hijmans, R. J. et al. 2005. Very high resolution interpolated climate surfaces for global land areas. - Int. J. Climatol. 25: 1965-1978.
Hines, H. M. 2008. Historical biogeography, divergence times, and diversification patterns of bumble bees (Hymenoptera: Apidae: *Bombus*). – *Syst. Biol.* 57: 58–75.

Ito, M. 1987. Geographic variation of an eastern Asian bumble bee *Bombus diversus* in some morphometric characters (Hymenoptera, Apidae) – Kontyu 55: 188–201.

Ito, M. and Sakagami, S. F. 1980. The bumblebee fauna of the Kurile Islands (Hymenoptera: Apidae). – *Low Temp. Sci. B* 38: 23–51.

Johnson, D. H. 1980. The comparison of usage and availability measurements for evaluating resource preference. – *Ecology* 61: 65–71.

Kawakita, A. et al. 2004. Phylogeny, historical biogeography, and character evolution in bumble bees (*Bombus*: Apidae) based on simultaneous analysis of three nuclear gene sequences. – *Mol. Phylogenet. Evol.* 31: 799–804.

Kraus, F. B. et al. 2009. Male flight distance and population substructure in the bumblebee *Bombus terrestris*. – *J. Anim. Ecol.* 78: 247–252.

Lamm, K. S. and Redelings, B. D. 2009. Reconstructing ancestral ranges in historical biogeography: properties and prospects. – *J. Syst. Evol.* 47: 369–382.

Lecocq, T. et al. 2013. Scent of a break-up: phylogeography and reproductive trait divergences in the red-tailed bumblebee (*Bombus lapidarius*). – *BMC Evol. Biol.* 13: 263.

Lecocq, T. et al. 2016. Conservation genetics of European bees: new insights from the continental scale. – *Conserv. Genet.* doi: 10.1007/s10592-016-0917-3

Lehmkuhl, F. and Owen, L. A. 2005. Late Quaternary glaciation of Tibet and the bordering mountains: a review. – *Boreas* 34: 87–100.

Lepais, O. et al. 2010. Estimation of bumblebee queen dispersal distances using sibship reconstruction method. – *Mol. Ecol.* 19: 819–831.

Løken, A. 1973. Studies on Scandinavian bumble bees (Hymenoptera, Apidae). – *Norsk Entomol. Tidsskr.* 20: 1–218.

Macdonald, M. 2001. The colonisation of northern Scotland by *Bombus terrestris* (L.) and *B. lapidarius* (L.) (Hym., Apidae), with comments on the possible role of climate change. – *Entomol. Monthly Mag.* 137: 1–13.

Macfarlane, R. P. and Gurr, L. 1995. Distribution of bumble bees in New Zealand. – *N. Z. Entomol.* 18: 29–36.

Markwick, P. 2006. Paul's palaeo pages/Early Oligocene 33 Ma, Palaeogene (Tertiary). – <www.palaeogeography.net/palaeogeography/...maps/pg_oligocene1/pjmg1oli1.html>.

McCluskey, A. 2012. Let sleeping queens lie. – *Buzzword Nov.* 2012: 5.

Miao, Y.-F. et al. 2012. What controlled mid-late Miocene long-term aridification in central Asia? Global cooling or Tibetan plateau uplift: a review. – *Earth-Sci. Rev.* 112: 155–172.

Mikkola, K. 1984. Migration of wasp and bumblebee queens across the Gulf of Finland (Hymenoptera: Vespidae and Apidae). – *Notulae Entomol.* 64: 125–128.

Morlon, H. et al. 2011. Reconciling molecular phylogenies with the fossil record. – *Proc. Natl Acad. Sci. USA* 108: 16327–16332.

Neter, J. et al. 1990. Applied linear statistical models. – *Irwin*.

Panfilov, D. V. 1957. On the geographical distribution of bumblebees (*Bombus*) in China. – *Acta Geogr. Sinica* 23: 221–239.

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

Pennell, M. W. et al. 2014. Geiger v2.0: an expanded suite of methods for fitting macroevolutionary models to phylogenetic trees. – *Bioinformatics* 30: 2216–2218.

Peters, G. and Panfilov, D. V. 1968. Hummeln (*Bombus*) und Schmarotzerhummeln (*Psithyrus*) Ergebnesse der zoologischen Forschungen von Dr. Z. Kaszab in der Mongolei. – *Reisbachia* 11: 117–182.

Prys-Jones, O. E. and Corbet, S. A. 1987. Bumblebees. – *Cambridge Univ. Press*.

Rasmont, P. 1988. Monographie écologique et zoogéographique des bourdons de France et de Belgique (Hymenoptera, Apidae, Bombinae). – *Faculté des Sciences agronomiques de l'Etat, p. 309 + Ixi*.

Rasmont, P. and Flagothier, D. 1996. Biogéographie et choix floraux des bourdons (Hymenoptera, Apidae) de la Turquie. Rapport preliminaire 1995-1996. NATO-OTAN TU-Pollination Project. – *Univ. de Mons-Hainaut, Mons, Belgium*.

Ray, N. and Adams, J. P. 2001. A GIS-based vegetation map of the world at the last glacial maximum (25,000–15,000 bp). – *Internet Archaeol.* 11: 1–44.

Ree, R. H. and Smith, S. A. 2008. Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis. – *Syst. Biol.* 57: 4–14.

Ree, R. H. and Sanmartin, I. 2009. Prospects and challenges for parametric models in historical biogeographical inference. – *J. Biogeogr.* 36: 1211–1220.

Renner, S. S. 2016. Available data point to a 4-km-high Tibetan Plateau by 40 Ma, but 100 molecular-clock papers have linked supposed recent uplift to young node ages. – *J. Biogeogr.* 43: 1479–1487.

Revell, L. J. 2012. Phytools: an R package for phylogenetic comparative biology (and other things). – *Methods Ecol. Evol.* 3: 217–223.

Ronquist, F. 1996. DIVA version 1.1 user's manual. – <www.ib.usp.br/~silvionihei/pdf/Ronquist%201996_DIVA%20user%20manual.pdf>.

Sakagami, S. F. 1976. Specific differences in the bionomic characters of bumblebees. A comparative review. – *J. Faculty Sci., Hokkaido Univ. (Zoology)* 20: 390–447.

Sanmartin, I. and Meseguer, A. S. 2016. Extinction in phylogenetics and biogeography: from timetrees to patterns of biotic assemblage. – *Front. Genet.* 7: 1–17.

Schmid, S. M. et al. 2004. Tectonic map and overall architecture of the Alpine orogen. – *Eclogae Geol. Helv.* 97: 93–117.

Schmid-Hempel, P. et al. 2007. Invasion success of the bumblebee, *Bombus terrestris*, despite a drastic genetic bottleneck. – *Heredity* 99: 414–422.

Schmid-Hempel, R. and Schmid-Hempel, P. 2000. Female mating frequencies in *Bombus* spp. from central Europe. – *Insectes Soc.* 47: 36–41.

Schmid-Hempel, R. et al. 2013. The invasion of southern South America by imported bumblebees and associated parasites. – *J. Anim. Ecol.* 2013: 1–15.

Scotese, C. R. 2003. Paleomap Project/Middle Miocene 14 Ma. – <www.scotese.com/>.

Searle, M. P. 1995. The rise and fall of Tibet. – *Nature* 374: 17–18.

Skorikov, A. S. 1923. Palaeartic bumblebees. Part I. General biology (including zoogeography)]. – *Izvestiya Severnoi oblastnoi stantsii zashchity rastenii ot vreditelei* 4 (1922): 1–160.

Skorikov, A. S. 1931. Die Hummelfauna Turkestans und ihre Beziehungen zur zentralasiatischen Fauna (Hymenoptera,
Bombidae). – In: Lindholm, V. A. (ed.), Abhandlungen der Pamir-Expedition 1928. Academy of Sciences of the USSR, pp. 175–247.

Smissen van der, J. and Rasmont, P. 1999. *Bombus semenoviellus* Skorikov 1910, eine für Westeuropa neue Hummelart (Hymenoptera: *Bombus*, *Cullumanobombus*). – Bembix 13: 21–24.

Tkalcu, B. 1974. Ergebnisse der 1. und 2. mongolisch-tschechoslowakischen entomologisch-botanischen Expedition in der Mongolei. Nr. 29: Hymenoptera, Apoidea, Bombinae. – Sborník faunistickych prací Entomologického oddelení Národního musea v Praze 15: 25–57.

Valencia-Barrera, R. et al. 2002. Bioclimatic indices as a tool in pollen forecasting. – Int. J. Biometeorol. 46: 171–175.

Wang, C.-S. et al. 2008. Constraints on the early uplift history of the Tibetan plateau. – Proc. Natl Acad. Sci. USA 105: 4987–4992.

Wang, L. et al. 2011. Phylogeography of the Sino-Hyalayan fern *Lepisorus clathratus* on "the roof of the world". – PLoS One 6: e25896.

Widmer, A. and Schmid-Hempel, P. 1999. The population genetic structure of a large temperate pollinator species, *Bombus pascuorum* (Scopoli) (Hymenoptera: Apidae). – Mol. Ecol. 8: 387–398.

Williams, P. H. 1985. A preliminary cladistic investigation of relationships among the bumble bees (Hymenoptera, Apidae). – Syst. Entomol. 10: 239–255.

Williams, P. H. 1991. The bumble bees of the Kashmir Himalaya (Hymenoptera: Apidae, Bombini). – Bull. Br. Mus. 60: 1–204.

Williams, P. H. 1996. Mapping variations in the strength and breadth of biogeographic transition zones using species turnover. – Proc. R. Soc. B 263: 579–588.

Williams, P. H. 1998. An annotated checklist of bumble bees with an analysis of patterns of description (Hymenoptera: Apidae, Bombini). – Bull. Nat. Hist. Mus. 67: 79–152, updated at <www.nhm.ac.uk/bombus/> accessed 2017.

Williams, P. H. et al. 2008. A simplified subgeneric classification of the bumblebees (genus *Bombus*). – *Apidologie* 39: 46–74.

Williams, P. H. et al. 2009. The bumblebees of Sichuan (Hymenoptera: Apidae, Bombini). – Syst. Biodivers. 7: 101–190.

Williams, P. H. et al. 2011. The bumblebees of the subgenus *Subterraneobombus*: integrating evidence from morphology and

DNA barcodes (Hymenoptera, Apidae, *Bombus*). – Zool. J. Linn. Soc. 163: 813–862.

Williams, P. H. et al. 2014. Bumble bees of North America. An identification guide. – Princeton Univ. Press.

Williams, P. H. et al. 2015a. Bumblebees, climate and glaciers across the Tibetan plateau (Apidae: *Bombus* Latreille). – Syst. Biodivers. 13: 164–181.

Williams, P. H. et al. 2015b. Genes suggest ancestral colour polymorphisms are shared across morphologically cryptic species in arctic bumblebees. – PLoS One 10: e0144544.

Williams, P. H. et al. 2016. Early-diverging bumblebees from across the roof of the world: the high-mountain subgenus *Mendacibombus* revised from species' gene coalescents and morphology (Hymenoptera, Apidae). – Zootaxa 4204: 1–72.

Yin, A. 2010. Cenozoic tectonic evolution of Asia: a preliminary synthesis. – *Tectonophysics* 488: 293–325.

Yin, A. et al. 2008. Cenozoic tectonic evolution of the Qaidam basin and its surrounding regions (part 3): structural geology, sedimentation, and regional tectonic reconstruction. – Geol. Soc. Am. Bull. 120: 847–876.

Yu, G. et al. 1998a. Pollen-based biome reconstructions for China at 0 and 6000 years. – J. Biogeogr. 25: 1055–1069.

Yu, G. et al. 1998b. Pollen-based reconstruction of vegetation patterns of China in mid-Holocene. – Science in China (D) 41: 130–136.

Yu, Y. et al. 2010. S-DIVA (statistical dispersal-vicariance analysis): a tool for inferring biogeographic histories. – Mol. Phylogenet. Evol. 56: 848–850.

Yu, Y. et al. 2014. RASP (reconstruct ancestral state in phylogenies) 3.0. – <http://mnh.scu.edu.cn/soft/blog/RASP>.

Zachos, J. C. et al. 2001. Trends, rhythms, and aberrations in global climate 65 Ma to Present. – *Science* 292: 686–693.

Zachos, J. C. et al. 2008. An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics. – *Nature* 451: 279–283.

Zhang, L. et al. 2013. Vegetation greenness trend (2000 to 2009) and the climate controls in the Qinghai-Tibetan Plateau. – J. Appl. Remote Sens. 7, doi: 10.1117/1.JRS.7.073572

Supplementary material (Appendix ECOG-03074 at <www.ecography.org/appendix/ecog-03074>). Appendix 1.