Genetic, morphological and acoustic differentiation of African trident bats (Rhinonycteridae: Triaenops)

Abstract

Rhinonycteridae (trident bats) are a small Palaeotropical family of insectivorous bats allied to Hipposideridae. Their taxonomy has been in a state of flux. Here, we use mitochondrial and nuclear sequences to evaluate species relationships, confirming the monophyly of both Triaenops and Paratriaenops. Although most Triaenops afer specimens are recovered as a group, mitochondrial analyses strongly support some Kenyan individuals as members of Triaenops persicus. Analyses of four nuclear introns (ACOX2, COPS7A, RODGI and STAT5A) strongly support the mitochondrial topology. Morphometric analysis of the skull, external morphology and echolocation calls confirm that the Triaenops from the Rift Valley in Kenya (Nakuru, Baringo and Pokot counties) are distinct from typical T. afer in coastal (Kilifi and Kwale counties) or interior (Laikipia and Makueni counties) colonies. We interpret these analyses to indicate that two species of Triaenops occur in East Africa: T. afer in coastal regions along the Indian Ocean and in the highlands of central Kenya and Ethiopia, and T. persicus in the Rift Valley of Kenya. Although they appear widely disjunct from Middle Eastern populations, Kenyan T. persicus might be more widely distributed in the Rift Valley; they are somewhat differentiated from Middle Eastern populations in terms of both cranial morphology and vocalizations.

bioacoustics, cranial, DNA, geographical variation, phylogeny, species delimitation

INTRODUCTION

Rhinonycteridae is a newly recognized family of Palaeotropical insectivorous bats. Long grouped with the family Hipposideridae, the ‘trident bats’ (Armstrong et al., 2016) were nevertheless distinctive enough to be separated by name in earlier classifications (e.g. subtribe Rhinonycterina; Koopman, 1994). However, an extensive genetic analysis indicated that the Rhinonycteridae split from the Hipposideridae during the Eocene (~39 Mya; Foley et al., 2015), shortly after their divergence from Rhinolophidae and approximately coeval with the splits of other bat families. Given that the taxonomy of the group is central to the principal subject of this study, we review it briefly here.

Taxonomic history

The trident bats include the extant genera Cloeotis Thomas, 1901, Paratriaenops Benda & Vallo, 2009, Rhinonicteris Gray, 1847 and Triaenops Dobson, 1871, in addition to the fossil genera †Archerops Hand & Kirsch, 2003, †Brachipposideros Sigé, 1968, †Brevipalatus Hand, 2005 and †Xenorhinos Hand, 1998 (Wilson et al., 2016). Rhinonicteris aurantia (J.E. Gray, 1845) is endemic to northern Australia, Paratriaenops is endemic to Madagascar [Paratriaenops auritus (Grandidier, 1912) and Paratriaenops furculus (Trouessart, 1906)] and the Seychelles [Paratriaenops pauliani (Goodman & Ranivo, 2008)], and the sole species of Cloeotis is restricted to eastern portions of sub-Saharan Africa, ranging across scattered localities from Kenya to South Africa (Monadjem et al., 2010; Benda, 2019).

Triaenops is the most species-rich genus of trident bats, and it has the largest geographical range and the most complicated taxonomic history. Discovery of the Middle Eastern species Triaenops persicus Dobson, 1871 (type locality, near Shiraz, Iran) was quickly followed by the naming of the continental African species Triaenops aferPeters, 1877 (type locality, Mombasa, Kenya). Until the mid-20^th century, these Palaearctic and Afrotropical taxa, respectively, were distinguished as distinct species. Dorst (1948) provided a detailed diagnosis to distinguish T. afer from T. persicus on the basis of morphology. His treatment of the two as distinct was later challenged by Harrison (1964), whose wider sampling of the Palaearctic range of T. persicus led him to regard them as conspecific. Conspecificity of the two was later affirmed by Hayman & Hill (1971) and Harrison & Bates (1991), meaning that all African, Arabian and south-west Asian populations of Triaenops were assigned to T. persicus: T. p. persicus in south-west Asia and the Arabian Peninsula, T. p. afer in eastern Africa and T. p. majusculus Aellen & Brosset, 1968 in the Republic of Congo and Angola (Simmons, 2005).

The Indian Ocean species allocated to Triaenops were revised by Ranivo & Goodman (2006), who concluded that Triaenops auritus G. Grandidier, 1912 is distinct from Triaenops furcula Trouessart, 1906. Goodman & Ranivo (2008) named the trident bat of the Aldabra Atoll in the western Seychelles (Triaenops pauliani), and Goodman & Ranivo (2009) clarified the nomenclatural status of Malagasy Triaenops, including Triaenops humblotiMilne-Edwards, 1881 and Triaenops rufusMilne-Edwards, 1881. They proposed the replacement name Triaenops menamena Goodman & Ranivo, 2009 for T. rufus, which had been based mistakenly on specimens from Yemen and not on the Malagasy endemic to which it was intended.

Benda & Vallo (2009) revised all Triaenops species s.l. using morphological and molecular characters. They described the new genus Paratriaenops for species endemic to Madagascar (Paratriaenops auritus and Paratriaenops furcula) and the Seychelles (Paratriaenops pauliani). However, the remaining Madagascar endemic (T. menamena) clustered loosely with specimens from mainland Africa and the Middle East and was retained in Triaenops. Their multivariate analyses of craniodental morphology indicated that T. persicus was monotypic, with a range extending from Pakistan to Yemen. They also discovered and described a new diminutive species from the Arabian Peninsula, Triaenops parvus Benda & Vallo, 2009. All mainland African populations of Triaenops were allocated to T. afer, whereas T. menamena and †Triaenops goodmani Samonds, 2007 (known from late Pleistocene fossils; Samonds, 2007) are endemic to Madagascar. Although Triaenops afer majusculus Aellen & Brossett, 1968 (type locality, Grotte de Meya-Nzouari, Kouilou, Republic of Congo) had been used to distinguish the widely disjunct population in Congo and Angola (Happold, 2013), Benda & Vallo (2009) found no basis for distinguishing it from T. afer, which they also regarded as monotypic. However, the genetic analyses used to guide their taxonomic decisions were limited to a 731 bp sequence of mitochondrial DNA (Cytb), making their conclusions provisional.

The following questions remain. With broader geographical and genetic sampling, do T. afer, T. menamena, T. parvus and T. persicus collectively constitute a monophyletic group? Is each recognized species also monophyletic? In particular, do all African Triaenops represent T. afer? Do distinctive morphologies and vocalizations accompany genetic differentiation in this group? To answer these questions, we undertook analyses of molecular, morphological and vocalization data from a broader range of Palaearctic, Afrotropical and Malagasy Triaenops.

MATERIAL AND METHODS

Selection of taxa for genetic analyses

The genetic dataset is based on 83 rhinonycterid individuals. We generated original genetic data from 65 individuals collected at 22 georeferenced localities and supplemented them with 18 mitochondrial sequences from nine localities downloaded from GenBank (Fig. 1; Supporting Information, Appendix S1). All individuals were sequenced for cytochrome b (Cytb) as an initial assessment of genetic diversity. The bats newly sequenced in this study were obtained over several decades, during the course of chiropteran surveys across sub-Saharan Africa and Madagascar, with geographically dense sampling in East Africa. Initial assignment of East African specimens to species was based on meristic, mensural and qualitative characters published in the bat keys of Patterson & Webala (2012). Collection protocols followed mammal guidelines for the use of wild mammals in research (Sikes & the Animal Care and Use Committee of the American Society of Mammalogists, 2016), and the most recent accessions were approved under Insitutional Animal Care and Use Committee #2012-003 of the Field Museum of Natural History (Chicago, IL, USA). The remaining African rhinonycterid, Cloeotis percivali Thomas, 1901, was included for context, but only Cytb records were available on GenBank. Lacking nuclear intron data, we draw no conclusions from its placement and do not discuss Cloeotis in this paper. See the Supporting Information (Appendix S1) for voucher numbers and institutions, locality data and GenBank accession numbers.

