Genetic Introgression and Morphological Variation in Naked-Back Bats (Chiroptera: Mormoopidae: Pteronotus Species) along Their Contact Zone in Central America

: Two sibling bare-backed bat species ( Pteronotus fulvus and P. gymnonotus ) have been traditionally differentiated by their size. However, intermediate specimens between the two species have been found in sympatric populations along southern Mexico and it has been suggested that they may be the outcome of a hybridization process between the two species. We used one mitochondrial (COI), three nuclear markers (PRKCL, STAT5A and RAG2) and 13 microsatellites to explore the evolutionary relationships between these two species and elucidate whether the intermediate morphotypes correspond to hybrid individuals. These markers have been analyzed in sympatric and allopatric populations of the two species plus the closely related species Pteronotus davyi . We conﬁrmed the species-level differentiation of the three lineages ( P. fulvus , P. davyi and P. gymnonotus ), but the phylogenetic hypotheses suggested by the nuclear and mitochondrial markers were discordant. We conﬁrm that the discordance between markers is due to genetic introgression through the mitochondrial capture of P. fulvus in P. gymnonotus populations. Such introgression was found in all P. gymnonotus specimens across its sympatric distribution range (Mexico to Costa Rica) and is related to expansion/retraction species distribution pulses associated with changes in forest distribution during the Quaternary climate cycles. Microsatellite analyses showed contemporary genetic contact between the two sympatric species and 3.0% of the samples studied were identiﬁed as hybrids. In conclusion, we found a historical and asymmetric genetic introgression (through mitochondrial capture) of P. fulvus into P. gymnonotus in Mexico and Central America and a limited contemporary gene exchange between the two species. However, no relationship was found between hybridization and the intermediate-sized specimens from southern Mexico, which might likely result from a clinal variation with latitude. These results conﬁrm the need for caution when using forearm size to identify these species in the ﬁeld and when differentiating them in the laboratory based on mitochondrial DNA alone.


Introduction
Molecular methods have contributed to understanding evolutionary patterns and obtaining more accurate inferences about the relationships between lineages, bringing in major changes in the recognition and delimitation of biological species. These methods have been particularly useful in groups in which criteria used for species delimitation do not correspond to clearly discernable morphological differences, which have traditionally been the main tools for recognizing and setting the boundaries between species [1,2].
Elucidating evolutionary history is basic for establishing taxonomic boundaries between taxa, understanding speciation processes and estimating the actual diversity of evolutionary lineages to formulate suitable conservation strategies [3]. Comparing phylogenetic patterns across lineages through comparative studies of genetic diversity, gene flow, phylogenetic relationships, diversification events and others [4] helps to elucidate whether taxa share a common history.
The use of molecular tools has led to an increased number of known species but also to clump species unto just one [3,5,6]. However, analyses of different molecular markers frequently lead to different taxonomic units or contrasting genealogical relationships between the taxa studied [7,8]. Such discordance between markers may be due either to an incomplete lineage sorting of ancestral polymorphism or to genetic introgression, among other processes [9,10], which pose additional challenges for phylogenetic studies [11].
Genetic introgression is an exchange of genetic material (hybridization) in which part of the gene pool of one lineage is mixed and becomes fixed in another evolutionary lineage [10]. This process is particularly interesting in systematic studies, given its confounding effect in phylogenetic inferences and its potential role in speciation processes, particularly through adaptive radiation [3]. For a genetic introgression to become fixed in a population, it is necessary that the hybrid offspring backcross with one of its parent species, thus leading to the permanent incorporation of the new genes from another lineage [10]. When genetic material is shared, the direction and degree of hybridization, as well as its evolutionary consequences, generally depend on factors such as the degree of overlap of the geographic ranges of the parent species, compatibility of genomes, breeding strategies, dispersal, social structure patterns and selective pressures, among others [12,13].
Hybrids found at the limits of the distribution ranges of parent species are particularly interesting for the interaction of intrinsic and extrinsic elements that determine the temporal and geographic distribution of the species [14]. Contact zones where morphological or genetic characteristics vary gradually are called clinal zones [15]. A higher proportion of hybrid individuals can be found in clinal zones mainly because they have a selective advantage over the parent species [15,16]. Clines pose additional challenges to the delimitation of species, particularly in species that show phenotypic plasticity and in which their differentiation is based upon morphological diagnosis [17].
Bats are an interesting model for studying introgression processes due to their high vagility, wide-ranging distribution and phenotypic plasticity. Species delimitation is frequently challenging: while some species can be readily recognized based on morphometric measurements [18,19], many others have been delimited only through molecular techniques due to the morphological conservatism found in species complexes occurring in several families such as Hipposideridae, Rhinolophidae, Miniopteridae, Vespertilionidae and others. [20][21][22] In addition, the use of more integrative approaches like IOTUs have recently been suggested [23].
The family Mormoopidae includes insectivorous bats of Neotropical distribution ranging from southern Texas to Brazil, including the Antilles [24,25]. The family has a complex taxonomic history that has led to repeated re-descriptions of taxa and taxonomic rearrangements [24][25][26]. The distributional boundaries between some species are still unclear, particularly for mainland species. The family had been traditionally considered to include two genera: Mormoops Leach, 1821 [27] with three species and Pteronotus Gray 1838 [28] with seven species subdivided into three subgenera [24,29]. Recently, however, a new taxonomic arrangement was proposed for Pteronotus [30][31][32][33], which elevates most of the subspecies to the species level, bringing the total number of species to 16 subdivided into the subgenera Pteronotus, Chilonycteris and Phyllodia [31]. The subgenus Pteronotus is characterized by wing membranes fused on the back, giving them a distinctive bare-back appearance. The main morphological difference between species refers to size, particularly forearm length (FA) [24]. Pavan and Marroig [31] recognized one large species, P. gymnonotus (FA ≥ 49 mm) and two visibly smaller species, P. fulvus and P. davyi (FA < 49 mm). The smaller species overlap in forearm length, but various studies have confirmed that, despite being considered as co-specific until recently, they do correspond to clearly differentiated evolutionary lineages at the species level [24,34,35]. The three species show very close phylogenetic relationships, having differentiated 2.5 million years (My) ago [31].
There is a sympatric zone between P. fulvus and P. gymnonotus in southeast Mexico, which also coincides with the distribution limit of P. gymnonotus (Figure 1b,d). Individuals with forearm length (FA = 46.1, 46.5 and 48.3 mm) close to both the upper limit of P. fulvus and the lower limit of P. gymnonotus [36,37] have been recorded in this area. These intermediate-sized individuals have been proposed as another P. davyi subspecies (cf. P. fulvus) [36,37]. The size variation observed in the subgenus Pteronotus has also been explained as part of a latitudinal gradient [24]. Pavan and Marroig [31] suggested a possible introgression of mitochondrial DNA (mtDNA) from P. fulvus into P. gymnonotus. Although these lineages represent phylogenetically close species [31,35] and share refuges (caves) [38], the occurrence of (current or past) hybridization between them in the sympatric zone has not been directly evaluated yet. Based on an intensive sampling conducted in allopatric and sympatric localities of P. fulvus and P. gymnonotus, this study pursued the following objectives: (1) evaluate whether forearm length and latitude are correlated; (2) discern the taxonomic composition of the subgenus Pteronotus; (3) determine whether hybridization between P. fulvus and P. gymnonotus has occurred and (4) confirm whether the intermediate-sized morphotypes correspond to hybrids between P. fulvus and P. gymnonotus.

