doi: 10.1111/jeb.12477

Early diversification trend and Asian origin for extent bat lineages W. YU*†, Y. WU* & G. YANG† *College of Life Sciences, Guangzhou University, Guangzhou, China †Jiangsu Key Laboratory for Biodiversity and Biotechnology, College of Life Sciences, Nanjing Normal University, Nanjing, China

Keywords:

Abstract

ancient biogeographic history; Asian origin; early diversification trend; extant bat lineage; palaeotemperature.

Bats are a unique mammalian group, which belong to one of the largest and most diverse mammalian radiations, but their early diversification is still poorly understood, and conflicting hypotheses have emerged regarding their biogeographic history. Understanding their diversification is crucial for untangling the enigmatic evolutionary history of bats. In this study, we elucidated the rate of diversification and the biogeographic history of extant bat lineages using genus-level chronograms. The results suggest that a rapid adaptive radiation persisted from the emergence of crown bats until the Early Eocene Climatic Optimum, whereas there was a major deceleration in diversification around 35–49 Ma. There was a positive association between changes in the palaeotemperature and the net diversification rate until 35 Ma, which suggests that the palaeotemperature may have played an important role in the regulation of ecological opportunities. By contrast, there were unexpectedly higher diversification rates around 25–35 Ma during a period characterized by intense and long-lasting global cooling, which implies that intrinsic innovations or adaptations may have released some lineages from the intense selective pressures associated with these severe conditions. Our reconstruction of the ancestral distribution suggests an Asian origin for bats, thereby indicating that the current panglobal but disjunct distribution pattern of extant bats may be related to events involving seriate cross-continental dispersal and local extinction, as well as the influence of geological events and the expansion and contraction of megathermal rainforests during the Tertiary.

Introduction Bats are a unique group of mammals, which are capable of self-powered flight and sophisticated echolocation. Bats comprise 205 extant genera and approximately 1100 species, and they represent one of the largest and most ecologically and morphologically diverse mammalian radiations (Jones et al., 2005; Simmons, 2005). With the exception of Antarctica, bats are found in almost all of the terrestrial habitats throughout the world (Simmons, 2005), where they play important roles in ecosystem functioning (Simmons & Conway, 2003). Large-scale molecular studies have resulted in major upheavals in chiropteran systematics Correspondence: Guang Yang, Jiangsu Key Laboratory for Biodiversity and Biotechnology, College of Life Sciences, Nanjing Normal University, Nanjing 210046, China. Tel./fax: +86 25 85891163; e-mail: [email protected]

in recent years (Eick et al., 2005; Teeling et al., 2005; Miller-Butterworth et al., 2007; Meredith et al., 2011), but the sparse fossil records (Teeling et al., 2005; Eiting & Gunnell, 2009) and lack of comprehensive phylogenetic studies [but see (Eick et al., 2005; Jones et al., 2005; Teeling et al., 2005)] have resulted in an evolutionary history that is poorly understood with respect to their early diversification dynamics and biogeographic history. Untangling chiropteran macroevolution has long been regarded as a thorny subject and Van Valen (1979, p. 103) stated: ‘One may hypothesize that bats did originate, but it is harder to go beyond this’. Identifying the evolutionary diversification rate of organisms is a crucial step in untangling their evolutionary history, where the ultimate goal is identifying the factors that might potentially affect the tempo of diversification (increases or decreases in the net diversification rate). Based on a phylogenetic supertree of bats, Jones et al. (2005) identified extreme heterogeneities of

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

1

2

W. YU ET AL.

diversification among lineages where the largest diversification rate shifts occurred at 30–50 Ma. Meredith et al. (2011) modeled the diversification trends of extant mammals using a phylogeny that incorporated almost all of the family level clades. However, neither of these studies explicitly modeled the diversification dynamics of bats over time. Surprisingly, the most notable key innovation of bats, that is flight (Simmons et al., 2008), does not appear to have triggered their diversification (Yu et al., 2012). Instead, a considerable level of heterogeneity was present within the crown bats group. The lack of an overall increase in the net diversification rate also contradicts the idea that flight triggered the transition into novel habitats during adaptive radiation (Simpson, 1949, 1953; Schluter, 2000; Yoder et al., 2010), which suggests that more nuanced drivers affected bat diversification. For example, the profound effects of temperature on ecosystem functioning and biodiversity (Mayhew et al., 2008), as well as the possibility of a close link between temperature and the energetic costs for bats (Stevens, 2004, 2006; Buckley et al., 2010), suggest that alternative drivers affected chiropteran diversification. Another unresolved and controversial issue in the macroevolution of bats concerns their biogeographic history (Springer et al., 2011). This is because the biogeographic history of bats cannot be inferred from the fossil record alone due to the fragmentary nature of the palaeontological records (Teeling et al., 2005; Eiting & Gunnell, 2009). The widespread oldest bats from the Early Eocene (Onychonycteris finneyi, about 52.5 Ma) (Smith et al., 2007; Simmons et al., 2008) as well as the large gap between the oldest record and the age estimated for crown bats based on molecular dating (around 58–71 Ma) (Eick et al., 2005; Jones et al., 2005; Miller et al., 2005; Teeling et al., 2005; Teeling, 2009; Meredith et al., 2011) have also made it difficult to determine their geographic origins (Gunnell & Simmons, 2005; Springer et al., 2011). Inferring the ancient biogeographic pattern and origin via phylogenetic reconstruction may be the only method available at present, but previous studies have obtained inconsistent results using this approach (Eick et al., 2005; Teeling et al., 2005; Springer et al., 2011). For example, based on a comprehensive molecular scaffold, Teeling et al. (2005) suggested Laurasia or Northern America as the ancestral area of bats using a parsimony method. Using dispersal–vicariance analysis, Eick et al. (2005) proposed an African origin for extant bats; however, Springer et al. (2011) inferred inclusive areas for crown bats [these reconstructions included at least five areas (Africa, Asia, Europe, North America and South America) and four areas (Africa, Asia, North America and South America)] using identical data and settings. In particular, all previous studies have inferred cross-Atlantic dispersal events between Africa and South America during chiropteran biogeographic his-

tory (i.e. emballonurids and noctilionoids) (Eick et al., 2005; Teeling et al., 2005). However, the parsimony and event-based methods utilized above have been criticized because they have various flaws, the most important of which is their intrinsic inability to incorporate past geographic and climatic information, thereby allowing the possible oversimplification and misrepresentation of the biogeographic pattern (Ree et al., 2005; Ree & Smith, 2008; Springer et al., 2011). Therefore, examining the biogeographic history of chiropterans using a more sophisticated model-based method (the dispersal-extinction cladogenesis model) (Ree & Smith, 2008) may facilitate the production of a more refined reconstruction of the ancient range evolution and geographic speciation. In this study, based on a fossil-calibrated phylogeny of nearly all extant chiropteran genera, as well as maximum likelihood-based analyses of the temporal shifts in the net diversification rates (Rabosky, 2006a, b; Stadler, 2011a) and geographic range evolution (Ree & Smith, 2008), we aimed to describe the early diversification pattern and ancestral biogeographic history of bats.

Materials and methods Phylogeny reconstruction and molecular dating Due to the unbalanced availability of sequences from different chiropteran genera, five mitochondrial and one nuclear gene segments were selected to construct the phylogenies, which included 173 of 205 currently recognized chiropteran genera and three outgroup taxa (Tables S1 and S2 in Supporting Information). This nonrandom dense sampling strategy provided a highly representative data set that covered more than 84% of the chiropteran genera from all the major lineages within Chiroptera, although the percentage of missing taxa at the species level was ca 85% (Table S1). The GenBank accession numbers are listed in Table S2. The protein-coding genes were partitioned by gene and were further divided into first + second vs. third codon sites, following Shapiro et al. (2006). We inferred the phylogeny using a Bayesian approach in MrBayes 3.1.2 (Ronquist & Huelsenbeck, 2003) where each data partition was assigned a best substitution model using the Bayesian information criterion scores determined by Modelgenerator (Keane et al., 2006), and the monophyly of two well-established suborders, that is Yinpterchiroptera and Yangochiroptera (Van den Bussche & Hoofer, 2004; Eick et al., 2005; Simmons, 2005; Teeling et al., 2005; Miller-Butterworth et al., 2007; Meredith et al., 2011), was constrained to improve the behaviour of the MCMC sampler (we refer to these topologies as ‘inferred chronograms’ or ‘ic’). To evaluate the influence of alternative higher level phylogenetic relationships on our diversification and biogeographic analyses, we also constrained the phylogenetic position

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

Diversification trend and Asian origin of bat

of Myzopodidae and the monophyly of Noctilionidae and Emballonuridae, according to the results reported by Eick et al. (2005) and Meredith et al. (2011) (we refer to the topologies that match these constraints as ‘constrained chronograms’ or ‘cc’). We used the penalized likelihood method implemented in r8s 1.7 (Sanderson, 2003) to estimate the divergence times with the palaeontological constraints, according to Teeling et al. (2005) and Meredith et al. (2011) (Table S3). We repeated the phylogenetic analyses in MrBayes using a fixed topology that matched the two maximum clade-credibility trees to estimate the uncertainty of the divergence dates. We obtained a set of 1000 divergence date estimates for the nodes in the phylogeny with standard deviations of the estimated ages. To consider the major influence of the maximal age constraints on divergence time estimations (Marjanovi
c & Laurin, 2007), we performed sensitivity analyses that involved various combinations of the upper (maximal) age constraints. The data matrix and chronograms generated in this study are deposited in Dryad (doi:10.5061/ dryad.r717q). Diversification analyses To investigate the early diversification rate, we defined an arbitrary but conservative threshold that included 85% stem nodes (approximately before 25 Ma) where the lineage-through-time (LTT) plot was not significantly influenced up to this level by the nonsampled taxa. Up to the threshold, the diversification analysis was performed in R using Laser (Rabosky, 2006a), TreePar (Stadler, 2011a) and TreeSim (Stadler, 2010, 2011b). A profile likelihood analysis was conducted to investigate the likelihood of diversification rate shifts over time. Subsequently, using the method of Meredith et al. (2011) with truncated chronograms, we fitted a candidate set of diversification models to two alternative chronograms under two time frames (25–63 and 35–63 Ma) and assessed their significance using the likelihood-ratio test (using Laser and TreePar) (Rabosky, 2006a, b; Stadler, 2011a). A likelihood analysis using TreePar (Stadler, 2011a) was also applied to the data sets derived from the sensitivity analyses to examine the effects of different combinations of maximal age constraints on the diversification pattern. The c statistic (Pybus & Harvey, 2000) was also calculated for three different time frames (49–63, 35–63, and 25–63 Ma) to reflect the changes in the diversification rate. Because of the lack of significance of the preferred Yule3rate models for the analyses that covered the period of 25–63 Ma using Laser, we simulated a candidate set of trees with different diversification rates (2000 trees for each condition) to evaluate the probability of Type II errors when determining the diversification rate with two shift time points using TreeSim in R (Stadler, 2011b).

