OUP user menu

Population Genetic Analysis of Myzopoda (Chiroptera: Myzopodidae) in Madagascar

Amy L. Russell , Steven M. Goodman , Isabella Fiorentino , Anne D. Yoder
DOI: http://dx.doi.org/10.1644/07-MAMM-A-044.1 209-221 First published online: 19 February 2008


The chiropteran family Myzopodidae is endemic to Madagascar and is characterized by several unique morphologies, such as sessile adhesive discs on the thumb and sole. A new species, Myzopoda schliemanni, was recently described from western Madagascar that is morphologically distinct and geographically disjunct from the eastern species, M. aurita, the only other member of this family. Geographic variation within Myzopoda has only recently been studied at the morphological level and has never been addressed at the genetic level. We used a combination of phylogenetic, coalescent, and population genetic analyses to characterize the speciation history of Myzopoda and to clarify current and former patterns of gene flow within and between Myzopoda. Mitochondrial DNA sequences were used to determine whether genetic data support the morphologically distinct species M. schliemanni, to infer the distribution of the common ancestor of extant Myzopoda, to estimate effective population sizes (Ne) and levels of migration between species, and to determine patterns of population structure within species. Phylogenetic and network analyses revealed the existence of 4 well-supported clades in Myzopoda, but could not resolve relationships among those clades. Divergent haplotypes within species may result from either recent gene flow between the 2 species or more likely from incomplete lineage sorting. Multiple coalescent-based methodologies produced concordant estimates of Ne for Myzopoda, but conflicting signals for migration between the species, probably reflecting differences in the underlying models used by the methods. We found significant genetic structure within M. aurita, but no correlation with geography. This pattern may result from recent gene flow facilitated by expansion of Ravenala stands, an important day-roost tree for Myzopoda, associated with anthropogenic deforestation and the opening up of new habitat for members of this genus.

Key words
  • Chiroptera
  • hybridization
  • incomplete lineage sorting
  • Madagascar
  • mitochondrial DNA
  • Myzopoda
  • Ravenala

Madagascar is renowned for its extreme levels of diversity and endemicity. Goodman et al. (2003) note an astonishing 100% endemicity for native terrestrial mammalian species, and similarly high levels exist for reptiles (95%—Raxworthy 2003) and amphibians (99%—Glaw and Vences 2003). Additionally, the numbers of species within the represented groups can be quite high. For example, more than 15% of living primate species only occur on Madagascar, an area that represents less than 0.4% of the total land surface on the Earth (Karanth et al. 2005).

The chiropteran fauna of Madagascar stands in contrast to terrestrial mammals with regard to high levels of diversity and endemicity. The island's bat fauna is not notably speciose, especially when compared to other Old World tropical islands, with 37 species as of early 2007. Comparing Madagascar with New Guinea and Borneo, Eger and Mitchell (2003) found that Madagascar, with 74–80% the land area of the other islands, respectively, had only 37–41% of the documented bat species. Recent field efforts and systematic studies of Malagasy bats are increasing the recognized diversity levels for the island, and 7 species new to science have been described since the study by Eger and Mitchell published in 2003 (Bates et al. 2006; Goodman and Cardiff 2004; Goodman et al. 2005, 2006a, 2006b, 2007a, 2007b). Although not as speciose as comparably sized Old World tropical islands, the level of endemicity observed among Malagasy Chiroptera (70%) far exceeds that of New Guinea (16%) or Borneo (6.5%). One of the more extraordinary examples from Madagascar is the endemic bat family, the Myzopodidae, which until lately was thought to be monotypic. Recently, Goodman et al. (2007a) described a new species, Myzopoda schliemanni, from the lowland areas of central-western Madagascar. This new species was distinguished from the eastern M. aurita based on a number of cranial, morphological, and pelage coloration characteristics. Myzopoda has long been considered an especially intriguing genus because of its remarkable morphological adaptations, its evolutionary uniqueness, its level of endemicity, and its supposed rarity. It possesses disk-shaped structures on its thumbs and feet, which are used to adhere to broad, smooth leaves such as those of Ravenala madagascariensis (family Strelitziaceae-Göpfert and Wasserthal 1995; Schliemann and Maas 1978), as well as cave walls (Kofoky et al. 2006). These structures and the associated behavior are similar to that of the neotropical genus

Thyropter a, although histological examinations of the discs indicate that they are convergent (Schliemann 1970, 1971).

Long considered a member of the superfamily Vespertilio-noidea, recent molecular work has instead supported the Myzopodidae's position as the basal member of the super-family Noctilionoidea (Teeling et al. 2005). Its closest extant relative is the New Zealand endemic family Mystacinidae, and fossil-calibrated molecular phylogenies place the origin of Myzopodidae at 52 (range 57–46) million years ago (mya), in the early Eocene (Teeling et al. 2005). Most likely due to the long evolutionary isolation of Myzopodidae, the place of origin of this family has been unresolved (Hoofer et al. 2003; Teeling et al. 2005). A fossil humerus from the early Pleistocene of Olduvai is described by Butler (1978:66) as agreeing “… very closely, except for its larger size, with that of Myzopoda aurita.” Beyond this short description and reference to an unpublished specimen, we are unaware of any recent evaluation of this material. Schliemann and Goodman (2003) expressed some doubt regarding the evolution of the Myzopodidae on the African mainland, and instead suggested that the family evolved in situ on Madagascar. Although the place of origin of Myzopodidae is obscure, the recent description of M. schliemanni, which is ecologically and geographically disjunct from M. aurita, presents questions regarding the mode of speciation and patterns of gene flow among and within living species of Myzopoda.