Materials and Methods
We captured a total of 247 bats (187 P. fulvus and 60 P. gymnonotus) at 14 different localities in Mexico, Nicaragua, Costa Rica and Brazil. The specimens were captured between 2004 and 2019 using harp-traps and mist nets ( Figure 1, Table 1). The bats were morphologically identified using the taxonomic keys by Medellín et al. [18] and Adams [39], using the forearm length as diagnostic characteristic and following the taxonomy proposed by Pavan and Marroig [31]. For each, specimens were measured the forearm, total length, weight and sexed recorded and a wing membrane biopsy of two mm diameter was collected with a biopsy punch (Fray Products Corp., Buffalo, New York, NY, USA). Most bats were released at the capture site, except for a few individuals that were kept as voucher specimens and deposited at the Universidad Autónoma Metropolitana, Iztapalapa Campus, under catalog entries RLW130901Pda11-RLW23072018Pgy2. Tissue samples were preserved in 70% ethanol. Ethical protocols set by the American Society of Mammalogists [40] and the División de Ciencias Biológicas y de la Salud (Division of Biological and Health Sciences) of Universidad Autónoma Metropolitana-Iztapalapa [41] were followed during the capture and handling of bats. Bat capture was carried out under the collection permits (see below) In addition, our own sampling, the genetic study, includes all the sequences available in GenBank (Table S1).

Correlation between Forearm Length and Latitude
The correlation between forearm length and latitude was evaluated separately for P. fulvus and P. gymnonotus using the Spearman's rank correlation coefficient (rs) and the data were tested for normality using the D'Agostino Omnibus test NCSS v. 11 [42]. This analysis included data from the specimens captured and additional data reported in the literature [43][44][45][46][47][48][49][50][51]. The final dataset included a total of 201 female and 387 male P. fulvus and 39 female and 50 male P. gymnonotus, from localities in Mexico, Guatemala, Nicaragua, Costa Rica, Panama, Venezuela and Brazil.

DNA Extraction, Sequencing, Editing and Alignment
Genomic DNA was extracted using the standard saline method [56]. Three nuclear DNA (nDNA) and one mitochondrial (mtDNA) gene fragments were amplified with Polymerase Chain Reaction (PCR). The nDNA fragments were (a) 422 bp of the protein kinase C iota intron (PRKC1), (b) 475 bp of the signal transducer and activator of transcription 5A intron (STAT5A) and (c) 717 bp of the recombination activating gene 2 (RAG2). For the mtDNA, a 607 bp fragment of the gene Cytochrome Oxidase subunit I (COI) was amplified. The primers and conditions described by [57][58][59] were used for each of these fragments (Table S2). The amplification products were sequenced on an ABI PRISM 370X equipment and the sequences were manually edited and aligned using the Clustal W algorithm implemented in Geneious v. 5.6.4 [60,61] and uploaded to GenBank (Table S1).

Species Tree
Seventy-four specimens (26 P. fulvus, 43 P. gymnonotus and 5 P. davyi) from localities representative of the distribution range of each species were selected to confirm and delimit taxonomic groups in the lineages studied of the subgenus Pteronotus. A species tree was constructed from the nDNA (RAG2, PRKC1 and STAT5A) and mtDNA (COI) data using the software *BEAST 2.4.8 [64] as available in the CIPRES webpage [65]; the resulting tree was visualized using the software FigTree 1.1.2 [66]. Three runs of 30,000,000 generations each were performed with the Bayesian and Maximum Likelihood phylogenetic inferences as initial hypotheses. A Yule speciation model was used with a 0.8 exponential mean (yule.birthRate), uncorrelated lognormal distribution, strict molecular clock and a substitution rate for each partition (ucld.mean) of 0.01 substitutions per site per million years (s/s/my) for PRKC1 and STAT5A, 0.005 s/s/my for RAG2 and 0.025 s/s/my for COI, as recommended by Pavan and Marroig [67]. Sequences from the species P. macleayii and P. quadridens (Table S1) were used as external groups.