3

One challenge in testing for temporal variation in diversification rates in this study is that our sampling procedure is both incomplete and phylogenetically nonrandom. To overcome these problems, a novel method for handling nonrandom incomplete sampling was also conducted using TreeSim in R (Cusimano et al., 2012). The maximum likelihood speciation and extinction rates that considered information about the nonsampling density (here 85.0%) were calculated using the TreePar package in R (Stadler, 2011a). Two different times of missing speciation events were set as 25–0 Ma, according to our threshold, and 35–0 Ma, which represented a more conservative threshold. For each threshold, the missing branching times were simulated 1000 times, and the simulated times were then added to the empirical branching times, thereby yielding 1000 complete data sets. Likelihood analysis (using TreePar) was then applied to these complete data sets, where their significance levels were determined using the likelihood-ratio test (Stadler, 2011a). Note that extinction was excluded from the model because it is difficult to obtain meaningful separate estimates of extinction and speciation parameters under discrete-shift models as the confidence intervals of the extinction rates are extremely large. To identify potential shifts in the diversification rate of lineages within our chronograms, particularly during 25–35 Ma, we applied the MEDUSA method (Alfaro et al., 2009) in R using package Geiger (Harmon et al., 2008) to our ‘diversity trees,’ which were derived from two maximum clade-credibility chronograms and the refined species richness data (Table S4). MEDUSA method begins by estimating birth and death values and an Akaike’s information criterion (AIC) score for a model with no diversification shifts in the chronogram. The method then fits models of increasing complexity by incorporating a branch where rates of diversification shift, with an additional birth and death value calculated for the clade where the shift point occurred. If the new model has an AIC score that is lower than the previous model by an AIC cut-off value, then the new model is retained. This step-wise procedure continues adding additional shift points until the AIC threshold criterion is no longer met. Ancestral area reconstructions using maximum likelihood The distribution ranges of bats were categorized into six areas based on the present and past separations of major landmasses (Fig. 1). The models of dispersal route availability were developed based on the geographic history and island chains, climatic data and related ancestral area reconstruction (AAR) studies (Scotese, 2001; Morley, 2003; Clayton et al., 2009) (Fig. 1). The distribution assignments for each genus were derived from ‘Mammal Species of the World’ (Simmons, 2005),

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

W. YU ET AL.

4

Model C1: Constrained

Model C2: Constrained

Dispersal between Africa and

Dispersal between Africa and

South America not viable

South America viable

0–70 Ma S S N E F A U

N E F A 1.00 0.10 0.10 0.10 1.00 0.10 1.00 1.00 1.00 1.00

E A

N

0–70 Ma U 1.00 0.10 0.10 0.10 1.00

S S N E F A U

N E F A U 1.00 0.10 1.00 0.10 1.00 1.00 0.10 1.00 0.10 1.00 1.00 0.10 1.00 0.10 1.00

F S

U

Model TV1: Temporally Variable Dispersal between Africa and South America not viable

0–5 Ma S S N E F A U

N E F A 1.00 0.10 0.10 0.10 0.10 0.10 0.10 0.75 1.00 0.75

5–30 Ma U 0.10 0.10 0.10 0.10 0.75

S S N E F A U

N E F A 0.75 0.10 0.10 0.10 0.10 0.10 0.10 0.75 0.75 0.75

45–70 Ma

30–45 Ma U 0.10 0.10 0.10 0.10 0.75

S S N E F A U

N E F A 0.75 0.10 0.10 0.10 0.75 0.10 0.50 0.50 0.50 0.50

U 0.10 0.10 0.10 0.10 0.50

S

N E F A 0.25 0.10 0.10 0.10 0.75 0.10 0.50 0.25 0.25 0.25

S

N E F A 0.25 0.10 0.25 0.10 0.75 0.10 0.50 0.25 0.25 0.25

S N E F A U

U 0.50 0.10 0.10 0.10 0.10

Model TV2: Temporally Variable Dispersal between Africa and South America viable

S N E F A U

N E F A 1.00 0.10 0.10 0.10 0.10 0.10 0.10 0.75 1.00 0.75

30–45 Ma

5–30 Ma

0–5 Ma S

U 0.10 0.10 0.10 0.10 0.75

S S N E F A U

N E F A 0.75 0.10 0.10 0.10 0.10 0.10 0.10 0.75 0.75 0.75

U 0.10 0.10 0.10 0.10 0.75

S S N E F A U

N E F A 0.75 0.10 0.10 0.10 0.75 0.10 0.50 0.50 0.50 0.50

45–70 Ma U 0.10 0.10 0.10 0.10 0.50

S N E F A U

U 0.50 0.10 0.10 0.10 0.10

Fig. 1 Areas assigned to the genera included in this study and the parameters of the four alternative biogeographic models examined. A, Asia (West of Ural mountains); E, Europe (East of Ural mountains); F, Africa and Madagascar; U, Australia/New Zealand/East of Wallace’s Line; N, Central/North America; S, South America. Probabilities of relative dispersal: 0.1, low; 0.25, medium–low; 0.5, medium; 0.75, medium–high; and 1, high. The different probabilities of dispersal among the four time frames are highlighted by the grey cells.

which were refined substantially to one or several ancient areas when palaeontological and/or phylogenic references were available (Table S4; see Data S1 for more detail). Given the potential effect of the convergent evolution of chiropteran morphology on AARs, as well as methodological or time constraints when using fossil data, we did not employ the method utilized by Teeling et al. (2005) and Simmons et al. (2008) to determine the positions of extinct fossil lineages. Instead, we applied a recently developed likelihood method for reconstructing ancestral areas (the dispersal–extinction– cladogenesis model of Lagrange) (Ree et al., 2005; Ree & Smith, 2008) to the maximum clade-credibility chronograms and each 1000 time-calibrated posterior tree under four different models. We set two constrained models (C1 and C2) where dispersal events were possible between all adjacent areas during the overall periods, whereas two other temporally variable models (TV1 and TV2) incorporated an explicit dispersal model in four different time frames (Fig. 1). Note that the direct dispersal route between South America and Africa was only allowed in models C2 and TV2. It was not clear how to compare the likelihood scores for models with changes in the specified ancestral ranges because the models were not nested, and with fewer

ancestral ranges, the additional parameters were fixed at a boundary value of zero (Ree et al., 2005; Ree & Smith, 2008). Therefore, the likelihood scores obtained using different models were compared directly.

Results Phylogeny and molecular dating In this study, the phylogeny of nearly all extant chiropteran genera (approximately 84%) was reconstructed using dense taxonomic sampling (Table S1). Five mitochondrial and one nuclear gene segment were used to construct a approximately 5.6 kb sequence matrix (Table S2). The inferred phylogenetic relationships above the family level were consistent with previous studies (Teeling et al., 2005; Miller-Butterworth et al., 2007), where monophyly was supported for four superfamilies with posterior probability values of 1.0, although the genus-level relationships within some families, such as Molossidae, Hipposideridae and Megadermatidae, were not well resolved (posterior probability < 0.85) (Fig. S1). To evaluate the influence of alternative relationships among families upon our diversification and biogeographic analysis, we also con-

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

Diversification trend and Asian origin of bat

strained topology matching using an alternative family level phylogeny (Eick et al., 2005; Meredith et al., 2011). It should be noted that the inference of early diversification dynamics is highly dependent on the deeper nodes, most of which were well-supported in our chronograms, and the diversification trend can tolerate a degree of uncertainty in the topology (Rabosky, 2006b). The estimates obtained from the inferred and constrained topologies produced similar ages and standard deviations. The divergence time estimates and the sensitivity analyses were comparable with those obtained using the maximum clade-credibility chronograms (Table 1). Moreover, the divergence date ranges derived using various combinations of the maximal age constraints were compatible with the standard deviations from the posterior time trees (Table 1). Although our studies yielded slightly older values, the estimated higher level divergence times were broadly consistent with previously reported family level chronograms (Eick et al., 2005; Teeling et al., 2005; Miller-

5

Butterworth et al., 2007; Meredith et al., 2011). The members of the suborder Yangochiroptera shared a common ancestor at 61
0.8 Ma, whereas Yinpterchiroptera originated at about 58
0.8 Ma (Table 1). At least 11 of the 19 families originated at approximately 50 Ma (Figs 2a and S1 Table 1). Temporal pattern of diversification The chiropteran diversification patterns from two chronograms were examined using LTT plots. Three phases characterized by distinct diversification rates were identified until approximately 25 Ma, that is a high net diversification rate persisted until approximately 49 Ma (first orange circle and patch in Figs 2b and 3a), followed by a dramatically decreased rate at about 35–49 Ma (blue circle and patch in Figs 2b and 3a), with a subsequent acceleration in the rate from approximately 35 Ma onwards (second orange circle and patch in Fig. 2b and 3a).