Over the past decade, members of this genus have been captured at a number of localities on the island. The vast majority of these records are of M. aurita in lowland eastern areas, in or near marshlands and rain-forest settings, and a handful of reports of M. schliemanni in lowland western sites with locally mesic conditions (Goodman et al. 2007a). Myzopoda appears to be well adapted to anthropogenically disturbed habitat. This adaptation may derive from Myzopoda' s association with R. madagascariensis, a common species in lowland disturbed habitats in the east and more mesic zones in the west (Blanc et al. 2003). This connection was 1st observed by Harry Hoogstraal in 1947, as an unpublished record of Myzopoda roosting on a Ravenala leaf (Schliemann and Maas 1978). The importance of this plant–bat co-occurrence can be inferred from the markedly similar upper elevation limits of both species in the eastern portion of the island: 1,000–1,100 m for Ravenala and less than 1,000 m for M. aurita (Goodman et al. 2007a).

In their species description, Goodman et al. (2007a) present 2 hypotheses for the divergence of M. aurita and M. schliemanni, both predicated on the assumption that M. schliemanni is derived from M. aurita. The 1st scenario involves a westward dispersal directly across the Central Highlands, perhaps at the approximate latitude of Antananarivo. This was considered an unlikely scenario, because the elevation in this portion of the island, which is now largely deforested, is notably greater than 1,000 m, and Ravenala is largely absent for a span of several hundred kilometers across the natural flora of the Central Highlands. However, a cataclysmic dispersal of several individuals might have followed this route, perhaps via cyclonic winds. A 2nd scenario involves westward stepping-stone dispersal across a lower-elevation region between the Central and Northern Highlands (Fig. 1). The presence of areas less than 1,000 m in elevation and a patchy distribution of Ravenala make this route more plausible than the Central Highland hypothesis. Further, the known distribution of M. schliemanni is in the region to the north and south of Mahajanga, the zone directly to the west in the gap between the Central and Northern Highlands. To address the question of the speciation history of Myzopoda and to understand former or current patterns of gene flow within and between M. aurita and M. schliemanni, as they have been defined on morphological ground, we have used mitochondrial DNA (mtDNA) sequence data from several populations of both species. We use a combination of phylogenetic, coalescent, and population genetic analyses to determine whether genetic data support the recognition of M. schliemanni based on morphological grounds, to assess whether the most recent common ancestor (MRCA) of M. aurita and M. schliemanni was distributed in eastern or western Madagascar, to estimate effective population sizes of and levels of migration among species, and to determine patterns of population structure within species. This work represents the first population-level genetic study of Myzopoda and is an important step toward understanding the evolution of this ancient and long-isolated family.

Fig. 1

Map of Madagascar, with sampling localities for Myzopoda. The general limits of the Central and Northern Highlands are shown.

Materials and Methods

Sample collection.—Specimen collection (see Table 1) followed protocols detailed by Goodman et al. (2007a), which followed guidelines of the American Society of Mammalogists (Gannon et al. 2007). Total genomic DNA was isolated using a DNeasy DNA isolation kit (Qiagen, Valencia, California), and stored in the provided elution buffer. Approximately 400 base pairs (bp) of the mitochondrial D-loop was amplified using the primers P(mt) and F(mt) (Wilkinson and Chapman 1991). The complete mitochondrial cytochrome-b gene (Cytb) was amplified using the primers pairs L14724/H15506 and L15171/H15915 (Irwin et al. 1991; Yoder et al. 1996). Polymerase chain reactions were performed in 50-μl reaction volumes containing 2.25 mM MgCl2, 0.25 mM deoxynucleo-side triphosphates, 2.5 U Taq polymerase (Promega, Madison, Wisconsin), 5 μl Promega 10X buffer, 5 μl genomic DNA, and 20 pmol of each primer. The amplification involved an initial denaturation at 94°C for 3 min, followed by 40 cycles of 94°C for 45 s, 52°C for 45 s, and 72°C for 1 min, with a final elongation step at 72°C for 4 min. The polymerase chain reaction product was then purified using gel band excision (Gel Excision kit, Qiagen), polymerase chain reaction purification (PCR Purification kit; Qiagen), or enzymatic treatment (ExoSAP-IT; USB, Cleveland, Ohio).

View this table:
Table 1

Sample information for Myzopoda. SF = Station Forestière.

M. auritaToliaraTolagnaro, Sainte Luce124°49′S47°10′;E
M. auritaToliaraForêt d'Analalava124°13′S47°19′E
M. auritaFianarantsoaKianjavato521°22′S47°52′E
M. auritaToamasinaTampolo417°17′33”S49°24′30”E
M. auritaToamasinaFoulpointe117°41′S49°31′E
M. schliemanniMahajangaForêt d'Andranomanintsy2316°31′27”S44°29′23”E
M. schliemanniMahajangaSF d'Ampijoroa116°18′36”S46°48′36”E

