A total of 92 T. brasiliensis individuals and three Molossus molossus individuals (used as outgroup in some of our analyses) were obtained from passive rabies surveillance carried out by the Argentine National Administration of Laboratories (ANLIS), Servicio Nacional de Sanidad y Calidad Agroalimentaria (SENASA), Institute of Zoonosis Luis Pasteur (Autonomous City of Buenos Aires), Zoonosis Urbanas of Buenos Aires, Institute of Zoonosis of Córdoba, and from scientific collections, namely: Colección de Mamíferos Lillo, University of Tucumán (CML), the Mammal Collection of the Argentine Museum of Natural Sciences Bernardino Rivadavia (MACN), and the Mammal Collection of the Instituto Nacional de Limnología (INALI-MA). The sampling included individuals from 13 provinces of Argentina divided in six regions (Fig. 1; Table SS1): AMBA (Área Metropolitana de Buenos Aires), Norte (Jujuy, Salta and Tucumán), Mesopotamia (Entre Ríos and Santa Fe), Centro (Buenos Aires and Córdoba), Cuyo (La Rioja), Patagonia (Chubut, Neuquén, Río Negro and Tierra del Fuego). The individuals were fortuitously captured in both residential and periurban areas during all seasons of the years 2016, 2017 and 2019–2021. Species-level identification of the individuals was performed morphologically using taxonomic keys (Díaz et al. 2021b).

DNA was extracted from both the wing membrane (plagiopatagium) and muscle tissue samples of bats using the High Pure PCR Template Preparation Kit (Roche), following the manufacturer’s instructions. DNA purity was evaluated using a Nanodrop spectrophotometer (ThermoFisher) and DNA concentration was assessed with a Qubit fluorometer (ThermoFisher).

To obtain nuclear SNP data for population analysis, we employed double-digest restriction site-associated DNA sequencing (ddRADseq) technology using SphI and EcoRI endonuclease enzymes. Library preparation, paired-end sequencing (with a NovaSeq 6000), and basic bioinformatic analysis were performed by IGATech (Italy) following Peterson et al. (2012). The latter included applying sequence quality filters and demultiplexing reads with the process radtags tool included in Stacks v2.53 (Catchen et al. 2013), and assembly by mapping to a reference genome. The reference used was the Tadarida brasiliensis (TadBra_v1_BIUU) genome [GenBank: GCA_004025005.1], the only T. ­brasiliensis whole genome sequence available at the time, which was very fragmented (with 24.3k scaffolds and average mean scaffold length of 25,995 bp). Reads were mapped with BWA-MEM (Li and Durbin 2009), with a mapping quality threshold of > 4. The gstacks program included in Stacks v2.53 package (Catchen et al. 2013) was used for loci assembly and SNP calling. We then applied the following filters to the SNPs obtained with gstacks using the populations command (Catchen et al. 2013): R = 0.8 (to retain only loci that are represented in at least 80% of the samples), one SNP per locus, a min-maf = 0.05 (to retain SNPs whose least common allele had a frequency of at least 5% in the samples), max-obs-het = 0.80 (to process a nucleotide site at a locus with an observed heterozygosity of at most 80%), and a mind = 0.2 (to filter out individuals with more than 20% missing genotypes). Exact tests for differences between observed (H_O) and expected (H_E) heterozygosity, based on Hardy-Weinberg Equilibrium (HWE), and likelihood ratio tests for genotypic linkage disequilibrium (LD) were carried out using the PLINK v.1.09 (Chang et al. 2015).

We assessed the genetic structure of T. brasiliensis in Argentina based on the ddRadseq SNP data using two multivariate clustering methods: Principal Components Analysis (PCA), an unsupervised method that identifies components of the total variation that maximize the global variance in the dataset, and Discriminant Analysis of Principal Components (DAPC), which requires prior assignment of individuals into groups and maximizes differences between them. We also performed an Analysis of Molecular Variation (AMOVA) to estimate the proportion of the total variation found within and between predefined groups, which is a measure of population subdivision. To delimit groups for DAPC and AMOVA, we used, alternatively, season and geographic regions of collection. The PCA was carried out with PLINK v1.9 (Chang et al. 2015), while the other analyses were performed using the adegenet (Jombart and Ahmed 2011), vcfR (Knaus and Grünwald 2017) and vegan v.2.8 (Oksanen et al. 2025) R packages. Results were summarised in plots obtained with the ggplot2 (Kassambara 2020) package for R. These analyses, along with the others in this study that employed R packages, were run using R v.4.3 in the RStudio environment (R Core Team 2021).

