Our genetic dataset is based on 453 hipposiderid individuals, the vast majority being represented by museum vouchers. We generated original genetic data from 319 individuals collected at 102 georeferenced localities, and complemented them with 134 mitochondrial sequences from 90 localities downloaded from GenBank (we obtained new sequence data for five individuals with prior GenBank records; see Suppl. material 1: Figure S1 and Appendix I). All individuals were sequenced for Cytochrome-b (cyt-b) in order to maximize assessment of genetic diversity; however, redundant haplotypes were removed for subsequent phylogenetic analyses (see Appendix I for complete list of individuals sequenced). The bats newly sequenced for this study were obtained over several decades in the course of small mammal surveys across sub-Saharan Africa and Madagascar, with relatively dense sampling in East Africa. Initial assignment of East African individuals to species was determined using meristic, mensural, and qualitative characters published in the bat keys of Thorn et al. (2009) and Patterson and Webala (2012). Collection methods followed mammal guidelines for the use of wild mammals in research and education (Sikes and the Animal Care and Use Committee of the American Society of Mammalogists 2016) and the most recent collections were approved under Field Museum of Natural History’s IACUC #2012-003. Only GenBank records for cyt-b were available for records of the Arabian-North African hipposiderid Asellia, which was included for context in the phylogenetic analyses. Lacking information from nuclear introns, we draw no firm conclusions from their placement and do not discuss Asellia in this paper (see Benda et al. 2011; Bray and Benda 2016).

Appendix I contains the institutions and voucher numbers, GenBank accession numbers, and locality information for our samples. The fact that museum voucher specimens were used wherever possible for the genetic analyses permits the genetic analysis to serve as a foundation for integrative taxonomic analyses of dental, cranial, and skeletal variation, using the same specimens. To avoid adding to current taxonomic confusion, we take a conservative approach in assigning names to clades in our analyses. Where a clade’s taxonomic identity was ambiguous or unknown, we referred to it simply as a numbered clade. Integrative taxonomic diagnoses of the various clades supported by our analyses will be necessary to determine which, if any, existing names may apply to them. However, to relate our results to those of earlier studies of African Hipposideros (Vallo et al. 2008, 2011; Monadjem et al. 2013), we cross-referenced specimens used in two or more analyses to equate the various non-binomial names that have been applied to these cryptic lineages.

Genomic DNA from preserved tissue samples was extracted using the Wizard SV 96 Genomic DNA Purification System (Promega Corporation, WI, USA). Fresh specimens were sequenced for mitochondrial cytochrome-b (cyt-b), using the primer pair LGL 765F and LGL 766R (Bickham et al. 1995; Bickham et al. 2004), and four unlinked autosomal nuclear introns: ACOX2 intron 3, COPS7A intron 4, ROGDI intron 7 (Salicini et al. 2011), and STAT5B (Matthee et al. 2001) for hipposiderid specimens and the sister group Triaenops afer (Rhinonycteridae; see Table 1 for primer information). PCR amplification, thermocycler conditions, and sequencing were identical to Patterson et al. (2018) and Demos et al. (2018). Sequences were assembled and edited using GENEIOUS PRO v.11.1.5 (Biomatters Ltd). Sequence alignments were made using MUSCLE (Edgar 2004) with default settings in GENEIOUS. Protein coding data from cyt-b were translated to amino acids to determine codon positions and confirm the absence of premature stop codons, deletions, and insertions. Several gaps were incorporated in the nuclear intron alignments, but their positions were unambiguous.

