Taxonomic sampling.—In order to address our objective, we obtained tissue samples of specimens from various institutions or individuals representing 15 (94%) of 16 genera in the subfamily Molossinae (Appendix I). No tissues from Tomopeas ravus or Platymops setiger were available. To our knowledge, we have samples representing the largest collection of molossids and we were able to sample the type species for all but 4 genera. The following type species were included: Cheiromeles torquatus Horsfield, 1824, Cynomops abrasus (Temminck, 1827), Eumops perotis (Schinz, 1821), Molossops temminckii (Burmeister, 1854), Molossus molossus (Pallas, 1766), Mops mops (de Blainville, 1840), Mormopterus jugularis (Peters, 1865), Myopterus daubentonii Desmarest, 1820, Neoplatymops mattogrossensis (Vieira, 1942), Nyctinomops femorosaccus (Merriam, 1889), Sauromys petrophilus (Roberts, 1917), and Tadarida teniotis (Rafinesque, 1814). Additionally, we included samples from Vespertilionidae (Myotis chinensis, M. velifer, M. yumanensis, and Antrozous pallidus) and Natalidae (Natalus stramineus and Chilonatalus [Natalus] micropus) as outgroups based on their familial relationship to Molossidae (Miller-Butterworth et al. 2007; Teeling et al. 2002, 2005).

Generation of sequence data.—We isolated whole genomic DNA from either wing punches; frozen liver or muscle tissue; or liver tissue stored in ethanol, dimethylsulfoxide, or lysis buffer (Longmire et al. 1997). All DNA was isolated using the DNeasy Tissue Kit (QIAGEN Inc., Valencia, California) following manufacturer's protocol, and samples stored in ethanol or dimethylsulfoxide were soaked in phosphate-buffered saline prior to extraction.

We used polymerase chain reaction to amplify 1 mitochondrial locus, nicotinamide adenine dinucleotide dehydrogenase subunit 1 (ND1), and 3 nuclear loci, beta fibrinogen intron 7 (βFIB), recombination activating gene 2 (RAG2), and dentin matrix protein 1 exon 6 (DMP1) using the following primers. For ND1, we used ER65 5′-CCTCGATGTTGGATCAGG-3′ and ER66 5′-GTATGGGCCCGATAGCTT-3′ (Petit et al. 1999). We used BFibI7L 5′-TCCCCAGTAGTATCTGCCATTAGGGTT-3′ and BFibI7U 5′-GGAGAAAACAGGACAATGACAATTCAC-3′ (Prychitko and Moore 1997) to amplify βFIB, and RAG2-179F 5′-CAGTTTTCTCTAAGGAYTCCTGC-3′ and RAG2-1458R 5′-TTGCTATCTTCACATGCTCATTGC-3′ (Stadelmann et al. 2007) to amplify RAG2. Lastly, we designed Den4F 5′-AGACAAGGAGGAAACTCCAGACT-3′ specifically for molossids and used Den10R 5′-GTTGCTCTCTTGTGATTTGCTGC-3′ (Van Den Bussche et al. 2003) to amplify DMP1. Polymerase chain reactions consisted of 200–500 ng of DNA, 3 U of Taq polymerase, 0.16 µM of forward and reverse primer, 2 mM MgCl₂, 0.16 mM deoxynucleoside triphosphates, and 1X reaction buffer for a total volume of 12.5 µl. The same thermal profile was used for each fragment with varying annealing temperatures: initial denaturation of 92°C for 2 min, followed by 39 cycles of 94°C for 1 min, annealing for 1 min (ND1 48°C, βFIB 48°C, RAG2 61–63°C, and DMP1 55–57°C), and 72°C for 1 min, with a final extension of 72°C for 10–30 min.

The βFIB polymerase chain reaction fragments were sequenced early in this project via cloning. Polymerase chain reaction products were separated on 1% low-melt agarose gels and ligated directly into a 2.1 TopoTA vector following manufacturer's protocols except using only one-fourth of the reagent volumes (Invitrogen, Carlsbad, California). We purified plasmids containing polymerase chain reaction products following the protocol of the QIAprep Spin Miniprep Kit (QIAGEN Inc.). Cycle sequencing of both strands was performed with dye-labeled M13 primers using ThermoSequenase (USB, Cleveland, Ohio) following the manufacturer's protocol. Two clones from each sample were sequenced and collected on a LICOR NEN Global IR2 sequencing system (LI-COR, Inc., Lincoln, Nebraska) using ESeq version 3 software.

Polymerase chain reaction fragments of ND1, RAG2, and DMP1 were visualized with electrophoresis and were cleaned using Exosap-IT (USB) before being sequenced on a Beckman Coulter CEQ 8000 genetic analyzer using the GenomeLab-DTCS Quick Start Kit (Beckman Coulter Inc., Fullerton, California). We followed the manufacturer's protocol but used one-fourth volume reactions. All DNA sequences obtained have been deposited in GenBank (Appendix I).