Additionally, we carried out a model-based clustering analysis with ADMIXTURE v2.3.4 (Alexander et al. 2009), which employs a maximum likelihood approach to determine individual ancestries without prior group assignments. This method uses a cross-validation procedure that enables the user to identify the value of K (total number of ancestral populations) for which the model has greatest predictive accuracy, as determined by “holding out” data points (Do et al. 2014). The figures that summarize the results were created using the ggplot2 (Kassambara 2020) package. In a final effort to identify groups within the T. brasiliensis sample, we reconstructed a phylogenetic tree based on the neighbor-joining method and Kimura 2-parameter (K2P) distances, employing the three genotyped individuals of Molossus molossus as outgroups. To obtain the tree and bootstrap support values for its nodes, we used the adegenet (Jombart and Ahmed 2011) and ape (Paradis et al. 2002) R packages. The geographic distribution of genetic variation was evaluated using Mantel’s test (Mantel 1967), which tests for a correlation between the pairwise K2P distances and the log-transformed linear geographic distances between sampling sites (Smouse et al. 1986).

Genetic diversity in the SNP data was assessed with the number of polymorphic loci, number of alleles, expected and observed heterozygosity, nucleotide diversity (π) and inbreeding coefficients (F_IS) calculated with vcftools (Danecek et al. 2011). A kinship analysis was carried out with COLONY v.2.0.6 (Jones and Wang 2010). Contemporary effective population size (N_e) was estimated from multi-locus diploid genotypes using the excess-heterozygotes method implemented in NeEstimator v.2.1, assuming random mating (Do et al. 2014).

The mitochondrial gene cytochrome b (cyt b) was amplified from six of our samples of T. brasiliensis (from Entre Ríos, Jujuy, La Rioja, Río Negro, Santa Fe, and Tucumán; Fig. 1) using published primers Bat05A, Bat14A, and Bat-Ep (Martins et al. 2007; Caraballo et al. 2020). Cytochrome b sequences were obtained with an Applied Biosystems 3500 Series Genetic Analyzer according to standard protocols and were edited with Unipro UGENE 51.0 (Okonechnikov et al. 2012). Additionally, we retrieved from GenBank 52 sequences of T. brasiliensis, 31 of T. teniotis (Rafinesque, 1814), and 18 of T. latouchei (Thomas, 1920) (Table SS2). Among the T. brasiliensis sequences, 20 were from South America (Chile and Brazil), 17 were from the Bahamas, and 32 were from the USA. A multiple sequence alignment was obtained with MEGA X (Kumar et al. 2018).

Estimates of inter- and intraspecific genetic distances were obtained using the K2P distance model. We performed maximum likelihood (ML) tree searches with IQTREE v.3.2 (Minh et al. 2020). For that, the alignment was partitioned into codon positions, and the preferred evolutionary model for each partition (TNe+I, F81+F, and TIM2+F+G4) was determined prior to tree searches using the same program. In addition, we also implemented tree searches and dating under a Bayesian phylogenetic approach using BEAST v.2.5 (Bouckaert et al. 2019). A Markov Chain Monte Carlo (MCMC) analysis was run for 10 million generations sampling every 2000 generations. Tracer v1.6 (Rambaut et al. 2018) was used to determine that the run had reached stationarity and that a 10% burn-in was appropriate. All trees were visualized using FIGTREE v.1.4 (Rambaut 2018). In this analysis, we applied a strict clock with a substitution rate of 2% per million years (Nabholz et al. 2009) and a coalescent exponential population model tree prior. Additionally, we carried out an Analysis of Molecular Variation (AMOVA) to assess the distribution of the variance in the cyt b gene within and among T. brasiliensis populations/clades and calculated pairwise overall fixation indexes (F_ST). These analyses were done with the R package vegan v. 2.8 (Oksanen et al. 2025), using 10,000 permutations to evaluate the statistical significance of the results. Finally, to clarify relationships among mitochondrial haplotypes within T. brasiliensis, we constructed a haplotype network applying the TCS method, which is based on Templeton’s statistical parsimony, as implemented in the program PopART (Templeton et al. 1992; Clement et al. 2002; Leigh and Bryant 2015).