Figure 1.

Distribution of genetic samples used in this study: inverted triangles, Paratriaenops auritus; triangles, Paratriaenops furcula; open squares, Triaenops afer; asterisk, Triaenops afer*; black circles, Triaenops menamena; blue circles, Triaenops parvus; open stars, Triaenops persicus.

DNA extraction, polymerase chain reaction and sequencing

Genomic DNA was extracted from fresh tissues using the DNeasy Blood and Tissue Kit (Qiagen) and the Wizard SV 96 Genomic DNA Purification System (Promega Corporation, WI, USA). Specimens with frozen tissue available were sequenced for mitochondrial Cytb using the primer pair LGL 765F and LGL 766R (Bickham et al., 1995, 2004) and four unlinked autosomal nuclear introns: ACOX2 intron 3, COPS7A intron 4, ROGDI intron 7 (Salicini et al., 2011) and STAT5A (Matthee et al., 2001). Primer information is contained in the Supporting Information (Table S1). Owing to the unavailability of frozen tissues of topotypic Triaenops persicus, we selected two specimens prepared as study skins (dating to 1963) and excised ~4 mm × ~4 mm of uropatagium for ultraconserved element sequencing and assembly.

Genomic DNA was extracted using the phenol-chloroform extraction protocol presented in the Supplementary Data S2 of McDonough et al. (2018), optimized for maximizing DNA yield and quality from historical tissues suitable for high-throughput sequencing methods (e.g. ultraconserved element sequencing). Sample concentrations were increased using Amicon Ultra-4 columns with Ultracel 30 membranes (Millipore, Fischer) bringing the final volume of DNA extract to 35 µL with total genomic DNA between 1000 and 3000 ng. Samples were submitted to Rapid Genomics (Gainesville, FL, USA) for genomic library preparation using an ultraconserved element MYbaits probe set (MYcroarray) that targets ~2500 loci. Enriched, pooled libraries were sequenced on an Illumina HiSeq platform (150 bp, paired-end).

We subjected fasta format reads delivered by Rapid Genomics to quality control and trimming using BBDuk in the Geneious Prime platform (v.2020.1.2; Biomatters). Reads with base pairs < Q20 (Q20, base call accuracy of 99%) and total length < 30 bp were discarded. Cytb sequence data were isolated by mapping reads to a reference sequence of Triaenops from our Sanger-generated Cytb alignment in Geneious. Outgroups included the genera Macronycteris and Hipposideros in the sister family Hipposideridae. Polymerase chain reactions, thermal cycling settings and sequencing were identical to the settings described by Demos et al. (2018) and Patterson et al. (2018). Chromatograms were assembled and edited using Geneious Pro v.11.1.5 (Biomatters). All sequences for each locus were aligned using Muscle (Edgar, 2004) with default settings in Geneious. Protein-coding sequences of Cytb were translated to amino acids to determine codon positions and examined for any premature stop codons and frameshifts. Several gaps were incorporated in the nuclear intron alignments whose positions were unambiguous. We resolved nuclear DNA to haplotypes in Phase (Stephens et al., 2001) and set the probability threshold to 0.70 following Garrick et al., (2010). Phase output files were formatted and assembled using the seqphase online platform (Flot, 2010).