Sequence alignment and single-gene analyses.—Sequencher version 4.7 and version 4.8 (Gene Codes Corporation, Ann Arbor, Michigan) were used to align overlapping fragments, and ambiguities were refined by eye. Alignment parameters used for the fibrinogen intron were adjusted to allow large gaps and minimized the number of single-base insertions–deletions (indels). Coding regions were translated when applicable. We also downloaded sequences from GenBank (Tadarida insignis AB079813 for ND1, Mormopterus acetabulosus EU487473 and Myotis chinensis AY726646 for βFIB, N. stramineus AY141024 for RAG2, and C. micropus AY141883 and N. stramineus AY141884 for DMP1) to increase our sampling. The final alignment from each data set was imported in PAUP*, version 4.0b10 (Swofford 2003) for parsimony analysis. Each data set was analyzed separately in Modeltest version 3.7 (Posada and Crandall 1998) to determine the appropriate model of evolution for likelihood analyses. We also used MrBayes version 3.1.2 (Ronquist and Huelsenbeck 2003) with 2 simultaneous runs of 1–2 million generations and the GTR+G model was used for each gene. Convergence was measured by values of the average standard deviation of split frequencies (average standard deviation of split frequencies less than 0.01 was considered convergence of the 2 simultaneous runs). Trees were sampled every 100 generations. We used a burn-in to discard 25% of saved trees. Bootstrap analysis (100 pseudoreplicates) using maximum-likelihood (ML) criteria was used as an additional measure of branch significance (Felsenstein 1985; Posada and Crandall 1998; Swofford 2003). We considered significantly supported nodes as those with both ML bootstrap (BS) values >70 (Hillis and Bull 1993; Suzuki et al. 2002) and Bayesian posterior probabilities (BPPs) >0.95. Finally, average genetic distances were calculated in PAUP* using uncorrected p for the mitochondrial data set and the appropriate model of evolution for each nuclear data set. Additionally, we plotted uncorrected p distance versus ML distance to look for saturation in ND1. Results from the saturation curve prompted us to also analyze ND1 sequence partitioned by codon position using Bayesian methods and to analyze the translated amino acid sequence using parsimony in MEGA version 4 (Tamura et al. 2007).

Analysis of combined data set.—Data from all 4 genes were combined for a single analysis to provide a better estimate of the relationships among genera (Gadagkar et al. 2005). To compensate for different taxonomic sampling in some genes, we created composite sequences from 2 individuals of the same species when possible. We retained taxa in the alignment if represented by 2 or more gene sequences except for Myopterus daubentonii and Eumops perotis (both with only a single gene sequence), which were retained because they were type species for these genera. This approach was used because missing data have been shown to be less important than missing taxa in a phylogenetic analysis (de la Torre-Bárcena et al. 2009; Wiens 2003).

Incongruence length-differences tests (Farris et al. 1995) were conducted to determine congruence between the different genes (Cunningham 1997) using the partition homogeneity test in PAUP* 4.0b10 (Swofford 2003). Because there was no significant incongruence among data sets (P  =  0.93), we created and analyzed a combined, multilocus alignment. As a 2nd justification for combining all genes, we evaluated congruence by noticing the lack of any significantly supported nodes in 1 data set that conflicted with strongly supported nodes in another (de Queiroz 1993). The combined data set was analyzed using a partitioned Bayesian analysis (MrBayes version 3.1.2—Ronquist and Huelsenbeck 2003) with 2 simultaneous runs of 2 million generations and a burn-in equal to 25% of saved trees. The data were partitioned by gene and the GTR+G model was used for each (estimates of rate variation were unlinked). In addition, we ran a parsimony analysis on the combined data set with gaps as a 5th state. Nodal support for this analysis was assessed with parsimony bootstrapping (1,000 pseudoreplicates—Felsenstein 1985).

Estimating divergence time.—We used the combined data set to estimate divergence times. Prior to estimating divergence times, we performed a likelihood ratio test (Felsenstein 1981) using the combined data set to test the hypothesis of a strict global clock. We then used the ML tree from PAUP* in r8s version 1.7 (Sanderson 2003) to rescale the branch lengths on the phylogram to years instead of substitutions per site. The oldest fossil molossid, Wallia, was used to calibrate the analysis with a minimum of 42 million years ago (mya) as the deepest divergence date for all molossids (Czaplewski et al. 2003). This tree served as the starting tree for analysis in BEAST version 1.5 (Drummond and Rambaut 2007). We performed initial analyses to optimize parameters. Nodal dates were examined using a Yule species prior and a lognormal distribution, and we attempted the analysis using GTR+I+G and HKY+I+G models of sequence evolution. Based on preliminary results, our final analysis consisted of 3 separate runs of 20 million generations with 10% burn-in using only the HKY+I+G model. Results from the final 3 runs were combined in Tracer version 1.5 (Rambaut and Drummond 2007) to check for effective sample sizes (Drummond et al. 2002) greater than 200 for all parameters and stable convergence on a unimodal posterior. Time estimates were calculated based on the combination of log and tree files from BEAST in TreeAnnotator version 1.5.4 (Drummond and Rambaut 2007).