We investigated demographic history using both the nuclear genome SNPs and the cyt b gene. To analyze SNP data, we first inferred the site frequency spectrum using the vcfR (Knaus and Grünwald 2017) package and then used fastsimcoal2 (Excoffier et al. 2013) to estimate the likelihood of the observed site frequency spectrum under two demographic models. The first model, “constant population size”, included effective population size N_e as the only estimated parameter, while the second model, “constant-rate population growth or decline”, estimated the current population size, the population size prior to growth or decline, and the time at which growth or decline occurred. N_e parameters were selected from uniform distributions of 10–10⁷, whereas the prior for the event time was set between 1000 and 300,000. For each model, we performed 100 independent runs of one million simulations each and 30 expectation-maximization (ECM) cycles to estimate the parameters. The substitution rate was set to 0.02 substitutions per site per million years and point estimates in generations were converted to absolute time (years) assuming a generation time of two years, based on the previous studies (Russell et al. 2011). To compare the fit of the two models, we chose the run with the highest likelihood for each model and contrasted their likelihoods with the Akaike’s information criterion (AIC). Finally, demographic trajectories were plotted using ggplot2, representing the median and the 95% confidence intervals of the estimated effective population sizes over time.

Second, we examined demographic trends in T. brasi­liensis based on the cyt b sequence data using two alternative approaches: neutrality tests and Bayesian Skyline Plots. The neutrality tests Fu and Li’s F, and Tajima’s D were implemented in the software DnaSP v.5.0 (Librado and Rozas 2009). Bayesian Skyline Plots (BSP) were obtained with BEAST v.2.5 (Bouckaert et al. 2019), using parameters established in previous bat studies (Russell et al. 2015). The cyt b matrix was separated into codon positions using unlinked substitution models, which were the same as in the phylogenetic analysis. A strict molecular clock was applied, using the aforementioned substitution rate. We repeated the analysis using an alternative rate of 0.046 substitutions per site per million years, which has also been proposed for the cyt b gene in bats (Ruedi and Mayer 2001; Baird et al. 2017). We used default values for the operators and set the group size to 3 and population size to 10,000. We ran an MCMC for 20 million generations, sampling every 1000 generations, and discarding the first 20% as burn-in. The results were checked for convergence by monitoring ESS values in Tracer v1.6 (Rambaut et al. 2018). Finally, the estimates of N_f (the effective size of the female population) were derived from the BSP results assuming a generation time of two years (Russell et al. 2011). The BSP graphics were created with the R packages devtools (Wickham et al. 2022) and ggplot2 (Kassambara 2020).

Approximately 1 billion reads were obtained, with an average of ~10.8 ± 1.6 million reads per sample. The mean coverage per sample was 6.7x and the average saturation at 6x was 99.5%. The assembly generated a total of 41,046 RAD loci, of which 40,976 were polymorphic, containing 660,413 polymorphic sites. After applying the filters with gstacks, 30,679 variable sites were retained, while the three M. molossus samples and one T. b. brasiliensis were removed from the dataset due to excess of missing genotypes. The mind filter was ignored to obtain a SNP table including the outgroup specimens for the Neighbor-Joining analysis. Finally, after applying the LD filter, our working dataset consisted of 29,715 unlinked SNPs. The transition-transversion rate (Ts/Tv = 2.4) was within the expected values given the random nature of ddRADseq SNPs.