Sequence alignments used in this study have been made available on the FIGSHARE data repository (https://doi.org/10.6084/m9.figshare.12811721.v1). All newly generated sequences have been deposited in GenBank, with accession numbers MT777711–MT777842 (see also Supporting Information, Appendix S1).

Phylogenetics, haplotype networks and species delimitation

jModelTest 2 (Darriba et al., 2012) on CIPRES Science Gateway v.3.3 (Miller et al., 2010) was used to determine nucleotide substitution models that best fit the data using the Akaike information criterion corrected for small sample size (AICc) for separate alignments of the four nuclear introns. Partitionfinder2 (Lanfear et al., 2016) on CIPRES was used to determine the best-fitting partitioning scheme and nucleotide substitution models for both the Cytb alignment and the concatenated alignment of four nuclear introns using the AICc under the ‘greedy’ search algorithm. Uncorrected Cytb sequence divergences (p-distances) between and within species were calculated in MEGA X v.10.0.5 (Kumar et al., 2018). We inferred maximum likelihood (ML) phylogenies in IQ-TREE v.1.6.10 (Nguyen et al., 2015; Chernomor et al., 2016) on the CIPRES portal for Cytb and the concatenated intron alignment. We inferred Bayesian phylogenies for Cytb and the concatenated intron alignment in MrBayes v.3.2.7 (Ronquist et al., 2012) on the CIPRES portal using the same set of genes as the ML analyses. We ran the two independent tree searches 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 using default heating values and sampled every 1000^th generation. A conservative 20% burn-in was applied, and convergence for all parameters was assessed in TRACER v.1.7 (Rambaut et al., 2018). Majority-rule consensus trees were assembled for each Bayesian analysis.

We inferred a species tree in BEAST v.2.6 (Bouckaert et al., 2019) using the StarBEAST2 algorithm (Ogilvie et al., 2017). Species tree analyses were 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 model of population size. Four independent replicates were run with random starting seeds and chain lengths of 1 × 10⁸ generations, with parameters sampled every 5000 steps. For the StarBEAST2 analyses, evidence of convergence and model parameter posterior distribution stationarity were assessed based on effective sample size values > 200 and examination of trace files in TRACER v.1.7. The burn-in was set at 10%, and independent replicates were assembled using LogCombiner v.2.5.1 and TreeAnnotator v.2.5.1 (Bouckaert et al., 2019).

Haplotype networks for Cytb were inferred using the median-joining network algorithm in POPART v.1.7 (Leigh & Bryant, 2015). Separate analyses were carried out for T. afer + T. parvus + T. persicus and for T. menamena. Based on the well-supported clades obtained in the Cytb phylogenetic analyses and the availability of intron samples, a species delimitation scenario with five candidate species [T. afer, T. afer* (henceforth used to designate Rift Valley samples of Triaenops from Baringo, Nakuru and West Pokot counties in Kenya), T. menamena, Paratriaenops auritus and P. furcula] was tested. We inferred the evolutionary isolation of their gene pools using the phased nuclear DNA dataset (ACOX2, COPS7A, ROGDI and STAT5A; 17 individuals) for joint independent species delimitation and species tree estimation using a multispecies coalescent model in the software BPP v.3.3 (Yang & Rannala, 2014; Rannala & Yang, 2017). Species memberships for BPP were identical to the assignments of individuals to species in the species tree analyses. The validity of our assignment of individuals to species was tested with the guide-tree-free algorithm (A11) in BPP. Given that delimitation probability in BPP is sensitive to parameter selection (Leaché & Fujita, 2010; Yang, 2015), we evaluated two replicates for each of four different combinations of divergence depths and effective population sizes priors (τ and θ, respectively; see Supporting Information, Table S2). Two independent Markov chain Monte Carlo chains were run for 5 × 10⁴ generations. Burn-in was set at 20%, and samples were drawn every 50^th generation. Species were considered to be well supported when delimitation posterior probability (PP) estimates were ≥ 0.95 under all four prior combinations (Supporting Information, Table S2).

Craniodental and mandibular analyses

We analysed the morphology of 148 Triaenops specimens distributed in Iran, Ethiopia, Kenya, Tanzania and Madagascar (Supporting Information, Appendix S1). Specimens are housed in six natural history museums: Field Museum of Natural History, Chicago, IL, USA (FMNH); Natural History Museum of Geneva, Switzerland (MGNG); National Museum of Natural History, Paris, France (MNHN); National Museum Prague, Czech Republic (NMP); Royal Ontario Museum, Toronto, Ontario, Canada (ROM); and Natural History Museum, Berlin, Germany (ZMB).

Using previously published measurements, our analyses included the holotype and paratypes of T. a. majusculus (holotype, MNHN 1968-412; paratypes, MNHN MP 19-04-64-24, MNHN MP 19-06-64-04, MHNG MG 1074.41, MHNG MG 1074.42, MHNG MG 1074.43, MHNG MG 1074.44, MHNG MG 1074.45 and MHNG MG 1074.46). A complete list of specimens used in the various analyses and the person responsible for measurements of craniodental, external and acoustic variables is provided in the Supporting Information (Appendix S1).

Morphometric data were collected only from adults, corresponding to those specimens with completely erupted and partly worn dentitions. Twelve craniodental and mandibular linear measurements were recorded, following Velazco & Gardner (2012), using digital callipers with 0.01 mm resolution: GLS, greatest length of skull; CIL, condyloincisive length; CCL, condylocanine length; BB, braincase breadth; ZB, zygomatic breadth; PB, postorbital breadth; MSTW, mastoid width; MTRL, maxillary toothrow length; MLTRL, molariform toothrow length; M2M2, width at M2; DENL, dentary length; and MANDL, mandibular toothrow length. Measurements are defined in the Supporting Information (Table S3).

Collecting localities were georeferenced (see Supporting Information, Appendix S1). We divided the specimens into 12 groups, based on their geographical locations and performed an exploratory data analysis (Tukey, 1977) using the log₁₀-transformed and standardized data to check the normality of each measurement and to highlight potential outliers. The effect of sex, measurer and their possible interaction were evaluated using multivariate analysis of variance (MANOVA) whenever possible for the available sample size; models were chosen on the basis of Wilk’s lambda statistic (P < 0.05). When these components of variation had significant effects, subsequent analyses used the residuals of a general linear model treating the morphometric variables as dependent and the significant sources of variation as independent. In cases where no effect was detected, the subsequent analysis was performed using log₁₀-transformed and standardized data.

Genetic information was available for a subset of the specimens we analysed morphologically and was used to assign specimens to species. To classify specimens lacking genetic information (‘unknown samples’), we performed a linear discriminant analysis using the individuals classified by Cytb gene sequences as training sets in the cross-validation tests and considering the linear discriminant analysis posterior probabilities in the classification decision. Principal components analysis (PCA) was also used to visualize the distribution and overlap of the ‘unknown samples’ in the morphospace of genotyped specimens.

After establishing the categories of all specimens, we estimated means, ranges and standard errors of the variables for each category. We performed a PCA to explore which morphological variables contributed the most to the variation of each axis and to identify which variables were most informative to discriminate populations and species in the subsequent discriminant function analysis (DFA). A DFA was applied on each group to test whether the morphometric variables could discriminate individuals according to the a priori geographical hypothesis based on molecular results.

We investigated gaps in the pattern of morphological variation using the statistical approach proposed by Zapata & Jiménez (2012), which takes into account both the morphological traits (craniodental and mandibular measurements) and geographical locations of individual specimens. This approach allowed us to assess the strength of the evidence supporting a gap in morphological variation between two hypothesized species and to investigate whether a morphological discontinuity between two hypothesized species could be explained by an alternative hypothesis of geographical variation, rather than a boundary between species.

Two or more local maxima (or modes) in the distribution of morphological variation might support the hypothesis that there are two species in a geographical locality (Futuyma, 1998). Assuming that the bimodal or multimodal distributions from continuous traits do not result from polymorphisms, ontogenetic variation or phenotypic plasticity, we tested the hypothesis of species limits by assessing the number of modes in the distribution of morphological traits of two hypothesized species using the probability density function. The probability density was then evaluated along the ridgeline manifold (α), a curve that contains all critical points (minima, maxima and saddles) in addition to the ridges of the density. To assess how distinct two hypothesized species are in their morphologies, we established a frequency cut-off following Wiens & Servedio (2000) to examine overlap of ellipsoidal tolerance regions (proportions β). If the proportions covered by tolerance regions were smaller than the frequency cut-off, we considered the two hypothesized taxa sufficiently distinct to support the hypothesis of a species limit. In order to confront a species limits hypothesis against the alternative of geographical differentiation within a single species, we: (1) estimated a geographical distance matrix using geographical coordinates of the collection localities; (2) extracted the eigenvectors of a principal coordinates analysis (PCoord or PCoA); and (3) used them as explanatory variables in a redundancy analysis (RDA). In the RDA, the multivariate morphological measurements were used as response variables and the geographical coordinates of the specimens as independent ones. Statistically significant results provide support for the hypothesis of species limits and indicate that the hypothesis of geographical variation within a species cannot explain the morphological discontinuity under consideration. Pairwise comparisons were performed among the following hypothesized species: (1) T. afer* vs. T. afer (all remaining African samples); and (2) T. afer* vs. T. persicus.

Finally, using measurements published in the description of T. a. majusculus (Aellen & Brosset, 1968), we included its holotype and paratypes in our analysis. Incomplete data for this taxon caused us to reduce the original set of 12 morphological measurements to a subset of nine and perform a PCA to visualize the morphospace occupied by T. a. majusculus. All analyses and graphs were carried out using R v.3.5.0 (R Development Core Team, 2019).

External morphology and analysis

Five external variables were taken in the field, typically with millimetre rulers, at the time of collection: total length (TTL), tail length (TL), hind foot length (RHF), ear length (EL) and forearm length (FA). Bat body mass was also recorded in the field, typically with Pesola balances.

We analysed the external morphology of 122 adult Triaenops (54 females and 68 males). These included T. afer from Kilifi, Kwale, Laikipia and Taita-Taveta counties in Kenya, T. afer* from Baringo and Nakuru counties in Kenya and T. persicus from Iran (Supporting Information, Appendix S1). The effect of sexual dimorphism was investigated using a MANOVA on the log₁₀-transformed variables. Following the same procedures described above for skull morphometrics, we used a PCA to visualize the distribution and overlap of ‘unknown samples’ lacking genetic identification in the morphospace occupied by genotyped specimens and estimated the mean, range and standard error of each variable for each taxon. Where ANOVAs indicated a significant effect of taxon treated as a grouping variable, we used Tukey’s post hoc tests to determine which groups differ significantly from others (Day & Quinn, 1989).

Echolocation call recordings and analysis

We recorded echolocation calls shortly after capture from individual hand-held Triaenops using a hand-held ultrasound detector (Pettersson D1000X; Pettersson Elektronik AB, Uppsala, Sweden; 384 kHz sampling rate, 16 bit resolution). Triaenops use a high duty-cycle form of echolocation, dominated by a narrow frequency band (‘constant frequency’) to which their hearing is attuned (see Taylor et al., 2005; Webala et al., 2019). Flight is unnecessary for these bats to generate characteristic calls (Webala et al., 2019). For sound analysis, a customized 512-point fast Fourier transform was used with a Hanning window for both spectrograms and the power spectrum. Following Jung et al. (2014), we characterized echolocation calls by measuring the peak frequency or frequency with maximum energy (FME), maximum frequency (StartF) and minimum frequency (EndF) using Kaleidoscope v.3.1.4b (Wildlife Acoustics, USA). The mean of ten calls with the best signal-to-noise ratios was measured for each bat. Vocalizations of 100 Triaenops individuals were recorded by PWW from 2012 to 2016 in Kilifi, Kwale, Laikipia and Nakuru counties in Kenya (Supporting Information, Appendix S1). Analyses of echolocation calls of T. persicus relied on published information (Benda et al., 2012) that used comparable equipment and procedures.

We log₁₀-transformed the variables and investigated the effect of sex following the same multivariate procedures described above for morphometric data. Owing to significant sex differences, separate analyses were conducted on males and females. We estimated the means, ranges and standard errors of the three call variables for bats from each county. Finally, we used one-way ANOVAs to test whether bats from each county had equal call frequencies, followed by Tukey’s post hoc test to determine which groups differed significantly from others.

RESULTS

DNA sequence characteristics

The ML and Bayesian inference (BI) gene tree analyses are based on 84 Cytb sequences that range in length from 728 to 1140 bp with 96% coverage (for information on these sequences, see Supporting Information, Appendix S1). We also include short sequences (104 and 168 bp, respectively) obtained from two museum skins of T. persicus from Iran (Supplemental Information, Appendix S1). Given that they were too short to archive in GenBank or use in the substitution network analyses, these sequences were used only in the ML and Bayesian phylogenetic analyses and are included in the 86-sequence FIGSHARE alignment (https://doi.org/10.6084/m9.figshare.12811721.v1). The numbers of base pairs for the sequence alignments used in Bayesian species tree analyses are as follows: ACOX2 (N = 14), 442–453 bp; COPS7A (N = 17), 537–652 bp; ROGDI (N = 16), 347–394 bp; STAT5A (N = 16), 426–511 bp; and in the four-intron concatenated alignment (N = 18), 1434–1984 bp. Introns and their best-supported substitution models estimated by jModelTest 2 are as follows: ACOX2, K80+G; COPS7A, GTR; ROGDI, HKY; and STAT5A, HKY. The best-supported partitioned substitution models estimated in Partitionfinder2 are, for Cytb: TRNEF+I+R, codon position 1; HKY+I, codon position 2; and GTR+I+R, codon position 3; and for the concatenated four-intron alignment: ACOX2 and COPS7A, TVM+G; and ROGDI and STAT5A, K81UF. Uncorrected Cytb p-distances among Triaenops species averaged 0.07 and ranged from 0.064 to 0.076 between pairs. Paratriaenops auritus and P. furcula were separated by 0.045. Average within-species distances ranged from 0.000 to 0.007 across species (Table 1).

Table 1.

Uncorrected Cytb p-distances within (italic) and between (below diagonal) Rhinonycteridae clades

Species	Paratriaenops auritus	Paratriaenops furcula	Triaenops afer	Triaenops afer*	Triaenops menamena	Triaenops parvus	Triaenops persicus
Paratriaenops auritus	0
Paratriaenops furcula	0.045	0.002
Triaenops afer	0.204	0.197	0.007
Triaenops afer*	0.203	0.199	0.072	0.007
Triaenops menamena	0.206	0.2	0.064	0.066	0.007
Triaenops parvus	0.218	0.211	0.076	0.063	0.077	0.002
Triaenops persicus	0.203	0.196	0.068	0.011	0.067	0.062	0.001

Analyses of p-distances were conducted using MEGA v.10.1.7. This analysis involved 81 nucleotide sequences. All gaps/missing data were removed using the pairwise deletion option.

*This group of Rift Valley Triaenops differs markedy from typical T. afer.

Mitochondrial genetic analyses

The MrBayes Markov chain Monte Carlo analyses converge successfully, with all parameters achieving effective sample sizes > 200. The ML and BI phylogenies are identical at deeper nodes in their trees, with minor differences for shallower relationships; only the ML topology is shown (Fig. 2). Cloeotis percivali is sister to Triaenops + Paratriaenops, and both Triaenops and Paratriaenops are strongly supported as monophyletic. Within Triaenops, T. menamena is poorly supported as monophyletic and appears as sister to the three other Triaenops species. The relationships among T. afer, T. persicus and T. parvus are mostly unresolved. Only T. parvus is well supported as monophyletic. In contrast, Triaenops sequences from the Rift Valley in Kenya (which we have designated T. afer*) are strongly supported in a clade with T. persicus from Yemen; T. afer* and T. persicus are minimally divergent (0.011) from each other (Table 1). Notably, despite the shortness of their sequences, the two T. persicus from Iran are confidently recovered with these Kenyan and Yemeni specimens. No geographical structure is apparent among the remaining East African T. afer, which range from Ethiopia through Kenya to central Tanzania. Both P. auritus and P. furcula are strongly supported as monophyletic in the BI tree, although P. furcula is weakly supported in the ML tree.

Figure 2.

Maximum likelihood phylogeny of mitochondrial cytochrome b sequences of Rhinonycteridae. The phylogeny was inferred in IQ-TREE, and its topology was similar to the Bayesian phylogeny calculated in MrBayes. Filled red circles on nodes denote bootstrap values ≥ 70% and Bayesian posterior probabilities (PP) ≥ 0.95. Filled black circles indicate bootstrap ≥ 70% and PP < 0.95. Open circles indicate bootstrap < 70% and PP ≥ 0.95. Unmarked nodes indicate bootstrap < 70% and PP < 0.95.

The Cytb haplotype network for T. afer, T. afer*, T. parvus and T. persicus (Fig. 3) confirms the close relationship between T. afer* and T. persicus and the strong distinction between T. afer* and T. afer. None of these taxa shares haplotypes, although those of T. afer* and T. persicus differ by only a few base pairs. The absence of geographical structure among T. afer populations is remarkable. Individuals from Fikirini Cave (in Kwale County, Kenya) encompass most of the documented haplotypic diversity of the species (Fig. 3). The Yemeni samples of T. persicus are genetically closer to some T. afer* haplotypes than the latter are to each other. The T. menamena Cytb haplotype network also shows a lack of genetic structure, although sample sizes are small (Supporting Information, Fig. S1). As in T. afer, single colonies of T. menamena (e.g. Ankarana) encompass most of the documented haplotypic variation.

Figure 3.

Substitution network plot for Cytb inferred in POPART v.1.7 for Triaenops afer, Triaenops afer*, Triaenops parvus and Triaenops persicus.

Nuclear intron genetic analyses

Using the concatenated nuclear intron alignment, the MrBayes Markov chain Monte Carlo analysis converged successfully, with all parameters achieving effective sample sizes > 200. The ML and BI phylogenies were identical at deeper nodes in their trees, with minor differences at shallower relationships; only the ML topology is shown (Fig. 4). As in the Cytb phylogeny, Triaenops and Paratriaenops are strongly supported as both sisters and monophyletic. Although introns were lacking for Palaearctic T. persicus, T. afer, T. afer* and T. menamena are each recovered as monophyletic, T. afer more weakly so. Triaenops menamena is strongly recovered as sister to T. afer* + T. afer. However, the ML and BI phylogenies based on the concatenated nuclear intron alignment offer no support for the monophyly of either P. auritus or P. furcula.

Figure 4.

Bayesian phylogeny of Rhinonycteridae based on four nuclear introns. The phylogeny was inferred in MrBayes, and its topology closely resembled the maximum likelihood phylogeny calculated in IQ-TREE. Filled red circles on nodes denote bootstrap bootstrap values ≥ 70% and Bayesian posterior probabilities (PP) ≥ 0.95. Filled black circles indicate bootstrap ≥ 70% and PP < 0.95. Open circles indicate bootstrap < 70% and PP ≥ 0.95. Unmarked nodes indicate bootstrap < 70% and PP < 0.95.

The StarBEAST analysis converged successfully, with effective sample size values > 200 for all parameters. The resulting species tree supports the monophyly of the tested Rhinonycteridae and of both Triaenops and Paratriaenops (Supporting Information, Fig. S2). However, relationships among Triaenops species are not resolved. The three putative species of Triaenops are recovered in a near trichotomy, again with T. menamena sister to T. afer* + T. afer, but with weaker support. Apparent strong support for the sister relationship of P. auritus and P. furcula is an artefact of the a priori assignment of individuals to species in species tree analyses. The extremely short branch lengths of these two taxa do not support their status as evolutionarily independent, at least on the basis of the four introns used in this study.

Results from the replicated BPP analyses show strong support for all three putative Triaenops species tested. However, no support was inferred for P. auritus or P. furcula as independent species under any combination of priors. Instead, the BPP analyses offer some support for considering P. auritus + P. furcula as a single species (Table 2).

Table 2.

Species delimitation results based on the four-intron dataset for Afrotropical Rhinonycteridae with four different parameter sets (PS1–PS4)

Putative species	PS1	PS2	PS3	PS4
Triaenops afer	1	1	1	1
Triaenops afer*	1	1	1	1
Triaenops menamena	1	1	1	1
Paratriaenops auritus	0.057	0.346	0.032	0.003
Paratriaenops furcula	0.057	0.346	0.032	0.003
Paratriaenops auritus+ furcula	0.943	0.653	0.968	0.997

Bold type indicates that delimitation of the putative species was achieved with ≥ 0.95 posterior probability in BPP. See Material and Methods section for parameter details.

*Specimens previously assigned to Triaenops afer that are supported in this study as a widely disjunct population of Triaenops persicus.

Craniodental and mandibular analyses

Five morphological groupings were identified in the cross-validation tests used to classify specimens that lacked genetic identification: (1) T. afer from coastal and inland sites in Kenya (Kilifi, Kwale, Laikipia and Makueni counties), plus Tanzania and Ethiopia; (2) T. afer from the Tsavo region of Kenya (Taita-Taveta County); (3) T. afer* from Rift Valley sites in Kenya (Nakuru, Baringo and West Pokot counties); (4) T. persicus from Iran; and (5) T. menamena from Madagascar (Supporting Information, Fig. S3). The Supporting Information (Fig. S3) shows the distribution of specimens lacking Cytb information (unknown group in each plot represented in black) in the morphospace defined by genotyped specimens.

After cross-validation of non-genotyped specimens, a PCA performed on all samples showed that principal component (PC) 1 accounts for 78.14% of the variation, PC2 for 5.53% and PC3 for 4.38%, together accounting for > 88.05% of the variation (Supporting Information, Table S4; Fig. S4). Principal component 1 summarizes variation associated with size, both isometric and allometric, with all individual vectors in the same direction, as indicated by their negative scores (Supporting Information, Table S4). The variables that explain most of the variance in PC1 are mainly associated with cranial and mandibular lengths (CCL, GLS, MTRL, MANDL and DENL), whereas the length of the rostrum (MLTRL) and breadth of the braincase (BB) account for variance in PC2.

The DFA performed on all samples produced high classification rates (> 0.7) of specimens in each of the a priori groups except for the T. afer from the Taita-Taveta group (0.58) (Supporting Information, Table S5). The first two axes of the DFA explain 94.07% of the variation and clearly separate three groups (T. menamena in Madagascar, T. afer* and T. persicus in Iran). In contrast, typical T. afer and T. afer from Taita-Taveta overlap slightly in their discriminant function (DF) scores, suggesting higher similarity between these groups (Fig. 5). The MANOVA showed a significant effect of sexual dimorphism, but no effect of the measurer or its interaction with sex was observed in the data.

Figure 5.

Linear discriminant function axes (DF1 and DF2) of craniodental and mandibular variables performed for all Triaenops samples after cross-validation tests. The density plot for the first and the second linear discriminants (LD1 and LD2, respectively) are presented on the right.

With significant sexual dimorphism, variable means, ranges and standard errors are tabulated for each morphogroup and sex separately (Tables 3 and 4; Supporting Information, Figs S5–S8). Note that T. afer* females and males presented the largest mean values for all morphological variables (Supporting Information, Figs S5–S8). A PCA based on nine craniodental and mandibular variables that included the hypodigm of T. a. majusculus indicated overlap of majusculus with all other morphogroups, except with T. afer* (Supporting Information, Fig. S9). Analyses of the external measurements showed closely comparable results; these are presented in the Supporting Information (Supplementary Material—External Morphology).

Table 3.

Descriptive statistics of craniodental and mandibular measurements for Triaenops females (in millimetres)

Measurement	Statistics	T. afer	T. afer*	T. persicus	T. afer* + T. persicus	F3,45	F1,47
	Number of females	43	3	3	6
GLS	Mean	18.72A	20.37B	19.37C	19.87B	39.91	52.14
	Min–max	18.18–20.04	20.01–20.57	19.17–19.69	19.17–20.57
	SE	0.04	0.18	0.16	0.24
CIL	Mean	16.54A	17.85B	17.12C	17.48B	21.61	33.61
	Min–max	15.74–17.74	17.36–18.17	16.97–17.27	16.97–18.17
	SE	0.05	0.25	0.08	0.2
CCL	Mean	15.8A, C	17.3B	16.23C	16.76B	32.07	36.48
	Min–max	15.37–17.11	16.87–18.02	16.10–16.33	16.1–18.02
	SE	0.04	0.36	0.06	0.29
BB	Mean	7.1A	7.75B,C	7.45C	7.60B	23.52	40.43
	Min–max	6.80–7.36	7.38–8.31	7.38–7.58	7.38–8.31
	SE	0.02	0.28	0.06	0.14
ZB	Mean	8.58A, C	9.38B	8.73C	9.05B	35.25	31.11
	Min–max	8.25–8.85	9.21–9.53	8.71–8.78	8.71–9.53
	SE	0.02	0.09	0.02	0.15
PB	Mean	2.73A, C	3.12B	2.79C	2.95B	11.09	11.67
	Min–max	2.45–3.03	2.99–3.25	2.73–2.87	2.73–3.25
	SE	0.02	0.07	0.04	0.08
MSTW	Mean	7.46A	8.04B	7.68C	7.86B	34.89	45.36
	Min–max	7.21–7.74	7.99–8.09	7.56–7.84	7.56–8.09
	SE	0.01	0.02	0.08	0.09
MTRL	Mean	6.19A	7.22B	6.59C	6.85B	45.67	62.37
	Min–max	5.82–6.54	6.85–7.28	6.55–6.66	6.55–7.28
	SE	0.02	0.13	0.03	0.13
MLTRL	Mean	4.64A	5.43B,C	5.17C	5.30B	65.02	115.2
	Min–max	4.42–4.91	5.23–5.63	5.06–5.28	5.06–5.63
	SE	0.01	0.11	0.06	0.08
M2M2	Mean	6.36A	6.77B,C	6.64C	6.70B	18.43	35.13
	Min–max	6.14–6.67	6.55–6.92	6.49–6.86	6.49–6.92
	SE	0.01	0.11	0.11	0.07
DENL	Mean	11.49A	12.71B	11.92C	12.31B	64.83	68.23
	Min–max	11.23–11.80	12.34–12.90	11.55–12.21	11.55–12.90
	SE	0.02	0.18	0.19	0.21
MANDL	Mean	6.71A	7.41B,C	7.23C	7.32B	21.86	42.81
	Min–max	6.19–7.50	7.25–7.51	7.20–7.28	7.20–7.51
	SE	0.03	0.08	0.02	0.05

Results are shown for Triaenops afer (Kenya: Kilifi, Kwale, Laikipia and Taita-Taveta counties, plus Tanzania and Ethiopia), T. afer* (Kenya: Baringo), T. persicus (Iran) and T. afer* + T. persicus (combined T. afer* and T. persicus). ANOVAs showed significant differences between all comparisons (P < 0.0001). F_3,45 corresponds to F-values from comparisons between T. afer, T. afer* and T. persicus. F_1,47 corresponds to F-values from comparisons between T. afer and T. afer* + T. persicus. Superscripts alongside means indicate group membership assessed by Tukey’s post hoc tests.

Abbreviations: Min‒max, range for minimum and maximum values; SE, standard error of the mean. Abbreviations for measurements are defined in the Material and Methods section.

Table 4.

Descriptive statistics of craniodental and mandibular measurements for Triaenops males (in millimetres)

Measurement	Statistics	T. afer	T. afer*	T. persicus	T. afer* + T. persicus	F4,69	F1,63
	Number of males	48	11	6	17
GLS	Mean	19.8A	21.58B	20.68C	21.26B	83.19	114.7
	Min–max	19.20–20.57	20.93–22.33	19.82–21.47	19.82–22.33
	SE	0.05	0.11	0.28	0.15
CIL	Mean	17.68A	19.16B	18.33C	18.86B	67.94	94.6
	Min–max	16.67–18.45	18.46–19.56	17.77–19.04	17.77–19.56
	SE	0.05	0.1	0.21	0.13
CCL	Mean	16.88A	18.26B	17.44C	17.97B	79.22	101
	Min–max	16.29–17.71	17.46–18.68	17.00–17.80	17.00–18.68
	SE	0.04	0.11	0.14	0.12
BB	Mean	7.22A	7.87B	7.55C	7.76B	89.32	128.4
	Min–max	6.94–7.64	7.61–8.32	7.32–7.84	7.32–8.32
	SE	0.01	0.06	0.07	0.06
ZB	Mean	9.07A, C	9.8B	9.2C	9.59B	50.97	50.06
	Min–max	8.67–9.56	9.53–10.24	9.08–9.34	9.08–10.24
	SE	0.03	0.06	0.04	0.08
PB	Mean	2.8A, C	3.1B	2.84C	3.01B	16.3	18.55
	Min–max	2.46–3.20	2.99–3.29	2.64–3.03	2.64–3.29
	SE	0.02	0.02	0.06	0.04
MSTW	Mean	7.69A	8.31B	8.01C	8.21B	62.35	97.24
	Min–max	7.30–8.18	8.07–8.56	7.73–8.37	7.73–8.56
	SE	0.02	0.05	0.09	0.05
MTRL	Mean	6.64A	7.34B	7.06C	7.24B	78.44	131.3
	Min–max	6.31–7.13	7.07–7.62	6.80–7.19	6.80–7.62
	SE	0.02	0.05	0.06	0.05
MLTRL	Mean	4.91A	5.48B, C	5.36C	5.44B	53.28	104.6
	Min–max	4.61–5.37	5.14–5.73	5.18–5.49	5.14–5.73
	SE	0.02	0.06	0.05	0.04
M2M2	Mean	6.5A	7.23B	6.96C	7.13B	84.84	144.8
	Min–max	6.23–6.91	6.97–7.53	6.83–7.09	6.83–7.53
	SE	0.02	0.05	0.04	0.04
DENL	Mean	12.27A	13.25B	12.63C	13.03B	50.2	66.69
	Min–max	11.78–13.04	12.63–13.64	12.24– 12.92	12.24–13.64
	SE	0.04	0.09	0.11	0.1
MANDL	Mean	7.19A	7.95B	7.64C	7.84B	91.4	147.7
	Min–max	6.92–7.56	7.63–8.15	7.27–7.88	7.27–8.15
	SE	0.02	0.04	0.095	0.05

Triaenops afer (Kenya: Kilifi, Kwale, Makueni and Taita-Taveta counties, plus Tanzania and Ethiopia), T. afer* (Kenya: Nakuru, Baringo and West Pokot), T. persicus (Iran) and T. afer* + T. persicus (combined T. afer* and T. persicus). ANOVAs showed significant differences between all comparisons (P < 0.0001). F_4,69 corresponds to F-values from comparisons between T. afer, T. afer* and T. persicus. F_1,63 corresponds to F-values from comparisons between T. afer and T. afer* + T. persicus. Superscripts alongside means indicate group membership assessed by Tukey’s post hoc tests.

Abbreviations: Min–max, range for minimum and maximum values; SE, standard error. Abbreviations for measurements are defined in the Material and Methods section.

The analysis of gaps in the pattern of morphological variation is summarized in the Supporting Information (Figs S10, S11; Table S6). The plot of f^ˆ(X) along the ridgeline manifold for T. afer and T. afer* presents a bimodal distribution (Supporting Information, Fig. S10, left), suggesting that there is a morphological gap between these populations. The corresponding plots of the proportions covered by tolerance regions reveals that phenotypic overlap is smaller than the frequency cut-off (Supporting Information, Fig. S10, right); T. afer* is sufficiently distinct from T. afer to support its recognition as a distinct species. After this result, we tested the species limit hypothesis against an alternative of geographical variation within a species. We found that the pattern of morphological variation in the pairwise comparison between T. afer and T. afer* was explained by a dummy variable representing a species boundary and by its interaction with the second spatial eigenvector (Supporting Information, Table S6).

Although the bimodality in the plot of f^ˆ(X) for T. persicus and T. afer* suggests a morphological gap separating these populations (Supporting Information, Fig. S11, left), the overlap in the distributions of both taxa failed to support the hypothesis of a species distinction between them (Supporting Information, Fig. S11, right).

Echolocation calls

A PCA performed for the whole sample (corrected for sexual dimorphism) and considering the three call frequency variables indicated that the first two PCs account for 99.8% of the variation and established two main groups: the first group is represented by T. afer*, which presents the highest scores on PC1; the other group is composed of T. afer from Kilifi, Kwale and Laikipia counties, which are variable and broadly overlapping (Fig. 6). An ANOVA registered significant sexual dimorphism in call frequencies; therefore, we tested each sex separately for significant differences in log₁₀-transformed frequencies among populations. Variable means, ranges and standard errors for each county separated by sex are presented in Table 5.

Table 5.

Descriptive statistics for call frequencies (in kilohertz) of females and males Triaenops from Kenya, grouped by counties

Call frequencies	Statistics	T. afer Kilifi	T. afer Kwale	T. afer Laikipia	F2,43	P-value
	Number of females	23	17	8
FME	Mean	87.42A	87.27 A	84.09B	7.27	0.0019**
	Min–max	81.94–91.68	83.51–90.45	81.65–86.07
	SE	0.54	0.5	0.45
StartF	Mean	88.54A	88.36A	84.88B	8.42	0.000818***
	Min–max	82.63–92.69	84.21–91.70	82.84–86.95
	SE	0.54	0.54	0.42
EndF	Mean	76.33A,B	78.75A	74.35B	4.704	0.0142*
	Min–max	71.60–82.63	72.61–85.46	71.63–80.66
	SE	0.63	1.07	1.02
Call frequencies	Statistics	T. afer Kilifi	T. afer Kwale	T. afer* Nakuru	F	P-value
	Number of males	20	24	8
FME	Mean	75.97A	76.6A	69.68B	91.74	< 2 × 10−16***
	Min–max	73.07–78.49	74.48–79.32	69.44–70.08
	SE	0.36	0.26	0.08
StartF	Mean	76.74A	77.41A	70.42B	71.65	2.83 × 10−15***
	Min–max	73.90–79.21	74.78–80.41	69.97–70.92
	SE	0.41	0.3	0.11
EndF	Mean	68.4A	67.77A	66.27A	2.236	0.118
	Min–max	62.03–73.03	64.04–71.02	63.88–68.53
	SE	0.65	0.42	0.54

Abbreviations: EndF, minimum frequency; FME, peak frequency; Min‒max, range of minimum and maximum values; SE, standard error (measure of the dispersion of sample means around the population mean); StartF, maximum frequency.

Superscript letters alongside means indicate the results for Tukey post hoc tests.

Included are F statistics and associated probability (P) from ANOVAs testing for differences among groups: *P < 0.01, **P < 0.001 and ***P < 0.0001.

Figure 6.

Principal components analysis based on call frequencies of recorded Triaenops individuals (EndF, minimum frequency; FME, peak frequency; StartF, maximum frequency) and grouped by counties. Triaenops afer (Kilifi, Kwale and Laikipia) and Triaenops afer* (Nakuru). The first two principal components (PCs) account for 99.8% of the variation.

For both males and females, ANOVAs reject the null hypothesis that Triaenops call frequencies from the different counties are equal. All three ANOVAs for females indicate that localities differ significantly. Tukey’s post hoc test for peak frequency (FME) and maximum frequency (StartF) indicate that significant differences are found in pairwise comparisons between Laikipia and Kilifi and between Laikipia and Kwale (Table 5); no difference is evident between Kwale and Kilifi females in either variable. Female Triaenops from Laikipia County call at peak frequencies significantly lower than those of Kwale and Kilifi females.

For male Triaenops, ANOVAs comparing Kilifi, Kwale and Nakuru counties are significant for peak and maximum frequencies but not for end frequency (Table 5). Tukey’s post hoc tests identify Nakuru males as calling at lower frequencies than either Kwale or Kilifi males. Both Kwale and Kilifi bats represent typical T. afer, whereas those from Nakuru represent T. afer*.

DISCUSSION

Application of names to African Triaenops

We have dubbed Triaenops from the Rift Valley of Kenya T. afer* for the purposes of this analysis, because their genetics, morphology and vocalizations all distinguish them from other African Triaenops matching the description of T. afer. Genetically, these bats scarcely differ from T. persicus in Yemen, differing by only 0.011 in Cytb. They are securely recovered in a clade with T. persicus, to the exclusion of all other Triaenops species (Fig. 2). Fewer base-pair substitutions separate Yemeni T. persicus samples from Kenyan T. afer* than separate some individuals of T. afer* (Fig. 3). Both T. persicus and T. afer* are separated from T. afer by ~7% sequence divergence in Cytb (Table 1). Although our nuclear intron analyses lacked samples of Middle Eastern T. persicus and T. parvus, the concatenated alignment of four nuclear introns reproduced the mitochondrial topology in all respects, with both T. afer and T. afer* well supported and reciprocally monophyletic (Fig. 4). Morphologically, T. afer and T. afer* are clearly distinguishable, with T. afer* significantly larger in both craniodental and external measurements (Supporting Information, Figs S5–S8, external measurements). The analysis of gaps in morphological variation found evidence for a species boundary between T. afer and T. afer*, but failed to find one between T. afer* and T. persicus (Supporting Information, Figs 10, 11). In keeping with its larger size (Jacobs et al., 2007), T. afer* echolocates at lower frequencies than T. afer (Table 5; Fig. 7). Thus, genetics, morphology and echolocation calls all demonstrate the existence of two species of Triaenops in Kenya. Although T. afer and T. afer* have not yet been recorded in sympatry, they certainly occur within 115 km of one another (Lolldaiga-Gilgil;, Supporting Information, Appendix S1).

Figure 7.

Density plot (top) and boxplot (bottom) showing the distribution pattern of call frequencies of recorded Triaenops (EndF, minimum frequency; FME, peak frequency; StartF, maximum frequency) in males recorded in Kilifi and Kwale counties (Triaenops afer) and Nakuru (Triaenops afer*). The x- and y-axes scales are equal to facilitate comparisons.

What to call these species is more complicated. We have little doubt about our application of the name T. afer, which Peters (1877) described from Mombasa, Kenya. Our robust reference samples from Kilifi and Kwale counties were taken within 100 km of Mombasa, both to the north and to the south of that city, and in the same coastal forest habitat. Good samples allow us to characterize typical T. afer closely across the various data partitions. In fact, the substitution network (Fig. 3) shows that the haplotypic diversity of single trident bat colonies in this region (e.g. those in Fikirini and Makuruhu Caves) encompass practically all of the known haplotypic variation of the species. Both caves contain individuals differing more from each other in Cytb than they do from T. afer in Ethiopia or central Tanzania. The morphology of typical T. afer also characterizes bats as far away as Ethiopia and Tanzania but differs somewhat from bats in the adjacent hinterland of Taita-Taveta County (Fig. 5). The vocalizations of T. afer also show slight differences between coastal and hinterland samples (female comparisons in Table 5). In view of this local variation, it is interesting that T. afer from Central Africa overlap broadly with Kenyan populations of T. afer but not with T. afer* (Supporting Information, Fig. S9). We therefore concur with Benda & Vallo (2009) that T. a. majusculus Aellen & Brosset, 1968 is a synonym of T. afer and conclude that it cannot be an available name for T. afer*.

We are equally confident concerning our ability to characterize at least the morphology and mitochondrial genetics of T. persicus. The samples used in our morphological analysis are near-topotypes from the Shiraz region of Fars Province, Iran. They agree well in morphology with broader Palaearctic sampling of Triaenops by DeBlase (1980) and Benda & Vallo (2009). Benda & Vallo (2009) conducted multivariate morphological analyses of Middle Eastern Triaenops that incorporated all relevant type specimens. They showed that: (1) there were no significant differences other than size between Iranian, Emirati, Omani or Yemeni populations and that all were referable as T. persicus; and (2) only two species of Triaenops occur along the southern margins of the Arabian Peninsula, namely T. persicus and their newly named T. parvus. Supporting their interpretation, our short Cytb sequences from typical Iranian T. persicus are identical at three (FMNH 96673) and nine (FMNH 96674) segregating sites to full sequences of Yemeni T. persicus (and to Kenyan T. afer*) and include no private alleles. Unfortunately, the vocalizations of Iranian Triaenops have not been published, but Benda et al. (2012) reported echolocation calls of T. persicus from Oman with peak frequencies between 76.5 and 82.6 kHz. Comparisons with the calls of male and female Triaenops from Kenya (Table 5) show that Arabian T. persicus uses nearly the same range of frequencies as do male and female T. afer. Although our sample lacks female calls, the range of peak frequencies used by male T. afer* falls below the ranges of all described species of Triaenops (Table 5; Benda et al., 2012; Ramasindrazana et al., 2013).

Thus, the Rift Valley population of T. afer* appears close to T. persicus genetically and is confidently recovered with it in a clade to the exclusion of all other Triaenops. The two lineages cannot have been separated for long. Nevertheless, T. afer* differs significantly from T. persicus in both morphology and vocalizations. What can we make of these differences? The Rift Valley populations of T. afer* in Kenya lie 1800 km from the nearest referred populations of T. persicus in Aden and 2500 km from Fars Province, Iran, where our morphological reference sample originated. Both morphology and vocalizations are subject to geographical variation and on a local scale. Within a genetically homogeneous sample of T. afer, we documented significant differences in morphology between Taita-Taveta and Kwale (only 200 km distant) and in female vocalizations between Kilifi and Laikipia (~500 km). What sort of differences should we expect within a genetically homogeneous T. persicus over the far vaster distances involved there? The gap analysis found evidence for species-level distinctions between T. afer and T. afer* (Supporting Information, Fig. S10) but failed to find such evidence between T. afer* and T. persicus (Supporting Information, Fig. S11). We conclude that the bats we designated T. afer* in the Rift Valley of Kenya must be recognized as African members of T. persicus.

Given the small genetic differences between T. persicus populations in Yemen and Kenya, it is logical to consider their spatial and temporal separation and whether these populations are really remote disjuncts. Triaenops has been recorded from several intervening localities in Ethiopia (Largen et al., 1974), including the Awash National Park (08°54′N, 39°55′E), which is situated in the Rift Valley and only 500 km from the Arabian Peninsula. Triaenops has also been recorded from Djibouti (Pearch et al., 2001), which is separated from Yemen by the narrow Strait of Bab al-Mandab, only 25 km wide and a natural corridor for African–Arabian colonization events. Interestingly, Djibouti Triaenops have the greatest skull lengths, which fall within the range of Kenyan T. persicus but outside (larger than) the range of T. afer (Pearch et al., 2001). It is possible that T. persicus extends its distribution through the Rift Valley to Djibouti and recently or regularly crossed the strait separating the Arabian Peninsula and the Horn of Africa. This raises other interesting questions on colonization events. Triaenops persicus and T. parvus are sister species (Fig. 2) and are likely to have diverged from T. afer (known only from Africa) in the Middle East. This suggests that African T. persicus populations represent back colonizations of Africa. An interesting follow-up to our work would be to address the dates and routes of these colonization events. Our integrative approach, combining molecular, morphometric and echolocation data to document the relationships of these trident bats, provides a firmer foundation for further biogeographical and taxonomic studies.

Implications for Paratriaenops

The strong genetic differentiation of Paratriaenops and Triaenops evident in our concatenated intron and species tree analyses (Fig. 4; Supporting Information, Fig. S2) reinforces earlier distinctions of these taxa based on mitochondrial evidence (Russell et al., 2007, 2008) and on morphology (Benda & Vallo, 2009). Foley et al. (2015) dated the divergence of these taxa at 22 Mya; that analysis recovered both Rhinonycteris and Cloeotis in successive splits off the lineage leading to Triaenops. The paraphyly of Malagasy rhinonycterids clearly supports the conclusion of Russell et al. (2008) that Madagascar was colonized at least twice in the history of this group. However, this interpretation hinges on the phylogenetic positions of Rhinonycteris and Cloeotis, which were not included in our analysis. A sister relationship of Paratriaenops + Triaenops, given the well-supported position of T. menamena as sister to African plus Arabian Triaenops, could indicate that only one colonization of Madagascar was involved. In this scenario, a descendent of T. menamena could have colonized the African mainland and given rise to the clade of T. afer, T. persicus and T. parvus. Additional genetic sampling of Cloeotis, Rhinonicteris and the missing Triaenops species might help to distinguish these alternatives, but it seems likely that extinction has strongly shaped the extant diversity of the group.

Genetic and distributional data provide mixed support for the validity of P. auritus and P. furcula as separate species. Previous morphological analyses of specimens assigned to these two taxa found consistent differences (Ranivo & Goodman, 2006). The Cytb genetic distance between these species (4.5%; Table 1) and the well-supported monophyly of P. auritus and moderately supported monophyly of P. furcula (Fig. 2) could be argued to support their current taxonomic status (also see Russell et al., 2008). In stark opposition, gene tree analyses of four independent nuclear loci under both ML and BI models did not recover any genetic structure within or between the two species (Fig. 4). Instead, the relationships inferred between P. auritus and P. furcula are consistent with extensive ongoing or recent hybridization. The allopatric distributions of these species do not provide any support for their reproductive isolation. Additionally, mitochondrial isolation by distance cannot be ruled out as the mechanism responsible for the genetic distance and topological relationship between populations assigned to P. auritus and P. furcula in our genetic analyses. Before the step is taken to combine these species, for which P. auritus would be the junior synonym, further genetic sampling is needed. In the material used in the present study, the northernmost locality for P. furcula (Namoroka, FMNH 175783) is ~330 km south of the southernmost locality for P. auritus (Betsiaka, FMNH 179370–179373); additional data are needed from the intervening zone to determine the nature and level of genetic separation between these forms.

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article at the publisher’s web-site:

Appendix S1. Sampling of Rhinonycteridae across various data partitions. All samples are documented by museum voucher specimens. Accession numbers identify sequences downloaded from GenBank or accessioned to GenBank for this study MT777711–MT777842. Institutional acronyms are as follows: DM, Durban Natural Science Museum, Durban; FMNH, Field Museum of Natural History, Chicago; MHNG, Natural History Museum of Geneva; MNHN, Muséum national de Histoire naturelle, Paris; NMK, National Museums of Kenya, Nairobi; NMP, Národní Museum, Prague; ROM, Royal Ontario Museum, Toronto; ZMB, Museum fur Naturkunde, Berlin. Columns headed by CD, Ext and Vocal correspond to samples used in analyses of craniodental morphology, external morphology and vocalization, respectively. Initials identify the person responsible for measurement of variables: BP, Bruce Patterson; JA, J. P. Adam; JM, Jessica Mohlman; JW, John Williams; PW, Paul Webala; SB, Stefania Briones; SE, Street Expedition; TA, T. Archer; VA, V. Aellen and A. Brosset.

Table S1. Primer information with respective name, primer melting temperature (Tm), sequence and publication.

Table S2. Prior schemes (PS) used in pairwise BPP analyses. Prior distributions on τ represent two relative divergence depths (deep and shallow) and on θ represent two relative mutation rate scaled effective population sizes (large and small).

Table S3. Craniodental and mandibular variables adapted from Velazco & Gardner (2012).

Table S4. Results of principal components analysis based on 12 log₁₀-transformed craniodental and mandibular variables. Coefficients of the first three principal components are tabulated; these account for 88% of overall variation. The amount of variation retained by each principal component (eigenvalues) and their respective proportion of variances are also tabulated.

Table S5. Cross-validation based on the linear discriminant analysis of craniodental and mandibular variables.

Table S6. Results from redundancy analysis after accounting for spatial linear trend for Triaenops afer and T. afer* from Nakuru–Baringo–Pokot counties. The second column provides the spatial eigenvectors (SEV) in descending order according to their first and second eigenvalues, followed by the dummy variable (Dum) and its interactions with spatial eigenvectors. Bold font indicates the rows corresponding to statistically significant regression coefficients (P < 0.05). The third to sixth columns provide the respective degrees of freedom (d.f.), variance, the F statistic and the significance of the respective regression coefficient (P).

Figure S1. Substitution network plot for Cytb inferred in POPART v.1.7 for Triaenops menamena.

Figure S2. Species tree of Triaenops and Paratriaenops estimated in StarBEAST2 using the four nuclear intron dataset. Numbers adjacent to nodes indicate posterior probabilities. Terminal tips in the tree that are statistically well supported [posterior probability (PP ≥ 0.95)] from BPP are indicated by ‘*’ preceding the clade name, and terminal tips that had a PP < 0.95 are indicated by ‘?’ preceding the clade name.

Figure S3. Principal components analysis (PCA) for craniodental and mandibular measurements showing the distribution and overlap of specimens without genetic information (‘unknown’; convex hulls in black) in the morphospace occupied by specimens containing cytochrome b gene sequences (represented in different colours).

Figure S4. Principal components analysis (PCA) performed for the whole sample after cross-validation tests.

Figure S5. Density plot of the 12 craniodental and mandibular variables (in millimetres) for female Triaenops afer, T. afer* and Triaenops persicus.

Figure S6. Box plot of the 12 craniodental and mandibular variables (in millimetres) for female Triaenops afer, T. afer* and Triaenops persicus.

Figure S7. Density plot of the 12 craniodental and mandibular variables (in millimetres) for male Triaenops afer, T. afer* and Triaenops persicus.

Figure S8. Box plot of the 12 craniodental and mandibular variables (in millimetres) for male Triaenops afer, T. afer* and Triaenops persicus.

Figure S9. Principal components analysis including the holotype and paratypes of Triaenops afer majusculus, based on nine craniodental and mandibular variables.

Figure S10. Inference of morphological gaps among Triaenops afer (green line) and T. afer* (blue line). The bimodal distribution on the left corresponds to the estimated probability function (pdf), fˆ(X), evaluated along the ridgeline manifold (α). The plot on the right corresponds to the estimated proportion, β, covered by the tolerance regions. The tolerance region overlaps above the frequency cut-off (dotted line).

Figure S11. Inference of morphological gaps among Triaenops persicus (red line) and Triaenops from Nakuru, Baringo and Pokot counties (blue line). The bimodal distribution on the left corresponds to the estimated probability function (pdf), fˆ(X), evaluated along the ridgeline manifold (α). The plot on the right corresponds to the estimated proportion, β, covered by the tolerance regions. The tolerance region overlaps below the frequency cut-off (dotted line).

Supplementary Material—External Morphology. Genetic, morphological and acoustic differentiation of African trident bats (Rhinonycteridae: Triaenops). This Supplementary Material file contains results for external morphology of Triaenops, including tables and figures.

ACKNOWLEDGEMENTS

D.M.R., T.C.D. and B.D.P. contributed equally to the manuscript. We thank the late Bill Stanley, Beryl Makori, Michael Bartonjo, Holly Lutz, Ruth Makena and David Wechuli for assistance in the field, and Stefania Briones for skeletal preparation and some of the cranial measurements. We thank Burton Lim and Jacqui Miller (ROM) for access to specimens in their care. Field collections in eastern and southern Africa were funded by a variety of agencies in cooperation with the Field Museum, especially the JRS Biodiversity Foundation. Field Museum’s Council on Africa, Marshall Field III Fund and Barbara E. Brown Fund for Mammal Research were crucial to fieldwork and analyses, as was the support of Bud and Onnolee Trapp and Walt and Ellen Newsom. For the collection of material on Madagascar, we acknowledge the Mention Zoologie et Biodiversité Animale, Université d’Antananarivo; Madagascar National Parks (MNP); Direction de la Gestion des Ressources Naturelles Renouvelables et des Ecosystèmes; and the Direction Générale de l’Environnement et des Forêts. We also thank Ara Monadjem for his insightful review of the previous versions of the manuscript. We thank the two anonymous reviewers for their comments and suggestions to improve the manuscript. We thank Jeronymo Dalapicolla and Bárbara Costa for productive discussions on analysis of gaps in morphological variation. We thank the Associate Editor and the Handling Editor for the editorial assistance on our manuscript.

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.