Bats (Chiroptera: Noctilionoidea) Challenge a Recent Origin of Extant Neotropical Diversity

Abstract

The mechanisms underlying the high extant biodiversity in the Neotropics have been controversial since the 19th century. Support for the influence of period-specific changes on diversification often rests on detecting more speciation events during a particular period. The timing of speciation events may reflect the influence of incomplete taxon sampling, protracted speciation, and null processes of lineage accumulation. Here we assess the influence of these factors on the timing of speciation with new multilocus data for New World noctilionoid bats (Chiroptera: Noctilionoidea). Biogeographic analyses revealed the importance of the Neotropics in noctilionoid diversification, and the critical role of dispersal. We detected no shift in speciation rate associated with the Quaternary or pre-Quaternary periods, and instead found an increase in speciation linked to the evolution of the subfamily Stenodermatinae (∼18 Ma). Simulations modeling constant speciation and extinction rates for the phylogeny systematically showed more speciation events in the Quaternary. Since recording more divergence events in the Quaternary can result from lineage accumulation, the age of extant sister species cannot be interpreted as supporting higher speciation rates during this period. Instead, analyzing the factors that influence speciation requires modeling lineage-specific traits and environmental, spatial, and ecological drivers of speciation.

Fossils, geographic range evolution, Neotropics, Phyllostomidae, speciation

The Neotropics harbors the most diverse terrestrial ecosystems on Earth, and the mechanisms behind this high biodiversity have remained contentious since the 19th century ( Wallace 1853 ; Bates 1863 ). Quaternary climate change (last 2.58 myr) has been invoked for several decades to explain how populations in lowland New World forests became genetically isolated, generating many descendant species ( Haffer 1969 ; Hooghiemstra and van der Hammen 1998 ). Geomorphological ( Mayle 2004 ), paleoecological ( Porto et al. 2013 ), and biogeographic data ( Poelchau and Hamrick 2013 ), however, suggest other mechanisms explain speciation in the region. Vicariance events resulting from the tectonic activity and habitat change during the pre-Quaternary (mainly the Neogene, between 2.58 Ma and 23.03 Ma) have been proposed as an alternative mechanism to explain neotropical biodiversity ( Bush 1994 ; Hoorn et al. 2010 ). The Quaternary and pre-Quaternary hypotheses have been treated as mutually exclusive, despite a variety of other possible explanations (e.g., the gradient hypothesis, the riverine barrier hypothesis, or the disturbance–vicariance hypothesis; see references in Rull 2015 ). As a result, although a multiple-hypothesis approach has been recommended to study the origin of neotropical diversity ( Rull 2015 ), most of the debate has focused on the evidence supporting either temporal model.

The age of crown groups and divergence times between sister species are the main evidence used to support these temporal models, and have been used as proxies for speciation rates during the Neogene or the Quaternary. For example, analyses of divergence time for 131 pairs of sister species of neotropical butterflies were interpreted as supporting the Quaternary origin of species diversity in that neotropical clade ( Garzón-Orduña et al. 2014 ). In contrast, an earlier meta-analysis of crown group ages from 105 neotropical time-calibrated phylogenies found that cladogenesis unfolded more or less continuously in the Neotropics ( Rull 2008 ), and neither the Quaternary nor the preceding periods were overwhelmingly more important in speciation (see also Zink et al. 2004 ). Another meta-analysis using crown ages of different clades from northern South America concluded that most extant lineages in the region have a pre-Quaternary origin ( Hoorn et al. 2010 ).

Even as the debate on the influence on speciation of geological events during particular periods has been reignited, it has become clear that incomplete taxon sampling, protracted speciation ( Etienne and Rosindell 2012 ), and null processes of lineage accumulation strongly influence divergence dates. Crown age and species divergences will also reflect those processes, in addition to any shifts in speciation rates. Testing the effects of environmental, spatial, and ecological drivers of clade diversity then requires explicit modeling of speciation rates ( Rabosky 2014 ).