Analyses of population structure with ADMIXTURE showed that the Argentine samples of T. brasiliensis adjusted better to a K=1 (with 0.52 of cross validation error and Loglikelihood mean = -1952762.7), suggesting the existence of a single cluster, i.e., a single population widely distributed in the Argentine territory (Fig. S1). The first two principal components of the PCA, which together explained over 10.4% of the total molecular variance, did not separate samples into distinct groups, confirming the absence of genetic clustering by geographic origin or collection season (Fig. 2). Accordingly, the AMOVA failed to identify significant genetic differentiation between seasons (Φ_ST = 0.0001, p = 0.439) or between geographic regions (Φ_ST = 0.0001, p = 0.346). In both treatments, 99.99% of the total variance was found within populations. This lack of structure suggests the absence of major barriers to dispersal across the studied localities, allowing for unrestricted gene flow. These findings are consistent with the Neighbor-Joining clustering analysis, which did not identify statistically significant clades within Argentina (bootstrap < 50%; Fig. S2). The Mantel test revealed lack of spatial genetic structure, as genetic and geographic distances were uncorrelated (r = –0.045, p = 0.78; Fig. S3)

The average nucleotide diversity of the SNP dataset was 0.24 (SD ± 0.13) and the average expected heterozygosity (H_e) was 0.192 (SD ± 0.12). Random mating throughout the entire Argentine population is evidenced by the low values of the inbreeding coefficient (F_IS= 0.17, SD ± 0.029). The kinship analysis resulted in 92 clusters, indicating no consanguinity among the whole sample. Estimates of N_e based on heterozygote excess were very large: 2 trillion individuals with Colony, “infinite” with NeEstimator and ~2.7 million individuals with fastsimcoal2.

The phylogenetic trees based on the cyt b gene, obtained via ML and BI analyses, generally coincided in topology and recovered three monophyletic clades corresponding to T. brasiliensis, T. teniotis, and T. latouchei (Fig. S4). The T. brasiliensis samples were clustered into three main, well-supported clades. The first to split off was a highly supported North American clade (NA), composed of individuals from the USA (Florida and Arizona) and the Little Bahama Bank (Grand Bahama and Abaco islands). The second clade (GBB) included all individuals from Eleuthera and Long Island, both of which are part of the Great Bahama Bank, as well as one individual from the Grand Bahama Island (likely a recent migrant, see Speer et al. 2017). The third clade (SA) included all individuals from South America (Fig. 3A, B). The analysis revealed a sister relationship between the latter two clades, but with relatively low bootstrap support (74%).

Genetic distance between New and Old World Tadarida species ranged between 15% and 17.7% (Table 1). At the interpopulational level within T. brasiliensis (i.e., between the NA, SA, and GBB clades), the distances ranged between 2.8% and 6.1%, while the intraclade distances were less than or equal to than 1%.

The AMOVA results revealed significant genetic differentiation among the three T. brasiliensis clades (F = 111.2, p < 0.0001), with 74.8% of the total genetic variation attributed to differences among populations, and around 25.2% of the variation occurring within populations. The nucleotide diversity corresponding to the mitochondrial gene cyt b for the whole T. brasiliensis dataset was 0.036. Taking each T. brasiliensis main clade separately, both the nucleotide and haplotype diversity estimates for the NA clade, and the SA clade were similar to one another and higher than the estimates obtained for the GBB clade. Haplotype diversity and segregating sites were high (> 0.5) while nucleotide diversity was low (< 0.5) in both continental populations (Table 2).

The network of T. brasiliensis haplotypes mirrored the results of the phylogenetic analysis revealing the existence of three distinct mitochondrial haplogroups within the species, corresponding to the clades observed in the phylogenetic tree (Fig. 3C). The haplogroup including all 27 individuals from South America was represented by a total of 13 different haplotypes. The seven individuals collected in Argentina (Table SS2), shared 2 haplotypes with individuals from Chile and Brazil, and the remaining 4 haplotypes differed in at most four substitutions from other South American haplotypes. A second cluster, differing in fifteen substitutions from the SA group, included the 17 individuals from the GBB clade and exhibited only two haplotypes distinguished by a single thymine to cytosine substitution. The third haplogroup encompassed 20 haplotypes found in individuals from the USA and islands of the Little Bahama Bank. This cluster, corresponding to the NA clade, differed from the SA cluster by 30 substitutions.

