INTRODUCTION

The Neotropics include some of the most biodiverse terrestrial regions in the world (Smith et al., 2014). While Neotropical hotspots of biodiversity are important because they hold a high proportion of endemic species, nonhotspot parts of the Amazon basin have contributed many of the new mammal species described for this region (Ceballos and Ehrlich, 2009). Documenting this biodiversity is often challenging, particularly for genera whose many species are broadly distributed and span multiple habitats across ecological gradients. Untangling genetic and morphological variation across these gradients is fundamental for understanding the micro- and macroevolutionary processes that shape these speciose regions.

While species are the most basic biological unit (Burbrink et al., 2022; De Queiroz, 2007; Mayr, 1982), our ability to properly delimit them is often challenged by inadequate sampling and issues with the data and tools used to describe biological variation. Reliable species delimitation is particularly problematic when dealing with cryptic species that are nearly indistinguishable phenotypically, such as is the case with many bat taxa (e.g., Demos et al., 2020; Morales and Carstens, 2018). In a practical sense, these issues are sometimes magnified in unusual or iconic species that exhibit unique phenotypic features that can lead to overconfidence of identifications and may effectively conceal undocumented variation. DNA barcode initiatives have been proposed to aid in the discovery and identification of cryptic diversity (Hebert et al., 2003). However, barcode studies have been criticized due to their failure to test clear hypotheses, reliance on poorly performing phylogenetic distance methods, use of single uniparentally inherited markers, and excessive dependency on monophyly in tree topologies as a criterion for identification (Collins and Cruickshank, 2012). These difficulties can be overcome by integrating different sources of data and methods to increase the reliability for the recognition of a new taxon (Fujita et al., 2012; Padial et al., 2010). In light of this context, phylogeographic approaches combining environmental, genetic, and phenotypic data provide an opportunity to delimit species using multiple lines of evidence (Hickerson et al., 2010). Complementing these with single and multilocus species tree estimation methods (e.g., Yang and Rannala, 2014) and machine-learning classification models applied to phenotypic data can elucidate the underlying patterns of speciation and provide an integrative and replicable framework to define species limits.

Among mammals, bats (order Chiroptera) have received significant systematic attention and at least 500 new species have been delimited in the past 35 years (Simmons and Cirranello, 2024). Many of these discoveries resulted from investigations of widespread species complexes including cryptic species (e.g., Chaverri et al., 2016; Demos et al., 2020). In the Neotropics, substantial new species diversity has been recognized within “easy to identify” widespread species (e.g., Pteronotus parnellii) using integrative taxonomic methods including molecular and morphological data (Clare et al., 2013; Dávalos, 2006; de Thoisy et al., 2014; Pavan and Marroig, 2016). Key in such studies is the application of genetics and morphology to samples from across the geographic range of the complex.

Trachops Gray, 1847, is a monotypic genus of phyllostomid bat (Chiroptera, Phyllostomidae) with T. cirrhosus Spix, 1823, as the single currently recognized species (Burgin et al., 2024; Simmons and Cirranello, 2024). This bat is widely distributed, occupying areas of low-elevation forests (ca. <500 m) from southern Mexico to southeastern Brazil (Cramer et al., 2001; Medellín, 2019; Simmons, 2005). Three subspecies are recognized: Trachops cirrhosus coffini Goldman, 1925, in Central America from México to Panamá, T. c. ehrhardti Felten, 1956, in the Atlantic Forest of coastal Brazil, and T. c. cirrhosus Spix, 1823, occurring between the other two subspecies in multiple South American countries associated with Amazonia (Cramer et al., 2001; Medellín, 2019; Simmons, 2005).