Evolutionary Relationships
The same dataset used for constructing the species tree (see above) was used to reconstruct and compare the evolutionary history of P. fulvus, P. gymnonotus and P. davyi. nDNA and mtDNA phylogenies were constructed with Bayesian inference (BI) running for 10,000,000 generations with five Markov chains, using the software MrBayes v. 3.2 [68] and a Maximum likelihood (ML) analysis by means of a heuristic search with 1 000 bootstrap replicates with the software PAUP* v. 4.0b10 [69]. The evolutionary model most consistent with each locus was estimated with the software jModeltest v. 2.1.6, using the Akaike Information Criterion [70,71]. The models selected were HKY for PRKC1, GTR+G for STAT5A, TPM2uf+I for RAG2 and HKY+G for COI.
Divergence times of the main nodes of the mtDNA phylogeny were estimated with the software *BEAST 1.8.0 [72] as available in the CIPRES webpage [65]. Three runs of 30,000,000 generations each were performed with a Yule tree prior (yule.birthRate = 0.8), uncorrelated lognormal distribution, strict molecular clock and the same mutation rate used for the species tree; the results were visualized with the software FigTree 1.1.2 [66]. Haplotype networks for the different markers (PRCK1, STAT5A, RAG2 and COI) were built using the software Network v. 4.6.1.3 [73] to infer the relationships between individuals; loops were resolved as per the criteria recommended by Pfenninger and Posada [74]. Genetic distances between haplogroups were estimated with the software MEGA v. 5.0.5 [75] using the p-distance model for nDNA and the Kimura 2-parameter model (K2P) for mtDNA.

mtDNA Diversity and Structure
Genetic structure and variation in P. fulvus and P. gymnonotus were evaluated with the COI mtDNA marker. The sample size was increased to a total of 300 individuals (192 P. fulvus and 108 P. gymnonotus) from sympatric populations of Mexico, Guatemala, El Salvador, Honduras, Belize, Nicaragua and Costa Rica) and from allopatric populations from Mexico, Panama, Venezuela, Guyana, Suriname, Peru and Brazil (Table 1 and Table S1). For these populations, genetic variation and structure for the two species were estimated by means of pairwise F ST values between species and localities, using the software Arlequin v. 3.5.1.2; significance was tested using 10,000 permutations [76]. The spatial distribution of haplotypes per locality and species was visualized on maps. Relationships between haplotypes were determined through parsimony network analyses using the Median-Joining method with the software Networks v. 4.6.1.3 [73]. Genetic distances (K2P) between haplogroups were calculated using the software MEGA v. 5.0.5 [75].

Microsatellinte Data Analyses
A total of 199 individuals (155 P. fulvus from six localities in Mexico and 44 P. gymnonotus from two localities in Mexico) were genotyped with the 13 microsatellites selected. The presence and frequency of null alleles on each locus were estimated with the soft-ware MICROCHECKER v. 2.2 [77] to determine the potential impact of null alleles on test results. F ST values and genetic distances with and without ENA (estimation of null alleles) correction were estimated using the software FREENA [78]. Subsequently, the F ST and genetic distance matrices were compared with a Mann-Whitney U test to test for significant differences between them, using the software NCSS v. 11 [42]. Deviation from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) tested with the software Arlequín v. 3.5.1.2 [76] and applying the Bonferroni's sequential correction (p < 0.05) [79].

nDNA Genetic Population Structure
The genetic diversity of each species at each locality was evaluated by estimating the number of alleles (Na), number of private alleles (Np) and observed (H O ) and expected heterozygosity (H E ) with the software GenAlEx v. 6.3 [80] and allelic richness (AR) with the software FSTAT [81]. Genetic variation and differentiation were evaluated through F ST and R ST (gradual mutation model, SMM) values between species and pairwise R ST values between species and localities, using the software Arlequin v. 3.5.1.2; significance was tested after 10,000 permutations [76].
Genetic structure was evaluated and potential hybrids at nuclear level were identified with the software STRUCTURE v. 2.2 [82]; were performed 10 7 MCMC iterations, 500,000 as burn-in with ten replicates for each genetic cluster (K), from K = 2 through K = 9. The most probable K value was determined based on the value of ∆K [83] estimated with the software Structure Harvester [84] and analyzing together the results from the ten replicates of each K value with the software CLUMPP v. 1.1.2 [85]. The software NEWHYBRIDS v.1.0 [86] was also used to detect hybrids, with a burn-in of 5 × 10 4 iterations followed by 1,000,000 five-chain sweeps.
Hybrids identification was based on assignment methods using the software STRUC-TURE [82] and NEWHYBRIDS [86]. According to STRUCTURE, were considered as hybrids those individuals that were not showing an assignment value of q < 0.90 to any of the two Pteronotus species. They were identified as hybrids in NEWHYBRIDS [86], to those individuals showing q > 0.75, according to the criterion number two of Burgarella et al. [87]; inconsistencies between the two analyses were resolved in favor of STRUCTURE, as recommended by these authors.