The divergence dating analysis using the 2% per million years substitution rate indicates that the first split in the species occurred around 2.1 million years ago (mya) (95% confidence interval: 1.49–2.72), separating the NA population from the other two analyzed populations, which diverged from each other around 0.9 mya. Interestingly, both continental populations coalesced concomitantly at approximately 0.27 mya, while the GBB population had a more recent common ancestor (Fig. 4A). The divergence dates were 0.9 and 0.4 mya when the alternative substitution rate of 4.6% was applied instead (Fig. S6A).

The model comparison showed a better fit of the observed SNP data to the decline/growth model (AIC = 59,687), effectively capturing the observed excess of low-frequency variants (Fig. S5). Demographic parameter estimates under this model indicate a significant populational expansion in the Southern Cone populations of T. brasiliensis. The median ancestral effective population size was estimated at approximately 306,000 individuals (95% CI 295,467–313,600), before a nearly 9-fold increase to a current N_e of ca. 2,700,000 individuals (2,566,452–2,789,485; Fig. S5). Based on a generation time of one generation every two years, the estimated onset of the expansion occurred approximately 270,000 years ago (0.27 mya; 265,678–277,014), a time similar to the estimated age of the SA population based on the mitochondrial locus cyt b. Different demographic histories were revealed by the cyt b variation when neutrality tests were applied to the NA, SA and GBB populations (Table 2). The continental populations (NA and SA) showed an excess of rare haplotypes (negative values in the Tajima’s D and Fu and Li’s F statistics), although only the SA population had a significant value of Tajima’s D. In turn, the population of the GBB had non-significant results for the three neutrality tests, suggesting that the population is evolving in a manner consistent with neutral expectations.

Demographic history differences between populations were also revealed by the Bayesian skyline plots (Fig. 4). According to the BSP, the demography of the SA population is uncertain due to the large confidence intervals obtained in the estimates of the effective number of females (N_f; Fig. 4B). However, the mean and median N_f estimates suggested a fivefold population expansion between approximately 160,000 and 100,000 years ago (kya). This expansion slowed down after that, reaching a plateau around 75,000 years ago and the population has remained relatively stable since then. This result contrasts with the one obtained for the NA population, for which the evidence of an expansion is stronger. In this population, the 6-fold expansion would have been more recent, during the last 110,000 years with a continuing growth trend until present (Fig. 4C). Both continental populations had similar coalescence times (~220,000 years), which was slightly more recent than the estimates obtained in the divergence dating analysis. Conversely, the population of the GBB appears to have maintained a more stable N_f of approximately 10,000 individuals since its coalescence event that was dated to only 15,000 years ago (0.015 mya), without evidence of population expansion or reduction (Fig. 4D). All those results were based on the lower mutation rate (2% per million years), while the same analysis employing a higher rate (4.6% per million years) revealed similar patterns but with proportionally (about half) younger ages (Fig. S6).

Our analysis of T. brasiliensis from Argentina based on ddRAD SNPs did not detect significant population structure across the studied area. Although the estimated genetic diversity was high across the country, the genetic distance between samples was not correlated with geographic distances. These results suggest that the samples included in our analyses, collected in an area covering around 4800 km in a north-to-south axis and 1800 km from east to west, constitute a single panmictic population. This implies that there are high levels of gene flow throughout the Argentine territory. The same conclusions were reached by analyzing the mitochondrial gene cyt b, which showed the existence of shared haplotypes between individuals collected thousands of kilometers apart, expanding the distribution of the panmictic South American population through Chile and at least part of the Brazilian territory. A similar pattern had been observed in North American populations, where analyses of different molecular markers failed to detect structure, despite apparent differences in migratory behaviour (Russell et al. 2005). These findings indicate a high degree of genetic connectivity across populations, reinforcing the role of gene flow in shaping population dynamics as expected given the high dispersal capability of T. brasiliensis reflecting their impressive flying abilities (McCracken et al. 2008; McCracken et al. 2016). It has been shown that T. brasiliensis is capable of flying distances exceeding 160 km during nocturnal foraging trips, returning to its roost within the same night, and maintaining horizontal flight speeds of up to 44 m/s (McCracken et al. 2016).

