Cryptic diversity in Hipposideros commersoni sensu stricto (Chiroptera: Hipposideridae) in the western portion of Madagascar

The Commerson’s leaf-nosed bat, Hipposideros commersoni sensu stricto, is endemic to Madagascar and is relatively common in the western portion of the island, where it is found in areas, including forested zones, from sea level to 1325 m. A previous study on morphological patterns of geographic variation within the species highlighted the presence of two distinct morphotypes; larger individuals in the north portion of the island and smaller individuals in the south. The main aim of this study was to use a combination of craniodental morphology and molecular data (mitochondrial and nuclear) to test previous hypotheses based on morphology and clarify the evolutionary history of the species group. We sequenced mitochondrial and nuclear genes from Hipposideros commersoni obtained from the western portion of Madagascar, and compared them with other African species as outgroups. We analyzed the sequence data using Maximum Likelihood and Bayesian phylogenetic inference. Divergence dates were estimated using Bayesian molecular clock approach. Variation in craniodental variables was also assessed from sequenced individuals. The molecular analyses suggest that H. commersoni is not monophyletic, with strong support for the presence of several independently evolving lineages. Two individuals amongst those sequenced from Isalo (south central) and Itampolo (southwest) form a separate clade (Clade A), distinct from other H. commersoni, and sister to continental African H. vittatus and H. gigas. Within the H. commersoni clade, the molecular data support two geographically distributed clades; one from the south (Clade B) and the other from the north (Clade C), which diverged approximately 3.38 million years ago. Morphometric data were consistent with the molecular analyses, suggesting a north–south break within H. commersoni. However, at some localities, animals from both clades occurred in sympatry and these individuals could not be differentiated based on external and craniodental measurements. Using a combination of molecular and morphological characters, this study presents evidence of cryptic diversity in H. commersoni on Madagascar. Further fine-scale phylogeographic studies are needed to fully resolve the systematics of H. commersoni. This study highlights the utility of the combined approach in employing both morphological and molecular data to provide insights into the evolutionary history of Malagasy population currently assigned to H. commersoni.


Background
Members of the Family Hipposideridae, known as Old World leaf-nosed bats, are one of the most widespread and abundant groups of insectivorous bats and inhabit tropical and subtropical regions of Africa, the Middle East, Asia and Australia [1]. To a large extent, species within this genus have been defined based on their external and craniodental morphology. In a recent summary, 70 species of Hipposideros were recognized [1], subsequently numerous other taxa have been described (e.g. [2][3][4][5]) and the taxonomy of the group, predominantly at the species level, is far from resolved. As a tool to understand the evolutionary history of members of this genus, closely related species are often placed in morphological species groups (e.g., [1,6]). As currently delineated, the H. commersoni group includes the Afro-Malagasy taxa H. commersoni (É. Geoffroy, 1813), described from Madagascar, and H. thomensis (Bocage, 1891), H. gigas (Wagner, 1845) and H. vittatus (Peters, 1852), from continental Africa and offshore islands. The last-named three forms were previously considered subspecies of H. commersoni sensu lato, but were recently raised to species rank [1] based on reputed morphology and echolocation call differences [7][8][9]). As members of the H. commersoni group s.l. have to date not been the subject of a detailed phylogenetic study, it is unclear if these taxonomic changes reflect the evolutionary relationships within this portion of the genus or are examples of morphological convergence [10].
Hipposideros commersoni sensu stricto is a widespread endemic to Madagascar and can be found from sea level to 1325 m, generally in forested zones [11]. Its diet (mostly Coleoptera) and activity in western Madagascar may change seasonally and may be related to possible intra-island movements [12][13][14][15]. On the basis of current information, Malagasy populations of H. commersoni demonstrate considerable geographic variation in morphological measurements and certain patterns cannot be explained by simple clines [16].
In this study, we examine genetic and morphological variation in H. commersoni using samples obtained from different areas of Madagascar within and outside dry forest formations, specifically the western half of the island, to explore aspects of their phylogenetic history and to help resolve the systematic relationships of the different morphotypes recovered by Ranivo and Goodman [16].