Morphological Analysis
The putative hybrid specimens accounted for 9.5% of the bats captured in the localities LO, AB, Ta and LV, from the Mexican states of Colima, Tabasco and Oaxaca ( Figure 2).

Correlation between Forearm Length and Latitude
The forearm data are not normal and the showed significant differences between sexes. No significant correlation was found between latitude and forearm length in P. fulvus females (rs = 0.1321, p = 0.8950). However, a significant negative correlation was observed in P. fulvus males (rs = 0.3271, p = 0.0205). In contrast, a positive correlation was found in P. gymnonotus females and males (rs = 2.8038, p = 0.0073; rs = 4.0997, p = 0.002, respectively) ( Figure S1).

Molecular Analysis
The species tree supported three lineages: P. davyi, P. fulvus and P. gymnonotus and suggesting with moderate support that P. fulvus is more closely related to P. gymnonotus than to its sibling species P. davyi ( Figure 3).

Correlation between Forearm Length and Latitude
The forearm data are not normal and the showed significant differences between sexes. No significant correlation was found between latitude and forearm length in P. fulvus females (r s = 0.1321, p = 0.8950). However, a significant negative correlation was observed in P. fulvus males (r s = 0.3271, p = 0.0205). In contrast, a positive correlation was found in P. gymnonotus females and males (r s = 2.8038, p = 0.0073; r s = 4.0997, p = 0.002, respectively) ( Figure S1).

Species Tree
The species tree supported three lineages: P. davyi, P. fulvus and P. gymnonotus and suggesting with moderate support that P. fulvus is more closely related to P. gymnonotus than to its sibling species P. davyi ( Figure 3).

Evolutionary Relationships
The 74 sequences examined-corresponding to P. fulvus, P. davyi and P. gymnonotus-contained 18 haplotypes in the joint sequences of the PRKC1 and STAT5A introns (897 bp), nine haplotypes in RAG2 (717 bp) and 32 haplotypes in the 607 bp fragment of COI (Table S4).
The Bayesian inference and maximum likelihood criteria yielded similar phylogenies for both marker types. However, a marked discordance was observed between the hypotheses derived from each marker type (nuclear and mitochondrial). Nuclear topologies comprised two monophyletic clades, one including the species P. fulvus and P. davyi (Group I) and the other including P. gymnonotus (Group II), separated by a genetic distance (p-distance) of 0.8%. The haplotype network for the PRKC1 and STAT5A introns clearly showed unique haplotypes to each species, unlike the RAG2-based haplotype network, where all three species shared some haplotypes (Figure 4b).
On the other hand, the phylogenies based on the mitochondrial marker COI revealed two clusters (Figure 4c). A first cluster (Group I) included mainly P. fulvus haplotypes (n = 9), as well as haplotypes of P. gymnonotus (n = 6) from its northern distribution range (Mexico, Guatemala, Nicaragua and Costa Rica) and corresponding to the sympatric zone with P. fulvus. The second cluster (Group II) comprised two haplogroups: the first included P. davyi haplotypes (n = 4), one P. fulvus haplotype from Mexico (H1) and one haplotype (H2) shared by P. fulvus and P. gymnonotus from Mexico and Costa Rica. The second haplogroup-although weakly supported (Pp = 73/76)-included all P. gymnonotus haplotypes (n = 10) from the southern distribution range (Panama, Suriname, Peru and Brazil) of the species (Figure 4). The genetic distance between the two haplogroups was 3.7% (K2P). The inter-specific genetic distance was 6.8% between P. fulvus and P. davyi, 3.7% between P. fulvus and P. gymnonotus and 7.5% between P. davyi and P. gymnonotus (K2P).

Evolutionary Relationships
The 74 sequences examined-corresponding to P. fulvus, P. davyi and P. gymnonotuscontained 18 haplotypes in the joint sequences of the PRKC1 and STAT5A introns (897 bp), nine haplotypes in RAG2 (717 bp) and 32 haplotypes in the 607 bp fragment of COI (Table S4).
The Bayesian inference and maximum likelihood criteria yielded similar phylogenies for both marker types. However, a marked discordance was observed between the hypotheses derived from each marker type (nuclear and mitochondrial). Nuclear topologies comprised two monophyletic clades, one including the species P. fulvus and P. davyi (Group I) and the other including P. gymnonotus (Group II), separated by a genetic distance (p-distance) of 0.8%. The haplotype network for the PRKC1 and STAT5A introns clearly showed unique haplotypes to each species, unlike the RAG2-based haplotype network, where all three species shared some haplotypes (Figure 4b).
On the other hand, the phylogenies based on the mitochondrial marker COI revealed two clusters (Figure 4c). A first cluster (Group I) included mainly P. fulvus haplotypes (n = 9), as well as haplotypes of P. gymnonotus (n = 6) from its northern distribution range (Mexico, Guatemala, Nicaragua and Costa Rica) and corresponding to the sympatric zone with P. fulvus. The second cluster (Group II) comprised two haplogroups: the first included P. davyi haplotypes (n = 4), one P. fulvus haplotype from Mexico (H1) and one haplotype (H2) shared by P. fulvus and P. gymnonotus from Mexico and Costa Rica. The second haplogroup-although weakly supported (Pp = 73/76)-included all P. gymnonotus haplotypes (n = 10) from the southern distribution range (Panama, Suriname, Peru and Brazil) of the species (Figure 4). The genetic distance between the two haplogroups was 3.7% (K2P). The inter-specific genetic distance was 6.8% between P. fulvus and P. davyi, 3.7% between P. fulvus and P. gymnonotus and 7.5% between P. davyi and P. gymnonotus (K2P). P. gymnonotus was the only species showing a well-marked intra-specific structure with two lineages, one corresponding to the northern end (Mexico-Nicaragua) and the other to the southern end (Panama-Brazil) of its distribution range. The divergence time estimated between the southern lineage of P. gymnonotus and the other groups (P. davyi, P. fulvus and P. gymnonotus North) was 1.021 My (0.705-1.506); for P. fulvus/P. gymnonotus North vs. P. davyi, it was 0.833 My (0.565-1.177); and for P. fulvus vs. P. gymnonotus North, it was 0.243 My (0.127-0.458) ( Figure S2).