The purified fragment was sequenced from both directions using the BigDye terminator cycle sequencing kit version 3.1 (Applied Biosystems, Foster City, California) in a 20-μl reaction containing 2 ul of the BigDye solution, 1.2 μl ABI 5X sequencing buffer, 5 pmol of primer, and 100–200 ng of purified polymerase chain reaction product. The temperature profile used was according to the manufacturer's instructions. The sequencing reactions were cleaned of unincorporated nucleotides using either genCLEAN (Genetix, New Milton, Hampshire, United Kingdom) or sephadex (Princeton Separations, Adelphia, New Jersey) dye terminator removal plates, and analyzed on an ABI 3100 automated sequencer (Applied Biosystems). Each individual was sequenced 4–15 times per locus from multiple polymerase chain reactions to resolve any ambiguities present in single sequencing passes. We used Sequencher version 4.2 (Gene Codes Corp., Ann Arbor, Michigan) to assemble and edit consensus sequences for each individual, and aligned the data by eye using MacClade version 4.0 (Maddison and Maddison 2000). All sequences were submitted to GenBank (Cytb: accession numbers DQ178324DQ178334 and EF43218932213; D-loop: accession numbers DQ178335DQ178345 and EF432214EF432238).

Phylogenetic analyses.—We performed Bayesian and maximum-parsimony phylogenetic analyses on the Cytb data using MRBAYES version 3.0b4 (Huelsenbeck and Ronquist 2001) and PAUP* version 4.0b 10 (Swofford 2002), respectively. Homologous portions of whole mitochondrial genomes from Mystacina tuberculata (family Mystacinidae; GenBank accessions AY960981 and NC006925) were included as outgroups for the analyses. Maximum-parsimony analyses were performed using heuristic searches with tree-bisection-reconnection branch swapping, and setting the maximum number of saved trees to 1,000. Parsimony bootstrap analyses were executed using heuristic searches, tree-bisection-reconnection branch swapping, and 100 replicates of the random addition search option. Bayesian analyses were performed with flat priors, running 4 chains of 10 million generations each, with sampling every 500 generations. The chains were heated using the temperature scaling factor T = 0.2. We specified an HKY+F model with 4 rate categories for the gamma distribution, the model having been selected using MODELTEST version 3.07 (Posada and Crandall 1998). After examining the likelihood profile using TRACER version 1.1 (Rambaut and Drummond 2004), we discarded the first 4,000 trees as a burn-in and constructed a 50% majority consensus of the remaining 16,000 trees in PAUP* (Swofford 2002). We calculated maximum-likelihood branch lengths for the Bayesian analysis, using the 50% majority-rule consensus tree as a constraint. All trees were midpoint rooted.

Although previous phylogenetic analyses of Chiroptera have placed the New Zealand endemic Mystacinidae as the closest extant relative to the Myzopodidae (Teeling et al. 2005), the branches connecting these 2 families are quite long. We therefore performed further outgroup analyses to determine whether the long branch observed between the Myzopodidae and the Mystacinidae was obscuring relationships among basal clades recovered within the Myzopodidae. We simulated a random outgroup sequence in Seq-Gen (Rambaut and Gras sly 1997), using the same evolutionary model parameters observed in the data set for Myzopoda. We then performed a maximum-parsimony analysis as described above, using this random sequence as an outgroup.

We performed a series of 3 parsimony bootstrap analyses to determine the statistical rigor of hypothesized relationships among basal clades revealed in the phylogenetic analyses (Carstens et al. 2004; Ruedi et al. 1998). For each analysis, we first conducted a maximum-likelihood search in PAUP* version 4.0b 10 (Swofford 2002) for the optimal unconstrained tree, using the evolutionary model and parameters estimated with MODELTEST version 3.07 (Posada and Crandall 1998). We then conducted a constrained maximum-likelihood search using the least-resolved constraint tree consistent with the hypothesis being tested (e.g., a test of reciprocal monophyly between the 2 species of Myzopoda specified 2 internally unresolved clades containing morphologically defined specimens of M. aurita or M. schliemanni) and specifying the same model and parameter estimates as in the unconstrained search. To determine whether these model and parameter estimates were optimal under the given constraint tree, we then ran the data through MODELTEST, using the constraint tree as a backbone for computing the likelihood of the data under each model. We then conducted a 2nd constrained likelihood search using the model and parameter estimates specified in the constrained MODELTEST analysis. The difference in likelihood scores between this tree and the unconstrained tree was used as the test statistic in the parametric bootstrap. Where the constrained maximum-likelihood search yielded more than a single most likely tree, we used TreeEdit version 1.0 (Rambaut and Charleston 2002) to identify the trees with the fewest possible resolutions. Of these trees, 1 was arbitrarily chosen and was fully resolved in TreeEdit. We then used Seq-Gen (Rambaut and Grassly 1997) to simulate 100 DNA sequence data sets on this resolved maximum-likelihood constraint tree. These simulations utilized the same model and parameter estimates as the constrained maximum-likelihood search. We then used PAUP to analyze the simulated data, with and without constraints. The test statistic described above was compared with the distribution of the difference in likelihood scores between constrained and unconstrained searches for the simulated data to determine the significance of the constraint hypothesis.