Nicotinamide adenine dinucleotide dehydrogenase subunit I.—We analyzed 888 base pairs (bp) of the ND1 gene from 33 molossid taxa representing 15 genera (Table 1). Divergences (uncorrected p) among molossid genera ranged from 9.6% between Chaerephon and Mops to 27% between Cheiromeles and Sauromys. Even though there were many variable sites, only 56 were parsimony informative (Table 1). An ML tree was obtained (−ln  =  8,916.7756) using the GTR+I+G model of evolution with the following parameters: base frequencies  =  0.3727, 0.2927, 0.0904; 6 substitution types and the rate matrix  =  6.8800, 36.1327, 4.8845, 0.7107, 97.4996; gamma shape parameter  =  1.1593, proportion of invariant sites  =  0.5043. After 2 million generations in the Bayesian analysis, the average standard deviation of split frequencies was 0.008416. Topology of the ML tree was not identical to that of the Bayesian tree (Fig. 2), but there was not significant conflict. The Bayesian tree grouped Cheiromeles, Mormopterus, and Sauromys (but with poor BPP support values), whereas the ML tree placed each of these taxa on different branches with all BS values <50%.

Fig.2.

Phylogram based on Bayesian analysis of 888 base pairs of Nicotinamide adenine dinucleotide dehydrogenase 1 partitioned by codon position. Support values on branches are Bayesian posterior probabilities (BPPs) followed by maximum-likelihood (ML) bootstrap (BS) values. Solid black dots indicate nodes with BPP ≥0.95 and ML (GTR+I+G) BS support ≥70. Dashes indicate values less than 50.

Table 1.

Descriptive information about each data set. ND1  =  Nicotinamide adenine dinucleotide dehydrogenase subunit 1; βFIB  =  beta fibrinogen intron 7; DMP1  =  dentin matrix protein 1 exon 6; RAG2  =  recombination activating gene 2; ML  =  maximum likelihood.

Overall, the ND1 data significantly supported more recent divergences such as the monophyly of genera (Mops, Otomops, Eumops, Nyctinomops, and Molossus). The only intergeneric relationships supported were a Chaerephon–Mops clade and a Molossus–Promops clade. Chaerephon species and Tadarida species did not form monophyletic clades although these results were not significantly supported. A clade containing 7 New World genera (Cynomops, Molossops, Neoplatymops, Molossus, Promops, Nyctinomops, and Eumops) also was recovered by the ND1 analysis (Fig. 2).

Because saturation was present in ND1 (Fig. 3), amino acid sequences were analyzed using parsimony methods. A total of 75 of 296 amino acid positions were variable and 55 were parsimony informative. We recovered 206 shortest trees (consistency index [CI]  =  0.6667, retention index [RI]  =  0.7723) and calculated a majority rule consensus tree (not shown). Only 2 nodes had significant support based on bootstrap percentages—1 for the family Molossidae and the other joining the 2 Otomops species together.

Fig.3.

Uncorrected p distances plotted against maximum-likelihood distances for the Nicotinamide adenine dinucleotide dehydrogenase subunit 1 gene.

Beta fibrinogen intron 7.—We analyzed 792 bp of βFIB from 20 molossid taxa representing 10 genera (Table 1). Divergences (TVM+G) among molossid genera ranged from 2.08% between Sauromys and Tadarida to 9.20% between Nyctinomops and Mops. The highest divergence of 18% was seen in 2 pairwise comparisons (between Nyctinomops and Mops and the outgroup Myotis).

Numerous small indels were present as well as 1 large one. Unique indels were a 4-bp insertion in Mormopterus acetabulosus, a 4-bp insertion in Promops centralis, a deletion of 56 bp in Tadarida brasiliensis, and a 9-bp deletion in Eumops patagonicus. Shared indels were a 3-bp deletion in both Myotis species; a 3-bp insertion in Tadarida, Sauromys, Otomops, Mormopterus, and Myotis; a 19-bp deletion shared by 2 Nyctinomops species; and a 5-bp deletion shared by Eumops, Promops, and Molossus.

This gene had the highest percentage of parsimony-informative sites of the nuclear markers (Table 1). These values were calculated without considering indels as characters. An ML tree was obtained (−ln  =  3,185.94933) using the TVM+G model of evolution with the following parameters: base frequencies  =  0.3117, 0.2142, 0.1880; 6 substitution types and the rate matrix  =  0.6600, 2.5844, 0.2737, 0.6923, 2.5844; gamma shape parameter  =  1.3519. For the Bayesian analysis, after 1 million generations in 2 simultaneous runs, the average standard deviation of split frequencies was 0.006761. Topology of the ML tree differed slightly from that of the Bayesian tree (Fig. 4) in 2 places but neither was significantly supported. The Sauromys–Tadarida clade was united with the New World clade and Otomops clustered with the Mops–Chaerephon clade.

