All samples were collected in CAR in 2008. Fruit bats were determined using the species identification key published by Bergmans [9]. A total of 23 individuals of E. gambianus from four localities and 97 individuals of M. pusillus from five geographic localities were included in this study. The localities are depicted in Fig. 2 and detailed in Table S1 (see Supplementary Material). In agreement with previous phylogenetic results [12], three outgroup species were used to root the tree: Epomops franqueti, which is another species of the tribe Epomophorini; Myonycteris torquata, which is a member of the tribe Myonycterini; and Rousettus aegyptiacus, which is the most distant outgroup, as it belongs to Rousettinae, a different subfamily from Epomophorinae.

Total genomic DNA was extracted from muscle or patagium samples using QIAGEN DNeasy Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol.

Three genes were sequenced for this study: two belonging to the mitochondrial genome (COI and Cytb) and one from the nuclear genome (a non-coding fragment: FGB). Primers used for polymerase chain reaction (PCR) were: UTyrLA and C1L705 for COI, UGluMA and LThrCH for Cytb, and 5′-CCA-CAA-CRG-CAT-GTT-CTT-CAG-CAC-3′ and 5′-GTA-TCT-GCC-ATT-TGG-ATT-GGC-TGC-3′ for FGB [13,14]. Amplifications were done in 20 μl using 2 μl of Buffer 10X with MgCl₂, 0.8 μl of dNTP (6.6 mM), 0.12 μl of Taq DNA polymerase (2.5 U, Qiagen, Hilden, Germany) and 0.32 μl of the two primers at 10 μM. The standard PCR conditions were as follows: 4 min at 94 °C; 35 cycles of denaturation/annealing/extension with 30 s at 94 °C for denaturation, 45 s at 50 °C for annealing and 90 s at 72 °C for extension; and 10 min at 72 °C.

PCR products were purified using ExoSAP Kit (GE Healthcare, Buckinghamshire, UK) and then sequenced in both directions using an automated DNA Sequencer (Applied Biosystems 3100). These two last steps were performed by GENOSCOPE (Ivry-sur-Seine, France). Nuclear fragments that contained multiple heterozygous single-nucleotide polymorphisms (characterized by double peaks in the chromatogram) were cloned using a QIAGEN PCR Cloning^plus Kit in order to isolate the exact sequence of each allele in heterozygous individuals. Sequences were edited and assembled using Sequencher 4.7 (Gene Codes Corporation). Sequences generated for this study were deposited in the EMBL/DDBJ/Genbank database (accession numbers JF728368-JF728760; Table S1 for details).

Sequences were aligned by eye using BioEdit version 7.0.9 [15]. No gaps were included in the alignments of COI and Cytb. For FGB, two unambiguous gaps were inferred to align the most divergent outgroup species, R. aegyptiacus, with other taxa.

Phylogenetic analyses were performed using Neighbor Joining (NJ) and Bayesian methods. The three genes (COI: 685 nt; Cytb: 1140 nt; FGB: 700 nt) were analyzed separately to evaluate their own signal and to detect any incongruence. The mitochondrial markers were then combined (COI + Cytb). MODELTEST version 3.7 [16] was run to determine the best-fitting evolutionary model. According to AIC criterion, the selected model was HKY + G for COI and Cytb, TVM + G for COI + Cytb and HKY for FGB. PAUP* 4b10 [17] was used for the NJ analyses with the selected model. Nodal support was assessed by the bootstrap analysis (1000 replicates). Bayesian inferences were conducted with MrBayes v3.1.2 [18]. Posterior probabilities (PP) were calculated using four independent Markov Chains run for 10,000,000 Metropolis-coupled MCMC generations with tree sampling every 1000 generations. TRACER [19] was used to calculate the effective sample size statistic to make sure that a sufficient number of generations were run and to decide on the length of the burn-in period. Analyses were run twice independently to check for convergence of the results.