Our results are also in agreement with findings of previous studies based on microsatellites, which revealed lack of isolation by distance within regional populations of T. brasiliensis (Speer et al. 2017; Morales et al. 2018). A significant correlation between geographic and genetic distances, however, was found in a continental scale analysis that included samples spanning North to South America (Morales et al. 2018). These results imply that the great diversity of ecoregions, thermal floors and topography found in both South and North America do not represent barriers to dispersal and gene flow in T. brasi­liensis, similar to patterns observed in other insectivorous bat species with migratory behaviour (Baird et al. 2017; Chipps et al. 2020; Loureiro et al. 2020). Nevertheless, ocean currents between the Northwest and Northeast Providence Channels appear to affect dispersal since the GBB population is markedly different from the NA population despite the short distances that separate the Great Bahamas Bank from the Little Bahamas Bank, the latter forming part of the NA’s population range (Speer et al. 2017). Migratory behavior in T. brasiliensis is probably triggered by temperature, alongside the seasonal availability of food resources, such as the North American moth migration in autumn, a key resource to rebuild their own pre-migratory fat reserves (Krauel et al. 2018). However, differences in migratory behavior have been observed across geographic regions and between sexes. Banding experiments suggested the presence of sedentary or non-migratory subpopulations in parts of the United States (Cockrum 1969), yet no underlying genetic differentiation has been identified to account for these migratory phenotypes (Russell et al. 2005). Furthermore, it appears that females travel longer distances, whereas males more often form resident populations during winter or do not migrate as far south into Mexico (Russell et al. 2005). In South America, similar patterns of disparate migratory behaviour between sexes have been observed in Uruguayan populations (Botto Nuñez et al. 2018).

Studies of large maternal colonies in Argentina, Brazil, and Uruguay showed that the annual cycle of T. b. brasiliensis female activity is similar to that reported in North America: Arrival to the shelters in early- to mid-spring, births in late spring (November – December), and gradual migration in late summer or early autumn (March - May) until the shelters are empty or just inhabited by a small group of individuals that remain during the cold seasons (Romano et al. 1999; Botto Nuñez et al. 2018). There is no information about migratory routes in South America, although it has been suggested that individuals migrate towards the tropics when they leave the maternal colonies (Díaz et al. 2019), the same pattern observed in North America where T. b. mexicana migrate as far as 1800 km to lower latitudes (Villa and Cockrum 1962; Cockrum 1969; Glass 1982). Our results showed an agreement between nuclear and mitochondrial markers in indicating a lack of genetic structure, suggesting similar migratory behaviour in males and females, and absence of female philopatry. Long range migrations of T. brasiliensis colonies potentially affect the circulation of pathogens and play a key role in the epidemiology of zoonotic diseases. Continuous gene flow between populations can facilitate the spread of diseases, either by introducing new pathogens or by modifying the dynamics of those already present in a given environment. In particular, the high connectivity of T. brasiliensis populations in Argentina points to its potential to transport and spread zoonotic viruses such as rabies lyssavirus and alphacoronaviruses across the South American continent (Caraballo et al. 2022, 2024).

Our results agree with previous molecular studies where microsatellite markers showed that in both resident and migratory T. brasiliensis populations of North America, the Caribbean, and South America have substantial genetic variability and heterozygosity (Russell et al. 2005; Morales et al. 2018). The uncovered genetic variation we observed in our dataset led to very high estimates of effective population sizes. We presented the results considering a generation time of two years, as reported in previous studies (Russell et al. 2011). However, ecological data from different colonies indicate that females mature at 9 months and births occur annually, which suggests that the generation time could be 1 year, and in this case, N_e estimates would be twice as large than the ones presented herein (Romano et al. 1999; Díaz et al. 2019).

The elevated N_e estimated for the Argentine population align with those reported from census monitoring of colonies within the country and neighboring regions. For instance, the largest documented colony, located in the province of Tucumán, Argentina, had been estimated to consist of tens of millions of individuals (Díaz et al. 2021a). This colony was significantly reduced in size after the closing of some of the tunnels it used as roosting sites. Other known large colonies in Argentina include one at the Facultad de Derecho, Universidad Nacional de Rosario building (Santa Fe, Argentina), originally comprising approximately 64,000 bats roosting in its roof (Romano et al. 1999), Ex-Bodega Chilecito in La Rioja with ca. 2000 specimens (Gómez et al. 2024), and the La Calera colony in Córdoba, with an estimated 900,000 individuals roosting in abandoned mine tunnels, ranking as the second largest T. brasiliensis colony in South America (Boero et al. 2020). Additionally, in Porto Alegre (Rio Grande do Sul state, Brazil) and several regions of Uruguay, colonies of between 1000 and 10,000 individuals have been recorded in various buildings (Botto Nuñez et al. 2018).

