Genomic DNA was extracted from a set of 56 specimens belonging to the genus Molossus, sourced from various locations within Argentina (Table S1), and obtained through three distinct methods. A group of specimens was captured by fieldwork specialists during ecological studies of rabies virus, who followed a protocol involving the capture and subsequent release of the specimens. Wing membrane tissue samples were obtained using a biopsy punch with a diameter of 3 mm. Captures were carried out under the approval of the Ethics and Safety Committee of the Universidad Nacional del Litoral (Expte. FCV-0869428-17) and the Ministerio de Medio Ambiente de la Provincia de Santa Fe (Expte No. ­02,101-00, 181, 129-9, Resolution No. 093/2018). Guidelines of the American Society of Mammalogists (Sikes and the Animal Care and Use Committee of the American Society of Mammalogists 2016) were followed. The second group consisted of ethanol-preserved tissue samples, including muscle and patagium, obtained from museum vouchers deposited in the Colección de Vertebrados of the Instituto Nacional de Limnología (INALI-A), Colección de Mamíferos of the Museo Provincial de Ciencias Naturales “Florentino Ameghino” (MFA-ZV-M), Colección de Mamíferos of the Museo Provincial de Ciencias Naturales “Dr. Ángel Gallardo” (MG-ZV-M), and Colección de Mamíferos of the Museo Argentino de Ciencias Naturales “Bernardino Rivadavia” (MACN-Ma). The remaining samples including intestine, patagium, and muscle tissue were collected from bat specimens from the “Instituto de Zoonosis Luis Pasteur” (Buenos Aires city, Argentina) and personal collections.

One nuclear gene, the intron 7 and short partial sequences of exons 7 and 8 of the β-fibrinogen (FGB) and two mitochondrial genes, cytochrome c oxidase I (COI) and cytochrome b (cyt b), were sequenced totaling 2567 aligned base pairs (COI — 657 bp; cyt b — 1140 bp; FGB — 770 bp). Primers and PCR conditions for each gene are given in Table S2. Sequences were obtained with both primers for each locus. Samples that failed to produce sequences were discarded, remaining 47 successfully sequenced specimens (Table S1). Coding sequences were checked for the presence of premature stop codons to discard the amplification of pseudogenes.

To test the resolution of each marker we downloaded all available sequences annotated as Molossus in GenBank. Although it is immensely valuable, this public database is susceptible to having lots of errors in the taxonomy of the sequences available in it, especially in taxa with cryptic species and taxonomic instability such as Molossus. Morphological species misidentification and/or lack of taxonomic update are the two sources of error that abound in Molossus sequences in GenBank.

To overcome this issue we designed a strategy consisting in updating all previous names according to the current taxonomy, taking into account species distributions (Loureiro et al. 2020; Montani et al. 2021, 2023; IUCN 2022; Pavé et al. 2023; Olímpio et al. 2024). For example, M. rufus was recently subdivided into M. nigricans, M. rufus, and M. fluminensis, with non-overlapping distributions in Central America, northern South America, and southern South America, respectively (Fig. S1). Another example is that of M. milleri, which is distributed in the Cayman Islands and Cuba, and sequences uploaded previously to its description were annotated as M. molossus. A second part of this strategy was the recognition of misidentified specimens. This was done by scrutinizing majority clusters (composed of more than 50% sequences) of a given species in each single-locus phylogeny and revising those sequences that behaved as outliers in those clusters. Again, this revision was done by contemplating the updated taxonomy, identifying also which species have low levels of genetic differentiation, and also checking if the putative wrong assignment was done to a specimen captured within the distribution range of the majority cluster. Even after the correction of erroneously annotated sequences, we noted that several species failed to produce monophyletic clades, so we included sequences representing all these lineages in the multi-locus analysis.

In addition to the analysis conducted for each marker separately, phylogenies were inferred from the concatenated mitochondrial loci on one hand, and from the three combined loci on the other. The combined analysis of cyt b and COI genes allows the reconstruction of the same evolutionary history (Hurst and Jiggins 2005). Different phylogenetic patterns have been reported between mitochondrial and nuclear DNA (Spinks and Shaffer 2009), which highlights the need to compare both datasets independently (Loureiro et al. 2019). For the full dataset, we concatenated the three genes including, whenever possible, specimens with sequences available for the three loci. In cases when this was not achievable, we constructed chimeras ensuring these were formed by individuals of the same species, country, and lineage in the single-locus trees (Table S3). At least two individuals/chimeras per species and lineage were included to test their monophyly. Since the FGB analysis resulted in a less informative topology, we gave priority to specimens with sequences from both mitochondrial loci. If a specimen lacked FGB data, we completed the chimera with another specimen, following the previously mentioned criteria. The sequences of specimens/chimeras that lacked a specific locus were coded as missing data.