Morphological and molecular sampling
In this paper, reference to Hipposideros commersoni is restricted to Malagasy populations and, hence, sensu stricto. A total of 22 H. commersoni (20 females and two males) were included in the molecular portion of this study (Table 1). During the past two decades intensive bats survey were carried out by different researchers across Madagascar, but with a distinct bias towards the west, where there are often extensive cave systems used as day roost sites for Hipposideros. The collection of H. commersoni specimens and associated tissues were greatly biased in this context. This species is present in eastern part of Madagascar but only a few specimens were available. Morphological analyses were only conducted on the 20 females. These samples come from collections made over the past 15 years from 11 localities across the western portion of Madagascar (Fig. 1). All voucher specimens are cataloged in the Field Museum of Natural History (FMNH) or Université d' Antananarivo, Département de Biologie Animale (UADBA). Samples used in the molecular study included the aforementioned material, as well as additional tissue samples of H. vittatus (n = 7) and H. gigas (n = 1), two morphologically similar species [17], from the FMNH and the American Museum of Natural History (AMNH) collections (Table 1, Fig. 1).

DNA sequencing
Genomic DNA was isolated using the NucleoSpin® Tissue kit (Macherey-Nagel, Germany), following the manufacturers protocol for tissue samples. Two mitochondrial (mtDNA) and two nuclear intron (ncDNA) markers were amplified.
The cytochrome b gene (Cyt b) was amplified using two sets of nested primers. The primers L14724AG (5′-ATG ATA TGA AAA ACC ATC GTT G-3′; [4]) and H15915 (5′-TCT CCA TTT CTG GTT TAC AAG AC-3′; [18]) were used to amplify a 1200 bp segment of Cyt b. In specimens in which L14724AG and H15915 did not amplify, the primers JorF (5′-GAC CTT CCA ACT CCC TCA AGC AT-3′; designed for study) and H15553 (5′-TAG GCA AAT AGG AAA TAT CAT TCT GGT-3′; [18]) were used to amplify a smaller 700 bp segment. The hypervariable portion of the control region (CR) of the mitochondrial genome was amplified in all specimens as a single fragment using primers P (5′-TCC TAC CAT CAG CAC CCA AAG C-3′) and E (5′-CCT GAA GTA GGA ACC AGA TG-3′; [19]). The 16 th intron of the signal transducer and activator of transcription 5A (STAT) was amplified using previously published primers (bSTATa 5′-GAA GAA ACA TCA CAA GCC CC-3′, bSTATb 5′-AGA CCT CAT CCT TGG GCC-3′; [20,21]). The 5 th intron of the organic solute transporter subunit alpha gene (OSTA5) was amplified using the primers OSTA5F (5′-TGM WGG YCA TGG TGG AAG GCT TTG-3)′ and OSTA5R (5′-AGA TGC CRT CRG GGA YGA GRA ACA-3′; [22]). The STAT marker was used based on the work of Eick et al. [20], who found high levels of intraspecific divergence for this marker in 58 bat species. Igea et al. [22] identified the intron OSTA5 to be an adequate marker for analyses of species delimitation, gene flow and genetic differentiation within two bat species. Cycle sequencing was performed using BigDye Chemistry (Version 3.1, Applied Biosystems, USA), and products analyzed on a 3100 ABI automated sequencer (Applied Biosystems). All heterozygous sites in the ncDNA were coded using the IUB code. All sequences were first aligned using ClustalW [23] as implemented in BioEdit [24], and thereafter manually to optimize homology. All new sequences were deposited in GenBank (Table 1).