Table 1 Divergence times of higher level relationship (Myr) within Chiroptera, which were derived from the r8s analysis based on inferred and constrained topologies. Inferred chronogram

Constrained chronogram

Sensitivity analysis

Sensitivity analysis

Nodes

Age

SD

Mean age

Range

Chiroptera [1] Yinpterchiroptera [2] Pteropodidae [3] Rhinolophidae [4] Rhinolophidae + Hipposideridae [5] Hipposideridae [6] Megadermatidae + Craseonycteridae [8] Megadermatidae + Craseonycteridae + Rhinopomatidae [7] Megadermatidae [9] Yangochiroptera [10] Emballonuridae [11] Emballonuridae + Noctilionidae* Emballonuridae [12] Noctilionidae + Myzopodidae [14] Noctilionidae [15] Phyllostomidae + Mormoopidae [17] Thyropteridae + Furipteridae + Noctilionidae [20] Furipteridae + Noctilionidae [21] Mormoopidae [19] Phyllostomidae [18] Vespertilionidae + Noctilionidae + Myzopodidae [13] Vespertilionidae + Myzopodidae* Vespertilionidae [22] Natalidae [23] Molossidae [24] Molossidae + Vespertilionidae + Miniopteridae [25] Vespertilionidae + Miniopteridae [27] Vespertilionidae [26]

63.0 61.5 38.4 51.2 41.9 37.6 45.5 50.1 22.7 58.4 56.0

0.6 1.6 1.6 1.8 2.0 3.2 1.7 1.8 0.8 1.4

61.6 37.5 52.0 43.0 38.1 45.0 50.7 22.4 58.0 55.4

60.1 36.9 49.7 39.2 34.4 42.1 49.1 22.0 57.4 54.7

62.8 41.2 52.2 43.3 38.8 46.9 52.5 23.6 58.9 57.5

49.1 55.0 49.6 40.1 42.3 34.2 36.0 34.0 56.4

1.8 1.1 1.5 1.0 1.9 1.8 2.5 0.2 2.7

49.3 54.7 49.3 40.4 42.1 34.7 36.8 34.6 56.3

48.4 54.1 49.1 39.4 41.6 32.7 33.2 33.8 54.4

51.2 56.4 51.2 41.4 43.6 36.4 39.1 35.8 59.2

53.3 21.2 33.9 51.7 49.1 44.3

2.4 2.0 3.3 1.3 1.5 1.6

53.1 20.9 35.5 51.3 48.3 43.8

51.6 18.2 31.1 50.8 47.7 43.1

55.3 21.4 35.9 52.8 51.2 45.4

Age

SD

Mean age

Range

63.0 61.3 37.1 51.8 43.1 38.1 44.4 50.0 22.7 58.0 53.5 56.8 47.5

0.8 1.5 1.5 1.7 1.9 3.1 1.6 1.5 0.7 1.4 0.9 1.7

61.6 37.8 51.9 43.2 38.2 45.0 50.4 22.3 57.9 53.0 56.5 47.4

60.4 36.2 49.8 41.2 35.7 40.1 48.4 21.2 57.1 51.8 55.1 46.3

63.1 39.5 52.2 43.6 39.9 47.9 52.5 24.6 58.9 57.5 59.1 50.2

49.9 40.3 42.9 34.9 36.2 34.0

1.2 1.0 1.4 1.8 1.6 0.2

49.7 40.5 42.3 34.8 37.0 34.6

49.1 39.1 40.2 32.7 33.2 33.8

51.2 42.4 43.9 36.4 39.1 35.8

56.2 53.8 21.4 36.7 52.0 49.3 44.4

1.0 1.1 1.9 3.6 1.2 1.4 1.5

56.3 53.8 21.4 36.2 52.2 49.1 44.3

55.5 51.3 19.3 31.1 50.8 47.8 43.1

58.5 54.9 21.9 38.9 53.1 50.9 46.3

The numbers that follow the node names are the related nodes shown in Fig. 4. The asterisk indicates the distinct node when using the constrained topology. ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

W. YU ET AL.

(b)

Ve

Early Eocene Climatic Optimum (EECO)

6 5 4

8 °C 4 °C 0 °C

1

Speedup Speedup Paleocene

70

70

60

60

50

Slowdown –4 °C Eocene

50

40

30

40

30

Miocene

Oligocene

20

10

20

Plio.

10

0 (Ma)

0 (Ma)

Pte

rop

odi dae

Outgroups Melonycteris Eidolon Pteralopex Desmalo Acerod pex Pterop on Notop us Macr teris Syc oglossus Sphaonycteris Th erias Ae oopteru Ba thalop s Aliolionyctes ris Ha n y c O plon teris La topte ycte C tide rop ris Pe hiron ns us D nth ax M yac eto C eg op r P yn ae ter Do ten opterops us A b oc r H pr so hi us R ar o t ni rus Eo ou pyi eles a ny set ony ct tus ct er er is is

ct

ilio

ni da e

us ss lo s og s ru al op ho pus e g m p ro i s M po omo pte ter s E o c ri Ep icr ony cte e ne M iss ny en me L yo tim cti us M yc ny ph os N ra l o r Pa hinooside R ipp ia s H sell cu A sellis ps ris te A elo Co inonics Rh eoti ps Cloiaeno ma Tr inopo ycteris Rh seon a Cra gaderm a Me ioderm Card derma Macro ris Nycte llonura Emba Mosia pteryx Balantio Cyttarops Diclidurus Cormura Peropteryx Centronycteris Rhynchonycteris Saccopteryx Saccolaimus Taphozou s

No

V Pl ampCar a y o VaMes tyrrh rodllia Ch mp o p h inu es iro yre y l s U St ro de ss la en de rm a od rm a Py Phy erm a g ll a Sp ha Amoder ops ero e m n tr a Ceycte ida ris Ardnturio Arti ops EncEctophbeus histh ylla e S turnnes R Glyp hinophy ira h lla Microonycteris nycte Mac ris tus Mystaro ci na Myzopod a

Miniopteridae Furipteridae Noctilionidae Thyropteridae Mormoopidae

12 °C

3

sp

er

tili

Ba

O Pi ton p y Ny istr cte S Ne co ct ellu ris o t a Ny La romoec lus s e S ctic ph ic us Ty coto ein otisia lon re op y p s Ny Hyp cterei ns Ch cto su s a p g V lino hilu o Scoesperlobuss tom tilio an His Ia es Eptetiotus Lasio s nycteicus Ari ri N elulu s Glauc ycticeiuss on Rhogyecteris essa Antrozo us Scotophi lu Cistugos Miniopterus Furipterus Noctilio Thyroptera Pteronotus Mormoops Diaemus s Desmodu Diphylla Tracholups m yl Macroph derma PhylloMimon us ostom a Phyll hostomtia Lop Tona um r s y p Vamtopteruina o rh Chr ncho cterisa Lo llony hyll a p ll y Ph Ero yphyeris h c yct g a Brapton pha llus o e s L los phy terisis G ono nyc cter us M so ny isc ris Mu ero ron cte ura a o e y o n ChCho ylon An tali eris lla a t H Pl nyc p h y o o Li ch n Lo

a=0 a = 0.9

2

e ida on

Co

Natalus Chilonatalus Nyctiellus Tomopeass Cheiromele Eumopss Promopus Moloss s opteru Morm omopss n Nycti mop Cynodarida Ta M o p s n pho ere ops Cha ss ps lo Mo tomoina O ur lus M pha la ce ivou tis o s rpio er K Myopu tis Ha c o dis y llus Eu erimstre teris a P ra c m y r s Pa on de tu s Idi Eu leco inu us P o r h i u r lla n s te r y L a as rb

Mo

Log (surviving Chiroptera lineages)

Natalidae ae lossid

Earliest Oligocene Glacial Maximum (EOGM)

7

(a)

Paleotemperature (°C)

6

e n i d a ae taci did Mys yzopo M

Rhin oloph idae oside ridae R h inopo C e r rida a onu s ball Me ma d Em eo ae Nycte gadermnaycteridtia tidae e ridae Hipp

Fig. 2 Genus-level chronogram of Chiroptera (inferred chronogram) and the associations between their diversification patterns and major changes in global palaeotemperature. (a) Genus-level chronogram showing the phylogenetic relationships within extant Chiroptera (inferred chronogram). The major lineages are coloured as follow: red, Pteropodidae; green, Rhinolophoidae; brown, Emballonuroidae; blue, Noctilionoidae; and black, Vespertilionoidae. The branch lengths are proportional to time. Clades with unusual diversification rates are indicated by the red star-shaped symbols. The phases characterized by distinct net diversification rates are indicated by the orange and blue circles. (b) Semi-logarithmic lineage-through-time (LTT) plots (circle: empirical LTT plots; dashed line: 95% CI) and constant rate diversification simulation. The black curve indicates the plots from the inferred chronogram, and the purple curve represents the plots from the constrained chronogram. The dark grey curves were obtained from 1000 complete data sets under a threshold of 25 Ma, and the light grey curves were from 1000 data sets under a threshold of 35 Ma. The two orange patches indicate distinct periods characterized by a rapid increase in diversification, whereas the blue patch indicates a period when the diversification rates were reduced. The red smoothed curve represents the temporal variation in palaeotemperature (Zachos et al., 2001). The green curve represents the null distribution of the LTT plots under a relative extinction rate of 0.9, and the blue curve represents the curve under a relative extinction rate of 0.

These patterns of diversification were also verified using Kolmogorov–Smirnov goodness-of-fit tests, maximum likelihood-based analyses of the diversification rate, c statistic, power analysis and simulations to handle nonrandom sampling (Rabosky, 2006a, b; Stadler, 2011a, b; Cusimano et al., 2012), and most results were