mtDNA Diversity and Structure
For the 133 specimens morphologically identified as P. fulvus (69.3%) and 39 identified as P. gymnonotus (36.1%) a total of 63 haplotypes were found in the two species; 34 haplotypes were unique to P. fulvus, 25 to P. gymnonotus and four haplotypes (H1, H2, H7 and H11) were shared. (Table S5). P. fulvus showed several haplotypes shared localities across the Gulf of Mexico and the Pacific regions. The largest number of shared haplotypes was found in localities of southern Mexico, while specimens from Central America showed unique haplotypes ( Figure S3a). Haplotypes of the species P. gymnonotus grouped in two lineages with a clear geographic pattern, being the haplotype H2 recorded in almost all the localities of southern Mexico and Central America, while the haplotype H47 was found in most of the South American localities ( Figure S3b).
The mtDNA haplotype network showed three clearly differentiated haplogroups with ≥30 mutational steps and genetic distances ≥6% ( Figure 5). Haplogroup I included 43 haplotypes from 183 P. fulvus specimens from Mexico-Costa Rica and 47 P. gymnonotus specimens from Mexico, Guatemala, Nicaragua and Costa Rica (sympatric distribution) and one specimen from Panama. This suggesting a mitochondrial capture of all the P. gymnonotus from the northern distribution range. The first haplogroup also shows two main haplotypes in a star-like structure, suggesting population expansion. Haplogroup II included 17 haplotypes, all from the 52 P. gymnonotus specimens from the southern distribution range with samples from Panama-Brazil (allopatric distribution) and showing expansion. Finally, the haplogroup III included three haplotypes, one from a P. fulvus from The divergence time estimated between the southern lineage of P. gymnonotus and the other groups (P. davyi, P. fulvus and P. gymnonotus North) was 1.021 My (0.705-1.506); for P. fulvus/P. gymnonotus North vs. P. davyi, it was 0.833 My (0.565-1.177); and for P. fulvus vs. P. gymnonotus North, it was 0.243 My (0.127-0.458) ( Figure S2).

mtDNA Diversity and Structure
For the 133 specimens morphologically identified as P. fulvus (69.3%) and 39 identified as P. gymnonotus (36.1%) a total of 63 haplotypes were found in the two species; 34 haplotypes were unique to P. fulvus, 25 to P. gymnonotus and four haplotypes (H1, H2, H7 and H11) were shared. (Table S5). P. fulvus showed several haplotypes shared localities across the Gulf of Mexico and the Pacific regions. The largest number of shared haplotypes was found in localities of southern Mexico, while specimens from Central America showed unique haplotypes ( Figure S3a). Haplotypes of the species P. gymnonotus grouped in two lineages with a clear geographic pattern, being the haplotype H2 recorded in almost all the localities of southern Mexico and Central America, while the haplotype H47 was found in most of the South American localities ( Figure S3b).
The mtDNA haplotype network showed three clearly differentiated haplogroups with ≥30 mutational steps and genetic distances ≥6% ( Figure 5). Haplogroup I included 43 haplotypes from 183 P. fulvus specimens from Mexico-Costa Rica and 47 P. gymnonotus specimens from Mexico, Guatemala, Nicaragua and Costa Rica (sympatric distribution) and one specimen from Panama. This suggesting a mitochondrial capture of all the P. gymnonotus from the northern distribution range. The first haplogroup also shows two main haplotypes in a star-like structure, suggesting population expansion. Haplogroup II included 17 haplotypes, all from the 52 P. gymnonotus specimens from the southern distribution range with samples from Panama-Brazil (allopatric distribution) and showing expansion. Finally, the haplogroup III included three haplotypes, one from a P. fulvus from Costa Rica, another haplotype from P. gymnonotus from Mexico and a third haplotype shared by eight P. fulvus from Mexico and Costa Rica and seven P. gymnonotus from Mexico; this haplogroup was separated from the others by a genetic distance ≥6.8% ( Figure 5). Genetic differentiation was high between the two species, with F ST = 0.665 (p < 0.001) and pairwise F ST between localities of the different species showed higher values within P. fulvus than between P. gymnonotus (Table S7).
Diversity 2021, 13, x FOR PEER REVIEW 11 of 19 pairwise FST between localities of the different species showed higher values within P. fulvus than between P. gymnonotus (Table S7).

Microsatellinte Markers Checking
Null alleles were found in several loci in at least one locality. However, since the FST and genetic distance values calculated either with or without ENA correction did not change significantly (z < 0.0001, p = 1.0 and z = 0.4916, p = 0.6230, respectively) no locus had to be excluded from the analysis. Only eight loci (Pps1, Pps7, Pps10, Pps13, Pps14, Pps15, Pps16 and NAC8) showed evidence of HWE deviation in at least one locality, but with no discernable pattern. On the other hand, as LD was not detected consistently among the loci examined, no population or locus was ruled out.