Noctilionoid bats (superfamily Noctilionoidea), and in particular fruit-feeding phyllostomids, are among the most common and diverse lowland mammals in the Neotropics ( Fleming and Kress 2013 ). The New World noctilionoid families Thyropteridae, Furipteridae, Noctilionidae, Mormoopidae, and Phyllostomidae account for ∼20% of global species richness among bats ( Teeling et al. 2005 ). The past and current distributions of two other noctilionoid families—Myzopodidae from Madagascar, with fossils in Africa, and Mystacinidae from New Zealand, with fossils in Australia—suggest a Gondwanan origin for the superfamily and open new questions on when and how ancestral noctilionoids reached the Americas ( Gunnell et al. 2014 ). Many species of extant New World noctilionoids occur in northern South America (Supplementary Fig. S1; available on Dryad at http://dx.doi.org/10.5061/dryad.s533p ), and biogeographic studies in a few genera of phyllostomids have found a South American origin for most species, followed by dispersals to Central America and the Caribbean ( Velazco and Patterson 2008 , 2013 ). This suggests South America has been an important area for diversification of New World noctilionoids, if not their area of origin. These characteristics make noctilionoids a good study system to analyze the timing of neotropical taxonomic diversification.

Despite their mostly neotropical distribution, there are fossil and extant noctilionoids outside the Neotropics that suggest a long history in other areas and make it critical to estimate their biogeographic history. Early Oligocene fossils (30.8–33.3 Ma) trace the most ancient noctilionoid lineages to North America ( Czaplewski and Morgan 2012 ; Morgan and Czaplewski 2012 ). These records suggest an ancient North American distribution for the two most species-rich noctilionoid families: Mormoopidae and Phyllostomidae ( Baker et al. 2012 ; Arita et al. 2014 ). To date, however, formal biogeographic analyses have not been conducted and the mode and timing of events leading to the present overwhelmingly neotropical distributions of noctilionoids remain to be analyzed.

Despite much recent progress, phylogenetic relationships within the superfamily remain unclear, in particular the position of Thyropteridae relative to other neotropical families ( Eick et al. 2005 ; Miller-Butterworth et al. 2007 ; Meredith et al. 2011 ; Teeling et al. 2012 ). Within-family relationships in Phyllostomidae—one of the most ecologically diverse family of mammals—have not been estimated with all available lineages. Instead, evolutionary studies at different taxonomic levels and on independent clades have been conducted with mostly mitochondrial sequences ( Dávalos and Jansa 2004 ; Dávalos 2007 ; Porter et al. 2007 ; Hoffmann et al. 2008 ; Redondo et al. 2008 ; Rojas et al. 2011 , 2013 ; Dumont et al. 2012 ; Velazco and Patterson 2013 ; Velazco and Lim 2014 ), or encompassing different subfamilies or genera to represent all species ( Baker et al. 2000 , 2003 ; Datzmann et al. 2010 ).

Here, we infer the dated phylogeny of New World Noctilionoidea with the goal of inferring both shifts in speciation and biogeographic history for this clade. We analyzed the greatest number of extant lineages, loci, and fossil-based node calibrations to date. The resulting phylogenies were used to infer the ancestral distribution of key clades in the superfamily, and estimate the roles of dispersal and vicariance in the diversification of this clade ( Matzke 2013 ). Finally, we analyzed the tempo and mode of taxonomic diversification, and compared those results to analyses of the timing of divergence between extant sister species. These analyses revealed the inadequacy of divergence times in determining the role of climate and habitat change at particular periods on the speciation process.

Materials andMethods

Taxonomic and Molecular Sampling

Generic and species-level taxonomy for Noctilionoidea follows Simmons (2005) , updated following the recent taxonomic literature ( Hurtado and Pacheco 2014 ; Parlos et al. 2014 ; Tavares et al. 2014 ; Velazco and Patterson 2014 ). Subfamily and generic taxonomy within Phyllostomidae follows Baker et al. (2003) . We sampled 90% of the 70 genera and 76% of the 254 species described hitherto for the New World Noctilionoidea (Supplementary Table S1). The sample included the 11 subfamilies of Phyllostomidae and all neotropical families of the superfamily Noctilionoidea. In addition, we used the Old World noctilionoids Mystacina tuberculata and Myzopoda aurita as outgroups. This data set included sequences from one undescribed species of the genus Micronycteris ( Porter et al. 2007 ), two newly described species, and one undescribed species in the genus Sturnira ( Velazco and Patterson 2014 ), as well as lineages from populations previously assigned to Lophostoma silvicolum and Pteronotus parnellii rubiginosus that likely represent independently evolving populations ( Porter et al. 2003 ; Clare et al. 2013 ).

We generated sequences for seven nuclear loci (Supplementary Table S2; 661 sequences, 141 taxa): introns of thyrotropin beta chain ( thy ) and signal transducer and activator 5A ( stat5a ), autosomal exons of brain-derived neurotrophic factor ( bdnf ), titin 6 ( ttn6 ) and recombination-activating protein 2 ( rag2 ), X-chromosome exon ATPase-7A ( atp7a ), and the 3′-untranslated region or phospholipase C beta 4 ( plcb4 ) ( Dávalos et al. 2014 ). To ensure the greatest possible coverage for the 194 taxa, we also downloaded sequences from GenBank for all these loci, and for the mitochondrial genes cytochrome b ( cyt-b ), cytochrome oxidase I ( coxI ) and the ribosomal RNAs 12S, tRNA ^val , and 16S (Supplementary Table S1).

We used the Perl script transAlign v1.2 ( Bininda-Emonds 2005 ) to align all nuclear exons and protein-coding mitochondrial genes. MAFFT v.6.710b ( Katoh et al. 2005 ; Katoh and Toh 2008 ) was used to align the non-coding nuclear sequences (protocol E-INS-i) and the mitochondrial ribosomal RNA sequences (protocol Q-INS-i, which accounts for RNA secondary structure). The concatenated alignment comprised 9587 base pairs. We used PartitionFinder v.1.0.0 ( Lanfear et al. 2012 ) to select the optimal partitioning scheme and model of DNA evolution for individual partitions. The analysis was performed with the greedy algorithm applying the Bayesian Information Criterion to make the search for optimal partitioning schemes as efficient as possible. We applied the PartitionFinder results (Supplementary Table S3) in subsequent maximum likelihood and Bayesian phylogenetic analyses.

Phylogenetic Estimation

Maximum likelihood (ML) was used to infer an initial tree with 100 rapid bootstrap replicates on the optimal partition scheme using RAxML v.7.2.7 ( Stamatakis 2006 ). We then estimated preliminary node ages for the ML phylogeny using a semi-parametric method based on penalized likelihood ( Sanderson 2002 ), as implemented in the R package ape ( Paradis et al. 2004 ). Eight nodes were constrained using as minimum and maximum ages those of the chronostratigraphic epochs in which fossils of New World noctilionoids were found (see below). This time-calibrated phylogeny was then used as a starting point for Bayesian analyses. Two independent Metropolis-coupled Markov chain Monte Carlo (MCMCMC) searches ran for 40 million generations sampling every 4000 steps in BEAST v.1.7.5 ( Drummond et al. 2012 ). To calibrate these phylogenies, we used eight fossil calibration constraints distributed throughout the tree, with hard boundaries set to constrain the minimal age of the stratum where the fossil was found and prior distributions shown in Table 1 (see also Fig. 1 ). Prior constraints for the fossils followed a normal distribution while the hard boundaries were set using uniform priors. We set the Yule process as the tree branching prior, and per-branch rates of molecular evolution were drawn from a log-normal distribution ( Drummond et al. 2006 ). We repeated the analysis using a birth–death prior and the results were almost identical (median likelihood _Yule  = −161,971.48, median likelihood _{birth–death}  = −161,972.46). Here we report the results of the analysis using the simpler Yule prior. We also performed a Bayesian inference of the phylogeny of noctilionids excluding the third codon position of cyt-b . Except for the position of Lonchorhina , a perennially problematic genus in the phylogeny of Phyllostomidae, the topology and the divergence times were consistent with the original analysis (Supplementary Table S4). In BEAST analyses, convergence of the chains to the stationary distribution was confirmed using Tracer v.1.5. We discarded 10% of the trees as burn-in and remaining trees were combined. In the original BEAST analysis, estimates of autocorrelation time between posterior samples were used to subsample the set of posterior phylogenies. Finally, we obtained an uncorrelated posterior sample of 871 time-calibrated trees. ML and Bayesian analyses were performed in the CIPRES Science Gateway ( Miller et al. 2010 ).

Figure 1.

ML phylogeny of New World Noctilionoidea (log-likelihood = 157,714). Symbols at the nodes represent bootstrap values. Numbers in square brackets identify the fossil-based node calibrations described in Table 1 . Outgroups appear in gray font.

Table 1.

Fossil calibration points used to estimate divergence times in New World Noctilionoidea

Fossil	Fossil age	Lower boundary	Mean (SD) [95% CI]	Location in Figure 1	References
Stem mormoopid fossil	Oligocene (Whitneyan)	30.8	31.0 (3.0) [30.9–37.66]	[1]	Czaplewski and Morgan (2003)
Notonycteris	Miocene (Laventan)	11.8	12.2 (3.0) [11.9–18.81]	[3]	Czaplewski et al. (2003b) ; Dávalos et al. (2014)
Palynephyllum antimaster	Miocene (Laventan)	11.8	12.2 (3.0) [11.9–18.81]	[6]	Czaplewski et al. (2003b) ; Dávalos et al. (2014)
Desmodus archaeodaptes	Pleistocene (Irvingtonian)	0.3	1.58 (5.0) [0.494–12.43]	[2]	Silva (1979)
Phyllonycteris major	Quaternary	0.001	0.905 (5.0) [0.183–11.84]	[5]	Choate and Birney (1968)
Cubanycteris silvai	Quaternary	0.001	0.905 (5.0) [0.183–11.84]	[7]	Czaplewski et al. (2003a) ; Mancina and García-Rivera (2005)
Brachyphylla nana	Quaternary	0.001	0.905 (4.0) [0.149–8.682]	[4]	Teeling et al. (2005) ; Dávalos (2009)
Phyllops vetus	Quaternary	0.001	0.905 (1.0) [0.076–2.95]	[8]	Teeling et al. (2005) ; Dávalos (2007)

Notes: Low boundaries were set to constrain the minimal age of the stratum where the fossil was found. Lower boundary, 95% confidence interval (95% CI), mean and standard deviation (SD) of the prior constraints are in millions of years ago (Ma).

We used a Bayesian Concordance Approach (BCA) implemented in BUCKy v. 1.4.0 to account for discordance between gene trees ( Ané et al. 2007 ; Larget et al. 2010 ). Briefly, the distribution of trees from the same locus obtained using Bayesian inference is summarized. Next, these files from different loci are used to estimate the primary concordance tree. This new tree comprises all clades with concordance factors (CF) >50%, and other clades that are not excluded by concordance among tree samples ( Baum 2007 ). The CF of a clade is the proportion of loci that have the clade. Unlike other approaches that account for discordance in the phylogenies inferred from different genes, BCA does not model a particular cause of conflict among gene trees.

To estimate tree CFs, we performed two independent MCMCMC searches in MrBayes v3.2.4 ( Ronquist et al. 2012 ) for each locus. Each analysis was allowed to run for 10 million generations sampling from the chain every 1000 generations. Convergence was assessed using Tracer v.1.5 (burn-in = 25% of the samples), and by examining the average standard deviation of split frequencies among parallel chains. Results of the Bayesian analyses were combined and summarized. We summarized the posterior distributions of tree topologies generated by MrBayes for each gene and then estimated a primary concordance tree in BUCKy. We also evaluated the effect of choosing four prior values (0.1, 1, 10, and 100) for the $\alpha$ parameter (i.e., the a priori level of discordance among loci). For each value of $\alpha$ , we ran two independent sets of MCMCMC. Each run used 1 cold chain and 3 heated chains, 100 million generations, and a 10% burn-in. The primary concordance trees produced using different $\alpha$ values were identical. We provide results for $\alpha =0.1$ . In this tree, some nodes show poor resolution (i.e., CF <0.5) (Supplementary Fig. S2).

To elucidate whether this pattern reflects actual incongruence between the loci, or weak phylogenetic signal in one or more of them, we conducted tree distance analyses. If the genes have weak signal, the distances of some single-gene trees to the primary concordance tree will overlap with distances from the random trees. If each locus has a strong signal and incongruence is small, then the distances from the concatenated tree to the primary concordance tree will be lower than the distances from the single-gene trees. If each locus has a strong signal and incongruence is high, then the distances from some of the single-gene trees to the primary concordance tree will be lower than the distances from the concatenated tree, while trees from other loci will show much greater distance. First, we calculated the topological distance between the primary concordance tree and the trees inferred using each locus, as well as the concatenated trees using the dist.topo function of the ape R library. This distance captures single-gene signal because we used the posterior distribution of topologies that we previously inferred in MrBayes for each locus. Then, we simulated 15,000 trees using the function rmtree in ape and calculated the topological distances between these and the primary concordance tree. By default, BUCKy prunes all taxa that are missing from one or more loci. As a result, the primary concordance tree had a reduced number of species compared to the maximum credibility tree inferred in BEAST. The taxa in the primary concordance tree, however, represent a random sample from the phylogeny of noctilionoids and include representatives from all families and most subfamilies (Supplementary Fig. S2).

Biogeographic Analyses

We defined three biogeographic regions: Central and North America (north of the Isthmus of Panama), South America (south of the Isthmus of Panama), and the Antilles ( sensuGenoways et al. 1998 , 2010 , including Grenada and the Grenadines). This was based on evidence supporting these three regions as centers of diversification of New World noctilionoids ( Dávalos 2006 , 2007 ; Arita et al. 2014 ). Taxa were assigned to one or more of these regions based on their current distribution compiled by the IUCN Schipper et al. 2008 ; data downloaded on October 2013), and the literature since 2008 (Supplementary Table S5).