Fig.4.

Phylogram based on Bayesian analysis of 792 base pairs of beta fibrinogen intron 7. Support values on branches are Bayesian posterior probabilities (BPPs) followed by maximum-likelihood (ML) bootstrap (BS) values. Solid black dots indicate nodes with BPP ≥0.95 and ML (TVM+G) BS support ≥70. Dashes indicate values less than 50.

The analysis of βFIB sequences significantly supported a Chaerephon–Mops clade resulting in nonmonophyly of either genus (Fig. 4). A New World clade containing Eumops, Nyctinomops, Promops, and Molossus had significant support and within it Promops–Molossus were sister taxa. A clade containing Tadarida and Sauromys had significant support; Sauromys was sister to Tadarida aegyptiaca, and this clade formed a sister to T. brasiliensis. There was not strong support for the monophyly of Molossinae in this analysis (Fig. 4).

Dentin matrix protein 1 exon 6.—We analyzed 496 bp from DMP1 for 20 molossid taxa representing 13 genera (Table 1). Divergence (TrN+G) among molossid genera ranged from 1.98% between Chaerephon and Mops to 6.22% between Chaerephon and Neoplatymops. DMP1 had the lowest number of parsimony-informative sites for the nuclear genes (Table 1).

An ML tree was obtained (−ln  =  1,842.8861) using the TrN+G model of evolution with the following parameters: base frequencies  =  0.3477, 0.2244, 0.2666; 6 substitution types and the rate matrix  =  1.0000, 2.4714, 1.0000, 1.0000, 6.7544; gamma shape parameter  =  0.6615. After 1 million generations in 2 simultaneous Bayesian runs, the average standard deviation of split frequencies was 0.007991. Topology of the ML tree was similar to that of the Bayesian tree (Fig. 5), except that Tadarida fulminans formed a separate branch between Mormopterus and all other molossines, and Neoplatymops clustered with Cynomops and Molossops in the ML tree.

Fig.5.

Phylogram based on Bayesian analysis of 496 base pairs of dentin matrix protein 1. Support values on branches are Bayesian posterior probabilities (BPPs) followed by maximum-likelihood (ML) bootstrap (BS) values. Solid black dots indicate nodes with BPP ≥0.95 and ML (TrN+G) BS support ≥70. Dashes indicate values less than 50.

Analysis of DMP1 significantly supported a Chaerephon–Mops clade, but did not support monophyly of Chaerephon. A New World clade (Cynomops, Molossops, Neoplatymops, Eumops, Nyctinomops, Molossus, and Promops) also was recovered with significant support. Within the New World clade, monophyly of the genera Eumops and Nyctinomops was supported. The position of Mormopterus as the most basal lineage was significantly supported. No significant resolution among Tadarida, Otomops, and Sauromys was recovered, but Sauromys formed a sister relationship with T. brasiliensis (Fig. 5).

Recombination activating gene 2.—We analyzed 1,040 bp from RAG2 for 22 molossid taxa, representing 12 genera (Table 1). Cheiromeles had a deletion of 3 codons that was not present in the other taxa. Divergence (TIM+G) among molossid genera was lower than for the other nuclear genes and ranged from 1.2% between Tadarida and Sauromys to 4.97% between Cheiromeles and Mops. In RAG2, approximately 10% of the sites were parsimony informative (Table 1).

An ML tree was obtained (−ln  =  3,382.17706) using the TIM+G model of evolution with the following parameters: base frequencies  =  0.2971, 0.2030, 0.2219; 6 substitution types and the rate matrix  =  1.0000, 3.9057, 0.6423, 0.6423, 6.2332; gamma shape parameter  =  0.3134. In the Bayesian analysis, after 1 million generations in 2 simultaneous runs, the average standard deviation of split frequencies was 0.004860. Topology of the ML tree was similar to that of the Bayesian tree (Fig. 6) except that Otomops joined a clade of Sauromys and Tadarida but by a short, unsupported branch.

Fig.6.

Phylogram based on Bayesian analysis of 1,040 base pairs of recombination activating gene 2. Support values on branches are Bayesian posterior probabilities (BPPs) followed by maximum-likelihood (ML) bootstrap (BS) values. Solid black dots indicate nodes with BPP ≥0.95 and ML (TIM+G) BS support ≥70. Dashes indicate values less than 50.

Significant results from the Bayesian analysis included a Chaerephon–Mops clade, but monophyly of Chaerephon was not supported based on RAG2 data. A New World clade (Cynomops, Nyctinomops, Eumops, Molossus, Promops, and Neoplatymops) had significant support; however, the only well-supported intergeneric relationship recovered was a Promops–Molossus clade. Cheiromeles was the most basal member of the family with significant support. A clade containing 4 species of Tadarida and Sauromys petrophilus was recovered although support was low (Fig. 6).