nDNA Genetic Population Structure
Considering all the loci, a total of 168 alleles were identified in 155 P. fulvus bats from six localities and 111 alleles in 44 P. gymnonotus bats from two localities. A total of 23 private alleles (Np) were also identified: ten in P. fulvus and 13 in P. gymnonotus, all of them with a low frequency (0.010-0.145). HO values. For P. fulvus ranged from 0.523 to 0.612, HE from 0.644 to 0.668 and the average AR from 4.208 to 4.333; the corresponding ranges for P. gymnonotus were HO: 0.549-0.574, HE: 0.568-0.622, average AR: 3.331-3.804 (Table S7). The between-species genetic differentiation value, RST = 0.089, was significantly (p < 0.001) and the pairwise analysis of RST between localities of the different species showed higher intraspecific RST values (Table S7).
The results from the Bayesian analysis carried out with STRUCTURE supported the split of the individuals into two groups (K = 2, Figure 6) corresponding to the species P. fulvus and P. gymnonotus; additionally, six P. fulvus specimens from the Mexican states of Colima, Tabasco and Oaxaca were classified as hybrids. Following the criterion 2 of Burgarella et al. [78], NEWHYBRIDS identified three individuals as hybrid, one in P. fulvus and two in P. gymnonotus from the Mexican state of Tabasco. In addition, NEWHY-BRIDS left 11 P. fulvus individuals as unclassified, three of which can also be considered as hybrids based on the probability values in STRUCTURE. One other bat from Tabasco that was identified morphologically as P. gymnonotus was classified as P. fulvus by STRUCTURE and NEWHYBRIDS ( Figure 6). In summary, 92.96% and 96.98% of the 199 individuals analyzed with STRUCTURE and NEWHYBRIDS, respectively, were classified consistently with the species assignment based on morphological traits but according to

Microsatellinte Markers Checking
Null alleles were found in several loci in at least one locality. However, since the F ST and genetic distance values calculated either with or without ENA correction did not change significantly (z < 0.0001, p = 1.0 and z = 0.4916, p = 0.6230, respectively) no locus had to be excluded from the analysis. Only eight loci (Pps1, Pps7, Pps10, Pps13, Pps14, Pps15, Pps16 and NAC8) showed evidence of HWE deviation in at least one locality, but with no discernable pattern. On the other hand, as LD was not detected consistently among the loci examined, no population or locus was ruled out.

nDNA Genetic Population Structure
Considering all the loci, a total of 168 alleles were identified in 155 P. fulvus bats from six localities and 111 alleles in 44 P. gymnonotus bats from two localities. A total of 23 private alleles (Np) were also identified: ten in P. fulvus and 13 in P. gymnonotus, all of them with a low frequency (0.010-0.145). H O values. For P. fulvus ranged from 0.523 to 0.612, H E from 0.644 to 0.668 and the average AR from 4.208 to 4.333; the corresponding ranges for P. gymnonotus were H O : 0.549-0.574, H E : 0.568-0.622, average AR: 3.331-3.804 (Table S7). The between-species genetic differentiation value, R ST = 0.089, was significantly (p < 0.001) and the pairwise analysis of R ST between localities of the different species showed higher intraspecific R ST values (Table S7).
The results from the Bayesian analysis carried out with STRUCTURE supported the split of the individuals into two groups (K = 2, Figure 6) corresponding to the species P. fulvus and P. gymnonotus; additionally, six P. fulvus specimens from the Mexican states of Colima, Tabasco and Oaxaca were classified as hybrids. Following the criterion 2 of Burgarella et al. [78], NEWHYBRIDS identified three individuals as hybrid, one in P. fulvus and two in P. gymnonotus from the Mexican state of Tabasco. In addition, NEWHYBRIDS left 11 P. fulvus individuals as unclassified, three of which can also be considered as hybrids based on the probability values in STRUCTURE. One other bat from Tabasco that was identified morphologically as P. gymnonotus was classified as P. fulvus by STRUCTURE and NEWHYBRIDS ( Figure 6). In summary, 92.96% and 96.98% of the 199 individuals analyzed with STRUCTURE and NEWHYBRIDS, respectively, were classified consistently with the species assignment based on morphological traits but according to the combined results of the two analyses, six individuals (3.01%) were regarded as hybrids (Table S8). The intermediate-sized P. fulvus specimens (forearm length ≥46 mm) did not correspond with any of the individuals identified as hybrids, either by mitochondrial capture or based on microsatellites. This shows that there is no relationship between putative hybrids identified based on forearm length and those identified by the phylogenetic analysis or genotypic assignment from microsatellites.  Table 1.

Discussion
Morphometric differences have been traditionally used for differentiating species in the family Mormoopidae; particularly, forearm length is commonly used for distinguishing between species in the subgenus Pteronotus [24]. However, individuals with intermediate forearms have been recorded in the zone of sympatry of P. fulvus and P. gymnonotus, which have been either described as a separate subspecies or reported as the result of hybridization between the two species [36,37,88]. Our study shows that bats with intermediate forearms (46 mm ≤ forearm length < 49 mm) were classified as P. fulvus based on both the nuclear genes and the microsatellites. Apparently lacking a genetic basis, the clinal variation in size (measured in terms of FA) observed in both species may be duethe result of phenotypic plasticity processes associated with environmental gradients The is suggested by the significant correlation found between size and latitude in this study and also reported by Smith [24], who found that the body size of P. fulvus increases from north to south across its distribution range, while the size of P. gymnonotus decreases from north to south.
The recognition of three evolutionary lineages corresponding to the species P. davyi, P. fulvus and P. gymnonotus, as identified by Pavan and Marroig [31] and Clare et al. [30] is supported also by the results of the species tree. However, the geographic delimitation of these lineages in Central America is still unclear and requires further detailed sampling. Although it seems clear from our study that only P. fulvus and P. gymnonotus occur in sympatry in Mexico.
The phylogenetic relationships between these three lineages inferred from nuclear and mitochondrial markers were discordant, suggesting the action of genetic interactions between the three Pteronotus. The discordances in the relationships built from nuclear and mitochondrial markers might be attributed to an incomplete lineage sorting or a genetic * * * * * * * * * a.
The intermediate-sized P. fulvus specimens (forearm length ≥ 46 mm) did not correspond with any of the individuals identified as hybrids, either by mitochondrial capture or based on microsatellites. This shows that there is no relationship between putative hybrids identified based on forearm length and those identified by the phylogenetic analysis or genotypic assignment from microsatellites.