Our results confirm previous studies showing differentiation at the cyt b mitochondrial gene between continental populations of T. brasiliensis, revealing three main clades among the analyzed sequences: the GBB, SA, and NA populations (Speer et al. 2017). No haplotypes were shared among those three populations that were separated by significant genetic distances. A previous study based on SNPs found genetic structuring among these populations, accompanied by cranial morphometric differences (Morales et al. 2018). We estimated the age of the first split within what we now accept as T. brasiliensis in around 2 million years, when the NA population started to diverge from the ancestor of the GBB and SA populations. The sister relationship between the latter populations is in disagreement with results from a previous study (Speer et al. 2017). In contrast, we found no phylogenetic resolution within each of those three main groups, in accordance with previous studies based on both mitochondrial and nuclear markers, which showed that the different subspecies proposed for North America are not supported by molecular data (Russell and McCracken 2006; Morales et al. 2018).

The genetic distances in the cyt b locus estimated between T. brasiliensis and its congeners from the Old World (15–1 7.7%), and between the different populations of T. brasiliensis were unexpectedly high (2.8–6.1%). Previous studies have shown that intraspecific distances in the cyt b gene are generally lower than 2.5% in bats, while distances between species of the same genus range from ~3.5 to16% (Bradley and Baker 2001; Chambi Velasquez et al. 2024). The distance estimates obtained herein, together with the strict monophyly of the involved species and populations, and the AMOVA results for the T. brasiliensis populations, indicate the need for a taxonomic revision of the genus. This review should include more comprehensive sampling of populations and species as well as complementary sources of data, such as nuclear locus sequences and morphological data.

Genetic diversity in the cyt b locus was similar in the SA and the NA populations, while the GBB population was significantly less diverse, as expected for insular populations due to the isolation and limited territory of islands in comparison to continents (Frankham 1998). Interestingly, nucleotide diversity in the continental populations of T. brasiliensis was higher than that found in T. teniotis, a species with a similar biology and flying abilities that inhabits the Old World (Amorim et al. 2019). This difference was reflected in differences in the demographic histories recovered for the two species (see below).

The demographic analyses based on SNP data revealed an expansion of the T. brasiliensis population in the Southern Cone. Contrary to the expectation of a late post-Pleistocene expansion, the coalescent-based analysis — supported by low confidence intervals across 100 independent replicates — places the onset of the demographic expansion at approximately 0.27 mya. This period coincides with the end of the MIS 8 glacial cycle, suggesting that the species underwent massive population growth during the Middle to Late Pleistocene, reaching a current N_e of approximately 2.7 million individuals (Lisiecki and Raymo 2005; Otvos 2015). This result suggests that posterior climate oscillations did not have an important impact on the demography of T. brasiliensis in South America (although we did not include in the model comparison more complex models with several demographic events). Such resilience could be attributed to the species’ high dispersal ability and generalist nature that allowed it to buffer the impacts of glaciations. Nevertheless, we did not test for more complex models of population growth.

The results of the neutrality tests applied to the cyt b matrix further support population expansions in both T. brasiliensis continental populations. The Bayesian Skyline Plots (BSP) are also in accordance with those results, although for the SA population, the large confidence intervals in the N_f estimates makes the case for expansion less robust than for the NA population. The BSP provided time frames that placed the onset of the expansions around 0.12 mya and 0.17 mya for the NA and SA populations, respectively. Our results support the conclusions of Russell et al. (2011), suggesting that population expansions in continental T. brasiliensis were most likely driven by environmental causes - in particular climate changes caused by the cyclic Pleistocene glaciations. As generalist species, the main limitation to the distribution of T. brasiliensis during glacial periods were likely temperature and the associated availability of insect prey. Therefore, populations were probably reduced and restricted to habitable tropical and subtropical refugia. During interglacial periods, populations previously confined to refugia were able to colonize higher latitude areas of North and South America, resulting in a substantial range expansion and concomitant population growth (Russell et al. 2011). The migratory behavior and the great dispersal ability of T. brasiliensis facilitated rapid colonization of temperate areas when conditions permitted (Russell and McCracken 2006).