Analysis of combined data set.—Parsimony analysis of the combined data set when gaps were coded as missing resulted in 5 equally parsimonious trees of 2,700 steps. There were 616 parsimony-informative characters used. A majority rule consensus resulted in poor resolution (tree not shown). A 2nd analysis using gaps as a 5th character state resulted in 1 most parsimonious tree with 3,080 steps (CI  =  0.4747, homoplasy index [HI]  =  0.5253, RI  =  0.5098). There were 673 parsimony-informative sites. Significant bootstrap values are mapped on the Bayesian tree (Fig. 7) except for 2 patterns that conflicted with the Bayesian tree. Parsimony analysis supported Nyctinomops femorosaccus and N. macrotis as sister taxa (BS  =  93) without N. aurispinosus. Parsimony analysis also supported a clade containing all Eumops except for E. hansae (BS  =  91). Partitioned Bayesian analysis of the same concatenated data set produced a similar topology to the parsimony bootstrap consensus tree, but with additional resolution (Fig. 7). After 2 million generations in 2 simultaneous runs, the average standard deviation of split frequencies was 0.006035. The combined Bayesian analysis of all 4 genes suggested some relationships not seen in the individual gene trees, although with poor support values. For example, a clade containing all Tadarida (plus Sauromys), a clade of all Old World taxa, and the clustering of Myopterus with Chaerephon and Mops were suggested by the Bayesian analysis with low support (Fig. 7). All other significantly supported relationships also were seen in at least 1 of the individual gene analyses, such as support for the subfamily Molossinae, the basal position of Cheiromeles, and a New World clade.

Fig.7.

Phylogram based on partitioned Bayesian analysis of 4 loci in a combined data set. Support values on branches are Bayesian posterior probabilities (BPPs) followed by bootstrap (BS) values. Solid black dots indicate nodes with BPP ≥0.95 and BS support ≥70. Dashes indicate values less than 50. Circled BPP values are on nodes that were not present in the parsimony tree because of an unresolved position for Nyctinomops aurispinosus, Eumops hansae, and Cynomops abrasus that failed to unite them with other New World taxa.

Divergence time.—A strict molecular clock was rejected (2ΔL  =  108.0556; p  =  0.00) and a relaxed molecular clock was used to estimate divergence times. The divergence of molossids from outgroup taxa was estimated to be 63 (51.5–75.3) mya. The position of outgroups in this analysis was reversed compared to the Bayesian analysis with vespertilionids as the sister to Natalidae + Molossidae. Additionally, the order of divergence within the Old World clade differed from the Bayesian analysis slightly, but was not significantly supported (Figs. 7 and 8). A split between an Old World clade and a New World clade was dated at 28.7 (24.4–33.2) mya, with a series of rapid divergence events soon after this within each clade. Cheiromeles and Mormopterus diverged before this event and were dated at 43.5 (42.0–46.5) and 33.3 (27.8–39.1) mya, respectively. Error bars were much longer for the earliest divergences and shorter within Molossinae (Fig. 8).

Fig.8.

Chronogram of the phylogeny of Molossinae generated from BEAST analysis of the combined data set using the relaxed uncorrelated lognormal molecular clock. Mean divergence estimates indicated at corresponding nodes with ages in millions of years before present. Gray bars represent the 95% highest posterior density interval for each node.

Mops–Chaerephon clade.—Data from all 4 genes agreed that there was significant support for a Mops–Chaerephon clade resulting in paraphyly of Chaerephon in all cases. This pattern was also recovered by Gregorin (2000) and by Freeman's (1981:102) quantitative analysis, and was seen in a supertree analysis of Jones et al. (2005). In addition, Chaerephon + Mops was 1 of the few statistically significant nodes recovered by analysis of cytochrome-b sequence (Agnarsson et al. 2011) even though neither genus was monophyletic in that analysis. No consistent picture emerges as to the relationships among these Old World species, but Mops condylurus from Africa clusters with Mops leucostigma from Madagascar when present, and these 2 form a close relationship with African Chaerephon and not other Mops. These results are in agreement with Rosevear (1965), who proposed that M. condylurus should be recognized as Chaerephon based on morphology, and with Peterson et al. (1995), who noted the intermediate morphology of M. condylurus and M. leucostigma.

In addition, the βFIB analysis clustered species of Mops and Chaerephon that are geographically closely distributed; M. mops and C. plicatus from the Asian region clustered together and M. condylurus and C. pumilus from the African region clustered together. Certainly wider taxonomic sampling of the species in this clade would be valuable for determining the validity of these generic designations. Both genera, Chaerephon and Mops, were elevated from Tadarida by Freeman (1981), Koopman (1984), and Simmons (2005), but not by Corbet and Hill (1986), Peterson et al. (1995), or Legendre (1984a). However, our results confirm that Chaerephon is distinct from Tadarida.