We estimated the fit of models of geographic range evolution to the phylogenies and distribution of New World Noctilionoidea using BioGeoBEARS ( Matzke 2013 ). This R package uses model selection criteria to compare likelihood implementations of different models of geographic range evolution: Dispersal–Vicariance Analysis (DIVA) ( Ronquist 1997 ), Dispersal–Extinction–Cladogenesis (DEC) ( Ree 2005 ; Ree and Smith 2008 ), and a model that decouples range evolution from cladogenesis (BAYAREA) ( Landis et al. 2013 ). In these models, four types of events may occur at cladogenesis: (i) the ancestral range is copied to both species, (ii) a new species evolves in a subset of the ancestral range, (iii) vicariance, and (iv) a new species jumps to a range outside of the ancestral range ( Matzke 2013 , 2014 ). The latter is known as founder event speciation. When this parameter is added, three new models are generated: DIVA+J, DEC+J, and BAYAREA+J. We tested these six biogeography models and used the Akaike Information Criterion (AIC) and Akaike weights to interpret how well the different models fit the geographic distributions and the maximum clade credibility (MCC) phylogeny. We obtained the overall probabilities of the anagenetic and cladogenetic events conditional on the geographic distributions, the phylogeny, and the best-fitting model from 1000 stochastic biogeographic mappings. This allowed us to quantify the relative role of dispersal and vicariance at cladogenesis in the diversification of noctilionoids.

Analysis of Diversification Rates and Speciation Events

If the Quaternary played a central role in the diversification of noctilionoids, then the speciation rate should increase during this period relative to previous time intervals, and more speciation events should have taken place in the Quaternary than expected given a null process of constant speciation and extinction. To test these predictions, first we visualized the lineage-through time (LTT) plot (using code modified from the ltt function of the phytools R library, Revell 2012 ) as the cumulative number of lineages in the time-calibrated phylogeny on a logarithmic scale graphed against time when lineages arise ( Nee et al. 1992 ). An uptick of this curve in the Quaternary would provide support for the mechanisms explaining diversification during this period ( Moyle et al. 2009 ).