generally consistent (Please see the Supporting Information for the further details). The differences in the Akaike’s information criterion (DAIC) score for the two best Yule3rate models when using Laser (Rabosky, 2006a) were close to the critical value (alpha = 0.05) for the analysis that covered the period from 63 to

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

Diversification trend and Asian origin of bat

(b) 0.20 0.15

60

50

Chattian

Rupelian

Bartonian

Priabonian

40

–12.0 65

Speciation rate (sp Myr–1)

0.15 0.10

60

50

40

Chattian

Rupelian

Bartonian Priabonian

Lutetian

0.05

70

–8.0 55

45 35 (Ma)

25

(c)

0.20

0

–4.0

30 25

(Ma)

Ypresian

Net diversification rate (sp Myr–1)

70

Lutetian

Danian Selandian Thanetian

0

Ypresian

0.05

Log-likelihood

0

0.10

Danian Selandian Thanetian

Net diversification rate (sp Myr–1)

(a)

7

0.12 0.08 0.04 0 65

55

45 35 (Ma)

25

30 25

(Ma) Fig. 3 Average net diversification in each geological stage and profile analysis of the diversification rates in Chiroptera. The most recent 25 Myr of bat evolutionary history were excluded from the analyses. (a) Average net diversification during each geological stage from different time trees. The upper histogram in black is based on the inferred chronogram, and the lower grey histogram shows the results from the constrained chronogram. The phase prior to 49 Ma is divided into three phases of equal duration due to the limited lineagethrough-time (LTT) plots in this period. (b) Log-likelihood of temporal changes in chiropteran diversification rates. The results are based on a sliding windows analysis using a fixed window of 14 Myr. (c) Net diversification based on the time curve of extant Chiroptera, which was inferred using a fixed window of 14 Myr. The grey patches in panels (b) and (c) indicate the phase of 35–49 Ma with the maximum log-likelihood score. The black and grey trends in panels (b) and (c) indicate the diversification curves derived from the inferred and constrained chronograms, respectively.

25 Ma (ic: DAICYule3rate(63–25) 4 = 2.5, P = 0.13; cc: DAICYule3rate(63–25) 4 = 1.9, P = 0.17). However, the likelihood analysis for the period from 63 to 35 Ma using Laser (Rabosky, 2006a) (ic: DAICYule2rate(63–35) 2 = 4.7, P < 0.05; cc: DAICYule2rate(63–35) 2 = 4.4, P < 0.05) and the v2 value obtained from the likelihood-ratio test using TreePar (Stadler, 2011a) (ic: v2Yule3rateð6325Þ3 = 8.0, P < 0.05; v2Yule2rateð6335Þ1 = 6.2, P < 0.05; cc: v2Yule3rateð6325Þ3 = 8.2, P < 0.05; v2Yule2rateð6335Þ1 = 5.1, P < 0.05) confirmed the observed fluctuation in the chiropteran diversification trend (Table 2). The lack of statistical significance for the Yule3rate models using Laser may have been attributable to a bias towards a high likelihood score for the empirical diversification data with alternate high-low net diversification rates when fitting rate-constant models or to a high AIC penalty during the analysis of diversification patterns that included multirate shifts (Rabosky, 2006b). Indeed, the power analysis showed the probability of rejecting rate-constant models (when they were incorrect) that were similar to our chronograms was < 50% (Fig. S2). According to the sensitivity analyses, the Yule3rate model provided a better fit for

100% of the data sets. Moreover, most detected a similar temporal pattern, and the time points found in the diversification analyses using the maximum clade-credibility chronograms were also similar (approximately 94% of the inferred chronograms; approximately 82% of the constrained chronograms). Although only approximately 39% of the inferred chronograms and approximately 30% of the constrained chronograms could be rejected based on significant constant rate models (v2critical value 3 = 7.81, P < 0.05), most of the data sets were close to the critical value (Fig. S3). To overcome the problems caused by our overdispersed sampling procedure, a novel method called CorSiM was used for nonrandom incomplete sampling (Cusimano et al., 2012). The speciation and extinction rates estimated for the tree simulations using CorSiM were s = 0.09 and l = 0 for both the inferred and constrained chronograms under a threshold of 25 Ma, and s = 0.08 and l = 0 for both chronograms under a threshold of 35 Ma. The resulting complete data sets are visualized as LTT plots in Fig. 2b, and the results are summarized in Table S5). In both cases, TreePar

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

5

5

8

Birth–death with one time shift

Yule3rate

Birth–death with two time shifts

5

4.7/4.4 6.7/6.4 0.2/1.7 0.4/0.6 0.0*/0.0*

–/–

–/–

62.8/61.7 64.8/63.7 58.3/59.1 58.5/59.1 59.3/58.0

0.0NS/0.0NS

–/–

–/–

21.4/20.7

4.5/3.5 3.6/3.5

4.5/3.5

2.5/1.9 4.0/2.3

DAIC

25.8/24.2 25.0/22.8

25.9/24.2

23.9/22.2 25.3/22.9

AIC

k = 0.07/k = 0.07 k = 0.07; a = 0/k = 0.07; a = 0 k = 0.16; k = 28.25/k = 0.16; k = 27.1 k = 0.50; xp = 0.79/k = 0.47; xp = 0.77 k1 = 0.15; k2 = 0.04; Ts1 = 49/k1 = 0.14; k2 = 0.04; Ts1 = 49

k1 = 0.14; k2 = 0.04; k3 = 0.10; Ts1 = 49; Ts2 = 35/k1 = 0.13; k2 = 0.04; k3 = 0.10; Ts1 = 49; Ts2 = 35 –/–

6.2/5.1 6.2/5.1 –/– –/– 0.0*/0.0*

0.1/0.1

66.1/63.7

0.1/0.1

0.0*/0.0*

69.2/66.2 69.2/66.2 –/– –/– 66.1/63.7

178.7/176.5

178.6/176.5

2.8/1.7

–/– 6.2/5.1

–/– 181.7/179 180.0/177.3

–/–

–/–

k = 0.09; k = 1101343.7/k = 0.09; k = 1251182.8 k = 0.09; xp = 0.04/k = 0.08; xp = 0.01 k1 = 0.14; k2 = 0.08; Ts1 = 49/k1 = 0.07; k2 = 0.12; Ts1 = 28 –/–

8.6/8.6 8.0/8.2

2DlnL(v2)

182.9/180.8 182.7/180.5

lnL

k = 0.09/k = 0.09 k = 0.07; a = 0.30/k = 0.07; a = 0.40

Parameters

TreePar

k = 0.07/k = 0.07 k = 0.07; a = 0/k = 0.07; a = 0 –/– –/– k1 = 0.14; k2 = 0.04; Ts1 = 49/k1 = 0.14; k2 = 0.04; Ts1 = 49 k1 = 0.13; a1 = 0.02; k2 = 0.04; a2 = 0.01; Ts1 = 49/k1 = 0.12; a1 = 0.18; k2 = 0.04; a2 = 0.01; Ts1 = 49

–/– k1 = 0.07; k2 = 0.10; Ts1 = 35/k1 = 0.07; k2 = 0.10; Ts1 = 35 k1 = 0.12; a1 = 0.66; k2 = 0.03; a2 = 0.72; Ts1 = 49/k1 = 0.15; a1 = 0.47; k2 = 0.02; a2 = 0.84; Ts1 = 47 k1 = 0.14; k2 = 0.04; k3 = 0.10; Ts1 = 49; Ts2 = 35/k1 = 0.13; k2 = 0.04; k3 = 0.10; Ts1 = 49; Ts2 = 35 k1 = 0.10; a1 = 0.41; k2 = 0.04; a2 = 0.01; k3 = 0.10; a3 = 0.01; Ts1 = 49; Ts2 = 35/k1 = 0.15; a1 = 0.01; k2 = 0.03; a2 = 0.52; k3 = 0.08; a3 = 0.34; Ts1 = 49; Ts2 = 35

k = 0.09/k = 0.09 k = 0.07; a = 0.28/k = 0.07; a = 0.41 –/–

Parameters

0.27NS/0.33NS

2.19*/2.11*

0.65NS/1.03NS

c

Rate-constant models are denoted as RC, the time points of diversification rate shifts are denoted by Ts, the number of parameters in each model is given (NP), and the c statistic and the AIC or likelihood (lnL) scores of the models are also shown. Values on the left of the slash are based on the inferred chronogram (ic), where those on the right are based on the constrained chronogram (cc). Boldface indicates the preferred model, where the significance of the observed DAIC value was assessed by simulating 5000 chronograms using Laser, and the significance of the observed v2 value was assessed by a likelihood-ratio test using TreePar. NS Model not significant compared with the best constant rate model (P > 0.05). *Model significant compared with the best constant rate model (P < 0.05).

Pruning most recent 50 Ma evolutionary history

Birth–death with one time shift

1 2 2 2 3

2 3

DDX Yule2rate

Pruning most recent 35 Ma evolutionary history Pure birth (RC) Birth–death (RC) DDL DDX Yule2rate

2

1 2

Pruning most recent 25 Ma evolutionary history Pure birth (RC) Birth–death (RC)

DDL

NP

Model

Laser

Table 2 Maximum likelihood analysis of the diversification rate and c statistic estimates using the two chiropteran chronograms.

8 W. YU ET AL.

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

Diversification trend and Asian origin of bat