Coalescent and population genetic analyses.—For each locus, we used MODELTEST version 3.07 to determine the best-fitting model of evolution. Using the prescribed model and parameter values, we utilized Arlequin version 2.001 (Schneider et al. 2000) to test for genetic structure among geographic regions and sampling sites. Parsimony networks were constructed for the combined dataset using TCS version 1.18 (Clement et al. 2000). Pairwise genetic distances were calculated in PAUP* version 4.0b 10 among Cytb and D-loop haplotypes under the Tamura-Nei + I (Tamura and Nei 1993) and HKY+I (Hasegawa et al. 1985) models, respectively.

Using the combined mitochondrial data set, we estimated demographic parameters for the Myzopoda system utilizing the isolation with migration model under a changing population size as implemented in the IM software package (Hey 2005; Hey and Nielsen 2004). For a given pair of populations (in this case, M. aurita and M. schliemanni), IM will estimate the diversity parameter 6 (=Neμ) for each population and for their MRCA. Other parameters in the model include the divergence time t, directional migration rates mji and mji, and the population splitting parameter s. This last parameter specifies the proportion of the ancestral population that evolved into 1 of the daughter populations (in this case, M. aurita). Although this software permits the violation of the infinite sites model, it is more powerful when there are no such violations in the data. Therefore, to make our data compatible with this stricter model, we pruned by hand sites from the data set that possessed an insertion–deletion (indel). We further pruned the data using IMgc (Woerner et al. 2007) to exclude both sites and individuals to produce the most information-rich contiguous DNA sequence segment that passes the 4-gamete test. Although the 4-gamete test is typically used to detect recombination in nuclear data, this test also can be used to detect homoplasy in mtDNA; therefore, our application of IMgc to these mtDNA data yielded a data set that is more powerful, albeit slightly smaller, for estimating population parameters. The resulting data set for these analyses included 902 bp of Cytb sequence for 10 M. aurita and 22 M. schliemanni and 292 bp of D-loop sequence, analyzed together as a single locus, for 9 M. aurita and 21 M. schliemanni. Results were obtained from 3 independent replicate runs of IM, each run using 20 Markov chains of at least 8 × 106 steps with a burn-in of 105 steps. Geometric heating was implemented between the chains, with the heating parameters gl and g2 set to 0.9 and 0.8, respectively. Parameter conversions were calculated using the geometric mean of the mutation rates of Cytb (4 × 10−8 substitutions site−1 year−1Hulva et al. 2004) and the D-loop (1 × 10−5 substitutions locus−1 year−1Russell et al. 2005) and a generation time of 11 months (estimated; S. M. Goodman, pers. comm.). Analyses were assumed to have converged on the stationary posterior distributions of the parameters if the effective sample size for each parameter was greater than 50 and if 3 independent runs generated similar parameter estimates.

Demographic parameters were also estimated for each species using the program LAMARC version 2.0.2 (Kuhner et al. 2005). This analysis utilized the combined mitochondrial data set to estimate θ and directional migration rates between species under a Bayesian framework. We assumed an F84 model of evolution (Felsenstein 1984) with a transition : trans version ratio (κ) = 9. The analysis utilized a single long chain of 20,000 samples with a sampling interval of 50 steps and a burn-in of 2,000 samples. This long chain was divided into 3 subchains with differential heating (temperatures 1.0, 1.2, and 1.8) to search the parameter space more thoroughly. The analysis was repeated 3 times with a different random seed to verify the accuracy of the results via convergence.


Phylogenetics.—Maximum-parsimony and Bayesian phylo-genetic analyses of the Cytb data are in complete concordance regarding clades with more than 50% support (Fig. 2). Both analyses support the division of Myzopoda into 4 well-supported clades, here named clades A–D. These same 4 clades also are evident as well-differentiated clusters of related haplotypes in statistical parsimony network analyses of the combined Cytb and D-loop data (Fig. 3). Two of these clades (A and D) contain all individuals morphologically consistent and codistributed with M. schliemanni, and the other clades (B and C) contain all individuals of M. aurita. Although the relationships among these clades are not able to be resolved with any confidence in phylogenetic analyses using either DNA (Fig. 2) or translated amino acid characters (not shown), the network indicates that M. schliemanni is paraphyletic with respect to M. aurita, with the D clade more closely related to the C clade of M. aurita than to other M. schliemanni haplotypes.

Fig. 2

Phylogeny of Myzopoda. The topology shown is from a Bayesian analysis of cytochrome-b gene sequences, and support measures are given as Bayesian posterior probabilities above the branch and maximum-parsimony bootstrap values below the branch. Numbers in parentheses indicate the number of individuals represented by that haplotype, and further details on the specific specimens involved can be found in Fig. 3. Identified clades (A, B, C, and D) are assigned to species as follows: clades A and D are M. schliemanni, and clades B and C are M. aurita.

Fig. 3

Maximum-parsimony network of combined mitochondrial data for Myzopoda. Haplotypes are shown as circles with the individual labels noted inside. Missing inferred haplotypes are shown as small circles. Connections between haplotypes are shown with the number of mutations noted. Clusters of related haplotypes are grouped within bold squares, with the clade name A–D.

A phylogram using maximum-likelihood branch lengths illustrates the long-branch relationship between Myzopoda and the outgroup Mystacina. Although the branch lengths within Myzopoda are consequently very small relative to that with the outgroup, the Bayesian phylogeny (Fig. 4) does indicate some structure among the 4 clades identified above, with clade D being the most basal, although it is not statistically significant. A random outgroup provided some slight statistical support (bootstrap = 53%) for the placement of clade B as the most basal within Myzopoda (figure not shown).

