Previous understanding of the relationships among genera of bats in the family Molossidae was based largely on phenetic analyses of morphological data. Relationships among the genera of this family have not been tested with molecular data and, thus, the objective of this study was to construct a phylogeny of representative members of free-tailed bats using DNA sequence data from 1 mitochondrial locus (Nicotinamide adenine dinucleotide dehydrogenase subunit 1 [ND1]) and 3 nuclear loci (dentin matrix protein 1 exon 6 [DMP1], beta fibrinogen intron 7 [βFIB], and recombination activating gene 2 [RAG2]) for members of the subfamily Molossinae and outgroups from the families Vespertilionidae and Natalidae. Data for each gene were analyzed separately using maximum-likelihood and Bayesian methods and also analyzed in a single combined analysis of a total of 3,216 base pairs. Divergence times were estimated from the combined data set using BEAST analysis. Few intergeneric relationships were significantly supported by mitochondrial data; however, monophyly of most genera was supported. Nuclear results supported a Chaerephon–Mops clade; a New World clade consisting of Eumops, Molossus, Promops, Molossops (including Neoplatymops), Cynomops, and Nyctinomops; and a basal divergence for Cheiromeles. Divergence analysis suggested a Paleocene origin for the family and a split between molossids in the Old World and New World around 29 million years ago. Generally, relationships recovered in our analyses reflected biogeographic proximity of species and did not support the hypotheses of relationship proposed by morphological data.
The family Molossidae Gervais, 1856, is the 4th largest family of bats, containing approximately 100 species divided into 16 genera (Simmons 2005); 7 Old World genera (Mops, Chaerephon, Platymops, Sauromys, Cheiromeles, Otomops, and Myopterus), 7 New World genera (Promops, Molossus, Eumops, Nyctinomops, Molossops, Cynomops, and Tomopeas), and 2 genera with members found in both the Old World and New World (Mormopterus and Tadarida). Currently, 2 subfamilies are recognized (Simmons 1998, 2005)—Molossinae and Tomopeatinae. The subfamily Molossinae contains 15 genera; however, if Neoplatymops is treated as a separate genus (elevated from Molossops—Peterson 1965; Peterson et al. 1995; Willig and Jones 1985) there are 16 genera. The subfamily Tomopeatinae (Sudman et al. 1994) is monotypic and includes the enigmatic species Tomopeas ravus. The number of species in the family continues to grow; new Mops (Stanley 2008), Chaerephon (Goodman and Cardiff 2004), Mormopterus (Goodman et al. 2008; Reardon et al. 2008), Molossus (González-Ruiz et al. 2011), and 5 new species of Eumops (Baker et al. 2009; Eger 2007; Timm and Genoways 2004) have recently been described resulting in a total of 110 species in the family.
Many species in this family are similar in appearance and only can be externally distinguished by subtle anatomical features such as the extent of ear joining, shape of the antitragus, wrinkles on the lips, and the presence of rump bristles (Freeman 1981). They range in size from as large as 85- to 86-mm forearm in Eumops dabbenei and Cheiromeles torquatus (Corbet and Hill 1992; Eger 2007) to as small as 25- to 28-mm forearm in Mops spurrelli and Molossops temminckii (Freeman 1981). Molossid bats are found throughout tropical and subtropical regions of the world and their most obvious feature is a free tail that is not enclosed in a tail membrane like the tail of many other bats. They are distinguished from members of other bat families by numerous skeletal modifications, many of which are related to the narrow wing morphology (Miller 1907; Peterson et al. 1995).
Species of molossids are notably difficult to capture because of their tendency to fly high above the ground and thus are poorly represented in museum collections compared to other families of bats. For example, surveys in Peru (Hice et al. 2004), Brazil (Bernard and Fenton 2002), West Africa (Fahr and Ebigbo 2003), Myanmar (Struebig et al. 2005), and Ecuador (Lee et al. 2010) find that less than 1% of all mist-net captures are molossids even though they are documented from these areas. Furthermore, in a survey of 5 large collections of mammals (Royal Ontario Museum, Texas Tech University Museum, Museum of Vertebrate Zoology, United States National Museum, and the American Museum of Natural History), only 11.7% of the total specimens from the order Chiroptera were from the family Molossidae; whereas 40% were from the family Phyllostomidae and 23% were from the family Vespertilionidae.
Presumably as a result of their rarity in collections, few systematic analyses of the entire family Molossidae have been conducted and none take a molecular approach using modern methods of phylogenetic analysis. Freeman (1981) and Legendre (1984a) analyzed relationships in the family using morphological data before cladistic methods were widely applied to phylogenetic questions. External characteristics and skeletal measurements were analyzed by Freeman (1981) to generate several phenograms that depict the pattern of relationships for many of the genera in this family. Freeman's hypothesis (Fig. 1) described 2 major groups, which were Mormopterus (including Sauromys and Platymops), Myopterus, Cheiromeles, Molossops (including Neoplatymops and Cynomops); and Tadarida, Chaerephon, Mops, Otomops, Nyctinomops, Promops, Molossus, and Eumops. In addition, Legendre (1984a) examined dental morphology and divided the family into 3 subfamilies that were not congruent with Freeman's (1981) groups—Molossinae (Molossus, Eumops, Molossops, Cynomops, Neoplatymops, Myopterus, and Promops), Tadarinae (Tadarida [including Chaerephon and Mops], Mormopterus [including Platymops, Sauromys, and Micronomus], Nyctinomops, Otomops, and Rhizomops), and Cheiromelinae (Cheiromeles).
More recently, Gregorin (2000) conducted the 1st cladistic analysis of 112 morphological characters from 76 species of molossids. He recognized 2 subfamilies, Tomopeatinae and Molossinae, but within Molossinae the phylogenetic relationships (Fig. 1) and classification differ substantially from Freeman (1981) and Legendre (1984a). Within Molossinae, he recognized 2 tribes, Molossini (Molossus, Promops, Cheiromeles, and Myopterus) and Tadaridini (all other genera). Clearly there is a lack of consensus regarding the classification and evolutionary history of this family.
The evolutionary history of bat families such as Vespertilionidae (Hoofer and Van Den Bussche 2003; Roehrs et al. 2010), Phyllostomidae (Baker et al. 2003), Emballonuridae (Lim et al. 2008), Pteropodidae (Colgan and da Costa 2002), Natalidae (Dávalos 2005), and Mormoopidae (Lewis-Oritt et al. 2001) has been elucidated with molecular phylogenetic approaches using multiple genes, but molecular data for relationships within the family Molossidae are sparse. Sudman et al. (1994) used cytochrome-b sequences and protein electrophoretic data to demonstrate that the genus Tomopeas was a molossid instead of a member of the Vespertilionidae (Miller 1907) but their taxonomic sampling was limited to 6 molossid genera, all from the New World tropics. Although numerous phylogeographic studies have been conducted in recent years on molossids, especially for African species of Otomops (Lamb et al. 2006, 2008), Chaerephon (Jacobs et al. 2004; Ratrimomanarivo et al. 2009a; Taylor et al. 2009), Mormopterus (Ratrimomanarivo et al. 2009b), and Mops (Ratrimomanarivo et al. 2007), the higher-level relationships remain unclear. The objective of our study was to examine relationships among the 16 genera in the subfamily Molossinae using molecular data based on both nuclear and mitochondrial DNA sequences.
Materials and Methods
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 MgCl2, 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).
Overall, we obtained sequences from 60 individuals representing 35% (38 of 110) of the species in the subfamily Molossinae (sensu Simmons 2005; Appendix I). Although amplification of all genes was attempted for all samples by adjusting concentrations of polymerase chain reaction reagents and altering thermal profiles, some genes did not amplify for some taxa. Therefore, the resulting analyses were not identical with regard to their taxon sampling. In total, the 4 genes comprised 3,216 aligned nucleotides. The concatenated alignment included 30 ingroup and 4 outgroup taxa and 15 taxa were represented by composite sequences from 2 different individuals of the same species.
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%.
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.
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.
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.
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.
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.
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).
The evolutionary history of the molossid subfamily Molossinae recovered in our study based on DNA sequences is not consistent with existing morphological predictions (Freeman 1981; Gregorin 2000; Legendre 1984a). We did not recover groups that match Legendre's (1984a) proposed subfamilies or Gregorin's (2000) 2 tribes (Molossini and Tadaridini), nor did any of the molecular hypotheses from any gene mirror the 2 groups described by Freeman (1981). For instance, Cheiromeles formed the sister lineage to the remaining members of Molossinae and did not cluster with any of the taxa (Mormopterus, Molossops, Molossus, and Promops) as hypothesized based on cladistic analysis of morphology (Freeman 1981; Gregorin 2000; Fig. 1). However, our molecular phylogeny did appear to share some relationships recovered in Freeman's (1981:102) phenogram. Furthermore, the recovery of a New World clade in this study was incongruent with traditional taxonomic schemes; neither morphological arrangement predicted this grouping. Legendre's (1984a) subfamily Molossinae is similar to our New World clade, but he did not include Nyctinomops and did include Myopterus from the Old World in this subfamily.
As has been observed in other studies, geography seemed to better indicate phylogenetic relationship (such as the New World clade and Old World clade in Fig. 7) than did morphological traits (Hoofer and Van Den Bussche 2003; Ruedi and Mayer 2001). Disagreements with morphological hypotheses could be due to convergence in morphological characteristics that traditionally have been used to classify molossids. For example, Freeman (1981) cautioned that characters such as palatal emargination, separation of the ears, basisphenoid pits, and lip wrinkles might have more to do with functional morphology of eating beetles and might not be ideal for recovering phylogenetic patterns. Furthermore, the morphology associated with flat-headedness (Freeman 1981; Peterson 1965) appears to have evolved multiple times in crevice-dwelling molossids. The genera Sauromys, Platymops, Mormopterus, and Neoplatymops have a high coronoid process and a mandibular condyle elevated above the toothrow. These characteristics could have evolved to compensate for the leverage required when the temporal muscle was forced to extend laterally with the flattening of the skull (Freeman 1981) and might not be useful synapomorphies. As the phylogeny of this group is refined, a better understanding of character evolution should emerge.
In summary, although there was slightly different taxonomic sampling for each gene, information from 4 different loci corroborated several hypotheses of intergeneric relationship. ND1 data resolved some recent divergences, but did not provide much resolution for the deeper divergences due to saturation. Nevertheless, the topology based on ND1 analysis was not in conflict with nuclear hypotheses of relationship. Our analyses significantly support a Mops–Chaerephon clade, a New World molossid clade, and a close relationship between Sauromys and Tadarida. Additionally, analysis of the combined data set suggests an Old World clade exists.
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 , Gregorin , Legendre [1984a], and Simmons ) 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.
We extend our deepest gratitude to the following individuals and institutions for their contributions of tissue loans and specimen data, without which this project would not have been possible: F. Anwarali, R. Baker, and H. Garner (Natural Science Research Laboratory, Texas Tech University); P. Benda (National Museum of Prague); J. Bickham (Purdue University); A. Bishop and R. Dowler (Angelo State Natural History Collections, Angelo State University); A. Cleveland (California Baptist University); J. Eger (Royal Ontario Museum); G. Eick and D. Jacobs (University of Stellensbosch); J. Fahr (University of Ulm); B. Patterson, J. Phelps, S. Goodman, and W. Stanley (Field Museum); D. Dittmann and M. Hafner (Lousiana State University Museum of Zoology, Louisiana State University); A. Gardner and K. Helgen (National Museum of Natural History); T. Kunz (Boston University); T. Lee (Abilene Christian University Museum of Natural History, Abilene Christian University); S. McLaren (Carnegie Museum of Natural History); J. Patton (Museum of Vertebrate Zoology, University of California Berkeley); S. Rossiter (Queen Mary, University of London); P. Taylor (Durban Natural Science Museum, Durban, South Africa); R. Timm (University of Kansas Natural History Museum, University of Kansas); and R. Van Den Bussche (Oklahoma State University). We thank the Economic Planning Unit of the Malaysian Prime Minister's Department, Department of Wildlife and National Parks, Universiti Malaysia Sarawak, M. T. Abdullah, and the Sarawak Forest Corporation for permission to use samples in this study. We also acknowledge the following students for help in the laboratory: R. Dolman, G. Hernandez, and M. McDonough. Special thanks to J. Apodaca and R. Rodriguez for use of their sequences in this project; to J. Lack for help with BEAST analyses; and to P. Freeman, M. Dixon, and M. McDonough for their input. Finally, this project was funded by an Angelo State University Faculty Innovation Grant, a Tri-Beta Research Grant, a Texas Academy of Science Grant, and a Carr Student Research Grant.