(Stadler, 2011a) preferred the Yule3rate model for 100% of the data sets, and it supported diversification rate shifts at 35 and 49 Ma (Table S5). Moreover, 100% of the complete data set significantly rejected the rate-constant diversification models under a threshold of 25 Ma, and approximately 99% of the complete data set rejected the rate-constant diversification models under a threshold of 35 Ma (Table S5). These results are consistent with our previous diversification analyses based on truncated chronograms. These results demonstrate the nonconstant fluctuations in chiropteran diversification, and they also support the 3.5-fold decrease in the net diversification rate during 35–49 Ma (0.04 sp Myr1 for both chronograms) compared with the initial high rate (63–49 Ma) (0.14 and 0.13 sp Myr1 for the inferred and constrained chronograms, respectively), as well as the subsequent 2.5-fold increase in the net diversification rate (0.10 sp Myr1 for both chronograms) compared with the previous slow rate until 25 Ma (the temporal threshold established in this study) (Figs 2b and 3a; Table 2). Variation in the lineage diversification rate Using the comparative method MEDUSA (Alfaro et al., 2009), we found evidence for a relatively high background net diversification rate (k = 0.07 sp Myr1 for both the inferred and constrained chronograms) compared with most other vertebrate clades (Alfaro et al., 2009; Yu et al., 2012) and an extremely low, but unrealistic, relative extinction rate (a = 3.1 9 108 for both chronograms) (Fig. S1). We also identified similar clades with significant diversification rate shifts using two types of time tree (Figs 2a and S1). No branches or nodes exhibited exceptionally high levels of diversification during 25–35 Ma. It is important to determine whether the clade shift model provides a better fit to the data than the temporal shift model; however, the MEDUSA likelihood obtained for a completely sampled unpartitioned Yule or birth–death model was different from those generated by Laser or TreePar (which only use information derived from the branching times), thereby demonstrating that MEDUSA and non-MEDUSA models cannot be compared in a model selection framework (Harmon et al., 2008; Alfaro et al., 2009). Furthermore, considering the different principles and focuses of clade shift model (MEDUSA) and temporal shift model (Laser and TreePar) on diversification, we argued that the results from MEDUSA and nonMEDUSA were complementary rather than conflicting. Biogeographic analyses The inferred and constrained maximum clade-credibility chronograms, models TV1 and TV2, which incorporated dispersal routes and the dispersal possibilities

9

between continents at historical intervals, produced higher likelihood scores (over 20 log-likelihood units) than the unconstrained models (C1 and C2) (Tables S6 and S7). The results obtained using models TV1 and TV2 yielded similar likelihood scores and AARs. Four models based on different maximum clade-credibility chronograms consistently suggested an Asian origin for extant bats (Figs 4a,b and S4, and Tables S6 and S7). These results were consistent with most of the AARs generated using time-calibrated posterior trees (> 97% for each model using the inferred chronogram and > 61% for each model using the constrained chronogram) (Tables S6 and S7). With different chronograms, models TV1 and TV2 also explicitly suggested the Asian origin of two suborders (Figs 4a,b and S4, and Tables S6 and S7). The optimal AARs for model TV2 using the inferred chronogram required 137 dispersal events, 14 vicariance events and 16 local extinction events to explain the current cross-continental distribution of bat genera, whereas the constrained chronogram required 127 dispersal events, 12 vicariance events and 16 local extinction events. With model TV2, direct dispersal events between South America and Africa were not generated by both AARs.

Discussion Early diversification pattern of bats If the early chiropteran ancestors experienced a rapid radiation, it has been hypothesized that there should be a signature of a rapid diversification burst, followed by a decrease in diversification due to intensified competitive pressures and saturated ecological niches (Yoder et al., 2010). The rapid diversification burst we detected from the origin of crown bats until 49 Ma and the significant decrease in diversification during 49–35 Ma appear to agree with this hypothesis. However, the second striking increase in the diversification rate at 25– 35 Ma and the specific times when the decreases in the rate occurred in our empirical data (Figs 2b and 3a; Tables 2 and S3) indicate that density-dependent cladogenesis (Yoder et al., 2010) may not have been the only mechanism related to the significant slowdown in the net diversification rate. The high net diversification observed at 49–63 Ma (0.14 sp Myr1 for the inferred chronogram, 0.13 sp Myr1 for the constrained chronogram), which corresponded to approximately twice the average diversification rate (0.07 sp Myr1), indicates that a rapid radiation occurred during this period (Figs 2b and 3a; Tables 2 and S3). Importantly, the widespread distribution and the diverse morphology of primitive ‘Eochiroptera’ (Simmons & Geisler, 1998; Smith et al., 2007; Simmons et al., 2008) suggest increasing trait variation as well as increased exploitation of unsaturated ecologi-

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

10

W. YU ET AL.

(a)

(b)

Model

Node 1 Node 2 Node 3 Node 4 Node 5

C1

A|A (0.62) A|A (0.87) A|A (0.48) A|A (0.90) A|A (0.75) A|A (0.69) A|A (0.97) A|A (0.39) A|A (0.87) A|A (0.58) A|A (0.91) A|A (0.78) A|A (0.72) A|A (0.97)

TV2

A|A (0.63) A|A (0.87) A|A (0.48) A|A (0.89) A|A (0.40) A|A (0.87) A|A (0.58) A|A (0.91) A|A (0.63) A|A (0.91) A|A (0.74) A|A (0.96) A|A (0.39) A|A (0.90) A|A (0.43) A|A (0.97) A|A(0.61) A|A(0.92) A|A(0.75) A|A(0.96) A|A(0.40) A|A(0.89) A|A(0.43) A|A(0.97)

Model

Node 8

C2 TV1

Node 6

A|A (0.75) A|A (0.78) A|A (0.92) A|A (0.95) A|A(0.92) A|A(0.95)

A|A (0.69) A|A (0.72) A|A (0.76) A|A (0.83) A|A(0.76) A|A(0.83)

Node 7

A|A (0.97) A|A (0.97) A|A (0.99) A|A (0.99) A|A(0.99) A|A(0.99)

Current distribution 3

Pteropodidae

A

F A U

A

C2 TV1 TV2

A|A (0.90) A|FAU(0.16) A|A (0.67) A|A (0.87) A|A (0.88) A|FAU(0.17) A|A (0.35) A|A (0.61) A|A (0.89)A|FAU(0.16) A|A (0.67) A|A (0.88) A|A (0.88)A|FAU(0.17) A|A (0.37) A|A (0.62) A|A (0.97) A|A (0.50) A|A(0.67) A|A (0.94) A|A (0.99) A|A (0.50) A|A (0.41) A|A (0.77) A|A (0.97) A|A (0.51) A|A (0.63) A|A (0.93) A|A (0.99) A|A (0.50) A|A (0.42) A|A (0.78)

A|A (0.76) A|A (0.59) A|A (0.76) A|A (0.60) A|A (0.91) A|A (0.79) A|A (0.89) A|A (0.79)

A|A

(0.62)

-

A|A

2 A

(0.43)

-

A

5 4

A|A (0.63) A|A (0.45) A|A (0.61) A|A (0.62) A|A (0.57) A|A (0.58) -

A

C1 C2 TV1 TV2

A

A|N (0.46) U|S (0.44) A|N (0.48) A|N (0.45) A|N (0.63) A|N (0.47) A|N (0.59) A|N (0.50)

1

N|N(0.48) S|S (0.46) N|N(0.49) S|S (0.44) N|N (0.64) N|N (0.48) N|N (0.62) N|N (0.52)

S|S (0.44) S|S (0.63) S|S (0.42) S|S (0.61) S|S (0.51) S|S (0.43) S|S (0.49) S|S (0.40)

S|S (0.56) S|S (0.71) S|S (0.55) S|S (0.69) S|S (0.37) S|S (0.54) S|S (0.39) S|S (0.51)

A|A (0.31)

A|A(0.50) N|N(0.72) NA|N(0.63)A|A(0.63) A|A(0.60) A|A(0.86) A|A(0.61) N|N(0.68) NA|N(0.25)A|A(0.73) A|A(0.57) A|A(0.82)

-

A|A (0.93) A|A (0.88) A|A (0.92) A|A (0.81)

14 A

South America

N

North and Central America

A|A (0.63)

A|A(0.64)

A|A (0.50) A|A (0.53)

A|A (0.66) A|A (0.60)

18 N

N

16

A

S

13

N

19

S 20

S

A

21

S S

A

Africa (including Madagascar)

A

Asia and insular SE Asia

N

Europe

A 26

N S F A U

Myzopodidae Mystacinidae

U

Phyllostomidae

N S

Mormoopidae Thyropteridae Furipteridae Noctilionidae

N S N S S

Miniopteridae

F E A U

60

F

N S

Vespertilionidae N S F E A U

A

24

Eocene

Paleocene

70

Emballonuridae

27

A N 23

Australia (including New Zealand )

F A

N 17

S

-

A|A(0.33)

F

U

Nycteridae

A

25 A A N 22

E

A E A U

N

15 N

A

S

A A

A

A|A (0.49) N|N (0.72)NA|N(0.63) A|A (0.62) A|A (0.60) A|A (0.86) A|A (0.62) N|N (0.68)NA|N(0.24) A|A (0.73) A|A (0.57) A|A (0.82)

A|A (0.90) A|A (0.88) A|A (0.89) A|A (0.81)

A

12

A

C2

N|N (0.81) NA|N (0.84) A|A (0.47) N|N (0.85) NA|N (0.53) A|A (0.46) N|N (0.87) NA|N (0.82) A|A (0.46) N|N (0.81) NA|N (0.41) A|A (0.43)

A

10

C1

NA|N(0.45) NA|N(0.47) NA|N (0.45) NA|N (0.40)

A

A

Node 22 Node 23 Node24 Node 25Node 26 Node 27 Emb+Noc Myz+Ves

TV2

A

A

Model

TV1

9

E A U E A U

A

Node 15 Node 16 Node17 Node 18Node 19 Node 20 Node 21 S|S (0.42) N|N(0.51) N|N (0.59) S|S (0.53) S|S (0.40) N|N (0.44) S|S (0.44) N|N (0.52) N|N (0.60) S|S (0.51) S|S (0.38) N|N (0.46) S|S (0.68) N|N (0.72) N|N (0.79) S|S (0.49) N|N (0.52) N|N (0.64) S|S (0.65) N|N (0.69) N|N (0.78) S|S (0.53) N|N (0.56) N|N (0.69)

A A

8

A 11

Model

A

A

7

A