Sequence analyses
The four markers (CR, Cyt b, STAT and OSTA5) were analyzed separately and then combined into a single data set. Gaps were treated as missing data. In addition, the markers were concatenated and analyzed according to origin of marker (mtDNA or nucDNA). The number of variable sites, number of parsimony informative sites and nucleotide frequencies were estimated for each data matrix in MEGA 6 [25].
Phylogenetic reconstruction was performed using both maximum likelihood (ML) and Bayesian (Bayes) approaches using the programs Garli 2.0 [26] and MrBayes 3.2 [27], respectively. The most appropriate substitution model for each gene (CR -HKY + I + G, Cyt b-HKY + I, STAT -TIM1 + I, OSTA5 -TIM1ef + I) was selected using Table 1 Details of specimens included in the molecular analysis (n = 22, 20 females and 2 males). The Hipposideros commersoni specimens are all from Madagascar; more precise details on collection localities are presented in Appendix 1 the Akaike information criterion (AIC) as implemented in jModelTest [28,29]. For the concatenated data sets, partitioned analyses were conducted, with data partitioned by gene, with the parameters of nucleotide substitution models unlinked across partitions. Each ML analysis was initiated from a random starting tree, with nodal support assessed using 1000 bootstrap replicates. Two independent Bayes runs of 5 million generations each were performed; each run consisted of four Monte Carlo Markov chains (MCMC), with topologies sampled every 250 generations. The program Tracer 1.6 [30] was used to determine that the effective sample size (ESS) had reached > 200 for all parameters. A 50 % majority rule consensus tree was constructed using the CONSENSE program in the PHYLIP package [31]. In each simulation the first 20 % of generations were discarded as burn-in, after a pilot run to determine that this was sufficient to achieve stationarity.

Molecular dating
No Rhinolophidae or Hipposideridae fossils are known from before the middle Eocene, but fossils referable to both families are reported from the middle to late Eocene of Europe [32,33], including H. schlosseri from the late Eocene of France [34]. As fossil calibration points are not available for H. commersoni s.l., we expanded the taxonomic sampling used in the molecular clock analysis to allow the use of fossil calibration points. Cyt b sequences were downloaded from GenBank for six species of Hipposideros: H. armiger (DQ865345), H. pratti (EF544427), H.
Divergence dates between clades were estimated from the expanded Cyt b data set using an uncorrelated relaxed lognormal Bayesian molecular clock approach [39], as implemented in BEAST 2.1.3.0 [40]. The HKY + I substitution model was used, with the Yule speciation model as tree prior. As an alternative to fossil calibrated estimate of divergence times, an additional molecular clock analysis was conducted using a fixed mean substitution rate of 1.30 × 10 −8 subs/site/year [5,41]. This analysis was performed using the strict molecular clock model in BEAST. All other parameters were the same as in previous analysis. The MCMC chains were run for 30 million generations, with topologies and parameters logged every 1500 generations. Results were evaluated using Tracer v1.6 [30]. The Effective Sample Size (ESS) values were > 200 for all parameters, suggesting the MCMC chains had sufficiently converged [40]. After discarding the first 25 % of generations as burn-in, the maximum clade credibility tree was constructed using TreeAnnotator 1.7.4 (available in the BEAST package), and then visualized with FigTree 1.3.1 [42].

Morphological measurements
The following standard external measurements were taken from specimens collected in the field before their preparation using a millimeter ruler accurate to the nearest 0.5 mm: total length (TL), tail length (TAIL), hind foot length (HF) (not including claw), ear length (EAR) and forearm length (FA). Further, body mass (WT) was recorded in grams using a spring balance accurate to the nearest 0.5 g.
Cranial and dental measurements were obtained from cleaned skulls of voucher specimens using digital callipers accurate to the nearest 0.1 mm and following for the most part Freeman [43]: cranialgreatest skull length (SL), condyle-basal length (CBL), greatest zygomatic breadth (ZYGO), minimum interorbital width (IOW), greatest mastoid breadth (MAST), rostrum length (ROST), palatal length (PAL); and dentaltotal tooth row (C1-M3), upper molar row (UP MOL R), width at upper canines (C1-C1), width at upper posterior molars (M3-M3), height upper canine (UP CANIN), dentary length (DENT LEN), moment arm of temporal (MOM1 COR), total lower tooth row (I1-M3) and lower tooth row (LOWER TR). Only adult specimens were used in this study, as defined by the eruption of all permanent teeth (often showing some wear), the complete ossification of the basiosphenoid suture, and the development of the sagittal crest. All external and craniodental measurements used in the analyses were made by a single individual (SMG). The number of adult male H. commersoni available in the morphometric dataset was limited, and given there is evidence of sexual dimorphism in this species [16,44], males were excluded from the morphometrics analyses. Intact skulls from 20 adult females were included in this study from 13 localities spanning the latitudinal distribution of H. commersoni in western Madagascar.