Heterozygosity of nuclear genes complicates phylogenetic analyses because the two alleles from a heterozygous individual are not necessarily most closely related to each other. In practice, for many phylogenies based on diploid gene sequences, allelic sequences are not separated and heterozygous sites are treated as polymorphic. This procedure can make phylogenetic inference ambiguous, especially for closely related taxa [20]. However, the pairing of different alleles can provide information on interbreeding and may be useful in delimitating species. Therefore, the nuclear sequences of FGB were analyzed using the median-joining algorithm available in Network version 4.5.1 [21].

The pairwise genetic distances were calculated with PAUP* 4b10 [17] using the Kimura 2-parameter (K2P) correction. For each species, we calculated the number of polymorphic sites, number of haplotypes, nucleotide diversity (π) and haplotype diversity (h) [22,23] using DnaSP version 5.10.01 [24].

An “Isolation with Migration” (IM) model was used to test introgression rate between the two fruit bat species. We used a Bayesian MCMC method in the IMa version 2.0 coalescent program [25], which attempts to fit the data to a null model where an ancestral population bifurcates into two allopatric populations 1 and 2. This MCMC simulation estimates the joint posterior probability of six demographic parameters: population size of species 1 and 2 (N_e1, N_e2) and of the ancestral species (N_eA), migration (m₁, m₂) and divergence time (t). In our study, subscript 1 refers to the M. pusillus species and subscript 2 refers to the E. gambianus species. Although migration typically refers to movements of individuals between populations, the m parameter is better interpreted as the rate of mtDNA introgression for the purposes of this study. In this way, the m₁ parameter refers to introgression from E. gambianus into M. pusillus, whereas the m₂ parameter refers to introgression from M. pusillus into E. gambianus. As recommended in Hey & Nielsen [26], the inheritance scalars were set to 1 for the nuclear locus (FGB) and 0.25 for the mitochondrial sequences (COI + Cytb) and the HKY model was used. To assess convergence, we checked effective sample sizes throughout the run and compared results between three independent runs. Following the approach used by Won and Hey [27], we first ran IMa using large, flat priors for each parameter. Based on those runs, we defined narrower upper bounds that encompassed the full posterior distribution of each parameter. In the final runs, priors were set to 5 for t, 2 for m, 1000 for N_eA and 100 for N_e1 and N_e2. The input files were executed in IMa using 10 Metropolis-coupled chains with 100,000 steps for burn-in followed by 10 million steps for parameter estimation.

The 5′ region of COI (685 nt) and the complete Cytb gene (1140 nt) were sequenced for all samples available for E. gambianus (23 individuals) and M. pusillus (97 individuals). Mitochondrial genes showed high nucleotide and haplotype diversities for both species (Table 1), suggesting large effective size and stable populations (i.e., no recent bottleneck). Among COI sequences, a total of 60 haplotypes were identified: seven specific to E. gambianus, 49 specific to M. pusillus and four shared between the two species (Table S1). Among Cytb sequences, thirteen haplotypes were found for E. gambianus, 71 for M. pusillus and none were shared between the two taxa (Table S1). Using the alignment combining COI and Cytb sequences, we identified thirteen haplotypes for E. gambianus, 83 for M. pusillus and none were shared between the two taxa (Table S1).

The Bayesian tree reconstructed from the alignment combining COI and Cytb genes (123 taxa, 1825 nt) is shown in Fig. 3. All individuals of E. gambianus and M. pusillus fall into a robust clade (PP = 1), but neither species was found to be monophyletic. The trees obtained from separate analyses of COI and Cytb genes (Figs. S3 and S1) show very similar relationships, but the COI tree is less resolved than the Cytb tree, because COI has much less phylogenetic information content (165 and 89 parsimony informative sites for Cytb and COI, respectively). Moreover, the analysis including one complete Cytb sequence of Epomophorus wahlbergi (accession numbers DQ445706) revealed that the genus Epomophorus is paraphyletic, as E. wahlbergi occupies a more basal position with respect to the clade uniting E. gambianus and M. pusillus sequences (Fig. S2).