The onset of the expansions of T. brasiliensis populations, according to the BSP, roughly coincide with the end of the long and intense MIS6 glaciation which ended circa 0.130 mya (Lisiecki and Raymo 2005). It was followed by the Last Interglacial (LIG) period (MIS5e), a very warm period, and glaciation cycles that were shorter and not as severe as in the MIS6 period, with temperatures slowly dropping until the Last Glacial Maximum at around 0.2 mya (Lisiecki and Raymo 2005; Otvos 2015). Although fluctuations in environmental conditions during the Pleistocene glaciations affected the entire world, there is a notable difference between North America and South America: South America has a much larger tropical area as compared to North America, which could sustain larger populations of generalist species such as T. brasiliensis during the glaciation periods. This may explain the difference in the timing of population expansion in these two continents, which appears to have occurred later in North America. Unfortunately, little is known about the environmental conditions of most of South America during the Pleistocene (Gowan et al. 2021). Although a large portion of the South American continent likely did not experience abrupt fluctuations in temperature, part of today’s southern range of T. brasiliensis, the Chilean and Argentine Patagonia, was covered in snow in glacial periods (Lisiecki and Raymo 2005). In contrast, the Salto Ander Egg formation in northeastern Argentina (dated to ~0.12 mya) provides evidence of a warmer LIG climate, as indicated by a tropical fauna comprising species now restricted to lower latitudes (Ferrero et al. 2017).

A chronological discrepancy emerges when results obtained for the SA population with the two types of markers are compared, with younger expansion dates being inferred by mitochondrial DNA (~0.170–0.12 mya) in contrast with older dates obtained with nuclear SNPs (~0.27 mya). Such discrepancies are technically to be expected in phylogeographic studies, because the effective size of the mitochondrial genome is around a quarter of that of the nuclear genome, which accelerates coalescence and tends to reflect more recent demographic shifts (Moore 1995; Palumbi et al. 2001). In this sense, cyt b recovered population growth following the MIS6 glaciation and leading into the Last Interglacial, whereas SNPs preserved a deeper signal of the basal expansion occurring in the Middle Pleistocene (MIS 8). Furthermore, the higher resolution of the genomic data enabled us to reduce the statistical uncertainty of parameters estimated with the mitochondrial marker, and suggests that the resilience of T. brasiliensis has enabled it to endure multiple glacial cycles as a large, well-connected population in South America.

Interestingly, the results obtained for T. brasiliensis differ from those obtained for the European free-tailed bat, which appear to have expanded its population much more recently. Amorim et al. (2019) studied a sample of 120 individuals of T. teniotis collected across Europe. They analyzed microsatellite and mitochondrial loci using the ABC method and found that a population expansion occurred after the Last Glacial Maximum, which ended around 20,000 years ago. This expansion occurred apparently from two separate source populations that survived in different refugia, a common pattern found among European species (Amorim et al. 2019).

Finally, the potential role of anthropogenic factors in the demography of T. brasiliensis cannot be entirely dismissed. Although it did not promote the onset of population expansion, it may be contributing to the maintenance of current populations and facilitating the dispersal and establishment of new populations in different regions where resources are not naturally available (e.g., Patagonia). This expansion may have been facilitated by modern anthropogenic changes, including the proliferation of artificial roosting structures (Russo and Ancillotto 2015; Stevens et al. 2021) and large-scale agriculture that provides sufficient insect prey to sustain massive maternal colonies (Russell et al. 2011). In fact, the southern tip of Argentina was not included in previous distribution maps, although we have reports from the National Rabies Surveillance that several individuals (one of those included in our analysis) have been recently found in Ushuaia, Tierra del Fuego. Interestingly, this expansion of the range to higher latitudes has also been observed in North American populations (McCracken et al. 2018).