Second, we examined the tempo of speciation of noctilionoids using Bayesian Analysis of Macroevolutionary Mixtures (BAMM) ( Rabosky 2014 ). In support of a predominantly Quaternary origin for neotropical diversity, this analysis should detect either a general, phylogeny-wide trend toward higher diversification rates in the Quaternary, or a similar pattern among some subclades.

To account for incomplete taxon sampling, we specified sampling fractions by genus. We ran two independent BAMM analyses. Four reversible-jump MCMC ran for 11,000,000 generations and were sampled every 1000 generations. We discarded the first 1,000,000 generations as burn-in and analyzed the effective sample size to check for convergence. We calculated Bayes factors (BF) to compare alternative macroevolutionary models. A value of 3–12 for the BF represents positive evidence, a value greater than 12 is strong evidence, and BF >150 is very strong evidence ( Raftery 1996 ). To account for protracted speciation that can bias these analyses toward slowdowns in speciation closer to the present ( Etienne and Rosindell 2012 ), we removed the last 180,000 years from the phylogenies (during this interval no branching occurs) and repeated the BAMM analysis on the truncated tree.

Finally, we quantified the number and timing of speciation events between extant sister species for the Quaternary and pre-Quaternary periods on the posterior sample of 871 time-calibrated phylogenies (Supplementary R script). We used the birth–death model of cladogenesis to generate a distribution of expected sister-species divergences without the influence of any additional biotic or abiotic factors ( Nee et al. 1994 ; Morlon 2014 ). To evaluate whether the observed differences in number of events between the two time periods could arise through variation in this null process, we simulated phylogenies using birth–death models with a range of values for speciation (0.08, 0.14, 0.2, 0.24) and extinction (0, 0.01, 0.025, 0.05) rates using the pbtree function of phytools . These values lie within the observed distribution of rates that were estimated using BAMM and encompass great variation in parameters (i.e., three times the maximum value for speciation rates, and five times the minimum value for extinction rates). For each combination of parameters (i.e., 16 scenarios), we simulated 871 trees with the same number of tips and root age as the MCC phylogeny of noctilionoids. Then, we calculated the number of Quaternary and pre-Quaternary speciation events for each tree in each scenario. The empirical pattern of differences in number of speciation events between the Quaternary and pre-Quaternary time periods was then compared to the 16 simulated patterns using a randomization $t$ -test (9999 randomizations, $\alpha =0.05$ ). We repeated these calculations in a set of 100 scenarios to account for a broader range of variation in speciation (0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.45, 5, 0.75, 0.9) and extinction rates (0, 0.01, 0.02, 0.025, 0.035, 0.045, 0.05, 0.055, 0.075, 0.09) on the timing of speciation events.

Given that the phylogeny inferred in this study sampled 76% of noctilionoids, we also performed a sensitivity analysis to test for the effect of missing taxa on the timing of divergence events between extant sister species. We simulated 1000 time-calibrated trees using a birth–death model. The birth and death rates were taken from the mean values of speciation and extinction rates estimated in the BAMM analyses. The root age of the simulated trees was the same than that of the MCC tree, and the number of tips was equal to the 252 species described hitherto for New World Noctilionoidea. Then, we randomly pruned 5%, 10%, 15%, 20%, 25%, 30%, and 35% of the tips from each simulated tree. For each of the eight data sets, we quantified the number of divergence events between extant sister species for the Quaternary and the pre-Quaternary periods, to evaluate the impact of missing taxa on estimates of speciation events during these intervals.

Results

Phylogeny of New World Noctilionoidea

Tree distance analyses suggested the different loci had strong signal with relatively low incongruence, making the concatenation tree a good estimate of phylogeny (Tonini et al. 2015). Compared to random trees, single-gene trees were topologically very similar to the primary concordance tree estimated in BUCKy. Concatenation trees inferred with BEAST showed the lowest topological distance to the primary concordance tree (Supplementary Fig. S3), suggesting that the Bayesian inference of the phylogeny of noctilionoids was robust to discordance among single-gene trees.

This is the first species-level molecular phylogeny comprising all families of New World Noctilionoidea. Monophyly of this group, families represented by at least two species, and all subfamilies of Phyllostomidae received strong support (bootstrap support or $\text{BS}⩾98\%$ , posterior probability or $\text{PP}⩾0.97$ ) ( Fig. 1 ). The most recent common ancestor (MRCA) of New World noctilionoids was dated to the end of the Lutetian stage in the Eocene (42.1 Ma, 95% high probability density interval 37.6–47.0 Ma) ( Fig. 2 ). In contrast to previous studies, Thyropteridae was neither sister to Mystacinidae ( Eick et al. 2005 ), nor to Furipteridae and Noctilionidae ( Teeling et al. 2005 , 2012 ; Miller-Butterworth et al. 2007 ). Instead, it was corroborated as the sister group to all other ingroup families (BS = 98%, PP = 0.97) ( Meredith et al. 2011 ). Speciation among extant mormoopids (32.8, 30.8–36.2 Ma) and phyllostomids (30.3, 27.3–33.9 Ma) began in the Oligocene. Thyropteridae diversified during the Miocene (20.6, 15.7–25.7 Ma). The two species of Noctilionidae diverged at the beginning of the Pleistocene (2.6, 1.7–3.6 Ma), consistent with multi-allele coalescent-based analyses suggesting recent speciation in this genus ( Pavan 2013 ).

Figure 2.

Time-calibrated Bayesian phylogeny of New World Noctilionoidea. MCC tree summarizing 871 uncorrelated Bayesian posterior molecular phylogenies (log-likelihood = −161,972 ± 17). Color fill at the nodes depicts posterior probabilities (PP). Bars at the nodes represent the 95% highest posterior density (HPD) intervals of node ages. Colors of the epochs follow the International Chronostratigraphic Chart. Outgroups appear in gray font. The paraphyletic genus Lonchophylla appears in blue.

Within Phyllostomidae, our results corroborated Macrotinae ( Baker et al. 2003 ; Datzmann et al. 2010 ; Dumont et al. 2012 ), and rejected Desmodontinae ( Wetterer et al. 2000 ) as sister taxa to all remaining phyllostomids (BS = 100%, PP = 1). The position of Lonchorhininae was uncertain: it was neither sister to a clade including Phyllostominae ( Baker et al. 2003 ), nor sister to a clade excluding Phyllostominae ( Rojas et al. 2011 ; Dumont et al. 2012 ). Although morphological characters and cyt-b sequences support a clade of nectar-feeding phyllostomids ( Dávalos et al. 2012 , 2014 ; Dumont et al. 2012 ), Glossophaginae and Lonchophyllinae were not sister clades, corroborating results previously obtained with fewer taxa ( Baker et al. 2003 ; Datzmann et al. 2010 ; Rojas et al. 2011 ). Instead, Lonchophyllinae was sister to a clade comprising Carolliinae, Glyphonycterinae, Rhinophyllinae, and Stenodermatinae (BS = 93%, PP = 1) ( Baker et al. 2003 ; Datzmann et al. 2010 ).