Statistical analyses
Shapiro-Wilk's test and Levene's test were implemented to assess the assumptions of normality and equality of variances of variable characters in the dataset. Analysis of Variance (ANOVA) was carried out using post-hoc Tukey tests, to assess morphological and craniodental differences between the derived genetic clades.
Principal component analysis (PCA) was conducted separately on external and craniodental measurements to examine possible segregation of the different molecular clades, as well as geographic variation in H. commersoni. Further hierarchical cluster analysis was implemented using Ward's method on both measurement data sets to provide additional confirmation of the factor loadings obtained and to identify natural groupings among samples (Tables 2 and 3) [45,46]. Data were log-transformed to improve normality and homoscadisticity. All statistical analyses were carried out using SPSS (version 21.0, IBM SPSS Statistics).

DNA sequencing
The four genetic markers were successfully amplified for all 31 taxa included in the molecular portion of the study ( Table 1). The aligned sequence data for each marker included (Table 1): CR, 481 bp (114 variable sites); Cyt b, 705 bp (60 variable sites); STAT, 476 bp (six variable sites); and OSTA5, 676 bp (nine variable sites). The nucleotide composition and the levels of variation of the two marker systems (mtDNA vs nucDNA) differed ( Table 4). The mtDNA partition contained the highest number of variable characters (174 variable sites), while the ncDNA data was more conserved (15 variable sites). For the STAT gene, only eight unique haplotypes were identified. The haplotypic diversity for this dataset is high (h = 0.80), but the nucleotide diversity is low (π = 0.00274). For OSTA5 gene, 10 unique haplotypes were identified. Once again low levels of nucleotide variability were observed (π = 0.00264). As expected, CR contained the highest proportion of variable characters (24 % variable characters) followed by Cyt b (9 % variable characters).

Phylogenetic analysis
The phylogenetic analysis of each nuclear marker independently resulted in largely unresolved trees, which is not surprising given the few number of variable characters observed (Additional file 1: Figure S1 and S2; Table 4). The  mtDNA marker topologies were better resolved (Additional file 1: Figure S3 and S4). There was no significant (ML bootstrap > 70 %; Bayesian posterior probability >95 %) conflict among the topologies recovered by the independent analysis of the four molecular markers [47] and the molecular data was concatenated (2336 bp, 118 variable characters). The ML and Bayesian analysis of the concatenated data matrix (mtDNA + ncDNA; Fig. 1 Uncorrected pairwise sequence distances for the two mtDNA regions (CR and Cyt b) are presented in Table 5. Genetically, Clade A is as distant from the H. gigas-H. vittatus species pair (respectively 3.2 % and 2.9 %) as it is from other H. commersoni placed in Clades B and C, again highlighting the uniqueness of this lineage.

Molecular clock dating
The maximum clade probability tree (Fig. 2) Fig. 2). The estimated divergence times using the substitution rate calibrated molecular clock resulted in more recent divergence dates (Additional file 2: Figure S5). For example, molecular clock estimates suggest that Clade A diverged from its sister taxa 4.20 MYA (95 % HPD 1.99-13.73) and the two sister Clades B & C shared their last common ancestor 2.55 MYA (95 % HPD 1.15-7.89, Table 6). The 95 % HPD intervals for divergence events from both analyses (fossil calibrated and substitution rate calibrated) were broad and but did not show considerable overlap. From Taylor et al. [37], Patterns of sequence variability are presented for two mtDNA regions (CR and Cyt b), two nuclear introns (bSTAT and OSTA5) and the combined data matrix. The total number of nucleotide sites, variable and parsimony informative sites, as well nucleotide frequencies are given for each partition and the combined data matrix R. mossambicus and R. fumigatus, diverged 6.96 MYA, which is older than our estimate and R. ruwenzorii and R. maclaudi about 2.99 MYA, which is more recent than our estimated. We suggest that using a calibration point allowed BEAST to estimate a more realistic clock rate. The substitution rate of 1.0 in the fossil calibrated clock allowed the determination of relative rates for each recovered clade.