Fig. 4

Phylogram of Myzopoda, using Mystacina as an outgroup. The topology of the tree was derived by Bayesian analysis, and the branch lengths by maximum likelihood.

We performed a series of parsimony bootstrap analyses to further explore the statistical rigor supporting hypothetical relationships among clades within Myzopoda. Specifically, we tested a set of 3 hypotheses: reciprocal monophyly of Myzopoda species; monophyly of M. schliemanni; and monophyly of a clade containing the clades C and D named above, as suggested by the network analysis. Enforcing the monophyly of M. schliemanni increased the likelihood score of the Myzopoda tree by 0.90914, enforcing the monophyly of a combined CD clade increased the likelihood score by 1.49053, and enforcing the reciprocal monophyly of Myzopoda species increased the likelihood score by 1.64645. After comparing these differences to their respective null distributions, we could not reject any hypothesis of monophyly (P = 0.21 for monophyly of M. schliemanni; P = 0.14 for reciprocal monophyly of Myzopoda species; P = 0.08 for monophyly of a combined CD clade).

Distance analyses.—We calculated average pairwise genetic distances within and between morphologically defined species (M. aurita and M. schliemanni) and phylogenetically defined clades (A, B, C, and D) for both Cytb (Table 2) and the D-loop (Table 3). Both loci show a greater distance between species (Cytb: 2.384%; D-loop: 7.100%) than within species (Cytb: 1.426% and 1.245%; D-loop: 2.634% and 2.976%). Intraclade distances are, by definition, quite low, with ranges from 0.319% to 0.593% for Cytb and from 0.846% to 1.155% for the D-loop. Average interclade distances are roughly similar for Cytb sequences and range from 2.190% between clades A (M. schliemanni) and C (M. aurita) to 3.211% between clades A and D of M. schliemanni. However, average interclade distances at the D-loop locus seem to fall into 2 classes. Distances between clades B and C of M. aurita (= 4.290%) and between clades C and D (= 4.469%) are similar, and both are significantly less than other interclade distances (range from 7.236% between clades B and D to 8.301% between clades A and D).

View this table:
Table 2

Genetic distances (%) between cytochrome-fr haplotypes for Myzopoda from Madagascar. Distances are calculated under a Tamura-Nei model with invariant sites (I = 0.7757).a

USNM 4488860.094
FMNH 1876210.0940.190
UADP 2260.9000.9960.996
UADBA RJB2000.7910.8870.8870.094
USNM 5770651.9322.0452.0452.2862.168
FMNH 1792731.9352.0462.0462.5522.4210.384
FMNH 1654542.0462.1592.1592.6642.5340.4800.094
FMNH 1773271.7141.8251.8252.0611.9452.2882.2952.410
FMNH 1876031.6061.7181.7181.9451.8322.1722.1752.2900.291
FMNH 1876182.0512.1692.1692.4012.2862.4012.4052.5220.6780.577
FMNH 1876152.1612.2782.2782.5202.4012.5202.5282.6450.7830.6780.094
FMNH 1876162.3992.5202.5202.5292.6433.0123.0123.1380.7760.8770.6760.777
FMNH 1876172.2792.3992.3992.4082.5222.8862.8883.0120.6760.7770.5770.6780.094
FMNH 1876102.0622.1822.1822.4062.2942.6512.6452.7680.6760.5770.5770.6780.4860.392
FMNH 1876082.2572.3752.3752.8502.7273.0922.8512.7312.8532.7323.2233.3463.6113.4783.234
  • a FMNH = Field Museum of Natural History; USNM = National Museum of Natural History; UADP = Universit&#00E9; d'Antananarivo, D&#00E9;partement de Pal&#00E9;ontologie; UADBA = Universit&#00E9; d'Antananarivo, D&#00E9;partement de Biologie Animale.

View this table:
Table 3

Genetic distances (%) between the D-loop haplotypes for Myzopoda from Madagascar. Distances are calculated under a HKY model with invariant sites (I = 0.6896).a

FMNH 1876650.557
FMNH 1876670.8400.850
FMNH 1876200.2770.8401.125
FMNH 1876170.8400.8500.5591.125
FMNH 1876161.1351.1450.2771.4230.848
USNM 5770656.6127.3847.8536.9716.9728.318
FMNH 1792737.0607.8538.3447.4277.4248.8380.848
FMNH 1839906.4427.2247.6966.8016.7768.2164.0493.696
USNM 4488866.8757.6748.1647.2417.2158.7044.4094.0500.278
UADP 2266.2967.0657.5136.6506.6367.9954.7085.0881.7372.048
UADBA RJB2006.3897.1657.6266.7456.7268.1274.0314.3971.1461.4480.557
USNM 4489306.9027.7038.1987.2697.2388.7504.4184.0580.2780.5612.0501.449
FMNH 1876087.5098.3238.8397.8877.8699.3807.4407.0324.0404.4035.4594.7553.689
  • a FMNH = Field Museum of Natural History; USNM = National Museum of Natural History; UADP = Université  d'Antananarivo, Département de Paléontologie; UADBA = Université d'Antananarivo, Département de Biologie Animale.