Historical Biogeography

The DEC+J model received the highest support among the six biogeographic models ( Table 2 ). The most frequent modes of range inheritance in the biogeographic history of New World noctilionoids were range copying within the same large area with 49.2% of total events, and anagenetic dispersal with 29.4% of events ( Fig. 3 ). The extent of the three biogeographic regions defined in this study precludes interpreting these events as “sympatric” speciation, and we use “the same large area” to clarify this point. Range copying and range subset within the same large area represent 84.1% of the total number of events at cladogenesis, or 160 of the 191 nodes in the phylogeny. Founder events accounted for 9% and vicariance for 6.9% of cladogenetic events. Anagenetic dispersal was more frequent between the two continental regions, particularly from South America to Central and North America, than between the mainland and the Antilles. Dispersal at cladogenesis (i.e., founder event) was also more frequent from South America to Central and North America and, to a lesser extent, from South America to the Antilles, than in any other direction ( Fig. 4 , Supplementary Table S6).

Figure 3.

Frequency distributions of event counts from 1000 biogeographic stochastic mappings on the MCC phylogeny. Given the size of the continental size biogeographic areas analyzed, we use “within same large areas” instead of “sympatric” to describe cladogenetic events involving range copying and range subset.

Figure 4.

Biogeographic analysis of New World Noctilionoidea. Pie charts at nodes indicate marginal ML probabilities for the estimated ancestral areas, according to the DEC model with founder event speciation (DEC+J, all rates in events per million years): rate of anagenetic range expansion $(d)=0.025$ , rate of anagenetic range contraction $(e)=0.1\times {10}^{\text{-11}}$ , relative per-event weight of jump dispersal at cladogenesis $(j)=0.031$ , log-likelihood = −294. Genera are indicated with vertical (two or more species) or horizontal (monotypic genus) lines to the right of the current distribution of taxa. Clades discussed in the text are named. Fossils are added with bars indicating the geological horizon where the fossil was found. The ancestral areas of the clades to which fossils 1 and 2 belong were estimated as equivocal. This prevented us from using these two fossils in support of biogeographic analyses. Inset represents the mean values of dispersal events between any pair of the geographic areas. S: South America, C: Central and North America, A: Antilles. The thickness of the arrows is in correspondence with the values.

Table 2.

Models of ancestral range estimation of New World Noctilionoidea.

Model	ln $L$	$k$	$d$	$e$	$j$	AIC	$w$
DIVA	−310.8	2	0.033	0.1 × 10 −11	−	625.7	0.221 × 10 −11
DIVA+J	−303.8	3	0.028	0.1 × 10 −11	0.028	613.7	0.896 × 10 −4
DEC	−302.3	2	0.028	0.1 × 10 −11	−	608.6	0.001
DEC + J	−294.5	3	0.025	0.1 × 10 −11	0.031	595.0	0.999
BAYAREA	−371.6	2	0.014	0.038	−	747.2	0.899 × 10 −33
BAYAREA+J	−303.7	3	0.006	0.017	0.049	613.4	0.100 × 10 −3

Notes: Models include Dispersal–Extinction–Cladogenesis (DEC), Dispersal–Vicariance analysis (DIVA), and a model that assumes that no range evolution occurs at cladogenesis (BAYAREA), as well as these three models allowing for founder event speciation (+ J). ln $L$ : log-likelihood, $k$ : number of parameters fitted, $d$ : rate of anagenetic range expansion, $e$ : rate of anagenetic range contraction, $j$ : relative per-event weight of jump dispersal at cladogenesis, AIC: Akaike Information Criterion, $w$ : Akaike weights. Best-fit model is in bold .

The DEC+J model estimated the ancestral region of all New World noctilionoids in South America (marginal ML probability = 0.50) ( Fig. 4 ). This region also received high support (i.e., >0.5) in 120 of the 191 nodes in the phylogeny, making noctilionoids a suitable group to account for speciation timing in the Neotropics ( Fig. 4 ). A widespread distribution involving all three biogeographic areas was highly supported for the MRCA of Noctilionidae (0.85). The ancestral areas of Thyropteridae, the MRCA of Mormoopidae and Phyllostomidae, and Mormoopidae were equivocal. South America + Central and North America had the highest support for the ancestral region of the MRCA of Phyllostomidae (0.40), followed by Central and North America (0.23). The origin of most phyllostomid subfamilies was traced to South America, including Micronycterinae (0.65), Desmodontinae (0.64), Phyllostominae (0.93), Lonchorhininae (0.95), Glossophaginae (0.51), Lonchophyllinae (1.00), Glyphonycterinae (0.87), Rhinophyllinae (1.00), and Stenodermatinae (0.98). The ancestral area of the MRCA of Macrotinae was inferred to be Central and North America (0.60), while the MRCA of Carolliinae was inferred as widely distributed on the mainland (0.68 for South America + Central and North America). Within the phyllostomid subfamily Glossophaginae, the Antilles received high support (0.63) as the ancestral region of the MRCA of Brachyphylla , Erophylla , Phyllonycteris , Monophyllus , Leptonycteris , and Glossophaga . The MRCA of short-faced bats was inferred to include South America and possibly the Antilles (0.43 for South America, 0.33 for South America + Antilles).

The geographic distribution of six of the eight fossils we used to calibrate the phylogeny matches the estimated ancestral regions of the MRCA of the clades to which these fossils were assigned ( Fig. 4 ). This provides an independent line of evidence supporting the estimates of ancestral ranges from extant species.

Diversification Rates and Speciation Events

The LTT plot reveals that cladogenesis slowed down in the last 2.5 Ma in the history of noctilionoids ( Fig. 5 b), and this may suggest protracted speciation ( Rosindell et al. 2010 ). We found strong evidence for a two-regime macroevolutionary model compared to the zero rate shifts model (BF = 27.7, Supplementary Table S4). Other multiple-regime macroevolutionary models showed BFs higher than the two-regime model but were discarded because the evidence for each of those was less than positive relative to the two-regime model. Rate-through-time dynamics identified a decoupling in the rates of speciation for the subfamily Stenodermatinae (mean [90% confidence interval]: 0.227 [0.138–0.292] species/Ma) from the rates of speciation for other New World noctilionoids (0.157 [0.129–0.196] species/Ma) (Supplementary Fig. S5). The macroevolutionary cohort matrix analysis revealed that stenodermatines are grouped in a well-supported cohort for speciation, different from the rate of other noctilionoids (i.e., there is a high probability that stenodermatines share the same macroevolutionary regime) (Supplementary Fig. S6). Results were consistent after accounting for protracted speciation, and this suggests that uncertainty surrounding ongoing speciation does not influence the estimation of shifts in speciation rate in noctilionoids.