In Fig. 4, the comparisons between pairwise COI genetic distances show that DNA barcodes cannot be used for discriminating between E. gambianus and M. pusillus. Indeed, interspecific distances, which are comprised between 0% and 2.5%, are not higher than intraspecific distances found for E. gambianus (0–2%) and M. pusillus (0–2.5%). As a consequence, a query sequence cannot be correctly assigned to either E. gambianus or M. pusillus.

The nuclear FGB gene (700 nt) was also sequenced for all samples available for E. gambianus and M. Pusillus. The nuclear gene showed low nucleotide and haplotype diversities for both species (Table 1). We found three alleles for E. gambianus and fourteen alleles for M. pusillus (Table S1). None of them was shared between the two species. We identified six heterozygous individuals for E. gambianus (26%) and seventeen for M. pusillus (17.5%) (Table 1). In contrast with mtDNA data, the comparisons of FGB sequences indicate a higher haplotype diversity for E. gambianus and a higher nucleotide diversity for M. pusillus (Table 1).

The allelic network of FGB is presented in Fig. 5. The three alleles of E. gambianus are closely related and separated by only one mutation. The fourteen alleles of M. pusillus can be divided into two linked groups: the first one forms a star-like topology with a main central allele (Mp1), from which ten alleles are derived by one mutational step (Mp2 to Mp9; Mp12; Mp13) and two alleles by two mutational steps (Mp10 and Mp11); the second group corresponds to allele Mp14. Interestingly, the allele Mp14 is found in only five heterozygous individuals of M. pusillus and it is always associated with the most frequent allele, Mp1. In addition, Mp14 is six mutational steps distant from other alleles of M. pusillus, but it differs by only two mutations from Eg1, the most common allele of E. gambianus.

The introgression (or migration) parameter curves obtained from the isolation with migration program (IMa version 2.0) are presented in Fig. 6. In the case of introgression of E. gambianus into M. pusillus, the curve has a peak near zero and the probability of the peak is very near to what the probability is at zero. In this case, we cannot reject an introgression rate of zero. However, in the opposite direction, i.e., introgression from M. pusillus into E. gambianus, the peak is far from zero and the estimated probability of zero introgression is zero; we can therefore reject a model of no introgression from M. pusillus into E. gambianus.

According to Bergmans [9], Epomophorus and Micropteropus can be ranged in the tribe Epomophorini with the three other genera Epomops, Hypsignathus and Nanonycteris. By analyzing mtDNA sequences of Cytb and 16S rRNA genes, Juste et al. [12] have concluded on a sister-group relationship between Epomophorus wahlbergi and Micropteropus pusillus, whereas Epomops franqueti was found to be more divergent. Other species of Epomophorini were however not included in the molecular analyses. Although our mtDNA results confirm that Epomophorus and Micropteropus are closely-related genera, it was highly unexpected to find a paraphyletic pattern for Epomophorus (Fig. S2) and a polyphyletic pattern for both E. gambianus and M. pusillus species (Fig. 3). Indeed, M. pusillus can be easily distinguished from E. gambianus and E. wahlbergi on the basis of morphological characteristics, such as body size (Fig. 1; weight: 20–35 g versus 56–155 g for E. gambianus and 54–125 g for E. wahlbergi; forearm length: 48–56 mm versus 75–100 mm for E. gambianus and 68–95 mm for E. wahlbergi) and the pattern of palatal ridges (Fig. 1; only five bridges for E. wahlbergi, see in Bergmans [9]). Since our analyses of the nuclear gene FGB agree with morphology, as both M. pusillus and E. gambianus possess distinct alleles (Fig. 5), we can conclude that mtDNA data are problematic. Three hypotheses can be advanced to explain the incongruence of mtDNA with nuclear and morphological data: the amplification of Numts, mtDNA introgression and incomplete lineage sorting of mtDNA haplotypes.