Morphometrics
As Clade A was only comprised of two individuals, statistical comparisons were only made between animals belonging to Clades B and C. Shapiro-Wilk's test (P > 0.05) and a visual inspection of their histograms showed that the variable characters in the dataset were normally distributed. Levene's test verified the equality of variances in the samples (P > 0.05). Analysis of variance (ANOVA) revealed significant variation in four of six external variables and 15 of 16 craniodental variables when taxa were sorted into Clades A, B or C following the results of molecular analyses (Table 7). Clade A is morphologically similar to Clade B, but is morphologically differentiated from Clade C ( Table 7). The larger Clade C bats had significantly greater total length, tail length, ear length and forearm length than specimens assigned to Clade B, yet there were no significant differences in hindfoot length and body mass between the clades (Table 7). In parallel, craniodental measurements were significantly larger in Clade C than Clade B, except minimum interorbital width ( Table 7).
The first two unrotated principal components (PC1 and PC2) explained 67.6 % of the total variance in external measurements (Fig. 3a) and 88.2 % of total variance in cranialdental morphology (Fig. 3b). PCA plots of external and      (Fig. 3). In the case of external measurements, PC2 showed an inverse relationship between tail length and hind foot lengthbats that had high loadings on PC2 had a relatively long tail but short hind foot, whereas bats that had low loadings on PC2 had a relatively short tail but long hind foot (Table 8). In the case of craniodental variables, PC2 indicated skull robustness ( Table 9) -bats that loaded high on PC2 had a relatively larger interorbital width than bats that loaded low on PC2. Two major clusters were recovered from the dendrograms produced by the hierarchical cluster analyses of external and craniodental variables, supporting the PCA analysis. The first cluster, recovered in both dendrograms, included all eight individuals from northern Madagascar assigned to Clade C and two animals genetically assigned to Clade B (FMNH 221308 from Ankarana and FMNH 175777 from Namoroka). The second cluster contained the smaller southern individuals from Clades A and B (Fig. 4), which confirms that the two genetically divergent animals belonging to Clade A are morphologically similar to animals from Clade B.