Figure 5.

Cladogenesis in the phylogeny of New World Noctilionoidea. a) Number of Quaternary and pre-Quaternary speciation events of extant species across the uncorrelated posterior sample of 871 time-calibrated trees of noctilionoids. b) LTT plot for the posterior sample of time-calibrated trees (gray). The LTT plot of the MCC tree is in black. The MCC tree is shown in the background along with estimated speciation rates along its branches. c) Number of Quaternary and pre-Quaternary speciation events of extant species across simulated distributions of trees under a birth–death model, using different values of speciation ( $b$ ) and extinction ( $d$ ) rates.

The median number (first–third quartiles) of speciation events between extant sister species in New World noctilionoids was 33 (31–35) in the pre-Quaternary and 31 (29–33) in the Quaternary ( Fig. 5 a). This is consistent with the observed slowdown in cladogenesis in the last 2.5 Ma ( Fig. 5 b). Estimated species-level divergence for the Quaternary (mean ± standard deviation) averaged 1.53 ± 0.66 Ma and ranged from 0.02 to 2.59 Ma. Estimated timing of speciation events between extant sister species in the pre-Quaternary averaged 6.66 ± 4.33 Ma and ranged from 2.59 to 25.72 Ma.

The number of speciation events in New World noctilionoids in the pre-Quaternary was on average significantly higher than in the Quaternary ( $t=15.57$ , $P<0.0001$ ) ( Fig. 5 a). The opposite pattern was observed in each of the scenarios of both the small set ( $n=16$ , Fig. 5 c), and the expanded set of simulations ( $n=100$ , Supplementary Fig. S7). The empirical pattern of differences in number of speciation events between the pre-Quaternary and Quaternary differed significantly from all 16 main simulated patterns (Supplementary Fig. S8). Sensitivity analyses show that the difference in number of speciation events between the Quaternary and pre-Quaternary decreases as the proportion of missing taxa increases. This follows from a disproportionate decrease in the number of speciation events in the Quaternary. Hence, the null model fails to generate a pattern similar to that observed, even when as many as 35% of the taxa are missing (Supplementary Fig. S9).

Discussion

Biogeographic analyses of the phylogeny of New World Noctilionoidea show that dispersal has been about five times more important than vicariance in the history of these bats. The timing of sister-species speciation in this diverse group of neotropical vertebrates was slightly higher in the pre-Quaternary than in the Quaternary period, in contrast with simulations based on a null model that consistently recover the opposite pattern. These simulations reveal that more divergence events between sister species in the Quaternary cannot be interpreted as evidence to test Pleistocene-related climate change as a mechanism driving diversification. This is because more divergence events in the Quaternary can result from lineage accumulation through a birth–death process. Further, the pattern of more speciation events in the Quaternary recurs even as the proportion of missing species increases. Diversification analyses support a clade-dependent shift in diversification rates at the base of a clade of mostly fruit-eating bats. Lineage-specific traits, then, have been more important than paleoclimatic and landscape changes in the diversification of New World Noctilionoidea. We review the systematics and biogeography of noctilionoids, before discussing in detail the general implications of our findings for elucidating speciation timing and testing different mechanisms of diversification on neotropical biota.

Evolution of New World Noctilionoid Families and Subfamilies

By estimating the phylogeny of New World Noctilionoidea, we provide the most detailed evolutionary history for one of the five major clades within Chiroptera. The new phylogenies inform systematic relationships among families and subfamilies of great ecological importance in the Neotropics. First, these trees help resolve a decade-long controversy over the position of Thyropteridae or disc-winged bats as sister to all other New World families of Noctilionoidea ( Eick et al. 2005 ; Teeling et al. 2005 , 2012 ; Miller-Butterworth et al. 2007 ; Meredith et al. 2011 ). Second, we resolve relationships within Mormoopidae, a small family of bats whose species—40% of which are island endemics—are the dominant mammalian insectivores of the West Indies ( Simmons 2005 ; Dávalos 2006 ). Finally, as the first multilocus species-level phylogenies of Phyllostomidae based on mostly nuclear sequences, these analyses generally confirm the relationships among species within genera that had been previously sampled using primarily mitochondrial genes and a small number of nuclear loci ( Dávalos and Jansa 2004 ; Dávalos 2007 ; Porter et al. 2007 ; Hoffmann et al. 2008 ; Redondo et al. 2008 ; Velazco and Lim 2014 ; Velazco and Patterson 2014 ).

However, aspects of the century-long debate on the taxonomy and systematics of the family remain open (see Wetterer et al. 2000 ). In particular, the phylogenetic position of Lonchorhininae remains unresolved (see Baker et al. 2003 ; Rojas et al. 2011 ; Dumont et al. 2012 ). Lonchorhina bats are morphologically distinct, so that relationships inferred based on morphology are also uncertain ( Wetterer et al. 2000 ; Dávalos et al. 2012 ). Additionally, our in-depth sampling highlights the need for future systematic revision of Hsunycteris and Lonchophylla ( Parlos et al. 2014 ), and species complexes such as L. silvicolum ( Medellín 1989 ; Porter et al. 2003 ), whose named populations, centralis and laephotis , are not sister taxa.

Mechanisms of Change in Geographic Range and Biogeographic History

The model that best explains the biogeography of New World Noctilionoidea, the DEC+J, implies within-continent cladogenesis and anagenetic dispersal as the most important events in the history of the group (see Fig. 3 ). Within-continent speciation is consistent with ecological differentiation as a mechanism of diversification among noctilionoids ( Dumont et al. 2012 ; Rojas et al. 2012 ). Future biogeographic analyses of areas within continents can test the extent of sympatric cladogenesis that cannot be captured by the large geographic scope of these analyses.

Founder event speciation is particularly important in the historical biogeography of clades currently distributed on oceanic archipelagos ( Matzke 2014 ). Here, we find support for founder speciation events at a continental scale. We found dispersal was more frequent than vicariance at cladogenesis. This suggests that dispersal played a more prominent role in diversification of neotropical bats than landscape reconfigurations, a result that is consistent with diversification in another group of flying vertebrates: neotropical birds ( Smith et al. 2014 ). Assuming noctilionoids first reached the Neotropics through South America (see Fig. 4 and below), and given the “splendid isolation” of this continent until ∼3.5 Ma ( Simpson 1980 ; Coates and Stallard 2013 ), ∼17 of the 191 nodes of the phylogeny for which founder event speciation events are inferred between these continents correspond to oceanic dispersal of ancestral bat populations. This primarily involves new lineages evolving from South American ancestors in Central and North America, or the Antilles (Supplementary Table S6).