New World molossid clade.—Significant support for 7 genera in a New World clade was recovered by 3 nuclear genes and suggested by the mitochondrial gene tree (BPP  =  0.89) including Eumops, Molossus, Promops, Molossops, Cynomops, Nyctinomops, and Neoplaytmops. This arrangement conflicts with Freeman (1981; Fig. 1) because of the inclusion of Molossops (and Neoplatymops) in the New World clade. She placed these genera with Myopterus, Cheiromeles, and Mormopterus. It also conflicts with Gregorin (2000; Fig. 1) because New World genera were represented in each of his 3 proposed tribes or subtribes and there was no clade similar to the New World clade that we recovered in our study. As expected based on previous work (Freeman 1981), the New World clade did not include T. brasiliensis—the only New World member of this genus. Tadarida last shared a common ancestor with the New World clade almost 29 mya (Fig. 8).

A sister relationship recovered between Molossops and Cynomops was congruent with Peters et al. (2002). Neoplatymops (considered a subgenus of Molossops by Freeman [1981], Gregorin [2000], Legendre [1984a], and Simmons [2005]) did not form a monophyletic clade with Molossops in our analysis, supporting its generic status (Willig and Jones 1985). We note that another species elevated from Molossops, Cabreramops aequatorianus (Ibáñez, 1981), has been recognized by Eger (2007), but it was not included in our study.

There is strong support for monophyly of Eumops and Nyctinomops except by RAG2 data. Comparison to the quantitative analysis of Freeman (1981:102) reveals substantial convergence in the morphology of Eumops species as evinced by the nonmonophyly of the genus. For example, E. perotis clustered with Otomops and large Tadarida based on similar morphotypes that are influenced by loud, low-frequency echolocation, adaptation to dry environments, and similar diets and foraging behavior (Jones and Rydell 2003; Rydell and Yalden 1997). The inclusion of Cynomops within Nyctinomops by RAG2 is unexpected and probably the result of the lack of variable characters in this fairly conserved exon. Other significantly supported nodes within the New World clade were the Molossus–Promops clade. Freeman (1981), Gregorin (2000), and Sudman et al. (1994) also found they were sister taxa. The divergence between these 2 genera (Molossus and Promops) appears to be one of the most recent in the family (14.2 mya) and is even more recent than the oldest divergence within the genus Eumops (19.2 mya). A 5 bp deletion in the βFIB intron was shared by Eumops, Promops, and Molossus, suggesting a close relationship, but there was not significant branch support.

Tadarida and Sauromys.—Our analyses support Freeman's (1981) and Gregorin's (2000) hypothesis that the closest relative to the common Mexican free-tailed bat (T. brasiliensis) is not other North American free-tailed bats, but instead is T. aegyptiaca from Africa. In addition, the close relationship between Tadarida and Sauromys seen in 3 of our analyses (βFIB, DMP1, and RAG2) has been suggested by other authors (Freeman 1981; Harrison 1962; Peterson 1965). Peterson et al. (1995) stated that Sauromys was more closely related to Mormopterus than to any other genus. Similarly, Gregorin's (2000) morphological analysis found that Mormopterus was sister to Platymops (including Sauromys), but our results did not support this view. These 2 species of Tadarida (T. brasiliensis and T. aegyptiaca) joined with the other Tadarida species in some of our analyses, forming a poorly supported monophyletic clade in the combined Bayesian tree (Fig. 7).

Legendre (1984a) proposed the genus Rhizomops for a single species, T. brasiliensis, based on dental morphology. This view was not widely accepted (Owen et al. 1990) and is not warranted based on our results. Based on this phylogenetic analysis of molecular data, it appears that T. brasiliensis, T. aegyptiaca, and S. petrophilus form a monophyletic clade distinct from other Tadarida species. Interestingly, this relationship also was recovered in a quantitative analysis of morphology by Freeman (1981:102). Harrison (1962) also thought that Sauromys was closely related to T. aegyptiaca and, in fact, proposed that it was a transition between T. aegyptiaca and Platymops. Samples from P. setiger will be important for testing this hypothesis. Sauromys (as well as Platymops) was classified as Mormopterus by some authors (Freeman 1981; Koopman 1993, 1994; Legendre 1984a), but not by others (Corbet and Hill 1986; Gregorin 2000; Peterson et al. 1995; Simmons 2005). The level of genetic divergence and the topology recovered in this study regarding Mormopterus and Sauromys does not support the inclusion of Sauromys in the genus Mormopterus. To maintain monophyly of the genus, we propose that Sauromys is more appropriately recognized as a member of the genus Tadarida having split from its sister taxon in this analysis, T. aegyptiaca, approximately 6.5 mya. Inclusion of P. setiger and additional Tadarida species in future studies will be important for determining the extent to which this genus should be applied.