Sequence alignments used in this study have been deposited on the FIGSHARE data repository (https://doi.org/10.6084/m9.figshare.11936250). All newly generated sequences were deposited in GenBank with accession numbers MT149315–MT149893 (see also Appendix I).

jMODELTEST2 (Darriba et al. 2012) on CIPRES Science Gateway v. 3.3 (Miller et al. 2010) was used to determine the sequence substitution models that best fit the data using the Bayesian Information Criterion (BIC) for cyt-b and the four nuclear introns. PARTITIONFINDER2 (Lanfear et al. 2016) on CIPRES was used to determine the sequence substitution models for the concatenated alignment of four nuclear introns using the Bayesian Information Criterion (BIC) with the ‘greedy’ search algorithm. Uncorrected sequence divergences (p-distances) between and within species/clades were calculated for cyt-b using MEGA X v. 10.0.5 (Kumar et al. 2018). Maximum-likelihood (ML) analyses were performed using the program IQ-TREE v. 1.6.10 (Chernomor et al. 2016; Nguyen et al. 2015) on the CIPRES portal for separate gene trees (cyt-b, ACOX2, COPS7A, ROGDI, and STAT5B) and a concatenated alignment, partitioned by gene, using the four nuclear introns. As in Hillis and Bull (1993), nodes supported by bootstrap values (BP) ≥ 70% were considered strongly supported. Gene tree analyses under a Bayesian Inference (BI) framework were inferred in MRBAYES v. 3.2.7 (Ronquist et al. 2012) on the CIPRES portal for the same set of genes as the ML analyses. Two independent runs were conducted in MrBayes, and nucleotide substitution models were unlinked across partitions for each nuclear locus in the concatenated alignment. Four Markov chains were run for 1 × 10⁸ generations for individual gene trees, and 2 × 10⁷ generations for the concatenated analysis, using default heating values and sampled every 1000^th generation. A conservative 25% burn-in was applied and stationarity of the MRBAYES results was assessed in Tracer v. 1.7 (Rambaut et al. 2018). Majority-rule consensus trees were constructed for each Bayesian analysis. Following Erixon et al. (2003), nodes supported by posterior probabilities (PP) ≥0.95 were considered strongly supported.

Haplotype networks for cyt-b were inferred using the median-joining network algorithm in PopArt v. 1.7 (Leigh and Bryant 2015). Separate analyses were carried out for the following clades, each consisting of four subclades: (1) Doryhina (D. camerunensis, D. cf. camerunensis, D. cyclops1, and D. cyclops2); (2) Macronycteris (M. commersoni, M. cryptovalorona, M. gigas, and M. vittata); (3) Hipposideros caffer1–4; (4) Hipposideros caffer5–8; and (5) Hipposideros ruber1–4.

Hipposiderid taxa included in the species tree analyses were assigned to either species or numbered clades based on clade support in the ML and BI gene-tree analyses of the cyt-b dataset. This in turn identified populations to be used as ‘candidate species’ in a coalescent-based species-tree approach implemented in StarBEAST2 (Ogilvie et al. 2017), an extension of BEAST v. 2.5.1 (Drummond et al. 2012; Bouckaert et al. 2014). Species tree analysis was conducted using the four nuclear intron alignments. Substitution, clock, and tree models were unlinked across all loci. A lognormal relaxed-clock model was applied to each locus under a Yule tree prior and a linear with constant root population size model. Four independent replicates were run with random starting seeds, and chain lengths of 1 × 10⁸ generations and parameters were sampled every 5,000 steps. For the StarBEAST2 analyses, evidence of convergence and stationarity of posterior distributions of model parameters was assessed based on ESS values >200 and examination of trace files in Tracer v. 1.7. The burn-in was set at 10% and separate runs were assembled using LOGCOMBINER v. 2.5.1 and TREEANNOTATOR v. 2.5.1 (Rambaut et al. 2018).

Based on the well supported clades obtained in the cyt-b gene tree analyses and available intron samples, a lineage delimitation scenario with 18 candidate species was tested. We inferred the evolutionary isolation of their gene pools using the phased nuclear DNA dataset (ACOX2, COPS7A, ROGDI, and STAT5A; 104 individuals) for joint independent lineage delimitation and species-tree estimation evaluated under the multi-species coalescent model using the program BPP v. 3.3 (Yang and Rannala 2014; Rannala and Yang 2017). This analysis was carried out to guide future investigations of the species status of evolutionarily isolated lineages inferred here. Supported lineages will be examined using an integrative species taxonomic approach, including morphological, morphometric, and acoustic characters, as well as ectoparasite associations and distributional data. Species/clade memberships for BPP were identical to individuals assigned to lineages in the species tree analyses. The validity of our assignment of individuals to populations was tested using the guide-tree-free algorithm (A11) in BPP. Because the probability of delimitation by BPP is sensitive to selected parameters (Leaché and Fujita 2010; Yang 2015), we evaluated two independent runs for each of four different combinations of divergence depth and effective population sizes priors (τ and θ, respectively; Table 2). Two independent MCMC chains were run for 5 × 10⁴ generations. The burn-in was 20% and samples drawn every 50^th generation. In total, eight BPP runs were carried out using four phased nuclear intron alignments. Lineages were considered to be statistically well supported when the delimitation posterior probabilities generated were ≥0.95 under all four combinations of priors.

The three Afrotropical hipposiderid genera differ substantially in terms of their internal genetic differentiation. Clades of Hipposideros are separated by cyt-b p-distances averaging 9.7% (2.5–16.1%), whereas Doryrhina clades average p-distances of 4.8% (3.0–5.7%) and Macronycteris clades 2.7% (2.6–2.9%). Distance values for these genera tend to fall at the lower end of values obtained with similar sampling intensity for species-ranked clades in other Afrotropical bat genera: 2.5% for Otomops (Patterson et al. 2018), 9.3% for Miniopterus (Demos et al. 2020), 10% for Scotophilus and Rhinolophus (Demos et al. 2018, 2019a), 13.5% for Myotis (Patterson et al. 2019), and 17% for Nycteris (Demos et al. 2019b). Fewer cyt-b substitutions on average for these hipposiderids does not limit support for individual clades, and because distances do not approach those characteristic of substitutional saturation, the cyt-b tree recovers much of the deeper phylogenetic structure evident with nuclear intron sequences (compare Figs 2, 5).

Both cyt-b and intron analyses securely recovered Doryrhina, Macronycteris, and Hipposideros as monophyletic. Doryrhina + Macronycteris are sister to the remaining hipposiderids. However, only the cyt-b analysis included the hipposiderid genera Aselliscus, Coelops, and Asellia alongside Hipposideros. That analysis recovered all four genera as monophyletic with strong support. Aselliscus and Coelops were recovered as sister to Hipposideros, with Asellia joining later, but these relationships lacked confident support.

Using a supermatrix approach on exemplars of 46 species of hipposiderids, Amador et al. (2018) found Hipposideros sensu stricto to be paraphyletic. They recovered a mostly Asian group of Hipposideros as sister to two subclades, Coelops + Aselliscus and Asellia + African hipposiderids excluding H. jonesi, which was recovered with the Asian taxa. Paraphyly in this molecular analysis echoed earlier indications of Hipposideros paraphyly from morphology (Bogdanowicz and Owen 1998; Hand and Kirsch 1998, 2003). In another supermatrix analysis of exemplars belonging to 49 hipposiderid species, Shi and Rabosky (2015) failed to recover Macronycteris as monophyletic; M. commersoni was sister to all remaining hipposiderids, but strangely it did not group with M. gigas. When the anomalous position of M. commersoni in their tree is ignored, their topology is highly similar to that of Figure 2, except that Asellia (Aselliscus, Coelops) become the sister of Hipposideros (Macronycteris, Doryrhina), rather than sister of just Hipposideros. Using both mitochondrial and nuclear loci, Lavery et al. (2014) found that 17 species of Asian, Oceanian and Australasian Hipposideros were monophyletic with respect to the genera Aselliscus, Coelops, and Anthops. Clearly, missing data and missing taxa compromise all of these phylogenetic appraisals, so that the question of hipposiderid and Hipposideros monophyly remains open. However, subject to its sampling limitations, there is clear support in our analyses of monophyly for Doryrhina, Macronycteris, and Hipposideros as we apply these names.

Despite employing different mitochondrial and nuclear loci and using different sets of taxa, the phylogeny recovered by Lavery et al. (2014) is largely congruent with that in Figure 5. Their earliest diverging species group of Hipposideros is the calcaratus group, not represented in our tree unless H. obscurus is a member (Table 1). Their next diverging unit is the diadema group, which is also positioned near the base of our tree. Their other two groups are paired: the galeritus group (which includes H. cervinus, indicating that this species is misclassified as calcaratus member) joined with the bicolor/ater group. In our intron analysis (Fig. 5), members of the larvatus and diadema groups join H. obscurus as sister to all remaining Hipposideros groups. The remainder form a trichotomy: H. coronatus, typically considered in the bicolor group; H. pygmaeus and H. cervinus, which are listed in different groups but were both considered members of the galeritus unit by Tate (1941); and the erstwhile bicolor group (sensu Hill 1963), which was subdivided into the ater subgroup (for Asian, Oceanian, and Australasian species) and the ruber subgroup (for Atrotropical ones) by Monadjem (2019).

The ater subgroup members included in our mitochondrial analysis (Fig. 2) form a well-supported clade consisting of H. bicolor, H. cineraceus, H. pomona, H. doriae, H. ater, H. khaokhouayensis, H. rotalis, H. halophyllus H. dyacorum, H. ridley, and H. durgadasi. This group is sister to all analyzed members of the ruber subgroup: the various clades allied with Hipposideros beatus, H. caffer, and H. ruber, as well as individuals of the Afrotropical species H. lamottei and H. fuliginosus. H. abae, which was previously considered in the speoris group (Simmons 2005; Murray et al. 2012), is clearly a member of the ruber group. Outside this pairing are the Asian species H. cervinus, H. coronatus, H. coxi, H. obscurus, and H. pygmaeus. Two Afrotropical species also lie outside the ruber + ater clade: H. jonesi and H. marisae, both thought to belong to the speoris group, appear as sisters in Figure 2A.

Parsimony, topological position, and the strong support of branching relationships in the mitochondrial and intron trees (Fig. 5; also Lavery et al. 2014) make it clear that the Afrotropical ruber group represents a comparatively recent colonization event from Asian ancestors–the ruber group is sister to the ater group and this pair has Asian sisters. However, although the basal dichotomy within Hipposideros includes an all Asian clade, lack of support for its sister(s) clouds the phylogenetic position of the H. jonesi-H. marisae clade–possibly sister to all sampled Hipposideros but more likely sandwiched between Asian clades. In any case, Figure 2 suggests that the H. jonesi-H. marisae clade resulted from an earlier African-Asian colonization event.

The lack of agreement in the phylogenetic position of H. diadema and H. larvatus between the concatenated intron tree (Fig. 5) and the species tree (Fig. 6) deserves comment, as both analyses were based on the same genetic dataset. The position of H. diadema-H. larvatus as sister to the ruber group (Fig. 6) runs counter to both our other genetic analyses (Figs 2, 5) and morphological assessments (Hill 1963; Murray et al. 2012; Table 1). This discrepancy is likely due to the generally weaker support for deep nodes within the tree; in the absence of saturation, this is often taken as evidence of rapid evolutionary radiations (e.g., Almeida et al. 2011). Lanier and Knowles (2014) used simulated data on deep phylogenies to show that species-tree methods do account for coalescent variance at deep nodes but that mutational variance among lineages poses the primary challenge for accurate reconstruction. In either case, vastly expanded genetic sampling via NGS techniques offers the most plausible avenues to clearer resolution.

However, the highly distinctive species H. megalotis belongs to its own species group (Table 1) and has not been included in any genetic analysis. Distributed in the Horn of Africa and the Arabian Peninsula, H. megalotis is the only hipposiderid with a fold of skin joining the base of the ear pinnae. Its uniquely specialized auditory system and derived dentition (e.g., loss of anterior premolars and enlargement of outer lower incisors), led Hill (1963) to regard it as a species that diverged early from the other groups of African Hipposideros. Including this species in future analyses would shed light on the group’s biogeography. Were there three colonizations of Africa by Asian groups of Hipposideros or could H. megalotis be sister to all Asian lineages of this genus? This information would greatly clarify ancestral geographic range inference.

The lineage delimitation analyses indicate that a number of hipposiderid lineages are either unnamed or unidentified, and also that a number of recognized species may not be genetically and evolutionarily independent.

Previous studies had indicated that both Hipposideros caffer (Vallo et al. 2008) and H. ruber (Vallo et al. 2011) appear to be complexes of cryptic species. The two are traditionally distinguished on the basis of size and pelage color, H. ruber being the larger and more brightly colored form, but this distinction is clouded by geographic variation in size and the presence of both reddish and gray-brown phases in both species. Our mitochondrial analyses identified four H. ruber lineages and eight H. caffer lineages in two distinct groupings among the sampled populations (Fig. 2). Four of the caffer lineages and three of the ruber clades were identified as putative species by the BPP analyses (Table 5). The large number of clades in East Africa is remarkable: Kenya and Tanzania each support four of the eight clades allied with H. caffer, and all but one of the eight clades known from throughout the continent occur in one or the other East African country. This undoubtedly reflects the region’s great landscape diversity, where West and Central African rainforests reach their eastern limit, southern savannas reach their northern limits, the Sahel reaches its southern limits, and all are riven by the African Rift Valley. It also is a product of our sampling intensity (see Suppl. material 1: Fig. S1).

Because some cyt-b sequences were used in multiple studies of this group, it is possible to relate our clade labels to those used by earlier studies (Table 6). Based on attributions made on morphological grounds by Vallo et al. (2008) and Monadjem et al. (2013), some well-supported but unnamed clades in our analysis can be identified. For instance, caffer1 has a distributional range and includes specimens previously identified as Hipposideros tephrus (Appendix I), while specimens of caffer4 come from near the type locality of H. caffer Sundevall, 1846, and may well represent that species. However, no samples confidently identified as H. ruber from the vicinity of its type locality have been sequenced, leaving the application of that name to clades in any of these trees purely conjectural. Applying formal names only after integrative taxonomic assessment is a responsible course as multispecies coalescent models like BPP can lead to over-splitting of species, especially when applied to geographically variable species complexes with parapatric distributions (Chambers and Hillis 2020).

Doryrhina is a poorly known genus characterized morphologically by the peculiar club-shaped processes on the central and posterior nose leaves. This trait is shared by the two recognized African species, D. cyclops and D. camerunensis, which differ chiefly in size (the latter is larger, with forearm lengths >75 mm). Although D. cyclops is considered to be monotypic, mitochondrial sequences clearly separate West African populations in Liberia and Senegal (cyclops1) from Central African populations in Gabon and Central African Republic (cyclops2), and these are substantially separated from D. camerunensis and a specimen referred to that species from Tanzania (Figs 2, 3). However, both the intron analysis (Fig. 5) and the species tree (Fig. 6) show little or no geographic structure. The BPP analyses confirm that none of the mitochondrial clades is behaving as an independent evolutionary lineage (Table 5). Geographic structure in mtDNA but continent-wide admixture in the nuclear genome could result from either male-biased dispersal with female philopatry or highly structured seasonal migrations, which are known in other hipposiderids. In any case, the genetic patterns of Doryrhina are hard to reconcile with its space-use behavior; individuals appear to have very small home ranges, on the order of a few hectares (Monadjem 2019). An integrative taxonomic review of the genus Doryrhina is needed to determine the validity of D. cyclops and D. camerunensis. It would also shed light on whether six Australo-Papuan species tentatively allocated to that genus (cf. Monadjem 2019) belong there or elsewhere. Tate (1941) had earlier allocated those species to the Australasian muscinus group, convergent on but separate from his Afrotropical cyclops group, but Hill (1963) later united these groups.

Our analysis included four of the five recognized species of Macronycteris, lacking only M. thomensis, which is endemic to São Tomé Island in the Gulf of Guinea. Two species, M. gigas and M. vittata, occur on the African mainland and two others, M. commersoni and M. cryptovalorona, occur on Madagascar. Macronycteris cryptovalorona was named only in 2016, on the basis of its strong genetic divergence from M. commersoni; it appears in Figure 2 as sister to all three remaining species of Macronycteris. Despite a search for diagnostic characters, Goodman et al. (2016) could not distinguish it morphologically from M. commersoni. Both species are known to occur in the same caves in south central and southwestern Madagascar (Goodman et al. 2016; Rakotoarivelo et al. 2019). On the other hand, M. vittata and M. gigas are distinguished typically on the basis of size and pelage color (cf. Monadjem 2019). They are also known to occur together in the same cave (Shimoni Cave in Kwale, Kenya; Webala et al. 2019), where they utilize echolocation calls with different peak frequencies: vittata at 64–70 kHz and gigas at 53.4–54.8 kHz. Both in Africa and on Madagascar, these pairs of taxa appear to act as distinct species, but the monophyly evident in the cyt-b sequences (Figs 2, 4) disappears in the nuclear intron analyses. BPP analyses fail to resolve any of the Macronycteris species, and none appear as monophyletic in the concatenated intron analyses.

Our results clearly underscore the importance of using multilocus datasets to evaluate phylogenetic and phylogeographic relationships at the genus and species level in mammals. Use of a single genetic system may lead to widely divergent conclusions regarding species identity and distribution. Toews and Brelsford (2012) reviewed cases of mito-nuclear discordance in animals generally. Fully 18% of the cases they reviewed had discordant patterns of mitochondrial and nuclear DNA. In most cases, such patterns are attributable to adaptive introgression of mtDNA, demographic disparities, and sex-biased asymmetries; in some cases they found evidence for hybrid zone movement or human agency. Discordant patterns of variation between mitochondrial and nuclear DNA have been reported in at least six other families of bats (Nesi et al. 2011; Furman et al. 2014; Naidoo et al. 2016; Hassanin et al. 2018; Demos et al. 2019a; Gürün et al. 2019). Gürün et al. (2019) implicated the role of sex-biased dispersal in causing such discordance, male dispersal spreading nuclear variation farther and faster than the movement of mitochondria. This may be a more general pattern in bats (see also Demos et al. 2019b). To understand the processes responsible for these discordant patterns of genome evolution, extensive genomic sampling and far fuller knowledge of natural history will be required.