Estimates of ancestral areas support a South American origin for New World noctilionoids (see Fig. 4 ). This result is consistent with a biogeographic scenario recently proposed by Gunnell et al. (2014) , in which the MRCA of neotropical noctilionoids reached South America through a series of dispersal events from Australia and across Antarctica. The analyses presented here include the greatest number of fossil-based calibrations of any previous phyllostomid phylogeny (see Baker et al. 2012 ), providing a strong foundation for inferring the timing of divergence events among noctilionoid bats. The estimated divergence time of New World noctilionoids and Mystacinidae (40.3–50.8 Ma), and the inferred age of the MRCA of New World noctilionoids (37.6–47.0 Ma) (see Fig. 2 ) suggest that noctilionoids reached the Neotropics during the Eocene from Australia and across Antarctica (see Hand et al. 2005 for mystacinid fossils from the late Oligocene and early Miocene in Australia).

Although geological analyses suggest that the landmasses of Australia and Antarctica were connected at ∼45 Ma ( White et al. 2013 ), and Antarctica and South America were connected until 41 Ma ( Scher and Martin 2006 ), no paleogeographic estimates are available to differentiate among land bridges, chains of islands, or over-sea dispersal (but see Gunnell et al. 2014 ). The paleoflora of Patagonia in South America and the Antarctic Peninsula suggests that a moist, frost-free corridor through Antarctica may have connected Australia and South America during the Paleocene–Eocene thermal maximum ( Zachos et al. 2001 , 2008 ; Reguero et al. 2002 ). If corroborated, this corridor may have facilitated the colonization of the New World by Old World noctilionoids.

There are two important gaps that need to be filled to test this scenario. The first is searching the fossil record of Australia, Antarctica, and southern South America for new noctilionoid fossils that would support dispersal through the proposed corridor between southern continents. The second is the inclusion of critical stem noctilionoid fossils from Florida in phylogenies. Relationships between these taxa and extant noctilionoids remain elusive ( Czaplewski and Morgan 2012 ; Morgan and Czaplewski 2012 ; Dávalos et al. 2014 ). The geographic distribution of these fossils, however, could upend the current narrative of noctilionoid dispersal via South America.

Biogeographic analyses point to South America, along with Central and North America, as the region of origin of Phyllostomidae (see Fig. 4 ). This result is also supported by analyses of geographic distributions pointing to South America and Central and North America as two important centers of diversification of phyllostomids ( Arita et al. 2014 ). The inferred ancestral distribution of phyllostomids was discontinuous, encompassing landmasses separated by oceanic barriers until the final closure of the Isthmus of Panama ( Coates and Stallard 2013 ). A possible connection between the two regions at 15 Ma is uninformative because the MRCA of phyllostomids (27–34 Ma) predates the proposed connection by at least 12 Ma ( Montes et al. 2015 ). The fossil record supports North America being part of the ancestral region of phyllostomids. The oldest fossil of Mormoopidae, the sister family to Phyllostomidae, is from the Oligocene of Florida (30.8–33.3 Ma) ( Czaplewski and Morgan 2003 ). Two newly described fossils from the same deposit and period (family Speonycteridae) appear to be closely related to mormoopids ( Czaplewski and Morgan 2012 ). Land exposure in the Caribbean between 48 and 38 Ma makes for another possible route of dispersal of early noctilionoids ( Iturralde-Vinent and MacPhee 1999 ). The South American MRCA of New World noctilionoids that dispersed from South America to southeastern North America dates back to this period (between 37 Ma and 47 Ma).

Given the current range of some phyllostomids, and if ancestral populations were as vagile as some extant lineages, a discontinuous ancestral distribution is plausible. Gene flow in two stenodermatines, Ardops nichollsi and Artibeus jamaicensis , has been maintained over time among islands of the Lesser Antilles that are separated by as much as 100 km ( Carstens et al. 2004 ; Larsen et al. 2010 ). The occurrence of shared haplotypes among the Greater Antilles for A. jamaicensis , and the smaller glossophagines Erophylla sezekorni and Erophylla bombifrons , indicates ancient migration between islands separated by at least 120 km of open waters ( Fleming et al. 2009 ; Muscarella et al. 2011 ). If it had the ability to disperse over water, the MRCA population of phyllostomids may have had a fragmented distribution. Then, the MRCA of phyllostomids split into the Central American lineage that evolved into the MRCA of Macrotinae and into the South American lineage that evolved into the MRCA of Karyovarians (see Fig. 4 ).

The Antilles have been proposed as an important area for the diversification of noctilionoids during the Miocene, with proposed insular origins for Mormoopidae, subtribe Stenodermatina (short-faced bats), and a clade within the subfamily Glossophaginae ( Dávalos and Jansa 2004 , 2006 , 2007 , 2009 ). Biogeographic analyses rejected an Antillean origin for both mormoopids and short-faced bats, and instead supported continent-to-island colonization in the Miocene followed by speciation on the islands (see Fig. 4 ). The Antilles was inferred as the ancestral region for the clade comprising all glossophagines that are currently restricted to the islands ( Brachyphylla , Erophylla , Phyllonycteris , and Monophyllus ) and the mainland genera Leptonycteris and Glossophaga (see Fig. 4 ). This clade represents the only example of island-to-mainland colonization in New World noctilionoids corroborated with new data. The timing of colonization is also consistent with Miocene sea-level minima as geological changes that facilitated range expansion first from the mainland to the islands, and then from the islands to the mainland ( Haq et al. 1987 ; Dávalos 2009 ).