Discussion
This study combines evidence from molecular (mtDNA and ncDNA) and morphological characters to provide support for the reciprocal monophyly of several independently evolving lineages within Hipposideros commersoni, occurring on Africa and nearby islands, as well on Madagascar. One of the principal questions addressed is the evolutionary history and systematic relationships of Malagasy populations currently assigned to H. commersoni, as well as African populations currently placed within the H. commersoni species group [1].
The results suggest that previous taxonomic treatments of the group underestimated species diversity of H. commersoni and that a cryptic species appears to be present. Although our geographic sampling did not cover the complete range of this species, specifically the eastern portion of the island, the results indicate nonmonophyly with respect to Madagascar of different recovered Malagasy clades.
Single locus or multilocus molecular data and/or morphological differences have been used previously to identify See materials and methods for definitions of variable acronyms. ns = not significant cryptic species diversity in southeast Asian [10,48,49] and African [50,51] hipposiderids. Based on an analysis of the H. larvatus species complex using mitochondrial control region markers, morphology and bioacoustics two forms were identified (H. khasiana and H. grandis) that were differentiated based on haplotypic structure and phonics [52].
In another example, H. khaokhouaensis from Laos is similar in general body size and shape to its sister species H. rotalis, but differs in aspects of the noseleaf, skull structure related to bioacoustics and echolocation frequency [4]. The phylogenetic relationships within the African H. ruber species complex were investigated [51] using Cyt b to determine the taxonomic status of two divergent genetic forms often found in sympatry in Senegal, which might represent cryptic species despite being morphologically indistinguishable. However, in this latter case, absence of nuclear gene flow between these two reputed forms remains to be investigated to demonstrate their reproductive isolation. All H. commersoni sequenced in the current study were from the western half or extreme north of Madagascar. The molecular analyses presented herein indicate that H. commersoni as currently diagnosed is not monophyletic with respect to Madagascar, and with strong support for the presence of divergent lineages. As two individuals from Isalo (FMNH 175970) and Itampolo (FMNH 184173) form a well-supported monophyletic group (Clade A), basal to African H. vittatus and H. gigas, and separate from the balance of Malagasy H. commersoni (Clades B and C), a single origin of this species complex on the island is not supported. The long branches separating these clades (ranging from 2.6 to 3.2 % uncorrected sequence divergence) indicate relatively deep independent evolutionary trajectories of several million years based on molecular clock inferences. Although Clade A was only significantly supported by the mitochondrial data, this is not surprising, given the relatively conservative nature of the nuclear markers   sequenced. The absence of haplotype sharing in OSTA5 gene for Clades B and C, however, does indicate at least some degree of genetic isolation between these two groups.
These results indicate that members of the H. commersoni species group do not represent a single widespread Afro-Malagasy taxon. Furthermore, the molecular data supports a certain level of divergence between Clades B (in the north) and C (in the south). Additional samples are needed A B Fig. 4 Hierarchical clustering dendrogram for a log-transformed external variables for 19 specimens and b log-transformed craniodental variables of 17 female Hipposideros commersoni. Collection locality information is assigned to each individual and color-coding is based on main lineages recovered by molecular data to ascertain whether the recovered phylogeny is an artefact of sampling or the result of isolation by distance. The sequence divergence of Clade A with respect to the balance of H. commersoni (Clades B and C) is comparable with that observed between African H. vittatus and H. gigas, and based on the molecular clock analysis, it is estimated that the Clade A lineage diverged from its sister taxa during the Miocene. Genetic divergence in mitochondrial genes varies widely among species. Avise [53] highlighted that due to the matrilineal nature of inheritance of mitochondrial genes, relatively deep divergences do not necessarily correspond to species boundaries. Further, significant nuclear gene flow may occur among divergent mitochondrial phylogroups. Using the published literature, Baker and Bradley [54] found an interval of 3.3-14.7 % uncorrected genetic distance between sister species of bats, and distances ranging from 0.6-2.3 % encompassing intraspecific variation. These values have been corroborated by recent studies of cryptic species of Asian Hipposideros, which show three different levels of interspecific divergences: 1) as low as 3.9 %, with supporting evidence from external and craniodental morphology, as well as bioacoustics [4]; 2) an intermediate level of 6.5 %, with corroborating evidence from bioacoustics [55]; and 3) as high as 13.4 %, with corroborating evidence from bioacoustics [52]. The sequence divergence values recovered in the current study separating Clade A from other H. commersoni (Clades B and C) suggest that previous taxonomic conclusions underestimated the species diversity of Malagasy bats currently classified as H. commersoni.
Within the portion of the phylogeny composed of most individuals assigned to H. commersoni, the molecular data support two largely geographically non-overlapping clades: a northern group (Clade C) with a relatively limited range and a southern group with a broader geographical distribution (Clade B). The molecular clock analyses indicate that these two clades diverged from one another approximately 3.38 MYA. Morphometric analyses are generally consistent with the molecular data, suggesting a north-south break between animals assigned to Clades B and C. The exception was in Ankarana (far north), where the two lineages cooccur but individuals from each clade could not be differentiated based on multivariate analyses of external and craniodental measurements (Fig. 4). Patterns of morphological variation were not uniform or falling along well-defined clines, such as latitude, and members of these two clades do not completely separate from one another. In term of genetics, based on currently available samples, Clades B and Clade C are differentiated, for example, the uncorrected Cyt b sequence divergence is 1.9 %. In the case of the two individuals falling within Clade A, they are genetically distinct from those in Clade B, but show no apparent morphological differentiation. Hence, we interpret this variation as some form of incipient speciation between animals assigned to Clades B and C.
Ramasindrazana et al. [44] have recently analyzed echolocation calls of animals referred to as H. commersoni captured in western Madagascar. They found latitudinal variation -animals from the north being larger and emitting lower call frequencies and those from the south smaller and emitting higher call frequencies. On average, females, referred to H. commersoni, from the north (Ankarana) deviate from the allometric relationship with lower resting frequency of echolocation calls than predicted from body size. These authors suggested that this pattern might be explained by either regional variation in bioacoustics, intra-island migratory movements or the presence of a cryptic species. The animals that deviated from the pattern were not sequenced in this current study and no further interpretation can be offered.