Nuclear mitochondrial pseudogenes or ‘Numts’ are segments of mtDNA translocated to the nuclear genome. These paralogous sequences are commonly found in animal genomes and can be accidentally amplified instead of the targeted mitochondrial genes [28]. The undetected presence of Numts in the analyses can result in erroneous interpretations of phylogenetic relationships, because Numts evolve under constraints that are different than those of mtDNA [29–31]. The nuclear and mtDNA genomes have different rates and patterns of mutations and, because Numts are not functional, they do not evolve under purifying selection, explaining why they readily accumulate stop codon and frameshift mutations [32].

Four major arguments suggest that the COI and Cytb sequences produced in the present study cannot be considered as being Numts: (1) as the primers used here for PCR and sequencing were specially designed to amplify mtDNA genes, all our chromatograms were perfectly readable, without double peaks; (2) there is a high identity (99.8%) between our Cytb sequences of M. pusillus and the partial sequence previously obtained by Juste et al. [12] using a different set of primers; (3) there is no stop codon or indel (insertion or deletion); and (4) the tree reconstructed from COI sequences is very similar to that obtained from Cytb sequences (Figs. S3 and S1), with several identical nodes robustly supported in both analyses (Fig. 3). Therefore, there is no doubt about the mitochondrial origin of our COI and Cytb sequences.

The occasional mating between distinct species can result in genetic introgression of the mtDNA genome or nuclear alleles of one species into the gene pool of another species. This is a cause of gene/species tree discordance. Mitochondrial introgression following hybridization has been widely inferred, especially when dealing with differences in genetic signal between mtDNA and nuclear markers [13,33–36], but documented cases in bats are very rare [37]. The mtDNA from one species may completely replace that of another species, without leaving any trace of nuclear introgression or morphological change [38]. In such cases, the mtDNA tree generally shows a typical topology, in which one of the two species is found to be monophyletic and the other paraphyletic. In our case, the mtDNA tree shows a more complex pattern, in which both species are polyphyletic. Assuming the hypothesis of mtDNA introgression implies therefore that interspecific breeding occurred frequently and at different periods of time. In addition, IMa analyses suggest that mtDNA introgression events occurred from M. pusillus into E. gambianus. According to this scenario, females of M. pusillus have been integrated several times into E. gambianus populations and all thirteen mtDNA haplotypes detected here for E. gambianus were transferred from M. pusillus.

Three main arguments concerning DNA variation, ecology and sexual behaviour render this scenario plausible.

Firstly, mtDNA variations between E. gambianus and M. pusillus are between 0 and 2.5% for COI and between 0 and 2.4% for Cytb. These low levels of variation are comparable with intraspecific divergence calculated for COI and Cytb sequences in several species of the family Pteropodidae (Tables S2 and S3). By contrast, the complete Cytb sequence available for E. wahlbergi is divergent from those of E. gambianus and M. pusillus (4.8–5.7%). These comparisons support therefore the fact that all mtDNA sequences produced here for E. gambianus and M. pusillus originated from only one species, namely M. pusillus.

Secondly, both species are found in sympatry from Senegal to CAR (Fig. 2), where they occupy the same type of habitat, i.e. the forest-savannah mosaic [9,39], which is a transition zone, or ecotone, between the tropical moist broadleaf forests and the drier savannahs. As a consequence, they were often collected in the same mist nets (personal observations) and were observed roosting in the same tree [40]. The geographic distribution of M. pusillus (Fig. 2) suggests, however, that it is more dependent on rainforests than E. gambianus. Hybrid zones often coincide with ecotones, as previously shown in African elephants [33] or mouse lemurs [41]. Astonishingly, the situation in African fruit bats is comparable with that of African elephants, as multiple events of mtDNA introgression have been inferred from the smallest forest species (M. pusillus/L. cyclotis) into the largest savannah species (E. gambianus/L. africana).