Trachops, infamously known for its fondness for frogs, is an iconic bat among phyllostomids because of its brown wooly pelage, large ears, and the many wartlike protrusions on its chin and snout that make it very recognizable (Cramer et al., 2001). Google Scholar ( https://scholar.google.com/, accessed November 2023) lists nearly 2500 publications have discussed the ecology and distribution of this taxon. Two studies (Clare et al., 2011; Ditchfield, 1996) have explored the genetic variability of Trachops and shown high genetic variation within the genus. Mitochondrial cytochrome b (Cyt-b) analyses of Trachops recovered five clades with molecular distances ranging from 5.5%–11% among individuals from Honduras in Central America and Brazil in South America (Ditchfield, 1996). At the time, these levels of genetic variation were considered 10 times greater than those found in other species of bats in the same family with similar geographic distributions (e.g., Artibeus lituratus and Carollia perspicillata; Ditchfield, 2000). An additional DNA barcode study suggested up to nine distinct mitochondrial cytochrome oxidase (COI) lineages within the genus with these distributed between Guyana and Mexico (Clare et al., 2011). However, patterns of variation in the nuclear Y chromosome intron (Dby) could not confirm these putative lineages as distinct (Clare, 2011). Both phylogeographic studies, nonetheless, included only limited sampling throughout the distribution and primarily relied on assessments of genetic distance of single genes.

Our goal was to provide better resolution to the phylogeography of Trachops and test for species limits throughout its range. Following previous studies, we tested the hypotheses that Trachops is a monotypic genus with high intraspecific genetic and phenotypic diversity. Given its broad geographic distribution, we predicted that landscape environmental features are correlated with the variation within the genus associated with the three known subspecies. To test this, we integrated multilocus genetics examined under phylogenetic and species delimitation methods, analyzed phenotypic data using machine-learning classification algorithms, and evaluated environmental features where each subspecies exists using an ecological niche modeling (ENM) framework.

MATERIALS AND METHODS

Genetic sampling and sequencing: Tissue samples from 160 individuals of Trachops were obtained on loan from 10 museum collections. This sampling included individuals referable to the three recognized subspecies collected at 82 localities ranging from Mexico to southeastern Brazil. DNA extractions were performed to amplify three mitochondrial and one nuclear genes (details in the online supplement:  https://doi.org/10.5531/sd.sp.68). Additional sequences were downloaded from Genbank ( https://www.ncbi.nlm.nih.gov/genbank/) adding 132 individuals for cytochrome c oxidase I (COI) and 1 for cytochrome b (Cyt-b). Genes from samples for which georeferencing was not available were not considered in the analyses. All gene sequences produced were deposited in GenBank (table ST2) and details for individual specimens and their sequences are available online ( https://doi.org/10.5281/zenodo.8247231).

The gene sequences were aligned in Geneious 9.0.5 (Biomatters) using the Muscle algorithm (Edgar, 2004). We examined the degree of nucleotide base saturation and calculated pairwise uncorrected p genetic distances for the cytochrome b gene (table ST3) in the package ape (Paradis, 2012) in R 4.1.2 (R Core Team, 2022). Three multigene DNA sequence alignments were produced: one that maximized the number of individuals sampled (N = 111) including 2341 bp of mitochondrial COI, Cyt-b, D-loop and the nuclear intron STAT5A; one with reduced sequence gaps among individuals (N = 44) but maintaining sampling from each unique geographic area including 2331 bp of the same genes; and one including a wider sampling breadth (N = 60) and reducing sequence gaps using 1915 bp of COI, Cyt-b, and STAT5A to estimate dates of divergence. Sequence gaps comprised 12.2% and 2.4% for the full and subsampled alignments, respectively, used in phylogenetic inference. All alignments are accessible online ( https://doi.org/10.5281/zenodo.8247231).

Phylogenetic analyses: Phylogenies were reconstructed using Bayesian inference (BI) and maximum likelihood (ML) for single genes and multigene sequence alignments. We used PartitionFinder 1.1 following a Bayesian information criterion (BIC; Lanfear et al., 2012) to select a model of nucleotide substitution on each genetic marker and in the combined alignments, respectively. Analyses resulted in HKY+G as the best-fit model used to produce the BI tree search. Bayesian phylogenies were estimated using MrBayes 3.2.6 (Ronquist et al., 2012) within CIPRES ( https://www.phylo.org) using the HKY+G model of nucleotide substitution over 10 million generations for 2 runs of 4 independent Markov chain Monte Carlo (MCMC). Chains were sampled every 1000 generations and a 25% burn-in was applied. Bayesian posterior probabilities were used to assess branch support. Maximum likelihood (ML) phylogenies were estimated using IQTree (Nguyen et al., 2015) in the Amarel Cluster of the Office of Advanced Research Computing at Rutgers University. We separately used the model finder plus (mpf) option (Kalyaanamoorthy et al., 2017) in IQTree to obtain the best model of nucleotide substitution on a partitioned scheme, which resulted in GTR+F+I+G4. This model was implemented for the full (N = 111) and subsampled (N = 44) alignments in a full tree search including 1000 ultrafast bootstrap replicates with 1000 SH-aLRT (approximate likelihood ratio test) to assess branch support (Hoang et al., 2018). Genbank sequences from Tonatia saurophila (Koopman and Williams, 1951) and Carollia perspicillata (Linnaeus, 1789) were used as out-groups in ML and BI. All phylogenies were examined in FigTree 1.4 ( http://tree.bio.ed.ac.uk/software/figtree/).

Molecular divergence dating: Estimates of divergence times were performed in Beast 2.5 (Bouckaert et al., 2014) using a combined sequence alignment composed of Cyt-b, COI and STAT5A from a subset of 60 individuals from all geographic areas. The mitochondrial d-loop sequences were removed from this analysis to reduce noise of the population-level differences in this marker. Sequences of five outgroup taxa (Anoura geoffroyi, Artibeus glaucus, Carollia perspicillata, Desmodus rotundus, and Tonatia saurophila) were used to provide calibration points. We selected a calibrated Yule model as a prior under an uncorrelated log-normal relaxed clock. All parameters were unlinked, and the tree topology was constrained maintaining monophyly in the calibration points and the node of the two most inclusive clades of Trachops following the topology previously estimated in the BI and ML analyses.

Five fossil calibration points were selected to estimate divergence times in the bat family Phyllostomidae (Datzmann et al., 2010; Ho and Phillips, 2009). The ages of the chosen calibration points were consistent with a more recent study of Noctilionoidea (Rojas et al., 2016). These calibrations refer to the most recent common ancestors (MRCA) of the subfamily lineages of (1) Desmodontinae (32 [28–36] Ma), (2) Phyllostominae (29 [25–32] Ma), (3) Glossophaginae (27 [24–31] Ma), (4) Stenodermatinae-Carolliinae (22 [18–27] Ma), and (5) between the genera Tonatia and Trachops (25 [20–28] Ma). The standard deviation for each node showed 95% high posterior density interval (HPD). Following da Silva Fonseca (2019), three independent MCMC runs with 100 million generations sampled every 10,000 were implemented. Results were examined in Tracer 1.7.2 to assess stationarity in likelihood scores and evaluating sampling parameters based on effective sample size (ESS) values >200 (Drummond and Rambaut, 2007). LogCombiner 2.5.0 was used to combine the trees of each run excluding 25% of trees as burn-in (Bouckaert et al., 2014). The consensus tree was obtained using TreeAnnotator (Bouckaert et al., 2014) and visualized in FigTree 1.4.

Genetic species delimitation: To test the hypothesis that Trachops consists of multiple species, we used the subsampled alignment of 44 individuals to evaluate species limits using three approaches. First, we implemented a single-locus species delimitation analysis in mPTP (Kapli et al., 2017). This analysis was run on the server  https://mcmc-mptp.h-its.org/mcmc and used the non-ultrametric Cyt-b tree produced in IQTree. MCMC runs followed 10 × 10⁶ generations sampled every 10 thousand after a 10% burn-in. Three separate analyses were set with different starting delimitation models: null model (i.e., all lineages considered as constituting one species), maximum likelihood model (i.e., MLE based delimitation), and random model (i.e., random delimitation). In each analysis, the option --multi was used to include the intraspecific differences among rates of coalescence with a minimum branch length of 0.001 and provide node support.

Second, we implemented two multilocus coalescent approaches for species delimitation; SVDQuartets implemented in *PAUP (Chifman and Kubatko, 2014; Swofford, 2002) and starBeast3 (Douglas et al., 2022) implemented in Beast 2.7 (Douglas et al., 2022). In both analyses, we tested two-species limits hypotheses: (1) the three recognized subspecific groups (i.e., T. c. cirrhosus, T. c. coffini, and T. c. ehrhardti), and (2) the seven lineages recovered by phylogenetic analyses (see Results). Only COI, Cyt-b, and STAT5A were used in these analyses. We separately reevaluated the two-species hypothesis (i.e., T. cirrhosus and T. ehrhardti) from da Silva Fonseca (2019). The highly variable mitochondrial D-loop is informative at the population level; thus, it was excluded from these analyses to improve convergence of the coalescent algorithms. SVDQuartets calculates quartet-based scores that correspond to phylogenetic divergence and assumes that each locus in the dataset has its own genealogy (Chifman and Kubatko, 2014). The best topology is estimated based on quartet scores, which are then used to produce the species tree. A total of 57,621 possible quartets of the subsampled sequence alignment were examined using an exhaustive search followed by 100 bootstrap replicates to estimate support values, which were mapped at each node of a consensus tree. We used starBeast3, a fully Bayesian method that implements the multispecies coalescent (MSC) to infer gene trees by sharing information between loci through the species tree (Szöllősi et al., 2015). This used the HKY+G best fit model of substitution estimated in PartitionFinder 1.1 as in the MrBayes analysis, under a species tree relaxed clock, and a Yule model tree prior run over 50 million generations. MCMC operator priors were adjusted based on recommended settings following 10 preliminary runs and analysis convergence was assessed by examining that effective sample size (ESS) values were >200 in Tracer 1.7.2 (Rambaut et al., 2018). A 10% burn-in was applied to the resulting 4001 trees, and a maximum credibility tree was reconstructed in TreeAnnotator with posterior probabilities mapped at each node.

Morphological data and analyses: A total of 1093 skulls of adult individuals of the three currently recognized subspecies were measured using calipers (table ST4). We examined and measured the holotypes of T. c. ehrhardti and T. c. coffini; unfortunately, the holotype of T. c. cirrhosus has been lost (Carter and Dolan, 1978; Kraft, 1983). Ten dimensions were measured to the nearest 0.01 mm: greatest length of skull (GLS), condyle-incisive length (CIL), width between canines (CB), zygomatic breadth (ZB), mastoid breadth (MB), interorbital width (POC), breadth of braincase (BBC), length of the maxillary toothrow (CM3), width between upper molars (M3M3), and skull height (SH). Some measurements could not be taken on some specimens due to damage. To minimize data gaps, we partitioned the morphological data by subspecies and geography, ensured that missing data did not exceed 40% of each measured dimension, and then imputed missing values using the R package mice (Penone et al., 2014; van Buuren and Groothuis-Oudshoorn, 2011). The final morphological dataset comprised all 10 dimensions from 818 specimens: T. c. cirrhosus (N = 598), T. c. coffini (N = 185), and T. c. ehrhardti (N = 35) and is archived in Zenodo ( https://doi.org/10.5281/zenodo.8247231).

Two sample t-tests were used to examine mean size differences between males and females across all morphological measurements taken in each subspecies of Trachops. We then tested two phenotypic hypotheses. First, that each subspecific lineage of Trachops cirrhosus (i.e., cirrhosus, coffini, and ehrhardti) is morphologically distinct, and second that the seven recovered phylogenetic lineages of Trachops correspond to morphologically distinct groups. The degree of phenotypic divergence among groups was examined under four supervised machine-learning classification algorithms: classification and regression trees (CART), linear discriminant analysis (LDA), random forest (RF), and support vector machine (SVM) in the R packages caret v.6 (Kuhn, 2022) and MASS v.7.3 (Venables and Ripley, 2002). Results from each model were summarized to select the best performing model based on classification accuracy (i.e., how well the classifier assigned each partition to the correct group). Classification models were trained using a fixed 75% random data partition and tested using the remaining 25% of the data implementing a k-fold cross validation approach of five replicates. We also evaluated whether the overall accuracy rate was statistically greater than the no-information rate (Kuhn and Johnson, 2013). Phenotypic species limits were then visualized in a two-dimensional plot of the first two discriminant axes.

Spatial data and ecological niche analyses: We examined areas of environmental suitability for Trachops using an ecological niche modeling (ENM) and niche decomposition approach based on 19 climatic variables from WorldClim at a resolution of 2.5 min (ca. 5 km²; Fick and Hijmans, 2017). Georeferenced observations for Trachops were obtained from VertNet ( http://vertnet.org/), literature searches, and unpublished data provided by researchers. As a criterion, only observations with voucher specimens deposited in natural history collections were considered, for a grand total of 1760 records. Duplicate localities were removed using the R package dismo (Hijmans et al., 2023), and then rarefied to a 5 km² spatial resolution to ensure spatial independence. The resulting dataset contained 510 unique occurrence localities and is available at Zenodo. Three biologically informed partitions were selected to produce ENM based on the currently recognized subspecies hypothesis: T. c. cirrhosus (N = 375), T. c. coffini (N = 85), and T. c. ehrhardti (N = 50), which also matched the three species limits tested using genetics and morphology. Furthermore, we used the genetic information produced herein to examine whether environmental suitability helps explain the differentiation within multiple groups of T. c. cirrhosus. These data were limited because partitions of the five lineages within T. c. cirrhosus show overlapping geographic ranges; thus, reliable identification of individual bats without genetics was not possible. Final datasets used in modeling T. c. cirrhosus clades included three of the five clades, i.e., east Amazonia (eam) blue clade, N = 35; Guiana Shield (gs) light blue clade, N = 36; and northwest Amazonia (nwam) yellow clade, N = 13. The low sampling of individuals from the orange clade (pi, Panama/Venezuela, N = 3) and from the green clade (swam, southwest Amazonia, N = 10) precluded any reliable niche estimation or analysis of occupied niche space (see fig. 1 for color and geographic reference).

ENM were estimated in MaxEnt 3.4.1 implemented in the package dismo (Hijmans et al., 2023; Phillips et al., 2006). Environmental background selection for each partition was carried out in the packages dismo and sf (Pebesma, 2018) by generating a buffer of 100 km around each occurrence locality and then sampling 5000 random points from within it. The home range of T. cirrhosus (sensu lato) is estimated to be approximately 5 km² and this bat is known to travel short distances of about 4 km to reach feeding sites (Bernard and Fenton, 2003; Jones et al., 2017; Kalko et al., 1999). This resolution, however, may not reflect the dispersal potential of individuals in a population distributed in a large area of continuous habitat. A conservative buffer of 100 km was obtained from the dispersal capacity attributed to Artibeus lituratus (Esbérard et al., 2017), a larger bat of the same family. This helped reduce the likelihood of underestimating the available background environments for Trachops. Each model was calibrated using a 75%–25% training-testing scheme on the regional background extent for T. c. cirrhosus (and individually for each of the five lineages), T. c. coffini, or T. c. ehrhardti, then projected to the entire range of the genus. Appropriate model parameterization was ensured by producing an ENM using default parameters of β = 1 and allowing the MaxEnt algorithm to converge into feature classes (fc, either linear (l), product (p), quadratic (q), hinge (h), threshold (t), or a combination of these). Then, species-specific model parameter tuning was conducted using the package ENMeval 2.0 (Kass et al., 2021) by testing all combinations of feature classes and a range of regularization values from 1–5. The results were compared to the default model, and best model parameters were chosen based on the corrected Akaike Information Criterion values (AICc).

FIGURE 1.

Phylogenetic relationships of 44 individuals of Trachops based on four sequence markers (see text). A. Map indicating the distribution of genetically sampled groups and indicating wet forest biomes in green. Color codes in each group correspond to the geographic localities across the distribution of the genus, namely coffini (red), cirrhosus as Panama Isthmus and Venezuela (orange), northwest Amazonia (yellow), southwest Amazonia (green), Guiana shield (light blue), east Amazonia (dark blue), and ehrhardti (purple). B. Maximum likelihood tree estimated in IQTree based on the concatenation of Cyt-b, COI, d-loop, and STAT5A sequences (Bayesian topology was identical and included in fig. SF1A). Values at nodes indicate ultrafast bootstrap from IQTree (left) and Bayesian posterior probabilities from MrBayes (right). Outgroups not shown. Names of terminal taxa include group designation, institution and catalog number of the specimen, and geographic locality. Letters a–m at nodes are simple labels used as aid to discuss the results.

Model performance for each group prediction was assessed using the area under the receiver operating characteristic curve (AUC) estimated in Maxent. AUC is a threshold-independent measure that varies from 0 to 1, where a score of 1 represents perfect discrimination and a score below 0.5 represents a model no better than random (Peterson et al., 2011). We considered an AUC score greater than 0.7 to represent good model predictions (Peterson et al., 2011). Given that AUC has been deemed unreliable for estimating performance of presence-background models (e.g., Lobo et al., 2008), we separately calculated the Boyce Index (BI) to assess model prediction in the R package ecospat (Di Cola et al., 2017; Hirzel et al., 2006). The BI uses a Spearman rank correlation coefficient, which values range from -1 to 1 (Hirzel et al., 2006). A positive BI value indicates that model predictions are consistent with the evaluation dataset, zero indicates random performance, and negative values indicate a poor match with the evaluation dataset (Hirzel et al., 2006). Final model projections used a minimum training presence threshold of 95% (MTPT) obtained from the model estimated by the Maxent cloglog output (Soto-Centeno and Steadman, 2015). This relaxed threshold provided a conservative prediction where training model visualizations contained at least 95% of locality points (i.e., a theoretical expectation of 5% omission rate of the training data (Waltari and Guralnick, 2009) to evaluate discontinuities in suitable habitat.

To quantify the variation in the occupied niches of T. c. cirrhosus, T. c. coffini, and T. c. ehrhardti, we used a centroid shift, overlap, unfilling, and expansion (COUE) scheme (Guisan et al., 2014). This framework uses an ordination approach to overcome biases associated with quantifying niche dynamics in geographic space, with differential sampling efforts and/or spatial resolution, implemented in ecospat v.3.2 (Broennimann et al., 2012; Di Cola et al., 2017: 201; Liu et al., 2020). The analysis was performed using all unique observations (N = 510), 19 WorldClim environmental variables, and elevation, and could not be performed for the five T. c. cirrhosus clades individually due to the paucity of data. We partitioned the spatial data based on the three subspecific groups, generated three pairwise comparisons (T. c. cirrhosus vs. T. c. coffini, T. c. cirrhosus vs. T. c. ehrhardti, and T. c. coffini vs. T. c. ehrhardti), and computed the degree of niche overlap in each case (i.e., Shoener's D and Warren's I; (Schoener, 1968; Warren et al., 2008, 2021). The value of niche overlap ranges from 0 to 1, where zero represents no overlap, and one represents complete overlap. Niche identity (i.e., niche equivalency; Warren et al., 2008) tests were performed in the R package ENMtools to assess whether the ecological niches of each subspecies of Trachops cirrhosus were significantly different or interchangeable from each other (Warren et al., 2008, 2021). This analysis consisted of pairwise comparisons of D (and I) overlap values to a null distribution based on 1000 random replicates (Di Cola et al., 2017; Warren et al., 2021). Identity was determined if the observed D (and I) overlap values were significantly different from the overlap values in the null distribution. We also performed 1000 replicates of symmetric background niche similarity tests in ENMtools to assess if the background of environmental niche space for each pairwise species comparison was more similar than expected by chance (Warren et al., 2008, 2021).

RESULTS

Phylogenetics and species limits: Pairwise uncorrected p genetic distances revealed high levels of variation among the seven different lineages within Trachops, with percent genetic distances ranging from 4.6% up to 10% (table ST3). Phylogenetic reconstructions for the complete (N = 111 individuals) and subsampled (N = 44 individuals) alignments resulted in identical topologies for both BI and ML analyses. For conciseness, we present results based on the subsampled alignment; results from the full alignment and for single-gene trees can be found in figure SF1. The BI and ML topologies for the multigene analyses resulted in log-likelihoods (-ln) of -9991.30 and -8643.632, respectively. Both congruently recovered seven monophyletic geographically structured clades with high support (fig. 1). A single T. c. ehrhardti lineage from the Atlantic Forest of southeastern Brazil was recovered as sister to the remaining groups. The divergence of T. c. ehrhardti was dated to the late Miocene ca. 7 Mya and separated from their closest geographic relatives in eastern Amazonia and northeastern Brazil by at least 2 million years (node a; table 1; fig. 1; fig. SF2). The most recently divergent lineage corresponded to T. c. coffini, which separated from northwestern South America and Panamá T. c. cirrhosus ca. 3 Mya (table 1; node j in fig. 1; fig. SF2). Five clades were found to comprise T. c. cirrhosus, with deeper divergences in eastern Amazonia of Brazil and more recent divergences in northwest South America and Panamá (i.e., 5.04 to 3.83 Mya; table 1; nodes b–j in fig. 1; fig. SF2). These five clades were monophyletic; nevertheless, some showed geographic contact zones at the core of the Amazon (i.e., east, northwest, and southwest Amazonia, and Guiana Shield; fig. 1).

The single-locus mPTP species tree analysis confirmed with strong support that Trachops could consist of up to seven species (i.e., each geographic group labeled in fig. 1; fig. SF3). The three independent methods used as starting delimitation in mPTP (i.e., null model, maximum likelihood, or random) all inferred species trees congruent with the seven geographically structured clades (fig. SF3). The hypothesis test of three species (i.e., T. cirrhosus, T. coffini, and T. ehrhardti) was recovered in SVDQuartets and starBeast3 both delimiting each lineage as a species with strong support (fig. 2A). For the test delimiting seven species (i.e., 7 species lineages recovered in IQTree, fig. 1), both methods recovered T. ehrhardti with strong support (fig. 2B, C). Notwithstanding, these analyses produced slightly different topologies than the Bayesian and maximum likelihood results and showed variable support for the clades of T. cirrhosus and T. coffini (fig. 2B, C). Specifically, the SVDQuartets seven-species tree analysis recovered T. cirrhosus from east Amazonia as a species with good support (bootstrap = 84). The starBeast3 seven-species tree analysis recovered T. cirrhosus from southwest Amazonia as a species with strong support (PP = 1) and T. coffini as a species with moderate support (PP = 0.89). Results from the reexamination of the two-species hypothesis recovered T. cirrhosus and T. ehrhardti as species with strong support confirming previous analyses (da Silva Fonseca, 2019; fig. SF4).

Phenotypic limits: Males were larger on average for nearly all measurements in T. c. cirrhosus and T. c. coffini (P <0.05), whereas no significant differences were observed between sexes in T. c. ehrhardti (P >0.05; table ST4). Among the four phenotypic classification algorithms tested, LDA and SVM resulted in similar accuracy, 0.91 (95% CI: 0.896, 0.920) and 0.90 (95% CI: 0.895, 0.914) respectively (fig. SF5). Given the higher accuracy of LDA, we discuss phenotypic results based on this model. In the first test of whether the lineages T. c. cirrhosus, T. c. coffini, and T. c. ehrhardti are phenotypically distinct, LDA achieved a percent discrimination of 96% on LD1 and 4% on LD2 (fig. 3). The accuracy of the LDA model was significantly better than the no information rate (0.738, P <0.005). The phenotypic classifier discriminated T. c. cirrhosus with 97% certainty and showed some phenotypic overlap with T. c. coffini (i.e., 3% cirrhosus incorrectly classified as coffini). The level of phenotypic discrimination in T. c. coffini was 90%, and 10% showed classification overlap with T. c. cirrhosus. In T. c. ehrhardti, 14% were correctly classified; this group had considerable phenotypic overlap with T. c. coffini (66%), and some classification overlap with T. c. cirrhosus (20%; fig. 3). The craniodental features discriminating Trachops into three morphological groups were associated with the length of the cranium, specifically length of the maxillary toothrow (CM3) and the greatest length of skull (GLS). Along LD1, T. c. cirrhosus was phenotypically distinguished from the other two groups. Along the same axis, T. c. coffini and T. c. ehrhardti are phenotypically convergent, despite their distant phylogenetic relationships and significant geographic separation (>3000 km). All three groups show a greater overlap along LD2 than LD1. We also examined the hypothesis that T. c. cirrhosus consists of five distinctive phenotypic groups congruent with the genetic results (i.e., nodes c, e, g, i, and k in fig. 1). Pairwise LDA comparisons of the first three linear discriminants showed significant phenotypic overlap among these groups (fig. SF6). The level of overlap observed was congruent with previous findings based on principal component analysis (da Silva Fonseca, 2019).

FIGURE 2.

Species tree analyses conducted in SVDQuartets and starBeast3 examining two hypotheses of delimitation in Trachops. Support values at nodes represent bootstrap or posterior probabilities where indicated. A. Three-species tree was recovered by both analyses with strong support (only starBeast3 tree shown). Node values represent bootstrap (left) and posterior probabilities (right). B. Seven-species tree from SVDQuartets recovered only ehrhardti as a species with strong support (i.e., bootstrap = 100%) and cirrhosus from east Amazonia was recovered with moderate bootstrap support (i.e., bootstrap >70%). C. Seven-species tree from starBeast3 recovered ehrhardti and cirrhosus from SW Amazonia with strong support (i.e., PP >0.98), and coffini and cirrhosus from Panamá with moderate support (i.e., PP = 0.89).

Ecological niche modeling and niche dynamics: Given that phylogenetic and phenotypic species-delimitation analyses favored the hypothesis of three distinct groups (see genetic and phenotypic delimitation above), we used these as the biologically relevant partitions of the spatial data, and therefore modeled them independently. Different species-specific parameters were used for calibrating individual ecological niche models (ENM; cirrhosus fc = h, rm = 1; coffini fc = lqh, rm = 2; ehrhardti fc = lqhp, rm = 2; table ST5). Each individually built ENM showed good performance for T. c. cirrhosus (AUC = 0.822; BI = 0.986), T. c. coffini (AUC = 0.849; BI = 0.981), and T. c. ehrhardti (AUC = 0.958; BI = 0.955). The ENM, calibrated with observations and background for each partition, broadly showed available suitable habitat when projected across the entire geographic range of Trachops (fig. 4). However, some patterns of differentiation were evident when examining the suitability values. Specifically, suitable habitat for T. c. cirrhosus decreased toward the southeast part of Brazil (fig. 4A). Areas of highly suitable habitat for T. c. coffini expand from the core of Mesoamerica into northwestern South America. Although Panamá, eastern Colombia, and southwest Venezuela appear as highly suitable, these areas are partly isolated from parts of Mesoamerica by low suitability in Nicaragua and in the Colombian Andes (fig. 4B). The suitable habitat available for T. c. ehrhardti is concentrated in the Atlantic Forest of southeastern Brazil and surrounded by areas of low suitability northwesterly into the dry diagonal (fig. 4C).

FIGURE 3.

Results from the machine-learning Linear Discriminant Analysis (LDA) classifier of phenotypic limits in Trachops. Overall accuracy of the model = 91.8% (95% CI 89.7%–93.6%). A. In biplot, percent group separation described for the first linear discriminant (LD1) was 96, and 4 for the second linear discriminant (LD2). Solid lines represent 68% data ellipses centered at the bivariate mean to visualize phenotypic differences among species. LD1 density values plotted for aid in visualization on the x-axis highlighting closer phenotypic affinities of coffini (red) and ehrhardti (purple) compared to cirrhosus (gray, all five geographic groups included). B. Map of the distribution of individual specimens measured for morphological analyses.

FIGURE 4.

Ecological niche models generated using three biologically informed partitions for Trachops: A. cirrhosus; B. coffini; and C. ehrhardti. An MTPT threshold was applied in model projections. Models were calibrated using each data partition and then projected onto the geographic range for the genus. Color gradients represent habitat suitability from high (red) to low (blue). Rectangular polygons in inset maps denote the geographic region from each data partition.

TABLE 1.

Divergence time estimation of clades within the bat genus Trachops estimated in Beast 2.5. Nodes a–m refer to labeled nodes in figure 1. Clades define the relationships separating the node at the most recent common ancestor at the divergence of each group. Abbreviations: cof = T. coffini; eam = eastern Amazonia; ehr = T. ehrhardti; gs = Guiana shield; nwam = northwest Amazonia; pi = Panama Isthmus; swam = southwest Amazonia. Main nodes of specieslike clades identified in bold font (see also fig. SF2).

Niche decomposition analyses revealed differential patterns of niche occupancy in pairwise comparisons. The comparison of occupied niche between T. c. cirrhosus and T. c. coffini showed the highest degree of overlap (Schoener's D = 0.38). In contrast, both pairwise comparisons of T. c. ehrhardti with either T. c. cirrhosus or T. c. coffini showed little overlap (Schoener's D = 0.12 and 0.18, respectively; table 2). Occupied niches in environmental space varied among pairwise comparisons (fig. SF7). The comparison of T. c. cirrhosus vs. T. c. coffini emphasized the high degree of niche overlap and no evidence of niche shift between these lineages (fig. SF7A). In contrast, comparisons of T. c. ehrhardti with either T. c. cirrhosus or T. c. coffini showed evidence that the areas occupied were characterized by niche shifts (fig. SF7B, C). Finally, the null hypothesis of niche identity was rejected in all pairwise species comparisons, revealing that niches were not identical and supporting the differences among the three groups (table 2; fig. SF8). The symmetric test of background similarity was rejected for both pairwise comparisons between T. c. cirrhosus with T. c. coffini and with T. c. ehrhardti (table 2; fig. SF8). In contrast, the test of background similarity could not be rejected for the comparison between T. c. coffini vs. T. c. ehrhardti (table 2; fig. SF8).

TABLE 2.

Pairwise group comparisons of niche overlap, identity (i.e., equivalency), and symmetric background similarity tests. Niche metrics of Schoener's D (Schoener, 1968) and Warren's I (Warren et al., 2008) values shown for each test. Boldface numbers of identity and background tests refer to significant differences P <0.05.

Model parameterization following ENMEVAL 2.0 resulted in different species-specific parameters for calibrating individual T. c. cirrhosus ENM (eam, fc = lq, rm = 1; gs: fc = h, rm = 2; nwam: fc = lqhp, rm = 2). ENM for two of the T. c. cirrhosus groups performed well, for east Amazonia (eam) AUC 0.829, BI = 0.937 and for Guiana Shield (gs) AUC 0.936, BI = 0.82. In contrast, T. c. cirrhosus for northwest Amazonia (nwam) show lower performance; AUC 0.747, BI = 0.742. Thus, results should be interpreted with caution. The T. c. cirrhosus clade from east Amazonia included the broadest distribution within the “cirrhosus” complex, spanning localities in the northeast Atlantic Forest, the Guiana Shield, the Dry Diagonal, and even central part of Amazonas state in Brazil along the Amazon River. Despite its broad distribution across a wide breadth of biomes, the ENM shows low suitability in areas where individuals from northwest Amazonia exist and in the southeast part of Atlantic Forest occupied by T. c. ehrhardti (fig. SF9 and fig. 4). The environmental suitability of T. c. cirrhosus from east Amazonia drops drastically along the Amazon River in central Amazonas state. Despite the low performance of the ENM of northwest Amazonia, we observed a similar pattern of a drastic decrease in suitability along the Amazon River (fig. SF9).

DISCUSSION

Deciphering species limits and their underlying variation is key for understanding the processes that shape poorly known yet speciose regions of the Neotropics. Species assemblies across the Neotropics span a combination of climatic, geographic, and topographic factors that help structure and maintain their diversity (Carvalho et al., 2023; Fluck et al., 2020). Biodiversity studies have shown that the most important predictors of hidden mammal species diversity are small body sizes and large geographic ranges (Parsons et al., 2022). In bats, widely distributed species showed variation that can be associated with the environmental features they occupy (Morales et al., 2016; Soto-Centeno and Simmons, 2022). Given the importance of the Neotropics as a center of biodiversity, documenting and delimiting taxa in this region can help better understand species richness and local diversity, which can aid in conservation initiatives. Herein, we studied the species limits of the iconic and widely distributed Neotropical bat Trachops and found evidence for the presence of at least three species in this genus.

Three groups of Trachops cirrhosus (i.e., cirrhosus, coffini, and ehrhardti) have traditionally been recognized as distinct at the subspecific level based on morphological evidence (Cramer et al., 2001; Medellín, 2019). Over the past two decades, two studies documented genetic variation within Trachops using mitochondrial DNA barcodes and found that this putative monotypic genus showed high genetic diversity that might be indicative of multiple species (Clare et al., 2011; Ditchfield, 2000). The presence of up to nine unique mitochondrial lineages with at least four found within the Guiana Shield documented by Clare et al. (2011) was quickly supplemented by discovery of the absence of reciprocally monophyletic clades in analyses of a single nuclear intron (Clare, 2011). However, due to the limitations of their data, both Ditchfield (2000) and Clare et al. (2011) relied on estimates of genetic distance and/or counts of clades identified within genes to document the presence of multiple independent specieslike lineages within Trachops. Our study included broader sampling, used assessments of biparentally inherited genes, morphology, and niche dynamics under a hypothesis-testing framework to examine species limits more comprehensively within the genus.

Relationships, divergence, and genetic species limits: Both Bayesian and maximum likelihood phylogenetic reconstructions produced identical topologies and showed the presence of seven strongly supported geographically structured clades of Trachops (fig. 1). Trachops originated likely in the early Miocene (Rojas et al., 2016). We estimated the divergence of Trachops from the outgroup Tonatia at ca. 25.55 Mya, which corresponds well with previous estimates in noctilionid bats (Rojas et al., 2016). Trachops underwent further diversification into the Pliocene beginning with a split of the clade corresponding to T. c. ehrhardti from other lineages about 7 Mya (table 1; fig. SF2). Five different clades referable to T. c. cirrhosus originated between 5–3 Mya and culminated in the Central American lineage corresponding to T. c. coffini diverging from T. c. cirrhosus at the Isthmus of Panamá about 2.96 Mya. This divergence age of T. c. coffini and the northernmost clade of T. c. cirrhosus coincides with the uplifting of the Isthmus (O'Dea et al., 2016), which may have facilitated the colonization and likely northward movement of these lineages. At the core of the Amazon, the east Amazonia (blue), southwest Amazonia (green), and northwest Amazonia (yellow) groups all formed monophyletic clades within T. c. cirrhosus (fig. 1). These clades also showed various contact zones, and some appear to be geographically delineated by major rivers in the Amazon Basin, as previously shown in other taxa (Musher et al., 2022; Smith et al., 2014). These rivers attained their modern configuration in the Neogene (i.e., 23–2.6 Mya). Specifically, the Amazon River separated T. c. cirrhosus of northwest Amazonia from east Amazonia and southwest Amazonia in a north to south arrangement. The Madeira and Purus rivers separated east Amazonia and southwest Amazonia. Finally, in contrast to four groups (Clare, 2011; Clare et al., 2011), we recovered two distinct clades associated with the Guiana Shield (light blue) and east Amazonia that are sympatric and showed a shallow divergence at ca. 5 Mya (table 1; fig. SF2). These two groups meet at the Raposa Serra do Sol mountain range in the northeastern part of the state of Roraima in Brazil.

We tested the hypothesis that Trachops is a monotypic genus that is widely distributed across the Neotropics using three coalescent-based species delimitation methods. The best-supported maximum likelihood phylogeny produced from an alignment of Cyt-b was used in the single-locus method mPTP, which recovered each geographically structured clade with high support, and suggested that Trachops could consist of up to seven species (fig. SF3). Further support of clade differences across all Trachops can be observed in the examination of percent sequence divergence of Cyt-b among clades, which ranged from 4.6% to 10% (table ST2). To examine the species limits of Trachops more objectively, we also used the multi-locus approaches SVDQuartets and starBeast3. In these analyses, we tested whether Trachops consisted of either three species (i.e., the three currently recognized subspecies: cirrhosus, coffini, and ehrhardti) or seven species (i.e., the seven monophyletic clades recovered from phylogenetic analyses). Both SVDQuartets and starBeast3 recovered the three species tree hypothesis with strong support across all nodes (fig. 2A). The three species in this phylogenetic arrangement are T. coffini, T. ehrhardti, and T. cirrhosus, the latter containing the five remaining clades (i.e., from north to south: Panamá Isthmus (orange), northwest Amazonia, southwest Amazonia, Guiana Shield, and east Amazonia; fig. 1). We used these new name combinations for the remaining discussion. While both analyses converged in similar topologies, most lineages in the seven species hypothesis were generally not strongly supported by either of the multilocus methods (fig. 2B, C). A highly supported lineage of T. ehrhardti was recovered by SVDQuartets and starBeast3. The latter method also showed the clade sister to T. ehrhardti with high support, and the lineage of T. coffini with moderate support, which provided additional confidence that these correspond to species level differences.

The hypothesis test of seven species of Trachops included five clades within the cirrhosus group. We observed conflicting evidence among results from mPTP, SVDQuartets, and Star-Beast3. While the five T. cirrhosus clades potentially represent a species complex based on the results of mPTP, we found no evidence from multilocus species delimitation analyses or morphology (see below) to support recognizing these clades as distinct species at this time. Previous studies investigating species limits have reported discordant results across some coalescent species delimitation approaches (e.g., reviewed in Carstens et al., 2013). The failure in the delimitation of the T. cirrhosus clades may have been due to lack of resolution in the genetic sequence data used and the small number of individuals sampled in each group. Simulation studies have suggested the use of 5–10 independent loci and at minimum 10 sampled individuals per group are necessary to ensure deep population-coalescent events are appropriately sampled (Carstens et al., 2013; Fujita et al., 2012). Our dataset did not meet these thresholds. Thus, additional efforts are needed to untangle further potential cryptic diversity within T. cirrhosus as we recognize it here. Genomic studies using ultraconserved elements (UCE), or restriction site associated DNA (RAD-seq), could provide the necessary resolution without the need of such high individual sample sizes to examine in detail both the relationships and the factors that help maintain diversity in Amazonian T. cirrhosus. Such studies are beyond the scope of the current study but hold promise for future research. We expect that additional data and analyses may confirm that there are multiple species in what we are here recognizing as a single species, T. cirrhosus.

Phenotypic delimitation and variation: The one-, three-, and seven-species hypotheses were also evaluated using an independent dataset of morphological measurements. Size variation in Trachops was first documented by Felten (1956a, 1956b) and Goldman (1925), who partitioned the three known subspecific groups based on external characteristics. No other study has evaluated the potential of multiple phenotypic groups within Trachops. A list of external and craniodental measurements was compiled by Cramer et al. (2001) and used to emphasize differences. However, character descriptions presented as measures of central tendency, as shown in Cramer et al. (2001), are of limited use because these cannot appropriately characterize distinct phenotypic distributions (Cadena et al., 2018). Based on a linear discriminant analysis of 10 craniodental measurements from 818 specimens across the entire distribution, we found support for patterns of variation in Trachops corresponding to three species (i.e., T. cirrhosus, T. coffini, and T. ehrhardti; fig. 3). Both T. coffini and T. ehrhardti are phenotypically distinct despite some morphological overlap with T. cirrhosus. Given their broadly disjunct geographic distributions, the morphological overlap between T. coffini and T. ehrhardti could likely represent a case of conserved phenotypes between these groups. This is plausible considering the earlier divergence of T. ehrhardti as a sister group to T. cirrhosus and T. coffini (fig. 1; table 1; fig. SF2). Individuals of T. cirrhosus that had similar phenotypic characteristics to T. coffini or T. ehrhardti occur in the Isthmus of Panamá and northwest and southwest Amazonia. This observation emphasizes the distinctiveness of T. ehrhardti, which is geographically disjunct, and confirms that T. coffini is still phenotypically distinct despite the closer phylogenetic affinities with geographically close T. cirrhosus.

Trachops is undoubtedly an iconic species with several key characteristics, such as large ears, brown dense wooly pelage, and wartlike protuberances on the rostrum, that make it unmistakable. Males of T. cirrhosus and T. coffini are significantly larger than females. In contrast, no differences were found between sexes in T. ehrhardti. From the LDA analysis, T. coffini and T. ehrhardti were generally smaller in size than T. cirrhosus for nine of the 10 measurements (i.e., excluding only the breadth of braincase, BBC). The size variation within T. cirrhosus, as defined herein, is broad and some individuals overlap in size with T. coffini and T. ehrhardti. The distinction among T. cirrhosus, T. coffini, and T. ehrhardti is evidenced in craniofacial features that seem associated with biogeographic breaks (da Silva Fonseca, 2019). Specifically, compared with T. cirrhosus, T. coffini can be distinguished by a steeper angle of the rostrum, a more pronounced talonid basin on the p1, and a lingual displacement of a larger p3 (fig. SF10, and da Silva Fonseca, 2019: fig. 3). Our description of the latter character provides additional diagnostic information to the original description by Goldman (1925), who observed a “vestigial lower premolar relatively large” in T. coffini. These characters, combined with a notch in the cutting edge of the upper incisors and a forearm <60 mm, differentiate T. coffini from T. cirrhosus (Goldman, 1925). Geographically, T. ehrhardti exists in the Atlantic Forest, south of the state of Espírito Santo, Brazil. We note the occurrence of T. ehrhardti (referred to T. c. ehrhardti, Goldman, 1925) in Bolivia as erroneous (see Jones and Carter, 1976). Several cranial characters differentiate T. ehrhardti from the widely variable T. cirrhosus. A comparison of T. ehrhardti with geographically adjacent T. cirrhosus from northeastern Brazil shows significant differences in mean measurements of the condyle-incisive length (CIL) and zygomatic breadth (ZB) as presented in the original description (Felten, 1956a). Measurements of skull length, particularly CIL and greatest length of skull (GLS) suggest that T. ehrhardti is generally smaller than T. cirrhosus, with the overlapping T. cirrhosus individuals being nearly exclusively females, which range in the smaller size of the latter species.

We observed reciprocal monophyly among five clades of T. cirrhosus that diverged between 5 to 3 Mya, yet some of them showed overlapping geographic ranges (fig. 1A). In other vertebrate systems, genetic structure across heterogeneous environments can be influenced by local adaptation and phenotypic disparities that are often associated with geographic distance (Myers et al., 2019; Seeholzer and Brumfield, 2018). Notably, we observed wide phenotypic variability in T. cirrhosus, and we could not use morphology to discriminate among the monophyletic genetic lineages associated with geography (fig. SF6). The data collected in this study, either genetic or morphological, does not provide enough resolution to determine whether the multiple genetic clades recovered in T. cirrhosus could constitute species. Further studies using genomic data are needed to understand whether phenotypes within this group are adaptively conserved or genetic divergences are driven by diversification over local environmental gradients, as well as to determine the role of landscape features in maintaining this phenotypic diversity.

Ecological differentiation: The ENM and niche-deconstruction analyses provided additional context in support of recognition of three species in Trachops (fig. 4, supplementary figures S7, S8). When niche models were projected across geographic space, we observed marked declines in habitat suitability at the geographic limits of T. cirrhosus, T. coffini, and T. ehrhardti (fig. 4). When modeling T. cirrhosus in geographic space (fig. 4) we observed low to zero habitat suitability in the southeast of Brazil where T. ehrhardti occurs. A similar result was observed when modeling T. ehrhardti; habitat suitability decreased drastically in areas outside the Brazilian Atlantic Forest. Previous studies have shown that Trachops in general uses small foraging areas, exhibits low mobility, and prefers the dense understory of unfragmented habitats, which in part has been attributed to its frog-hunting behavior (Halczok et al., 2018; Kalko et al., 1999). Thus, it has been proposed that Trachops could serve as an indicator species of habitat degradation (Kalko et al., 1999). Generally, a considerable part of the niche occupied by T. cirrhosus contains distinct environmental areas beyond the occupied or available niche of T. ehrhardti (i.e., shaded gray area outside the purple dotted lines in fig. SF7B). The pairwise comparisons of occupied niches of T. ehrhardti were characterized by a niche shift (fig. SF7B, C), had low overlap, and were not deemed identical (table 2) to niches of either T. cirrhosus or T. coffini (fig. SF7A). The current floristic differences between the southern and northern Atlantic Forest and the biotic and abiotic interactions involved may represent ecological barriers separating T. ehrhardti. We hypothesize that the distinction of T. ehrhardti could be driven by a steep gradient of unsuitable habitat beyond the southeastern Atlantic Forest that creates a barrier to dispersal.

The pairwise comparison between T. cirrhosus and T. coffini was characterized by higher niche overlap and no evidence of niche shift, which was expected under the hypothesis of niche conservatism (Wiens and Graham, 2005) and due to the close relationship of the two taxa (table 2; fig. SF7). Nevertheless, the niche identity and background similarity hypotheses were both rejected, suggesting that the niches occupied by T. cirrhosus and T. coffini are not identical and their available environmental background was significantly more similar (or different) than expected given their ranges (table 2; fig. SF8). The small home range of Trachops and its tendency to avoid fragmented forest patches (Halczok et al., 2018) likely contributed to the observed differentiation in niche occupancy between T. cirrhosus and T. coffini.

The background similarity test was not significant between T. coffini and T. ehrhardti and their occupied niches are characterized by a niche shift (table 2; figs. S7, S8). Areas of niche space occupied by T. coffini are beyond the niche available and occupied by T. ehrhardti (i.e., the shaded red area outside the dotted purple line in fig. SF7C). This suggested that these two species occupy distinct ecological niches. Taken in sum, pairwise comparisons of T. cirrhosus, T. coffini, and T. ehrhardti show distinctions in occupied niche space that support the recognition of three species. However, we emphasize that we measured only environmental features related to climate. In light of its sensitivity to habitat fragmentation, additional landscape-level studies at finer resolution could help determine whether natural or artificial (i.e., driven by deforestation) geographic barriers maintain the patterns of differentiation in Trachops.

Finally, the ecological niche models examining three of the five clades within T. cirrhosus (as recognized herein) showed evidence supporting the distinctiveness of these clades. The blue T. cirrhosus clade from east Amazonia (eam) included the broadest distribution within the cirrhosus complex. This blue clade spanned localities in the northeast Atlantic Forest, the Guiana Shield, the Dry Diagonal, and even the central part of Amazonas state in Brazil along the Amazon River. Despite its broad distribution across a wide breadth of biomes, the ENM for the east Amazonia clade shows low suitability in areas where individuals from yellow northwest Amazonia (nwam) exist and in the southeast part of Atlantic Forest occupied by T. ehrhardti (fig. SF9, fig. 4). The environmental suitability of T. cirrhosus from east Amazonia drops drastically along the Amazon River in central Amazonas state. Despite the lower performance of the ENM of the northwest Amazonia, we observed a similar pattern of a drastic decrease in suitability along the Amazon River (fig. SF9). This rapid change in environmental suitability between the east and northwest Amazonia clades suggests that the Amazon River may serve as a strong biogeographic barrier between these lineages. These two clades (i.e., blue and yellow) have a divergence of 1.74 million years (i.e., node b = 5.04 and node h =3.30 in table 1 and fig. SF2) and a 6.7% genetic distance in Cyt-b. The low dispersal ability of Trachops (Halczok et al., 2018; Kalko et al., 1999) might make it improbable for individuals to cross a large biogeographic barrier such as the Amazon River, which ranges between 3–10 km wide with flood plains of up to 50 km wide in the state of Amazonas (Trigg et al., 2012). In sum, these observations support the hypothesis that the Amazon River could help maintain isolation between the east and northwest Amazonia clades. ENM for the light blue Guiana Shield (gs) clade of T. cirrhosus showed high environmental suitability in the north areas of the Guiana Shield extending toward the coast. Suitability dropped quickly in habitats present in Venezuela where the northwest Amazonia clade exists, but also in the southern parts of Guyana and Suriname where the east Amazonia clade showed greater suitability. Given that Trachops is generally known from elevations lower than 500 m (Medellín, 2019). These results suggest the presence of another potential biogeographic barrier at the Raposa Serra do Sol in the northeastern part of the state of Roraima in Brazil. This mountain range extends from southeastern Venezuela eastward to Brazil and Guyana and has an average elevation of 487 m (65–2798 m). Future work to test whether these biogeographic barriers effectively preclude gene flow among clades of T. cirrhosus would require additional sampling and an integrative approach of present day and historical demographics, and associations with the geographic and geological features that help maintain differentiation.

Taxonomic implications: The integrative analyses we presented herein provided evidence that the bat genus Trachops is not monotypic. Our findings showed congruence of multiple methods demonstrating that three groups of Trachops (i.e., cirrhosus, coffini, and ehrhardti) are phylogenetically distinct, exhibit morphological disparities, and occupy distinctive ecological niches. Based on this evidence, we propose the current three subspecies of Trachops should be recognized at the species level. Name combinations in this taxonomic arrangement include Trachops cirrhosus Spix, 1823, ranging from Costa Rica through the Isthmus of Panamá into the northern countries of South America and ranging across east, northwest, and southwest Amazonia, the Guiana Shield, and northeast Brazil (see fig. 1), Trachops coffini Goldman, 1925 ranging from Costa Rica north into southern Mexico, and Trachops ehrhardti in the Atlantic Forest throughout southeast Brazil (see fig. 1). The holotype for Trachops coffini, Goldman, 1925, is from “Guyo,” Peten, Guatemala (National Museum of Natural History, USNM 244266). Trachops ehrhardti was originally described by Felten (1956), with a holotype from Santa Catarina, Brazil (Senckenberg Museum Frankfurt, Germany, SMF11716).

Unfortunately, the holotype of Trachops cirrhosus has been lost (Carter and Dolan, 1978; Kraft, 1983). This is significant because of the apparent paraphyletic nature of T. cirrhosus as we recognize it here. Although recognizing paraphyletic entities at the species level is not ideal, the reality of evolution means that some such entities exist in nature due to peripheral speciation, and the rejection of all such taxa outright is not justified (Crisp and Chandler, 1996; Hörandl and Stuessy, 2010). While we recognize the potential problems associated with taxonomic recognition of paraphyletic entities (be they at the species or subspecies level, e.g., Burbrink et al., 2022), in the case of T. cirrhosus we argue that this is the best solution at the current time given the data available. In the future, we anticipate that more comprehensive analyses of the clades detected within the species T. cirrhosus as defined here will be necessary to determine whether this taxon represents a species complex, and if so, which lineages represent distinct species. Notably, studies including more specimens, more comprehensive sampling across its geographic range, additional genetic markers, and additional phenotypic comparisons may result in recognition of multiple species within what we are here recognizing as a single species. Future revisionary work within T. cirrhosus should include designation of a neotype (with associated genetic material and sequences) from Amazonian Brazil where the original type specimen is thought to have originated (see Husson, 1962, who restricted the type locality to the state of Pará, Brazil; see also Vanzolini, 1981). We do not take that step here because we believe that it is better reserved for a complete revision of T. cirrhosus sensu stricto that addresses identity of its junior synonyms (angusticeps Gervais, 1856 and fuliginosus Gray, 1865).

We recognize that additional work is necessary to dissect the multiple clades detected in Trachops but emphasize that this should not come at the cost of accepting morphologically diagnosable and genetically distinct groups. The recognition of three species in Trachops allows us to better define the variation within the genus, and has conservation and protection benefits, particularly for Neotropical regions that are threatened by high rates of deforestation. For example, the current IUCN Red List status of Least Concern for Trachops is based on the simplification that all “T. cirrhosus” (sensu lato) represent a single, variable, widespread species. Nevertheless, the negative effect of human driven habitat degradation combined with the low mobility of this species led to the reclassification of T. c. coffini in Mexico as Federally Threatened (Medellín, 2019). We expect that future reviews of this taxon by the IUCN would produce similar results now that it will be recognized as a distinct species. Trachops ehrhardti from the Atlantic Forest will also be in urgent need for conservation evaluation. Despite efforts to completely stop illegal deforestation, estimates show that close to 13% of the native Atlantic Forest of Brazil remains, and only 9% of that is protected (De Praga Baião et al., 2023). Deforestation and fragmentation are attributed to high demand for urban and large-scale agricultural land. Yet it is unknown how this level of habitat loss and fragmentation has affected forest-dependent bats. We hope that the current taxonomic changes proposed herein for Trachops spark new interest in evaluating conservation status of these bats and their use as bioindicators throughout their distribution.