Coalescent analyses.—Analyses using the isolation with migration model under a changing population size resulted in a demographic history with well-supported estimates for most parameters. Multiple runs converged on similar estimates for all parameters, and results are presented in Table 4. Modes of posterior probability estimates of 6 are approximately equal between M. auritaA = 9.239; 90% highest posterior density [HPD] = 0.722, 56.992) and M. schliemanniS = 9.370; 90% HPD = 0.415, 24.420; Fig. 5a). Although estimates of θ were significantly greater for the common ancestor of M. schliemanni and M. aurita than for either of its contemporary daughter species (θAS = 110.677; 90% HPD = 59.806, 227.703), significant error is attached to our estimate of this parameter, and the mitochondrial data utilized here seem unable to localize on an informative posterior distribution. These analyses yielded a distinct peak for the divergence time parameter t at 1.397 (90% HPD = 0.283, 16.890); however, there was some difficulty in identifying the upper limit of the posterior distribution (Fig. 5b). The data support an estimate of 0.0005 for the population splitting parameter s, although the 90% HPD covers a rather large interval from 0.0005 to 0.9995 (Fig. 5c). Furthermore, estimates of migration under the IM model are near zero in both directions (mAS = 0.01; 90% HPD = 0.01, 1.557; mSA = 0.01; 90% HPD = 0.01, 1.063; Fig. 5d).

View this table:
Table 4

Demographic parameter estimates from coalescent analyses for Myzopoda from Madagascar. Directional migration estimates (mij are given as migration rates per generation from species i into species j.a

M. schliemanniM. auritaMRCA
θNeθNeθNet (years)smASmSA
Run 18.626123,9419.695139,297112.1211,610,94084,7370.00052x10‼72x1(T7
Run 210.429149,8417.642109,799117.8241,692,87981,5790.00052x1(T72x10"7
Run 39.056130,10810.380149,131102.0841,466,72654,2110.00052x10∼72x1(T7
Run 10.00759,8110.01299,89630.55231.454
Run 20.00867,0400.01198,39545.07040.839
Run 30.00757,5080.012101,87135.42330.764
  • a MRCA = most recent common ancestor; t = divergence time parameter; s = population splitting parameter; A = Myzopoda aurita; S = Myzopoda schliemanni.

Fig. 5

Posterior probability densities for demographic parameters from IM. Results for each parameter are shown from 3 independent runs, a) Posterior distribution for θ. Estimates for Myzopoda aurita are shown in black, those for M. schliemanni are shown in dark gray, and those for the most recent common ancestor are shown in light gray, b) Posterior distributions for divergence time, c) Posterior distributions for the population splitting parameter s. d) Posterior distributions for migration rates. Estimates for mAS are shown in black, and those for mSA are shown in gray.

Analyses utilizing the LAMARC software implemented a less parameter-rich model, estimating only θ and migration for the contemporary species (Table 4; Fig. 6). A unimodal posterior distribution with narrow peaks was obtained for both 0 parameters (Fig. 6a). Maximum posterior estimates for M. auritaA = 0.012; 90% confidence interval [90% CI] = 0.006, 0.023) were greater than those for M. schliemanniS = 0.007, 90% CI = 0.004, 0.011), although the 90% CIs overlapped substantially, qualitatively matching the results of the IM analyses. Under the LAMARC model, examination of the data indicates significant amounts of migration in both directions (mAS = 37.015; 90% CI = 0.052, 127.225; mSA = 34.353; 90% CI = 0.085, 110.456; Fig. 6b). The posterior distributions of the migration parameters were very similar in each direction. A more detailed model adding population growth parameters for each species also was implemented, but this was abandoned because the data were not able to converge on an informative posterior distribution for the growth parameters (results not shown).

Fig. 6

Posterior probability densities for demographic parameters from LAM ARC. Results for each parameter are shown from 3 independent runs, a) Posterior distribution for θ. Estimates for Myzopoda aurita are shown in black, and those for M. schliemanni are shown in gray, b) Posterior distributions for migration rates. Estimates for mAS are shown in black, and those for mSA are shown in gray.

Population genetic analyses.—Mantel tests and analyses of molecular variance (AMOVA—Excoffier et al. 1992) were used to explore geographic structuring between and within species of Myzopoda. A Mantel test revealed no significant correlation between geographic and genetic distances (P = 0.539). Significant structure was revealed in a 2-level AMOVA treating each species as a single population (ϕST = 0.404, P = 0.008). Because M. schliemanni was sampled from only 2 locations, 1 of which was represented by a single individual, we were not able to investigate patterns of genetic structure within that species. However, M. aurita was sampled from Analalava and Tolagnaro in the southeast, from Tampolo and Foulpointe in the northeast, and from Kianjavato in central-eastern Madagascar (Fig. 1). An AMOVA among these 3 geographic groups revealed no structure among the groups (ϕCT = − 0.096, P = 0.591), but significant structure among populations within the groups (ϕsc = 0.979, P < 0.0001) and within populations (ϕST = 0.976, P = 0.002).


Genetic support for M. schliemanni.—Phylogenetic and network analyses largely support the existence of a differentiated gene pool of Myzopoda in central-western Madagascar. Average intraspecific diversity measures for Cytb are 1.25% for M. schliemanni and 1.43% forM. aurita, whereas interspecific measures average 2.38%. This interspecific measure approaches that observed between other chiropteran sister species at this locus (X̄ = 6.83%, range = 2.50–16.42%—Bradley and Baker 2001).