Discussion
Morphometric differences have been traditionally used for differentiating species in the family Mormoopidae; particularly, forearm length is commonly used for distinguishing between species in the subgenus Pteronotus [24]. However, individuals with intermediate forearms have been recorded in the zone of sympatry of P. fulvus and P. gymnonotus, which have been either described as a separate subspecies or reported as the result of hybridization between the two species [36,37,88]. Our study shows that bats with intermediate forearms (46 mm ≤ forearm length < 49 mm) were classified as P. fulvus based on both the nuclear genes and the microsatellites. Apparently lacking a genetic basis, the clinal variation in size (measured in terms of FA) observed in both species may be due -the result of phenotypic plasticity processes associated with environmental gradients The is suggested by the significant correlation found between size and latitude in this study and also reported by Smith [24], who found that the body size of P. fulvus increases from north to south across its distribution range, while the size of P. gymnonotus decreases from north to south.
The recognition of three evolutionary lineages corresponding to the species P. davyi, P. fulvus and P. gymnonotus, as identified by Pavan and Marroig [31] and Clare et al. [30] is supported also by the results of the species tree. However, the geographic delimitation of these lineages in Central America is still unclear and requires further detailed sampling. Although it seems clear from our study that only P. fulvus and P. gymnonotus occur in sympatry in Mexico.
The phylogenetic relationships between these three lineages inferred from nuclear and mitochondrial markers were discordant, suggesting the action of genetic interactions between the three Pteronotus. The discordances in the relationships built from nuclear and mitochondrial markers might be attributed to an incomplete lineage sorting or a genetic introgression; since both processes yield similar genetic patterns, it is difficult to distinguish them [9,89]. The fact that we found that P. fulvus, P. gymnonotus and P. davyi share mtDNA haplotypes in a small clade that shows no discernable geographic structure (Costa Rica, Oaxaca and Tabasco in Mexico and the Lesser Antilles) suggests that this clade represents an ancestral polymorphism conserved in these lineages due to incomplete sorting since their separation 0.833 My ago. On the other hand, all the P. gymnonotus specimens from Mexico-Costa Rica are genetically closer to P. fulvus than to the South American P. gymnonotus according to the mtDNA-based analyses. This suggests that P. gymnonotus is polyphyletic with respect to this marker, as opposed to the monophyletic groping observed in all P. gymnonotus with the nuclear markers. This explains why the genetic distance (based on the mtDNA marker, COI) intraspecific are higher than interspecific.
Pteronotus gymnonotus and P. fulvus from or near the sympatric zone (Mexico-Panama) have very similar or identical mitochondrial haplotypes. Both the topologies and the network suggest that these two species have had several contacts and finally, completing a mtDNA capture process, as Pavan [88] suggested based in one specimen from Guatemala. The results of this study show that bats from the northern population of P. gymnonotus have replaced their mtDNA with that of the resident species P. fulvus. According to our dating, this process could have occurred during the interglacial periods, when this species could advance northward following the formation of new tropical environments along Central America and the Pacific and Atlantic coasts of Mexico. This mitochondrial capture is also shown by the paired F ST between localities of the same species, with higher F ST values between some localities of P. fulvus than between P. fulvus and P. gymnonotus. The fact that two of the three P. gymnonotus bats from Panama showed the mtDNA corresponding to the South American lineage while the other showed the captured mtDNA suggests that the boundary of this genetic exchange is located around this latitude in Central America.
The genetic exchange observed in the mtDNA was not detected at the nuclear level, where we found a total differentiation between the sequences of the two species, as expected [90] and as observed in other vertebrates [91]. This differential pattern in mtDNA and nDNA has been extensively documented and explained by various mechanisms [6,[91][92][93][94][95]. Overall, mitochondrial replacement is a hybridization mechanism widely recorded in mammals [96] and particularly in bats, including species of similar size (e.g., Myotis myotis and M. blythii [97]) and, as in this case, species differing in size (e.g., Eptesicus serotinus and E. nilsonii [98,99]).
Theoretical models predict this type of exchange to be usually asymmetric and characteristic of species that have undergone abrupt expansion [90]. The asymmetry is mainly explained by demographic differences between the resident (donor) and the colonizing species (undergoing demographic expansion), which usually receives and fixes the new genetic material (Hubbs' or neutral model) [92,100]. These demographic differences are still observed in Mexico, where the sympatric caves typically harbor large populations of P. fulvus and only few individuals of P. gymnonotus [55]. This genetic capture process does not necessarily involve selective forces, although it seems they may be relevant in some cases [96,101].
In fact, an alternative or supplementary scenario cannot be ruled out, in which selective forces might have favored the movement of local "advantageous" genes from P. fulvus to P. gymnonotus to improve the adaptation to the local environment [102,103]. Some studies have suggested adaptive selection in mitochondrial genes in mammals [104,105], recognizing that mitochondrial metabolism is highly sensitive to environmental conditions, especially variations in temperature [101]. According to our dating (1.021-0.243 Ma), these lineages diverged recently in the Pliocene-Pleistocene ( Figure S2), periods when drastic changes in the distribution of tropical forests in Central America and the north of South America took place following the successive glacial-interglacial cycles and fluctuations in sea level that led to the alternating expansion and retraction of the fauna associated with those biomes [106]. These environmental changes, driven by glacial-interglacial cycles, have affected the patterns of gene flow between populations and created similar genetic patterns in markedly different species [106,107] and may have promoted the diversification of the subgenus Pteronotus in this region [67].
The values of genetic differentiation (F ST = 0.150 and R ST = 0.089) between P. gymnonotus and P. fulvus are generally higher than the differentiation reported for other bat species [97,108,109], even within the same genus (e.g., Pteronotus alitonus and P. rubiginosus) [110]. These results confirm that, although both species currently, share shelters in the sympatric zone and their reproductive cycles partly overlap [39,[111][112][113], they remain differentiated as independent evolutionary lineages and as clearly distinct species usually recognizable morphologically by their different size [24,31]. As demonstrated by microsatellites, this does not imply that their gene pools are necessarily isolated. In fact, our results identified six hybrid individuals (3.01% of the total sample). This relatively low hybridization level may be related to differences in the phenology of their reproductive cycles, or to morphological or acoustic differentiation, as reported in other Pteronotus species that share shelters and where asymmetric hybridization has also been detected [101].

Conclusions
This study examines size variation and evaluated historical and current hybridization between two sympatric species. Phylogenetic analyses have suggested that both species were in contact and exchanged gene material on several occasions, resulting in a capture of P. fulvus mtDNA by all the P. gymnonotus individuals across the sympatric distribution range (Mexico-Costa Rica), related to pulses of population expansion/retraction in P. gymnonotus associated with Quaternary climatic changes. Such mitochondrial capture might have also been influenced by an evolutionary advantage; a hypothesis that deserves further research. Our microsatellite analyses revealed the existence of ongoing genetic contact between the two species in sympatry, identifying that 3.01% of the specimens sampled are hybrids. No relationship was found between hybridization and the intermediate bats from southern Mexico. However, forearm length was found to be correlated with latitudinal gradients, indicating a possible plasticity response to environmental gradients and the need to be cautious when using forearm length as a diagnostic between these two species, particularly in the sympatric zone.
The results from this study demonstrate that genetic exchange between sympatric, evolutionarily close bat species is an evolutionary phenomenon that may be more common than currently recognized. These results also confirm, once again, the importance of combining multiple molecular markers when studying the identity of lineages and reconstructing evolutionary relationships between species. Therefore, molecular identification results based only on mtDNA must be interpreted with caution; nuclear genetic data, as well as morphological and ecological evidence [7], need to be taken into consideration to ensure the accuracy of any systematic conclusions. Finally, and from a conservation standpoint, the introgression presented in this study does not seem to imply any problem to the genetic make-up and species identity for any of the two involved species. However, conservation concern may rise in relation to the Mexican lineage of P. gymnonotus which is much rare than the relatively common P. fulvus and which population clearly represents a different evolutionary lineage from the nominal species. Therefore, this Mexican lineage needs to be considered as a conservation unit by itself that probably deserves taxonomic recognition.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/d13050194/s1, Figure S1. Correlation of forearm sizes and latitude of P. fulvus and P. gymnonotus. Figure S2. Bayesian reconstruction with divergence time for nDNA+mtDNA in BEAST. Bars show the 95% interval for the high posterior density (HPD). Divergence time value above each branch, time in millions of years ago. Figure S3. Distribution haplotype for a. P. fulvus and b. P. gymnonotus. Table S1: Study localities, type of amplified marker (B = COI; C = PRKC1; D = STAT5A; F = RAG2; M = Microsatellites) and GenBank access number. Table S2. Primers used for the amplification of fragments of RAG2, intrones (PRKC1 and STAT5A) and COI. Table S3. New microsatellite loci [114][115][116][117]. Table S4. Haplotypes list from P. fulvus, P. davyi and P. gymnonotus for PRKC1 + STAT5A, RAG2 and COI. The number of haplotypes, the studied site and the number of individuals corresponding to each locality are shown in parenthesis. Table S5. Haplotype list for COI marker used for Pteronotus fulvus and P. gymnonotus. The number of haplotypes, the studied site and the number of individuals corresponding to each locality are shown in parenthesis. Table S6. Genetic diversity statistics of microsatellites for each study site of P. fulvus and P. gymnonotus. Samples number (N); summary and mean of the alleles per locus (Na); exclusive alleles (Np); observed heterozygosity (HO), expected heterozygosity (HE); summary and mean of allelic richness (AR). Study localities abbreviations correspond to the locations presented in Table S1. Table S7. Calculation of F ST up with COI and R ST down with microsatellites. *Significant. Table S8. Number of individuals identified as pure P. fulvus puro (Pfu); P. gymnonotus (Pgy); unclassified hybrids (SC); and hybrids according to STRUCTURE and NEWHYBRIDS with to criterion 2 of Burgarella et al. [87]. The percentage of the total number of individuals sampled by locality is included in parenthesis.