Thirdly, introgressive hybridization is known to occur more frequently between sister-species [42], for which prezygotic isolation (behavioural and mechanical) and postzygotic isolation (zygote mortality and hybrid inviability and sterility) have not been fully accomplished. Here, both species share very similar reproductive patterns. Females have two parturition periods during rainy seasons around April and September/October in West Africa. Furthermore, adult males exhibit similar sexual characteristics (white tufts of erectile hair grow within circular pocket-like folds of skin on the shoulders) and have the same typical courtship behaviour (with courting calls, the display of the white shoulder patches and wing flapping) [40,43].

However, three arguments dispute an introgression from M. pusillus into E. gambianus. Firstly, the male courting call is at a different frequency: 1750 and 2800 Hz for Epomophorus and M. pusillus, respectively [44]. Secondly, the two species have different male karyotypes: M. pusillus is XY₁Y₂, with 2n = 35 and a fundamental number (FN) of 64 [45], whereas E. gambianus is XO, with 2n = 35 and a fundamental number (FN) of 68 [46]. The karyotype available for E. gambianus was established from two specimens from Kenya and Zimbabwe, i.e. from two populations of E. gambianus crypturus. Since this subspecies has been elevated to full species status in recent classifications, i.e. Epomophorus crypturus [1,2], we can expect a different karyotype for specimens from CAR. Thirdly, and probably most importantly, their body size is very different, as adult males of E. gambianus are four times heavier than females of M. pusillus (Fig. 1; [39,47]). Such morphological differences may preclude copulation and cause mortality in pregnant M. pusillus females carrying an overly large hybrid foetus.

As an alternative to the hypothesis of mtDNA introgression, the phenomenon of incomplete lineage sorting, due to the retention of ancient polymorphisms, can be invoked to explain the polyphyletic pattern observed here for mtDNA haplotypes of E. gambianus and M. pusillus. If the mtDNA genome has not reached coalescence at the species level, some lineages from one species may be more closely related to a lineage in the other species than to other lineages within the same species. Incomplete lineage sorting can take place between sister-species, when the recent speciation event occurs before sorting is completed [48,49]. A recent speciation event between E. gambianus and M. pusillus is here supported by the fact that both mtDNA and nuDNA markers show very low interspecies distances (< 2.5% for mtDNA and < 1% for nuDNA).

The shorter the coalescence time, the less likely the gene will suffer from the ‘incomplete lineage sorting’ problem. In this regard, mtDNA has an advantage, because it coalesces in one-fourth the time of nuclear genes [48,49]. Therefore, if incomplete lineage sorting of mtDNA haplotypes is possibly the cause of the polyphyly, we should assume a faster time of coalescence for the FGB nuclear gene than for the mitochondrial genome.

A strong positive selection in the FGB gene may have led to selecting one or several particular allelic forms rapidly in one or both species before the coalescence of mtDNA. Actually, positive selection can be advanced for M. pusillus, because most FGB alleles detected in this species are highly derived with respect to those found in E. gambianus. Indeed, Mp1 to Mp13 alleles differ by five mutations from the ancestral allelic form shared by E. gambianus and M. pusillus (Fig. 5). In accordance with the hypothesis of positive selection, the primitive allele Mp14 is very rare in M. pusillus (2.6%) and it is always found in heterozygosity with a derived Mp1 allele, suggesting that at least one derived allele is necessary to ensure the function of FGB in M. pusillus.

The FGB gene codes for the beta polypeptide chains of the fibrinogen, a glycoprotein involved in blood coagulation. In humans, several polymorphisms in the FGB gene were found to affect the concentration of clotting factor and therefore the plasma viscosity [50,51] or to protect against myocardial infarction [52]. In addition, previous physiological studies on fruit bats have concluded that there is a linear correlation between body mass and heart pulse rate [53]. In other words, the smallest bats, which have the highest wing beat frequency, have the highest heart rate. Therefore, we can expect to have a higher selective pressure for FGB alleles in M. pusillus than in E. gambianus.