F Rhinolophidae Hipposideridae F Rhinopomatidae F Craseonycteridae A Megadermatidae F

A A 6

50

40

Oligocene

30

N

Miocene

20

10

Plio.

Molossidae

N S F E A U

Natalidae

N S

Plt.

C1

Node 9 Node10 Node 11Node 12 Node 13 Node 14

0 (Ma)

Fig. 4 Ancestral area reconstructions (AARs) of extant chiropteran families based on the genus-level AAR using the model of Lagrange. (a) AARs of extant chiropteran families under model TV2 using the inferred chronogram. AARs with the highest likelihood are shown as coloured boxes at each node. Each box indicates the ancestral ranges inherited by the daughter lineage originating from the node. (b) Summary of the relevant AARs for the nodes under different models and chronograms. The numbers in the table correspond to the circle tags in panel (a). Note that the areas on the left-hand and right-hand sides of the bar are inherited by the upper and lower branches, respectively. The upper AAR in each cell was reconstructed using the inferred chronogram, whereas the lower AAR was derived from the constrained chronogram. ‘Ves,’ ‘Emb’ and ‘Myz’ in panel (b) represent Vespertilionidae, Noctilionidae and Myzopodidae, respectively.

cal niches throughout the globe, thereby indicating the potential adaptive nature of this radiation (Simpson, 1949, 1953; Schluter, 2000). In theory, the occupation of ‘adaptive zones’ (Simpson, 1949, 1953; Schluter, 2000) via the evolution of key innovations and the dispersal into new habitats could have triggered a rapid adaptive radiation caused by an ecological release from competition and other selective pressures (Yoder et al., 2010). A recent palaeontological study showed that an Eocene bat, O. finneyi, was capable of powered flight but not echolocation, which implies that flight evolved before echolocation in chiropteran history (Simmons et al., 2008), thereby suggesting that the acquisition of flight may have triggered an adaptive radiation because it allowed ancient bats to access new ‘adaptive zones’ by improving their ability to exploit new habitats and food resources. In addition, given the presumed echolo-

cation capabilities of some Early Eocene bats (e.g. Icaronycteris, Archaeonycteris, Palaeochiropteryx and Hassianycteris) (Simmons & Geisler, 1998), the acquisition of echolocation may be another intrinsic factor that promoted the radiation by strengthening and refining their capacities for orientation and obstacle/prey detection. The high diversification rate observed during the Ypresian (Fig. 3a) also coincided with the warming trends in Tertiary and the ‘Early Eocene Climatic Optimum’ (EECO) (Zachos et al., 2001), which was a period characterized by a major increase in plant diversity and the expansion of the megathermal rainforest region, where the peaks of Tertiary insect diversity and abundance would have increased the habitat complexity (Wing & Sues, 1992; Wilf & Labandeira, 1999; Teeling et al., 2005; Morley, 2007). We suggest that these extrinsic factors would also have promoted the long-

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

Diversification trend and Asian origin of bat

term radiation of bats by facilitating their exploitation of new ecological opportunities. Considering the profound effects of temperature on ecosystem functioning and biodiversity, and the close link between bat ecology and external temperature (Stevens, 2004, 2006; Mayhew et al., 2008; Buckley et al., 2010), we also examined the association between the diversification rate and palaeotemperature trends to investigate the abiotic factors that drove chiropteran diversification. An association would be expected if the diversification trend of bats was influenced or determined by palaeotemperature changes. The initial rapid diversification burst (0.14 sp Myr1 for the inferred chronogram, 0.13 sp Myr1 for the constrained chronogram) and the subsequent decrease in the diversification rate until 35 Ma (0.04 sp Myr1 for the inferred and constrained chronograms), respectively, coincided with a warming trend from the Early Palaeocene until the EECO and a subsequent cooling phase initiated in the Middle Eocene that persisted until the Late Eocene (approximately 35 Ma) (Figs 2b and 3a; Tables 2 and S3), thereby suggesting that this extrinsic physical factor may have played an important role in chiropteran diversification during this period. The drastic decline in the palaeotemperature during the long cooling phase led to substantial climatic and environmental changes, including a contraction of the megathermal rainforest region and reductions in the diversity and abundance of plants and insects (Labandeira, 2005; Jaramillo et al., 2006; Wilf et al., 2006; Morley, 2007; Currano, 2009). These changes may have reduced ecological opportunities as well as intensified the competitive pressures and risk of extinction for bats, which could have decreased speciation and/or increased the extinction rate. Dramatic increase in the diversification rate in the Late Eocene–Oligocene The unexpected increase in net diversification (from 0.04 to 0.10 sp Myr1 in the inferred and constrained chronograms) in the Late Eocene and Oligocene (25–35 Ma) cannot be explained by the palaeotemperature (Figs 2b and 3a; Tables 2 and S3) because this period was characterized by the radical cooling events of the ‘Earliest Oligocene Glacial Maximum’ (approximately 33.5 Ma) (Katz et al., 2008) and persistent cool conditions (Zachos et al., 2001) (Fig. 2b). Furthermore, according to the crown bat chronograms for different ages, the time that increased diversification occurred overlapped with the cooling period during the Late Eocene–Oligocene (approximately 38–34 Ma) (Fig. S5). Indeed, palaeontological and phylogenetic studies of the diversification of plants and insects, and of megathermal plant diversity, link the large-scale extinctions and major floral and faunal turnovers observed from the Late Eocene until the Oligocene with the drastic palaeoclimatic shift (Labandeira, 2005; Jaramillo et al.,

11

2006; McKenna & Farrell, 2006; Morley, 2007). Therefore, what are the factors that might have contributed to increased diversification under these severe environmental conditions? Although no specific node/clade exhibited an exceptionally high diversification rate in the Late Eocene and Oligocene (25–35 Ma) (Figs 2 and S1), the overall increase in the diversification rate was probably attributable to the accumulation of lineages within three clades (Pteropodidae, Phyllostomidae + Mormoopidae and Vespertilionidae), which accounted for approximately 70% of the new lineages during this period (Figs 2a and 5). The lineages identified did not match perfectly with the lineages that exhibited significant diversification shifts during 30–50 Ma according to Jones et al. (2005), but it was difficult to evaluate the differences and specific contributions because the topologies, ages and methods were different. However, both sets of results emphasize the importance of the Middle and Late Eocene in shaping the current nonuniform biodiversity and differences in diversification among chiropteran lineages. We suggest that the rapid diversification under these severe environmental conditions may have been related to novel or refined innovations or adaptations in bats, which facilitated unusual radiations. For the phyllostomids + mormoopids and pteropodids, which are clades that are typically associated with tropical niches (Simmons, 2005), the selective pressures related to reduced prey availability and habitat resources may have intensified during this period because of the contraction of tropical areas and rapid global cooling trends (Labandeira, 2005; Jaramillo et al., 2006; McKenna & Farrell, 2006; Morley, 2007). Shifts in feeding habits to frugivory, nectarivory and omnivory to escape the constraints of insectivory (Simmons & Conway, 2003) may have facilitated the expansion of their ecological niches and opportunities (Simmons & Conway, 2003; Stevens, 2006; Buckley et al., 2010). The herbivorous skeletal adaptations of a presumed pteropodid fossil from the Late Eocene appear to confirm this type of feeding shift during their early history (Ducrocq et al., 1992, 1993). In addition, reconstructions of the ancestral feeding traits of phyllostomids suggest that adaptive changes occurred throughout phyllostomid history (Monteiro & Nogueira, 2011; Rojas et al., 2011; Dumont et al., 2012). In contrast to the subtropical and tropical distributions of phyllostomids and pteropodids, existing vespertilionids appear to be able to extend their habitats into temperate or even into subarctic zones (Altringham, 1996; Parker et al., 1997; Stevens, 2004; Buckley et al., 2010). Specific intrinsic reproductive, physiological and behavioural traits (e.g. torpor, hibernation, migration, delayed fertilization and multiple offspring) may be ‘key innovations’ that facilitated their radiation, which could have helped to overcome the energetic constraints associated with inhabiting areas beyond the warm megathermal

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

W. YU ET AL.

(b)

10

0

100 50

10

Molossidae Pteropodidae Rhinolophidea Emballonuridea Noctilionidea Vespertilionidae &Miniopteridae

1 60

50

40

30

20

10

Pteropodidae

20

10 5

Vespertilionidae

30

0

(Ma)

rainforest zone during this difficult period (Altringham, 1996; Stevens, 2004; Yu et al., 2012). Further research is required in this area in the future. Asian origin and ancient biogeographic history of bats Another unresolved issue in the macroevolution of bat concerns their geographic provenance (Eick et al., 2005; Teeling et al., 2005; Springer et al., 2011). Using the maximum likelihood AAR method (Ree et al., 2005; Ree & Smith, 2008), all of the models derived from maximum clade-credibility chronograms indicated an Asian provenance for extant bats, which was also confirmed by the AARs based on time-calibrated posterior trees (Figs 4a,b; Tables S6 and S7). These reconstructions are consistent with the traditional Laurasian origin hypothesis of bats (Cracraft, 1973; Teeling et al., 2005), but they also suggest an Asian origin, instead of Northern American or African sources (Eick et al., 2005; Teeling et al., 2005). Eocene fossil records of bats in Asia are rare compared with those from Northern America (Gunnell & Simmons, 2005), but the recent discovery of a fossil deposit from the oldest and only known bat fauna from the Early Eocene in India, which represents the most abundant and diverse bat

0 Pteropodidae

20

Vespertilionidae

30 (Ma)

Phyllostomidae&Mormoopidae

40

Molossidae

50

Percentage of new lineages

60

Molossidae

Molossidae Pteropodidae Rhinolophidea Emballonuridea Noctilionidea Vespertilionidae &Miniopteridae

5

Phyllostomidae&Mormoopidae

0

10

1

Number of surviving Chiroptera lineages

10 Rhinolophidea

50