Although fossil taxa were not included as terminals in the phylogeny used for biogeographic analyses, the fossil record independently corroborates the estimates of ancestral areas based on extant species. The fossils used to calibrate the phylogeny shown in Figure 4 were assigned to nodes based on clade-specific diagnostic characters ( Czaplewski and Morgan 2003 ) or recent character-based estimates of phylogeny ( Dávalos et al. 2014 ), and also account for the fossil locality and horizon ( Czaplewski and Morgan 2003 ). These fossils are the best available evidence of ancestral range in neotropical noctilionoids ( Morgan and Czaplewski 2012 ). The geographic locality of each fossil then represents an independent line of evidence for testing estimates of ancestral range. If the geographic location of a fossil is included in the range estimated for the MRCA of the clade in which each fossil belongs, then fossils do not reject the estimate of ancestral area. This is the case for six of the eight fossils ( Notonycteris , Brachyphylla nana , Phyllonycteris major , Palynephyllum antimaster , Cubanycteris silvai , and Phyllops vetus ), and seems unrelated to the age of the horizons in which the fossils were found (see Fig. 4 ). For the remaining two fossils (the stem mormoopid and Desmodus archaeodaptes ), the ancestral regions for the corresponding MRCAs were inferred as equivocal in biogeographic analyses. Given the broad distribution of these two taxa, we do not expect qualitative changes in results when these and other fossils are included as terminal taxa in the phylogeny (an approach that is under development in BioGeoBEARS; http://phylo.wikidot.com/fossil-data-in-biogeographical-analysis-in-biogeobears ). Extant Desmodontinae are widely distributed in North, Central, and South America, and Pleistocene-era fossils in the subfamily have been collected from these regions ( Morgan et al. 1988 ; Czaplewski et al. 2003a ), as well as the Antilles ( Koopman 1958 ; Suárez 2005 ). Fossil Mormoops of similar age have been uncovered on all three biogeographic regions as well ( Morgan and Czaplewski 2012 ). New approaches that integrate fossil data in biogeographic analysis will clarify the historical biogeography of New World Noctilionoidea, but are unlikely to compromise our main conclusions, since the most ancient ancestral areas are estimated to encompass both mainland regions.

New World Noctilionoids do not Support a Quaternary Origin for Extant Neotropical Biodiversity

Traditionally, South America has been viewed as the region of origin for the family Phyllostomidae ( Lim 2009 ). Our results do not support South America as the singular area of origin for the family, but they do show this continent has been much more important for subsequent phyllostomid diversification than for other New World noctilionoid families. South America was the ancestral region of 9 of the 11 subfamilies of Phyllostomidae, and for the most species-rich genera of this family (e.g., Sturnira , Platyrrhinus , and Artibeus ) (see Fig. 4 ). This supports noctilionoids as a useful clade to test hypotheses on the impact of the Neogene on species diversification in the Neotropics ( Hoorn et al. 2010 ) and the Quaternary ( Hooghiemstra and van der Hammen 1998 ).

Analyses that accounted for phylogenetic uncertainty in the evolutionary history of Noctilionoidea found more speciation among noctilionoids in the pre-Quaternary than in the Quaternary (see Fig. 5 a). Simulations of birth–death trees under different values of speciation and extinction (including a pure-birth Yule process), however, indicate that using the number of speciation events as a proxy for speciation rate is a flawed approach (see Fig. 5 c and Supplementary Figs. S7–S8).

In each of the simulated scenarios, the average number of speciation events between extant species during the Quaternary was always significantly higher than during the pre-Quaternary. This indicates that an empirical pattern in which the number of speciation events is lower during the pre-Quaternary than in the Quaternary can arise through lineage accumulation, without any special influence of the Quaternary on speciation. As long as the speciation rate exceeds the extinction rate, one can expect more lineages, and therefore a greater number of speciation events closer to the present, even when the speciation rate remains constant through time. For clades in which fewer speciation events are traced to the pre-Quaternary compared to the Quaternary (see some examples in Rull (2008) ), support for a Quaternary origin of species diversity in the Neotropics requires Pleistocene speciation rates to be effectively higher than pre-Pleistocene rates.

An alternative to comparing the number of speciation events is to model time-varying diversification rates ( Morlon et al. 2011 ; Stadler 2011 ; Rabosky 2014 ), or the relationship between environmental change and diversification rates ( Condamine et al. 2013 ). We detected a shift in speciation rate in the branch leading to the MRCA of Stenodermatinae, at least 15 myr before the start of the Quaternary (Supplementary Figs. S4–S6). Crucially, this shift is clade- and not time-dependent. Instead of finding an increase in speciation rates closer to the present, as expected if the Pleistocene had disproportionate impact on diversification, more complex models of diversification rates in New World noctilionoids identified more clade-specific shifts (Supplementary Fig. S6). Finally, these patterns were robust to protracted speciation. The evidence from modeling diversification rates is therefore compelling in rejecting any of the temporal patterns expected if Quaternary climate change had disproportionately influenced speciation among noctilionoids.

Cryptic species have recently been found in Noctilionidae ( Pavan 2013 ) and Mormoopidae ( Clare et al. 2013 ; Thoisy et al. 2014 ). We expect that new species will continue to be found, particularly in the species-rich Phyllostomidae ( Velazco and Lim 2014 ; Velazco and Patterson 2014 ), but also in other noctilionoid families ( Velazco et al. 2014 ). Most of the new species described in those studies are allopatric sister species. Even if the divergence times of the new pairs of extant species correspond to the Quaternary, this would have little effect on our general conclusions because the main underlying pattern is one of clade-specific shifts in diversification and many of the new discoveries correspond to the rapidly diversifying clade of frugivorous phyllostomids ( Velazco and Lim 2014 ; Velazco and Patterson 2014 ). Our results from BAMM analyses are consistent with those from a recent analysis of diversification dynamics for the entire order Chiroptera ( Shi and Rabosky 2015 ). Despite including fewer noctilionoid species than here (i.e., 176 vs. 194) and the broader taxonomic scope of that study (i.e., all 20 families of Chiroptera), they found a single diversification shift at the base of Stenodermatinae. Their result supports our conclusion that the Quaternary did not play a substantive role in diversification of neotropical noctilionoids.

Conclusions

Biogeographic analyses revealed an island-like distribution in the ancestor of phyllostomids as well as the great importance of dispersal in the evolution of noctilionoids. Despite evidence of a long history outside the Neotropics, this area is the most important to diversification of that family. Simulations showed that dating the divergences of current sister species cannot effectively account for the role of particular geological periods on the diversification of any clade, since significant patterns can emerge through a birth–death process of lineage accumulation. Explicit modeling of diversification rates is, therefore, indispensable to analyze the causes of neotropical biodiversity. Such explicit modeling reveals that, instead of being associated with particular geological epochs, changes in diversification in New World noctilionoids are clade-specific and likely associated with corresponding traits.

SupplementaryMaterial

Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.s533p .

Funding

This work was supported by the National Science Foundation (DEB-0949759 and DEB-1442142) to LMD. DR was supported by Foundation for Science and Technology, Portugal (www.fct.pt), fellowship SFRH/BPD/97707/2013. The Portuguese Foundation for Science and Technology supported CESAM RU (UID/AMB/50017) through national funds and FEDER funds, within the PT2020 Partnership Agreement.

Acknowledgments

We thank Laurel Yohe, James Herrera, and two anonymous reviewers for their helpful comments to a previous version of this paper. For tissue loans, LMD thanks Nancy B. Simmons and the staff at the Ambrose Monell Cryo-Collection at the American Museum of Natural History; Robert J. Baker, and the staff at the mammal tissue collection of the Museum of Texas Tech University; Chris Conroy at the Museum of Vertebrate Zoology of the University of California, Berkeley; Kris Helgen and the staff at the National Museum of Natural History; Burton Lim and Judith Eger at the Royal Ontario Museum; and Bruce Patterson and John D. Phelps at the Field Museum.

Author notes

Associate Editor: Robb Brumfield