The status of the genus Tadarida is complex. The type species for Tadarida is T. teniotis, but poor resolution in our analyses prevents any definite conclusions regarding the evolutionary position of this genus. According to Peterson et al. (1995), Tadarida and genera recently recognized as separate from Tadarida (such as Chaerephon and Mops) have not been adequately diagnosed and share many characteristics with Eumops. Our results seem to provide evidence of the distant relationship between Tadarida and Eumops, clarify the distinct lineage that is Chaerephon and Mops, and confirm the validity of the genus Nyctinomops that was elevated from Tadarida by Freeman (1981). Peterson et al. (1995) thought elevation of Mormopterus by Freeman (1981) was justified, and our analyses suggest agreement but we cannot adequately resolve the validity of Mormopterus without additional samples.

Position of Cheiromeles.—The basal divergence of Cheiromeles we observed in this study is consistent with Legendre's (1984a) taxonomic arrangement placing this taxon in its own subfamily. It is also consistent with the relationships proposed by the quantitative analysis of Freeman (1981:102) that placed Cheiromeles outside of even the vespertilionid outgroup. The 2 species of Cheiromeles (naked bats) found in Southeast Asia are large-bodied and display unusual morphological traits unlike any other molossid. They are nearly hairless and possess a subaxillary pouch used for wing storage, an opposable thumb, and, finally, a dense protrusion of stiff hairs on the 1st toe of the hind foot (Leong et al. 2009). Their gular sac produces an oily secretion that is applied to their naked wings and skin to maintain the leathery texture (Leong et al. 2009). In Freeman's (1981) monograph on the family Molossidae, she stated that Cheiromeles was a derived member of an ancient molossid lineage because it has the most derived dental configuration in the family. In her cladogram, Cheiromeles was closely related to Myopterus, Molossops, and Mormopterus but the synapomorphies were not consistent within these genera. Gregorin (2000) found a similar arrangement (Fig. 1) in his cladistic analysis of the family. Because we were not able to generate Myopterus sequence for all loci, the sister relationship between Cheiromeles and Myopterus proposed by Freeman (1981) and supported by the supertree analysis of Jones et al. (2002) could not be adequately tested. However, the other members of this morphologically similar group (Fig. 1) did not cluster with Cheiromeles in any of the molecular analyses.

Unresolved phylogenetic position of Myopterus and Otomops.—A New World origin for Myopterus was proposed based on Freeman's (1981:119) analysis but because these bats possess an unusual combination of morphological characters, it has been difficult to determine the closest relative of this Old World genus. Our analysis is not conclusive because of missing data in this taxon, but Myopterus falls within the Old World clade in our combined analysis (Fig. 8). This poorly supported result as well as the position of Otomops in the Old World clade should be considered inconclusive at this point.

Divergence patterns in the family Molossidae.—Divergence of Molossidae from Vespertilionidae and Miniopteridae was estimated between 54 and 43 mya by Miller-Butterworth et al. (2007). Jones et al. (2005) and Teeling et al. (2005) had similar estimates for the origin of the family. Our analysis suggests that the family Molossidae originated even earlier (Paleocene, 63 mya [95% highest posterior density; 51, 75] ± 12) and fossil evidence suggests that at least 1 lineage represented by the fossil Wallia scalopidens was present in North America 42–44 mya (Legendre 1985; Storer 1984). A fossil history for this 20-million-year gap is unaccounted for and is consistent with calculations by Teeling et al. (2005) that the fossil record underestimates the origin of bat lineages by an average of 73%. The 2nd oldest molossid fossil, Cuviermops, is from 39 mya in France. Both of these early fossils show that the molossid group was already diverse in the middle Eocene. However, because our divergence date for the origin of the family was earlier than expected based on previous studies (Miller-Butterworth et al. 2007; Teeling et al. 2005), we conducted another BEAST analysis using the same minimum constraint for divergence of the molossid lineage as these studies did (37 mya). The divergence of molossids from the outgroup was estimated at 44.2 mya (data not shown) when the younger constraint was used. Thus, it appears that the older estimate (63 mya) for the origin of the family was a result of using different fossil priors.

Both a Eurasian origin (Teeling et al. 2005) and an African origin (Eick et al. 2005; Lim 2009) have been proposed for the family Molossidae, but conclusions were equivocal in these studies. Examination of our data cannot reject either hypothesis. In any case, based on fossil discoveries, dispersal to North America must have occurred before mid-Eocene. Dispersal to North America could have occurred along the common trans-Atlantic route used by several groups in the Paleocene–Eocene but Beringia also would have been a possible route at this time (Sanmartin et al. 2001). After dispersal to North America, subsequent global cooling and drying trends in the mid- to late Eocene (Zachos et al. 2001) could have contributed to isolation or extinction of some of these early molossid lineages. Persistence of molossids in Southeast Asia is likely because the severity of global cooling was not as great here during this time and tropical habitats persisted (Tsubamoto et al. 2004). However, without sequence data from Tomopeas and Mormopterus from the New World it is difficult to fully understand the origins of the family and to distinguish among dispersal hypotheses.