20

Rhinolophidea

Number of surviving Chiroptera lineages

100

30

Emballonuridea

Percentage of new lineages

(a)

Emballonuridea

12

Fig. 5 Semi-logarithmic lineagethrough-time (LTT) plots for the six major clades and the percentages of new lineages during 25–35 Ma. Note that the two upper panels are based on the inferred chronogram, whereas the lower panels are derived from the constrained chronogram. (a) Semilogarithmic LTT plots for the six major bat lineages. (b) Percentage of new lineages in each clade during 25– 35 Ma.

fauna during this period (Smith et al., 2007), suggests a rich chiropteran diversity during the Eocene in Asia. It should be noted that the earliest known fossil bats, that is O. finneyi and Icaronycteris index, were found only in the Early Eocene beds of North America (Simmons et al., 2008) and Asia (Smith et al., 2007). By contrast, the North American origin (Teeling et al., 2005) might not be supported because I. sigei, which was previously found only in North America, was also found in the Asian deposit (Smith et al., 2007). The proposed African origin (Eick et al., 2005) also seems to be debatable because Springer et al. (2011) reconstructed an inclusive provenance using the identical data and settings. Moreover, the early fossil records from the Southern Hemisphere are also poor, with only three fragmentary specimens known from the Early Eocene (Ravel et al., 2011). With different chronograms, both the TV1 and TV2 models explicitly suggest an Asian origin for two suborders, which resembles the reconstruction of Teeling et al. (2005) (Fig. 4; Tables S6 and S7). Furthermore, the similarities of the likelihood scores and the ancestral reconstructions with models TV1 and TV2 (Fig. 4; Tables S6 and S7) indicate that allowing an early direct dispersal between Africa and Southern America (probability of 0.25) would not change the AAR optimizations

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

Diversification trend and Asian origin of bat

greatly. The main disjunct distribution patterns between New World and Old World groups within the Emballonuridae and Noctilionidae appear to be the outcomes of a series of cross-continental dispersal events from Asia rather than direct cross-ocean dispersal via vegetation rafts or cross-ocean stepping stones (Eick et al., 2005; Teeling et al., 2005) (Figs 4a,b and S4, and Tables S6 and S7). New World emballonurids descended from an Old World group that dispersed into North America via the North Atlantic or Beringia land bridges during the Late Eocene, which coincided with the decline in sea levels at high latitudes during that period (Miller et al., 2005) (Figs 4a,b and S4). This biogeographic scenario is consistent with the chronological patterns in the palaeontological record of emballonurids (Czaplewski et al., 2008; Gunnell et al., 2008). The AARs also suggest that the disjunct geographic distribution of the Noctilionidae, where the sister group of the Mystacinidae is restricted to New Zealand, whereas the remaining lineages are confined mostly to the Neotropics (due to the unresolved phylogenetic relationship of the Myzopodidae, this family was excluded from the superfamily Noctilionidae in the present study), resulted from a continental extinction event and a series of dispersal events from Asia to South America via North America approximately 48 Ma (Fig. 4a,b; Tables 1, S6 and S7). The continental extinction scenario for the Mystacinidae is supported by the continental fossil record (Hand et al., 2001). The provenance of the remaining tropical noctilionoids is ambiguous (node 16 in Fig. 4), but the time frame of the dispersal events between Northern and Southern America during the Middle Eocene coincides with the expansion and contraction of megathermal plants that occurred from the Late Palaeocene to the Middle Eocene (Morley, 2007), as well as the uplift of the Panama isthmus during the Middle Eocene (Scotese, 2001). These AARs provided new insight into the biogeographic history of bats. Given the abundant dispersal events in the AARs, cross-continental dispersal appears to have played a major role in shaping the current panglobal distribution of bats. The endemic and disjunct distribution pattern of extant bats may be attributable to their high dispersability, local extinction events and the influence of geological events, such as tectonic movements, as well as the expansion and contraction of megathermal rainforests during the Tertiary.

Acknowledgments We are particularly grateful to Prof. Michel Laurin, Dr. Gavin Thomas and many anonymous reviewers for their valuable comments on the early version of the manuscript. We also thank Junxiao Xu for helpful guidance of writing Python script, as well as Prof. Emma Teeling, Prof. Gabor Csobra, Zuo Chen and Tong Shen

13

for the helpful discussions. This study is supported by the National Science Fund for Distinguished Young Scholars (grant no. 31325025) to G.Y, National Natural Science Foundation of China (31300314, 30870116, 30670277, 31172045, and 31110103910), the Program for New Century Excellent Talents in University, the Ministry of Education of China (NCET-07-0445), the Research Innovation Program for College Graduates of Jiangsu Province (CXZZ11_0883) and the Priority Academic Program Development of Jiangsu Higher Education Institutions to G.Y.

References Alfaro, M.E., Santini, F., Brock, C., Alamillo, H., Dornburg, A., Rabosky, D.L. et al. 2009. Nine exceptional radiations plus high turnover explain species diversity in jawed vertebrates. Proc. Natl. Acad. Sci. USA 106: 13410–13414. Altringham, J. 1996. Bats Biology and Behaviour. Oxford University Press, New York. Buckley, L.B., Davies, T.J., Ackerly, D.D., Kraft, N.J.B., Harrison, S.P., Anacker, B.L. et al. 2010. Phylogeny, niche conservatism and the latitudinal diversity gradient in mammals. Philos. Trans. R. Soc. Lond. B Biol. Sci. 277: 2131–2138. Clayton, J.W., Soltis, P.S. & Soltis, D.E. 2009. Recent long-distance dispersal overshadows ancient biogeographical patterns in a pantropical angiosperm family (Simaroubaceae, Sapindales). Syst Biol. 58: 395–410. Cracraft, J. 1973. Vertebrate Evolution and Biogeography in the Old World Tropics: implications of Continental Drift and Palaeoclimatology. Academic Press, London. Currano, E.D. 2009. Patchiness and long-term change in early Eocene insect feeding damage. Paleobiology 35: 484–498. Cusimano, N., Stadler, T. & Renner, S.S. 2012. A new method for handling missing species in diversification analysis applicable to randomly or nonrandomly sampled phylogenies. Syst. Biol. 61: 785–792. Czaplewski, N.J., Morgan, G.S. & McLeod, S.A. 2008. Chiroptera. Cambridge University Press, Cambridge. Ducrocq, S., Buffetaut, E., Buffetaut-Tong, H., Jaeger, J.-J., Jongkanjanasoontorn, Y. & Suteethorn, V. 1992. First fossil flying lemur: a dermopteran from the late Eocene of Thailand. Palaeontology 35: 373–380. Ducrocq, S., Jaeger, J.J. & Sig
e, B. 1993. Un m
egachiroptere dans l’eocene sup
erieur de Tha€ılande—incidence dans la discussion phylog
enique du groupe [A megabat in the Late Eocene of Thailand: phylogenetic implications in the discussion of group]. Neues. Jahrb. Geol. Pal€aontol. Monatschr. 9: 561–575. Dumont, E.R., Davalos, L.M., Goldberg, A., Santana, S.E., Rex, K. & Voigt, C.C. 2012. Morphological innovation, diversification and invasion of a new adaptive zone. Philos. Trans. R. Soc. Lond. B Biol. Sci. 279: 1797–1805. Eick, G.N., Jacobs, D.S. & Matthee, C.A. 2005. A nuclear DNA phylogenetic perspective on the evolution of echolocation and historical biogeography of extant bats (Chiroptera). Mol. Biol. Evol. 22: 1869–1886. Eiting, T. & Gunnell, G. 2009. Global completeness of the bat fossil record. J. Mamm. Evol. 16: 151–173. Gunnell, G. & Simmons, N. 2005. Fossil evidence and the origin of bats. J. Mamm. Evol. 12: 209–246.

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

14

W. YU ET AL.

Gunnell, G.F., Simons, E.L. & Seiffert, E.R. 2008. New bats (Mammalia: Chiroptera) from the late Eocene and early Oligocene, Fayum Depression, Egypt. J. Vertebr. Paleontol. 28: 1–11. Hand, S.L., Archer, M. & Godthelp, H. 2001. New Miocene lcarops material (Microchiroptera: Mystacinidae) from Australia, with a revised diagnosis of the genus. Mem. Assoc. Australas. Palaeontol. 25: 139–146. Harmon, L.J., Weir, J.T., Brock, C.D., Glor, R.E. & Challenger, W. 2008. GEIGER: investigating evolutionary radiations. Bioinformatics 24: 129–131. Jaramillo, C., Rueda, M.J. & Mora, G. 2006. Cenozoic plant diversity in the neotropics. Science 311: 1893–1896. Jones, K.E., Bininda-Emonds, O.R.P. & Gittleman, J.L. 2005. Bats, clocks, and rocks: diversification patterns in Chiroptera. Evolution 59: 2243–2255. Katz, M.E., Miller, K.G., Wright, J.D., Wade, B.S., Browning, J.V., Cramer, B.S. et al. 2008. Stepwise transition from the Eocene greenhouse to the Oligocene icehouse. Nat. Geosci. 1: 329–334. Keane, T., Creevey, C., Pentony, M., Naughton, T. & Mclnerney, J. 2006. Assessment of methods for amino acid matrix selection and their use on empirical data shows that ad hoc assumptions for choice of matrix are not justified. BMC Evol. Biol. 6: 29. Labandeira, C.C. 2005. The fossil record of insect extinction: new approaches and future directions. Am. Entomol. 51: 14–29. Marjanovi
c, D. & Laurin, M. 2007. Fossils, molecules, divergence times, and the origin of lissamphibians. Syst. Biol. 56: 369–388. Mayhew, P.J., Jenkins, G.B. & Benton, T.G. 2008. A long-term association between global temperature and biodiversity, origination and extinction in the fossil record. Philos. Trans. R. Soc. Lond. B Biol. Sci. 275: 47–53. McKenna, D.D. & Farrell, B.D. 2006. Tropical forests are both evolutionary cradles and museums of leaf beetle diversity. Proc. Natl. Acad. Sci. USA 103: 10947–10951. Meredith, R.W., Janecka, J.E., Gatesy, J., Ryder, O.A., Fisher, C.A., Teeling, E.C. et al. 2011. Impacts of the cretaceous terrestrial revolution and KPg extinction on mammal diversification. Science 334: 521–524. Miller, K.G., Kominz, M.A., Browning, J.V., Wright, J.D., Mountain, G.S., Katz, M.E. et al. 2005. The phanerozoic record of global sea-Level change. Science 310: 1293–1298. Miller-Butterworth, C.M., Murphy, W.J., O’Brien, S.J., Jacobs, D.S., Springer, M.S. & Teeling, E.C. 2007. A family matter: conclusive resolution of the taxonomic position of the long-fingered bats, Miniopterus. Mol. Biol. Evol. 24: 1553– 1561. Monteiro, L. & Nogueira, M. 2011. Evolutionary patterns and processes in the radiation of phyllostomid bats. BMC Evol. Biol. 11: 137. Morley, R.J. 2003. Interplate dispersal paths for megathermal angiosperms. Perspect. Plant Ecol. Syst. 6: 5–20. Morley, R.J. 2007. Creaceous and Tertiary climate changes and the past distribution of megathermal rainforests. In: Tropical Rainforest Response to Climatic Changes (M.B. Bush & J.R. Flenley, eds), pp. 1–31. Praxis Publishing, Chichester. Parker, D.I., Lawhead, B.E. & Cook, J.A. 1997. Distributional limits of bats in Alaska. Arct. Alp. Res. 50: 256–265. Pybus, O.G. & Harvey, P.H. 2000. Testing macro–evolutionary models using incomplete molecular phylogenies. Philos. Trans. R. Soc. Lond. B Biol. Sci. 267: 2267–2272.