Coalescent analyses place the time of divergence of these 2 species of Myzopoda at approximately 73,509 years ago (90% HPD = 14,912, 888,947 year ago; Table 4). This estimate is consistent among runs using different random seeds (Fig. 5b), and is generally coincident with the late Sangamon–Eowisconsin interglacial period of the Pleistocene (79,000 years ago). Pleistocene glacial periods are associated with colder and drier conditions in Madagascar, particularly in the high interior of the island (Burney et al. 2004). During such epochs, rain-forest formations associated with warm and wet microclimates would have been restricted to lowlands of the east and northwest (Burney 2003). During the interglacial periods, lower-lying humid forests would have extended upslope, which would have included the elevational expansion of the range of Ravenala, and potentially provided a corridor for the dispersal of Myzopoda across the Central Highlands or considerably facilitated movements across the zone between the Central and Northern Highlands.

Both phylogenetic and network analyses further support the division of each species into 2 major clades; these are referred to as clades A and D within M. schliemanni and clades B and C within M. aurita. Although phylogenetic analyses were not able to resolve basal relationships among these clades with any statistical rigor, a statistical parsimony network suggests a relationship discordant with the taxonomy in which clade D is more similar to M. aurita than to other M. schliemanni. Genetic distance analyses support this pattern, finding that clades A and D are the most differentiated at both Cytb (3.21%) and the D-loop (8.30%). Parsimony bootstrap analyses exploring hypothetical relationships among these clades were not able to reject any of the tested hypotheses of reciprocal monophyly between the species, monophyly of M. schliemanni alone, or monophyly of a combined CD clade.

We present 2 possible explanations for this divergent clade D, a single haplotype present in 2 individuals that are morphologically and distributionally consistent with M. schliemanni. If this haplotype was maintained from the common ancestor of M. schliemanni and M. aurita, its presence in the modern population of M. schliemanni would be a stochastic effect of incomplete lineage sorting. Alternatively, the presence of the D haplotype in M. schliemanni may be due to recent gene flow from M. aurita. These hypotheses might be tested using rapidly evolving nuclear markers such as microsatellites. If the mitochondrial pattern we observe here is due to incomplete lineage sorting, then we do not expect this stochastic force to extend over multiple loci. If the mitochondrial pattern is due to recent introgression from M. aurita, then we expect to see a similar pattern at nuclear loci. Analyses indicating a possible large ancestral population size (Fig. 5a) might suggest tentative support for the hypothesis of incomplete lineage sorting; however, a notably low genetic distance between clades C of M. aurita and D of M. schliemanni might indicate some support for a hypothesis of recent gene flow.

Center of origin for Myzopoda.Goodman et al. (2007a) presented hypotheses for the divergence of M. schliemanni from M. aurita, which proceeded from the assumption that the eastern species is the older of the 2, and that the origin of M. schliemanni involved some westward dispersal across the island. We used coalescent analyses under the isolation with migration model to test the assumption of Goodman et al. (2007a) regarding the pattern of speciation in Myzopoda. This model includes a population splitting parameter (s), which is defined here as the proportion of the MRCA that evolved into M. aurita. If the MRCA was distributed in eastern Madagascar, then the species divergence event would have involved a small proportion of the population dispersing west to found M. schliemanni, and the estimate of s should approach 1.0. Alternatively, if the MRCA was distributed in western Madagascar, then the estimate of s should approach 0.0. Our IM analyses result in an estimate of s of 0.0005 (Fig. 5c), which suggests an ancestral distribution for Myzopoda in western Madagascar and is in conflict with the consideration of Goodman et al. (2007a) of eastern Madagascar as the ancestral range of Myzopoda in Madagascar.

On the basis of several different types of data, Wells (2003) presented hypotheses of the stages of environmental change on Madagascar starting about 80 mya and associated evolution of the modern forest biomes. He concluded that the western dry deciduous forests, falling within the range of M. schliemanni, evolved after the mid-Paleocene (approximately 58 mya), whereas the more recent eastern humid forest, home to M. aurita, originated at the Eocene-Oligocene boundary (approximately 37 mya). Using a dispersal-vicariance analysis combined with a local molecular clock technique based on 3 different gene regions, Kress and Specht (2006) estimated that differentiation of the genus Ravenala occurred 55–49 mya. Hence, on the basis of these various analyses the divergence of the western dry forest biome took place before that of the eastern humid forests and the date of the evolution of Ravenala is remarkably concordant with Wells (2003) mid-Paleocene hypothesis for the origin of the western forests. Given the close association of Myzopoda with the distribution of Ravenala, which provides an important day-roost tree for this bat, the association of this bat's remarkable suckerlike structures may be closely tied to the evolution of this tree's large fleshy leaves (see below).