The next major event in the divergence of molossid groups is the split between the Old World and New World clades that took place around 28.7 mya (Fig. 8). One scenario that could describe this event is dispersal and subsequent radiation of a common ancestor of the New World clade from the Old World. This scenario also has been proposed for the divergence of Old World and New World emballonurid lineages (Teeling et al. 2005) around 30 mya via a “rafting” or “island hopping” event from Africa, as has been hypothesized for caviomorph rodents and primates (Flynn and Wyss 1998). Vicariance is another possible scenario. If lineages that were already in South America were separated from the Old World lineages due to cooling temperatures in the north, it would have reinforced their isolation. This hypothesis is consistent with the occurrence of the fossil of the now extinct Mormopterus faustoi from Brazil (Czaplewski et al. 2003; Legendre 1984b; Paulo Couto 1983) that shows molossid lineages were established in South America in the Oligocene. Around the Oligocene–Miocene boundary were a series of intermittent glaciations that subsequently could have promoted speciation and rapid divergence of lineages within the Old World and the New World clades (Zachos et al. 2001). In fact, there are at least 5 divergence events within each clade that occur close together in time and correspond to the Oligocene–Miocene boundary.

The close relationship of the North and South American species, T. brasiliensis, with Sauromys and T. aegyptiaca from Africa requires a biogeographic explanation. Morgan (1991) hypothesized a Neotropical origin for Tadarida and subsequent movement into North America at the Interchange in the late Pliocene. An alternate hypothesis (Czaplewski et al. 2003) supported by fossil evidence is that the genus Tadarida originated in Eurasia, immigrated to North America in the Miocene, and then moved into South America during the late Pliocene continental interchange. The 2nd hypothesis fits our topology and the divergence dates in our analysis for the separation of T. brasiliensis from the African sister taxa at 18 (14.1–21.9) mya. This lineage might have dispersed across the Bering land bridge but transoceanic dispersal from Africa cannot be ruled out, especially for the ancestor of a species that is known to be able to migrate hundreds of miles at very high altitudes and cover long foraging distances (Cockrum 1969; Williams et al. 1973). Similarly, it is not necessary to hypothesize arrival of this taxon in South America after the closing of the Panamanian isthmus. The presence of the hoary bat (Lasiurus cinereus) in Hawaii is further evidence that it is not necessary to invoke land bridges to explain the dispersal of bat lineages. Lim (2009) hypothesized 5 different dispersal events within Molossidae from Africa to South America. However, he acknowledged uncertainty with this hypothesis because of the lack of a comprehensive phylogenetic hypothesis for the family.

We recovered several species pairs with recent divergences (less than 2 mya [Fig. 8]). Three of these species pairs contained members from Madagascar and Africa (C. pumilus–leucogaster, M. condylurus–leucostigma, and Otomops martiensseni–madagascarensis). These results suggest at least 3 separate dispersal events by molossid bats to Madagascar and are consistent with previous work on these genera (Lamb et al. 2008; Ratrimomanarivo et al. 2009a).

Conclusions.—Our results generally are not congruent with traditional morphological hypotheses of relationship but are most similar to Legendre's (1984a) classification scheme. Simmons (1998) suggests that the 2 molossid subfamilies (Molossinae and Tomopeatinae) should be retained and an additional subdivision should be recognized at the tribal level.

We agree and therefore propose a subdivision of the subfamily Molossinae into 4 tribes: tribe Molossini containing the New World taxa (Molossus, Eumops, Molossops, Cynomops, Neoplatymops, Nyctinomops, and Promops), tribe Tadarini containing Old World taxa (Tadarida, Chaerephon, Mops, Platymops, Sauromys, Myopterus, and Otomops), tribe Cheiromelini, and tribe Mormopterini. Without additional data, it is impossible to know if subfamily Tomopeatinae is a valid separate lineage or alternatively is a member of Simmon's (2005) Molossinae. The only molecular analysis to include Tomopeas (Sudman et al. 1994) could not address this question because their study included only New World members of the family Molossidae. Until new material becomes available, this question will be difficult to answer.

An understanding of the evolutionary history of this group is essential for interpreting and understanding decades of comparative studies (cranial morphology, penile morphology, genetics, echolocation, ecology, reproductive physiology, and neurophysiology) that have been conducted in the absence of a rigorous phylogenetic framework (Arlettaz et al. 2000; Brown 1967; Gregorin 2003; Lim 2009; Reichart et al. 2010; Smith et al. 1986; Warner et al. 1974). Furthermore, a well-supported phylogeny is necessary for understanding factors affecting the prevalence, transmission, and history of diseases, as well as the emergence of new diseases such as Ebola virus that has been found to replicate in 2 molossid species (Swanepoel et al. 1996). We consider the phylogenetic hypothesis presented here to be the 1st step to understanding evolutionary patterns in this group. Undoubtedly, additional taxa and additional genes in future studies will clarify the hypotheses of relationship in Molossidae.