Rabosky, D.L. 2006a. LASER: a maximum likelihood toolkit for detecting temporal shifts in diversification rates from molecular phylogenies. Evol. Bioinform. Online 2006: 247–250. Rabosky, D.L. 2006b. Likelihood methods for detecting temporal shifts in diversification rates. Evolution 60: 1152–1164. Ravel, A., Marivaux, L., Tabuce, R., Adaci, M., Mahboubi, M., Mebrouk, F. et al. 2011. The oldest African bat from the early Eocene of El Kohol (Algeria). Naturwissenschaften 98: 397–405. Ree, R.H. & 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., Moore, B.R., Webb, C.O. & Donoghue, M.J. 2005. A likelihood framework for inferring the evolution of geographic range on phylogenetic trees. Evolution 59: 2299–2311.
Ferrero, V. & Navarro, L. 2011. When did Rojas, D., Vale, A., plants become important to leaf-nosed bats? Diversification of feeding habits in the family Phyllostomidae. Mol. Ecol. 20: 2217–2228. Ronquist, F. & Huelsenbeck, J.P. 2003. MrBayes 3: bayesian phylogenetic inference under mixed models. Bioinformatics 19: 1572–1574. Sanderson, M.J. 2003. r8s: inferring absolute rates of molecular evolution and divergence times in the absence of a molecular clock. Bioinformatics 19: 301–302. Schluter, D. 2000. The Ecology of Adaptive Radiation. Oxford University Press, Oxford, UK. Scotese, C.R. 2001. Atlas of earth history. PALEOMAP project. Arlington, TX. http://www.scotese.com. Shapiro, B., Rambaut, A. & Drummond, A.J. 2006. Choosing appropriate substitution models for the phylogenetic analysis of protein-coding sequences. Mol. Biol. Evol. 23: 7–9. Simmons, N.B. 2005. Order Chiroptera. In: Mammal Species of the World: A Taxonomic and Geographic Reference, 3rd edn (D.E. Wilson & D.M. Reeder, eds), pp. 312–529. The Johns Hopkins University Press, Baltimore, MD. Simmons, N.B. & Conway, T.M. 2003. Evolution of ecological diversity in bats. In: Bat Ecology (T.H. Kunz & M.B. Fenton, eds), pp. 493–535. University of Chicago Press, Chicago. Simmons, N.B. & Geisler, J.H. 1998. Phylogenetic relationships of Icaronycteris, Archaeonycteris, Hassianycteris, and Palaeochiropteryx to extant bat lineages, with comments on the evolution of echolocation and foraging strategies in Microchiroptera. Bull. Am. Mus. Nat. Hist. 235: 1–182. Simmons, N.B., Seymour, K.L., Habersetzer, J. & Gunnell, G.F. 2008. Primitive Early Eocene bat from Wyoming and the evolution of flight and echolocation. Nature 451: 818–821. Simpson, G.G. 1949. Tempo and Mode in Evolution. Columbia University Press, New York. Simpson, G.G. 1953. The Major Features of Evolution. Columbia University Press, New York. Smith, T., Rana, R.S., Missiaen, P., Rose, K.D., Sahni, A., Singh, H. et al. 2007. High bat (Chiroptera) diversity in the early Eocene of India. Naturwissenschaften 94: 1003–1009. Springer, M.S., Meredith, R.W., Janecka, J.E. & Murphy, W.J. 2011. The historical biogeography of Mammalia. Philos. Trans. R. Soc. Lond. B Biol. Sci. 366: 2478–2502. Stadler, T. 2010. TreeSim in R-Simulating trees under the birth-death model. http://cran.r-project.org/web/packages/ TreeSim/index.html. Stadler, T. 2011a. Mammalian phylogeny reveals recent diversification rate shifts. Proc. Natl. Acad. Sci. USA 108: 6187–6192.

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY

Diversification trend and Asian origin of bat

Stadler, T. 2011b. Simulating trees with a fixed number of extant species. Syst. Biol. 60: 676–684. Stevens, R.D. 2004. Untangling latitudinal richness gradients at higher taxonomic levels: familial perspectives on the diversity of New World bat communities. J. Biogeogr. 31: 665–674. Stevens, R.D. 2006. Historical processes enhance patterns of diversity along latitudinal gradients. Philos. Trans. R. Soc. Lond. B Biol. Sci. 273: 2283–2289. Teeling, E.C. 2009. Bats (Chiroptera). In: The Timetree of Life (S.B. Hedges & S. Kumar, eds), pp. 449–503. Oxford University Press, New York. Teeling, E.C., Springer, M.S., Madsen, O., Bates, P., O’Brien S, J. & Murphy, W.J. 2005. A molecular phylogeny for bats illuminates biogeography and the fossil record. Science 307: 580–584. Van den Bussche, R.A. & Hoofer, S.R. 2004. Phylogenetic relationships among recent chiropteran families and the importance of choosing appropriate out-group taxa. J. Mammal. 85: 321–330. Van Valen, L. 1979. The evolution of bats. Evol. Theory 4: 104– 121. Wilf, P. & Labandeira, C.C. 1999. Response of plant-insect associations to Paleocene-Eocene warming. Science 284: 2153–2156. Wilf, P., Labandeira, C.C., Johnson, K.R. & Ellis, B. 2006. Decoupled plant and insect diversity after the end-Cretaceous extinction. Science 313: 1112–1115. Wing, S.L. & Sues, H.-D. 1992. Terrestrial Ecosystems Through Time. University of Chicago Press, Chicago. Yoder, J.B., Clancey, E., Des Roches, S., Eastman, J.M., Gentry, L., Godsoe, W. et al. 2010. Ecological opportunity and the origin of adaptive radiations. J. Evol. Biol. 23: 1581–1596. Yu, W., Xu, J., Wu, Y. & Yang, G. 2012. A comparative study of mammalian diversification pattern. Int. J. Biol. Sci. 8: 486–497. Zachos, J., Pagani, M., Sloan, L., Thomas, E. & Billups, K. 2001. Trends, rhythms, and aberrations in global climate 65 Ma to present. Science 292: 686–693.

Supporting information Additional Supporting Information may be found in the online version of this article: Data S1 Additional details related to each procedure described in the Materials and Methods and the Result of temporal pattern of diversification.

15

Table S1 Total number of species and genera in each family, and the number (percentage) of samples used in this study. Table S2 GenBank accession numbers of taxa included in this study. Table S3 Fossils used for calibrating the molecular clock and related references. Table S4 Diversity and distribution of each genus used in this study. Table S5 Maximum likelihood analyses of the diversification rate based on the respective models and the complete CorSiM data sets. Table S6 Results of the biogeographic models analysed using the inferred chronogram, including likelihood scores (lnL), estimates of dispersal (D) and extinction (E) rates (events per Myr) under four different geographic range evolution models. Table S7 Results of biogeographic models analysed using constrained chronogram, including likelihood scores (lnL), estimates of dispersal (D) and extinction (E) rates (events per Myr) under four different models of geographic range evolution. Figure S1 Genus-level chronograms for Chiroptera, which were obtained using Bayesian and relaxed molecular clock approaches (inferred chronogram). Figure S2 Probability of rejecting the null hypothesis of a constant diversification rate when the birth rates shifted. Figure S3 Pie chart summarizing the statistics for the 2Dv2 values of the chronograms derived from the sensitivity analyses. Figure S4 Ancestral area reconstructions (AARs) of extant chiropteran genera under model TV2 with stratified dispersal probabilities between areas based on the inferred chronogram. Figure S5 Semi-logarithmic lineage-through-time (LTT) plots under different assumptions about the root age of Chiroptera. Data deposited at Dryad: doi:10.5061/dryad.r717q Received 23 December 2013; revised 9 August 2014; accepted 11 August 2014

ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY. J. EVOL. BIOL. doi: 10.1111/jeb.12477 JOURNAL OF EVOLUTIONARY BIOLOGY ª 2014 EUROPEAN SOCIETY FOR EVOLUTIONARY BIOLOGY