Effective population size of and migration between species of Myzopoda.—Proving relatively robust to differences in the demographic models used by IM and LAMARC, our estimates of Ne for both species of Myzopoda are roughly on the order of 105, with estimates for M. aurita being slightly greater than those for M. schliemanni. Under the isolation with migration model, we estimate the effective population size of M. aurita at 132,742 (90% HPD = 10,379, 818,853), whereas that of M. schliemanni is 134,630 (90% HPD = 5,961, 350,860). Under the LAMARC model, we estimate the effective population size to be 100,054 (90% HPD = 55,792, 197,508), whereas that of M. schliemanni is 61,453 (90% HPD = 36,113, 97,328). Despite the similarities in estimates of Ne for these species, their distributions are markedly different, with M. aurita having a much broader range than M. schliemanni (Fig. 1). If these species live in the same densities per unit area, these results imply that the range of M. schliemanni may be bigger than is currently documented. Alternatively, 1 of the species may have undergone a change in range, population size, or both (increase in M. aurita or decrease in M. schliemanni) that is not yet observable in the genetics of the species. Estimates of the effective population size of the MRCA are extremely large (Ne = 1,590,182; 90% HPD = 859,277, 3,271,590); however, these mitochondrial data do not appear to be informative regarding this parameter and the posterior distribution for 0MRCA is nearly flat (Fig. 5a). We place no confidence in this estimate and emphasize that we neither use nor advise the use of this estimate in drawing conclusions regarding the actual size of the ancestral population or patterns of change in population size.

The 2 coalescent analyses yield conflicting results regarding migration. Under the isolation with migration model, we esti mate migration to be negligible in both directions (Table 4; Fig. 5d). The LAMARC model yields substantial and nearly equal estimates of migration in both directions (mAS = 37.015, 90% CI = 0.052, 127.225; mSA = 34.353, 90% CI = 0.085, 110.456; Table 4; Fig. 6b). We consider that these conflicting results are due to differences in the models utilized by the 2 methods. Specifically, IM was implemented under a model that permits changes in the population sizes of the 2 daughter populations, whereas LAMARC assumes an equilibrium state between drift and gene flow. Thus, IM could accommodate some degree of incomplete lineage sorting, whereas LAMARC would attribute similarities between the 2 populations such as the presence of clade D in M. schliemanni to gene flow. It remains to be determined which of these 2 possibilities is more biologically accurate for Myzopoda.

Population structure within M. aurita. — Patterns of population structuring in M. aurita were explored using statistical parsimony network analyses (Fig. 3) and AMOVA. Although the network analyses reveal 2 significantly differentiated gene pools within this species, there appears to be no geographic association. Clade B contains haplotypes from the southeast and northeast, and clade C contains haplotypes from the southeastern, central eastern, and northeastern parts of the species' range. In both the southeast (Tolagnaro and Analalava) and northeast (Tampolo and Foulpointe), sampling localities that are very close geographically are quite distinct genetically (at Cytb, 2.17% sequence difference in the southeast and 1.93% in the northeast). This level of genetic diversity does not seem to extend to the local level, because multiple clades have not been detected within a sampling locality of M. aurita.

The patterns observed in the network analyses also are reflected in AMOVAs. We found significant structure within regional groups (northeast, central east, and southeast) and within populations in M. aurita. More variation was observed within groups than was expected, reflecting the discordance between geographic and genetic distances discussed previously. At the within-population level, we observed less variation than was expected. We consider this a spurious result of low sample size at this level of the analysis.

These results may indicate substantial movement of M. aurita about the landscape. Randrianandrianina et al. (2006) discuss the possibility of seasonal movements in this species. They note that acoustic surveys of Réserve Spéciale d'Analamazaotra in July–August and an inventory in September revealed the presence of M. aurita, whereas similar surveys in February detected no trace of the bat. Echolocation calls of this species are highly distinctive (Göpfert and Wasserthal 1995), and it is unlikely that they would be overlooked or mistaken in an acoustic survey.

Several studies have noted an association between Myzopoda and R. madagascariensis (Goodman et al. 2007a; Göpfert and Wasserthal 1995; Schliemann and Maas 1978). Similarly to the New World genus Thyroptera with the genera Heliconia (family Musaceae) and Calathea (family Marantaceae—Findley and Wilson 1974; Vonhof and Fenton 2004), Myzopoda uses the suckerlike appendages on its thumbs and soles to adhere to the broad, smooth leaves of Ravenala, which it uses as a day-roost. In much of eastern Madagascar, dense stands of Ravenala are present as a dominant pioneering species where vast areas of primary lowland forest have been cleared (Goodman et al. 2007a). These zones may act to increase the population size and facilitate the movement of Myzopoda across the landscape by providing day-roost sites. If so, Myzopoda may be among the few Malagasy vertebrate species to see any benefit from the rapid deforestation affecting the island.


We are grateful to the Direction des Eaux et Forêts and Association National pour la Gestion des Aires Protégées for issuing permits to conduct faunal surveys on Madagascar and Prof. O. Ramilijaona and Dr. D. Rakotondravony for assistance with numerous administrative details. Field research associated with this paper has been generously supported by World Wildlife Fund—Madagascar, Ellen Thorne Smith Fund of the Field Museum of Natural History, John D. and Catherine T. MacArthur Foundation, and the Volkswagen Foundation, and was assisted by F. Ratrimomanarivo and T. Tianarifidy. ALR was supported by a Gaylord Donnelley Postdoctoral Fellowship from Yale University for a portion of this study. Tissue samples were kindly made available by the National Museum of Natural History, Washington, D.C., and Project Ramanavy, and for access to this material, we are grateful to M. Carlton and R. Jenkins, respectively. We thank M. Cox and M. Vonhof for comments on the manuscript and L. Wilmé for Fig. 1 and the French résumé.

Literature Cited

View Abstract