A Bayesian phylogenetic analysis was conducted for each locus and the concatenated dataset including representatives of all known species of Molossus. We used sequences of Promops centralis and Eumops auripendulus as outgroups. Nucleotide substitution models were estimated using MrModeltest2 (Nylander 2004) under the Akaike Information Criterion (correcting by the number of taxa), being HKY+G for FGB and COI, and HKY+I+G for cyt b. In the multi-locus analyses, we inferred a substitution model separating each codon position being K80+I, JC and HKY+G for 1^st, 2^nd, and 3^rd positions of COI, and SYM+G, F81, and GTR+I+G for cyt b.

The phylogenetic analysis was performed in MrBayes version 3.2.7 (Ronquist et al. 2012), using computational resources from CCAD – Universidad Nacional de Córdoba (https://ccad.unc.edu.ar), which are part of SNCAD – MinCyT, República Argentina. For each single-locus analysis and for the concatenated mitochondrial (mtDNA) dataset, two independent runs for 1×10⁷ MCMC (Markov chain Monte Carlo) generations, sampling every 1,000 generations, were carried out. The multi-locus analysis including nuclear and mitochondrial genes (FGB, COI, cyt b) was run for 2×10⁷ MCMC generations, sampling every 1,000 generations. Convergence was assessed by analyzing the potential scale reduction factor (PSRF), and the average standard deviation of split frequencies (ASDSF). The burnin phase was set up in the generation that fulfilled PSRF values of 1.00–1.02 for all estimated parameters and standard deviations lower than 0.01. We also checked that the estimated sample sizes were >200. Trees were visualized with Figtree (http://tree.bio.ed.ac.uk/software/figtree). For the comparison of phylogenetic resolution, we considered a posterior probability (bpp) of 0.95, 0.75 and <0.95, and <0.75 as high, moderate, and low node support.

A median-joining network was constructed with PopART v1.7 (Leigh and Bryant 2015). The software necessitates that sequences have few or no undefined states. Since the nuclear marker (FGB) was absent in several individuals/chimeras, we opted to conduct the analysis solely considering the mtDNA loci. Additionally, individuals lacking sequences for either of the two mtDNA loci were excluded. Following this procedure, all lineages were represented, except for Molossus sp. 2. For this lineage, we created a chimera using a cyt b sequence and a COI sequence from two different individuals (Molossus sp. 2 Ecuador-A1 and Molossus sp. 2 Bolivia-A7, respectively).

To compare intra- and interspecific distances we generated a p distances table with the concatenated matrix including the three loci. Pairwise genetic distances were calculated with MEGA X (Tamura et al. 2021). Of the two sequences of M. currentium, one was considered as M. fluminensis (MH185138), while the other was treated as M. currentium (MH185139), based on the results of the phylogenies obtained in this study (see below).

We recorded the following data of each specimen analyzed: locality, age class (i.e., subadult or adult), sex, external and cranial measurements, and body mass, following Barquez et al. (1999). External and cranial measurements were obtained for 64 adult specimens of Molossus species from Argentina, including six M. currentium, 12 M. fluminensis, 13 M. melini, 20 M. molossus, and 13 specimens of the novel Molossus species (Molossus sp. nov.). Measurements were taken using a digital caliper, supplemented in some cases by data obtained from the existing literature (Thomas 1901; Cláudio et al. 2020; Montani et al. 2021; 2023; Pavé et al. 2023; Olímpio et al. 2024). We included five external [total length (ToL), tail length (TL), hindfoot length (HFL), ear length (EL), forearm length (FA)] and 14 cranial measurements [greatest length of skull with incisors (GLS w.i), greatest length of skull excluding incisors (GLS wo.i), condylobasal length (CBL), postorbital constriction (PC), braincase breadth (BB), zygomatic breadth (ZB), mastoideal breadth (MB), maxillary toothrow length (LMxT), palatal length (PL), width across canines (C-C), width across second molars (M2-M2), mandible length (LM), mandibular toothrow length (LMdT), height of the sagittal crest (SAR)] for each species. In some cases, information was recorded directly from labels of museum specimens (File S1).

Univariate (using one-way ANOVA or Kruskal–Wallis test with corresponding pairwise post hoc tests) and multivariate statistical analyses (employing principal component analysis—PCA—coupled with PERMANOVA and its corresponding pairwise post hoc tests) of cranial and external linear measurements were conducted in R Studio (RStudio Team 2020) utilizing R version 4.2.1 (R Core Team 2021). The following R packages were utilized: FactoMineR version 2.8 (Lê et al. 2008), ggplot2 version 3.4.3 (Wickham 2016), Vegan version 2.6–4 (Oksanen et al. 2008) and pairwiseAdonis (downloaded from the author’s github site: pmartinezarbizu/pairwiseAdonis/pairwiseAdonis). For shape based comparisons, an exploratory two dimensional geometric morphometric (2D GMM) study of nasal cavity perimeter of the three species that showed the lesser morphological differentiation (M. molossus, M. melini, and Molossus sp. nov.) was performed. The tpsDig2 version 2.32 (https://sbmorphometrics.org) was employed for manual definition of 13 landmarks and one semi-landmark on picture files (*.jpg) of rostral views. The R based geomorph package version 4.0.6 (https://cran.r-project.org/web/packages/geomorph/index.html) was used for Generalized Procrustes Analysis (GPA) and PCA. Detailed listings of the datasets employed for univariate, multivariate, and 2D GMM morphometry, as well as the statistical methods applied and the resulting test outcomes, are provided in File S1.

The trees obtained for all nuclear and mitochondrial loci confirmed the monophyly of the genus Molossus. However, significant discrepancies were observed in internal nodes, particularly between FGB and both mitochondrial genes (Figs S2–S4; File S2). While the nuclear gene fails to recover most species as monophyletic, it can enhance the resolution of another dataset, such as the mtDNA genes (Fisher-Reid and Wiens 2011). The mtDNA dataset produced a similar but better-resolved phylogeny compared to cyt b, which is the better-resolved single-locus phylogeny (Fig. S5; Files S2, S3).

The phylogeny inferred from the concatenated dataset depicts higher levels of phylogenetic resolution, with the most robustly supported internal nodes (609 parsimony informative sites of a total of 2567) (Fig. 1). While there are minor discrepancies between the mtDNA tree and the multilocus tree, these are slight differences. The only significant disparity involves the positioning of sequences of M. aztecus, indicating potential introgression events (discussed in File S2). Consistent with the single-locus and mtDNA phylogenies, M. fentoni, M. verrilli, M. milleri, and M. alvarezi form highly-supported basally-diverging groups. In addition, M. rufus (with high support) and M. bondae (with low support) are also monophyletic groups.

The clade of M. fluminensis has moderate support (bpp = 0.84) but is still recovered as monophyletic, and splits into two well-supported subclades. Within this group, there are two previously annotated M. rufus that were updated according to their distribution to M. fluminensis (Fig. S1) and a sequence annotated as M. currentium from Paraguay, a region where the distribution of both M. currentium and M. fluminensis overlap. Notably, the other sequence of M. currentium (also from Paraguay) clusters with M. nigricans, M. sinaloae, and M. rufus. These three species do not overlap with M. currentium, so the most reasonable explanation is that this second haplotype represents the valid M. currentium and the former is, in fact, M. fluminensis.

The species M. pretiosus and M. aztecus from Mexico form a highly supported clade with M. bondae, but fail to be reciprocally monophyletic, indicating that these species are genetically close. Interestingly, the three representatives of M. aztecus from Brazil form a monophyletic clade, sister to the former, although there are two distinct lineages (A and B), as evidenced by the long branch within this group. In contrast, in the mtDNA phylogeny, M. aztecus from Brazil is divided in two non-related clades, one related to M. fluminensis, and the other related to M. melini (Fig. S5; File S2).

The clade of M. melini includes several specimens previously identified as M. molossus, not only from the originally described localities but extending also to the Buenos Aires province. Interestingly, several specimens that fall in the clade of M. melini, do not have the typical ochraceous to orange fur coloration, but instead are chocolate to grayish brown with the venter paler, and that could be the reason for their initial identification as M. molossus (Fig. S7). Sister to this clade is the M. molossus clade which includes the subspecies M. molossus daulensis. This clade groups two of the three M. molossus lineages found with COI: the canonical lineage, and a lineage composed of sequences from Panama and Brazil (Fig. S3; File S2). We conclude that these two lineages are indeed M. molossus. However, as occurs with the cyt b- and COI-based phylogenies, there is another clade of specimens identified as M. molossus with moderate support (bpp = 0.81), that splits into well/moderately supported clades, one (bpp = 0.81) including specimens from Ecuador, Peru, Guyana, Suriname, and Bolivia, while the other is exclusive of Argentinian bats (bpp = 0.97). This last group should be treated as Molossus spp., consisting of Molossus sp. nov. (Argentinian clade) and Molossus sp. 2 (sister clade), since it is reciprocally monophyletic to the canonical lineage of M. molossus (which has a wider distribution and includes the subspecies M. molossus daulensis). It is noteworthy that, except for one sample, Molossus sp. 2 formed part of the M. molossus clade in the mtDNA analysis (Fig. S5) and in the full dataset it groups with Molossus sp. nov.

The genetic distance analysis illustrates the general panorama of Molossus (Fig. 1; Tables S4, S5). Intraspecific distances depict low values, ranging from 0 to 1.7% (mean 0.85%), the highest values belonging to M. molossus, M. pretiosus, M. coibensis, and M. fluminensis. The only exception is M. aztecus, depicting 2.5%, a remarkably high value for intraspecific variation (but see an extended discussion in File S2). As expected for this genus, there is a wide variation of interspecific distances, with a mean of 4.0%, but including extremely low values such as those found between M. pretiosus, M. nigricans, M. bondae, M. rufus, M. sinaloae, and M. currentium, all of which are <1.0%. The distance between M. molossus and Molossus sp. nov. is 2.8% which is indicative of interspecific variation. In contrast, the distances between M. molossus and Molossus sp. 2, and between Molossus sp. nov. and Molossus sp. 2 are 1.1 and 1.4, respectively, and will be discussed below. As a whole, Molossus appears to be a genetically stable group. The average genetic distance among morphologically recognized species using mtDNA markers was estimated to be approximately 3.3% (Loureiro et al., 2019), a finding consistent with the results of this study. These low levels of genetic variation are comparable to certain molossid bats from the Old World but are notably lower compared to others such as Cynomops, Mops, and Eumops, where cyt b average distances range from 8.6% to 13.4% (Ratrimomanarivo et al. 2007; Baker et al. 2009a; Moras et al. 2018). While interspecific distances in our study are notably lower compared to the >5% threshold suggested by Baker and Bradley (2006) for species delimitation, it is crucial to recognize that such thresholds serve as guidelines and can vary across different organisms (Baker et al. 2009b). These thresholds are not universally applicable and should be determined empirically for each taxonomic group, taking into account specific characteristics such as generation time, dispersal rates, and geographic range (Hebert et al. 2003). Finally, we believe that genetic distances should prompt taxonomic reassessments but should not be considered the sole criterion for species delimitation.

The mtDNA median-joining network illustrates significant differences among M. melini, M. molossus, and Molossus sp. nov. (Fig. S6). Molossus molossus exhibits more than 30 substitutions distinguishing it from the other two species, while M. melini and Molossus sp. nov. display twice this distance, indicative of their classification as distinct species. It is worth noting that certain taxa, such as M. fluminensis, M. coibensis, or the subspecies M. molossus daulensis, exhibit high levels of intraspecific divergence at the mitochondrial level. The chimera constructed using two sequences representing Molossus sp. 2 is closely associated with Molossus sp. nov. However, as will be discussed further, the fragmentary availability of sequences representing Molossus sp. 2 dictates a cautious approach in its treatment.

In order to contrast the above inferred phylogeny with morphological data, we undertook morphometric comparisons using a subset of genetically identified specimens from which we have enough confident linear metric external and cranial measurements. Molossus sp. nov. clearly differentiates from M. fluminensis by consistently exhibiting lower linear measurement values across all external and cranial variables evaluated (one way ­ANOVA, Fig. 2). However, Molossus sp. nov. tends to overlap with the other small species, being M. molossus the more distinguishable specially when considering the skull measurements (Fig. 2; File S1 [Table SM2] for the analysis with genetically defined specimens plus misidentified genetically specimens). Furthermore, between Molossus sp. nov. and the recently described M. melini only three out of fourteen linear metric cranial variables evaluated were statistically different (p<0.05): LMxT, PC, and LMdT (Fig. 2).

Global skull size variation between species, analyzed by PCA using a subset of cranial measurements of genetically confirmed specimens of M. fluminensis, M. melini, M. molossus, and Molossus sp. nov. (see File S1 for specimen details), showed that the convex hull enclosing the specimens of Molossus sp. nov. maps separated from M. molossus (p<0.001) and M. fluminensis (p<0.003) clusters and in a lesser extent with M. melini (p<0.05) (Fig. 3; File S1 [Table SM4]). In this analysis, the first two components summarize 93% of the variation. PC2 explained only 5.3% of the variation, being clearly dominated by PC and GLS wo.i variables (Fig. 3B). The PC1 explained the largest amount of variation because most cranial variables contribute to it with almost equal strength (Fig. 3B). Accordingly, the centroids of species groups are distributed along PC1 from the smallest (M. molossus, M. melini, and Molossus sp. nov.) on the left to the largest (M. fluminensis) on the right side (Fig. 3A).

Based on the genetic and morphometric evidence shown above, we describe the members of the lineage Molossus sp. nov. as a new species, as follows:

Family Molossidae Gervais, 1856

Genus Molossus É. Geoffroy, 1805