Evolution of Malagasy Hipposideros
A particularly striking result of the current analyses is the existence of a previously unrecognized clade of Malagasy H. commersoni (Clade A), estimated to have diverged from sister taxa (Malagasy H. commersoni Clades B and C, and African H. gigas and H. vittatus) during the late Miocene (5.81 MYA). Hipposideros gigas and H. vittatus are more closely related to H. commersoni Clades B and C, suggestive of two dispersal hypotheses. The first scenario is that Clade A and Clade B-C originated from two independent African mainland-to-Madagascar dispersal events, with Clade A arriving on the island during the Miocene (approximately 5.8 MYA) and Clade B-C more recently (approximately 3.38 MYA). A second hypothesis is that the H. commersoni group evolved on Madagascar and at some point after the end of the Miocene, a population related to Clade A crossed the Mozambique Channel and colonized the African continent leading to two recognized extant forms, H. vittatus and H. gigas, that are morphologically and karyologically similar [17,36]. Following this second hypothesis, speciation took place within the Malagasy population, giving rise to Clades B and C representing the most recent branch of this lineage.
Madagascar was cooler and drier during periods of Pleistocene glaciation, which lead to habitat shifts and forced some taxa to retreat into refugia [56][57][58], in different high mountain areas [59]. Expansion from refugia would have occurred during warmer periods. Hipposideros bones from relatively recent geological deposits are known from several sites on the island [60]. Subfossils from Tsimanampetsotsa, extreme southwest, identified as Hipposideros were slightly smaller than typical H. commersoni [61], which occur in this region today [11]. Samonds [62] conducted research on Hipposideros subfossils from Anjohibe Cave in the northwest of Madagascar, and the excavated fossils were dated between 10,000 and 80,000 years ago. Samonds [62] identified three morphological forms of Hipposideros from these deposits: 1) those fitting with extant H. commersoni; 2), H. besaoka, which was described as a new species, being larger and more robust than H. commersoni; and 3), Hipposideros sp. cf. H. commersoni, which appeared to have some dental differences from modern H. commersoni. Hipposideros besaoka and H. commersoni were sympatric and presumably living in the Anjohibe Cave during the same period, and they show a small amount of overlap in some dental measurements, not related to sexual dimorphism [62]. Our morphological and molecular data suggests parallel results in modern populations of H. commersoni, with Clade B and C known to occur in sympatry at one northern locality. This raises the intriguing possibility that one of the phylogenetic clades identified in this paper (Clade B or Clade C), might be referable to H. besaoka and, in this case, this species is not extinct. Further fine-scale phylogeographic studies using variable nuclear markers such as microsatellites are needed to clarify species boundaries and give a greater understanding of the processes underpinning the evolution of these taxa across Madagascar.

Geographically correlated population structure
Within H. commersoni (Clade B-C), the molecular data support two regionally associated clades: a small-bodied southern group with a broad geographical distribution (Clade B) and a large-bodied northern group (Clade C) with a relatively limited range. The molecular clock analyses indicate that these two clades diverged from one another approximately 3.38 MYA. Morphometric data are consistent with the molecular data, suggesting a north-south break in distribution. These two lineages are not completely allopatric. In Ankarana, sequenced individuals assigned to these two genetic clades could not be distinguished using external and craniodental measurements (Fig. 3). The morphometric data in the present study is consistent with conclusions of a previous study on geographic variation in morphology of this taxon in western Madagascar [16]. Specimens grouped into two distinct morphotypes, a larger morphotype found in northern Madagascar (from Analamerana to Ankarana and south to Bemaraha) and a smaller morphotype widely distributed in the south, from Isalo to Tsimanampetsotsa. Ranivo and Goodman [16] found that male H. commersoni do not show the same pattern and are largely homogenous in size across these zones.
At least three other Malagasy bat species, Paratriaenops furculus [63], Chaerephon leucogaster [64] and Myotis goudoti [65] show similar haplotypic segregation along a latitudinal gradient. However, the latitudinal distribution of different clades and the calculated expansion periods of the other species differ from late Pleistocene in M. goudoti to early Holocene in C. leucogaster, suggesting that no common historical process underlies the different demographic events between these taxa [64,65]. Ranivo and Goodman [16] found both H. commersoni morphotypes in Isalo. The morphologically divergent animals from Isalo included two specimens (FMNH 175973 and 175975) that were collected on the same day and at the same cave site as the Isalo specimen (FMNH 175970) analyzed in our molecular study, which falls into Clade A. This latter specimen morphologically aligns with the smaller southern individuals, while FMNH 175973 and 175975 are of the larger northern morphotype. This may indicate some form of intraisland movements.
In eastern Africa, seasonal fluctuations in abundance of prey utilized by large hipposiderids are pronounced, which can result in food shortages during the cool dry season. These shifts in the resource base have been invoked to explain local seasonal movement in H. vittatus/ H. gigas to areas with greater food abundance [7,66]. It is unclear if H. commersoni remains inactive in caves during times of resource shortage or if local populations migrate to other sites. Large hipposiderid bats have high wing loading and low to medium aspect ratios [67], which may favor relatively quick, long-distance movements, allowing certain populations to track food resources [68,69]. The colonization and speciation history of H. commersoni on Madagascar, as represented by a single species occurring on the island, is certainly more complex than currently understood. Further studies including increased spatial sampling and the use of additional molecular markers particularly faster evolving nuclear markers are needed to fully resolve the evolutionary history and associated systematics of the different clades occurring on Madagascar.

Conclusions
This study provides evidence, particularly from mitochondrial data, for the existence of at least two sympatricallyoccurring species of the genus Hipposideros on Madagascar. Absence of nuclear gene flow between groups remains to be established to verify their reproductive isolation, yet the lack of haplotype sharing in OSTA5 for Clades B and C indicates some degree of genetic isolation between these clades. Subfossil evidence indicates that in the recent geological past two species, H. commersoni and the presumed extinct H. besaoka, occurred in sympatry [62]. Given that we have recovered two genetically distinct lineages of H. commersoni (Clades B and C) living on occasion in sympatry, this might indicate that one of them is H. besaoka and, hence, still extant. A detailed morphological comparison of the type series of Samonds' [62] H. besaoka with modern H. commersoni represented in our data is needed to test this intriguing possibility and crucial before the description of a possible undescribed species. All listed specimens were included in the molecular analyses. and with the exception of those underlined, were also included in the morphometric craniodental analysis. For voucher specimens used in molecular analyses, clade assignments are given based on the supermatrix data presented in Fig. 2. Localities correspond with those presented in Fig. 1