Abstract
The Brazilian free-tailed bat (Tadarida brasiliensis) is a widely distributed Neotropical Molossidae species with significant ecological importance in pest control. Despite its broad range, the genetic diversity and population structure of South American populations remain poorly understood. This study assessed the genetic diversity, population structure, and demographic history of T. brasiliensis in southern Brazil (Capão do Leão) compared to available sequences from South, North, and Central American populations using COI and D-loop mitochondrial markers. Bayesian and coalescent-based analyses revealed deep divergence (~8.9 Mya) between northern and southern lineages, likely influenced by Andean uplift and Pleistocene climatic shifts. High genetic diversity was observed in southern populations, with distinct clades suggesting historical isolation. Population structure analyses confirmed significant differentiation between regions, with isolation-by-distance as a key driver. Demographic tests indicated post-glacial expansions (~0.5–0.12 Mya) in southern populations. These findings suggest that T. brasiliensis in South America comprises a single evolutionary lineage with complex demographic dynamics shaped by historical biogeographic barriers. The study highlights the need for further research on migratory connectivity and conservation strategies for this ecologically important species.
Key words
COI; D-loop; Demographic history; Gene flow; Phylogeography; Population genetics
INTRODUCTION
The Brazilian free-tailed bat, Tadarida brasiliensis (Chiroptera: Molossidae), is a Neotropical species of significant ecological and economic importance playing a crucial role in pest control, benefiting agriculture throughout its range (Kunz et al. 1995, Cleveland et al. 2006). This species is known for its impressive flight capacity, reaching high speeds and long distances (McCracken et al. 2016). This high ability is thought to be linked to the species’ broad geographic distribution and migratory behavior (Wilkins 1989, McCracken 2003, Russell et al. 2005a, Speer et al. 2017, Llaven-Macías et al. 2021). The distribution of T. brasiliensis extends from the United States to Argentina. It is widely distributed across the southern USA, Central America, the Caribbean Islands, and throughout South America (Wilkins 1989, Simmons 2005, Fabián & Gregorin 2007, Barquez et al. 2015, Rodríguez-San & Allendes 2016, Zegarra et al. 2020, Arévalo-Cortés et al. 2024). In South America, species distribution models have identified the Andean ecosystems on both sides of the Andes Mountain Range as the most suitable habitats for T. brasiliensis, while its presence in the Amazon region is considered unlikely (Escobar et al. 2015, Amaral et al. 2023). Interestingly, the Peruvian population exhibits a discontinuous distribution, suggesting that the Andes may act as a biogeographical barrier, limiting the species’ dispersal (Zegarra et al. 2020). Overall, most South American populations’ distribution and migratory patterns remain poorly understood (Botto-Nuñez et al. 2018, Boero et al. 2020).
The species inhabits several environments, including urban areas, forests, and deserts, often utilizing human-made structures (or also caves in the northern distribution) to form large and highly gregarious colonies (Schwartz 1955, Russell & McCracken 2006, Fabián & Gregorin 2007, Hristov et al. 2010, Zegarra et al. 2020). The formation of large colonies and the high dispersal ability (Kunz et al. 1995, Cleveland et al. 2006) underscores the importance of further studies on its genetics structure, phylogeography, and migratory behavior.
Tadarida brasiliensis demonstrates remarkable ecological and geographical adaptability (Wilkins 1989, Russell & McCracken 2006). Its broad distribution range and environmental adaptability have led to the recognition of nine distinct subspecies: T. b. brasiliensis, widely occurring in South America; T. b. cynocephala, occurring in the southeastern United States; T. b. mexicana, which occurs from the western United States to South Mexico; T. b. intermedia, occurring in Central America; and five subspecies occurring on the Caribbean islands (T. b. antillularum in Dominica; T. b. bahamensis in the Bahamas; T. b. constanzae in Haiti; T. b. murina in Jamaica; and T. b. muscula in Cuba) (Schwartz 1955, Russell & McCracken 2006, Morales et al. 2018). The differentiation among such subspecies is mostly based on skull morphology, genetic data (isoenzyme and D-loop control region sequences), geographic distribution, and migratory behavior (Schwartz 1955, Owen et al. 1990, Russell et al. 2005a, b, Russell & McCracken 2006). However, analysis of microsatellite data show gene flow between different populations, suggesting some level of connectivity across the species’ range (Morales et al. 2018). Thus, the genetic and morphometric differentiation of all populations, including the South America populations, has not been considered sufficient to be classified as independently evolving lineages or subspecies (Morales et al. 2018). Therefore, it is hypothesized that the phenotypic differences among T. brasiliensis subspecies are associated with regional climatic variations, leading to the recommendation that individuals be classified as either migratory or non-migratory (Morales et al. 2016).
Mitochondrial (CYT-b, COI, D-loop) and nuclear sequences (microsatellites) have been used to elucidate the phylogeographic patterns and historical demography of T. brasiliensis (Russell et al. 2005a, b, 2011, Clare et al. 2007, Morales et al. 2016, 2018, Speer et al. 2017). These studies have revealed high genetic diversity and significant genetic differentiation among populations, likely driven by geographic barriers and historical climatic events, as the transition between the Nearctic and Neotropical regions. These findings suggest that while T. brasiliensis is capable of extensive migration, local populations may experience limited gene exchange, leading to distinct genetic signatures (Russell & McCracken 2006, Morales et al. 2018). However, studies on the genetic diversity and structure of South American populations are scarce (Morales et al. 2018), limiting our understanding of the connectivity and migratory behavior of T. brasiliensis in the Southern Hemisphere.
In this study, we present a comprehensive assessment of genetic diversity, population structure, and demographic dynamics of T. brasiliensis across its distribution by examining a southern Brazilian population and comparing it with populations from South, Central and North America. First, to evaluate whether distinct T. brasiliensis species exist, we compared available COI sequences with those from our southern Brazilian population. We hypothesized that interpopulation genetic distance values would be less than 0.6%, consistent with the suggested barcode gap for bats (Clare et al. 2007). Using DNA barcoding methodologies, we assessed geographical clustering of sequences and examined COI genetic distances both within and between populations. Second, we tested three hypothesis of population diversity, structure and dynamics, as follow: (i) the northern and southern populations are genetically isolated due to geographical barriers and historical climatic events; (ii) the southern populations exhibit signatures of post-glacial expansion and long-term stability, contrasting with the panmictic northern populations; and (iii) the genetic diversity of the southern Brazilian population results from migratory connectivity with other southern populations. To achieve this, we analyzed both avaliable and new COI and D-loop data for T. brasiliensis, integrating phylogeographic, coalescent-based, and population genetic approaches. This provides the first comprehensive assessment of South American T. brasiliensis, addressing key gaps in our understanding of its evolutionary history and ecological resilience in south Brazil.
MATERIALS AND METHODS
Tissue and DNA extraction
For this study, 53 wing tissue samples, preserved in 90% ethanol at -20°C, were obtained from the tissue collection of the Laboratório de Aves e Mamíferos, Universidade Federal de Pelotas, Rio Grande do Sul, Brazil. The samples were collected from T. brasiliensis individuals captured at a colony in Capão do Leão, southern Brazil (31°48’3.820”S, 52°24’30.906”W; Figure 1) during 2014 and 2015, with proportional monthly representation of both sexes (Supplementary Material - Table SI). This colony inhabits the attic of an abandoned residential structure built in 1980 within a wetland-bordered area, occupying 238 m³. The population peaks at approximately 30,000 individuals during warmer months, and remains active year-round, despite reducing in size during the colder periods (Rui, personal observation). Bat collection was conducted under SISBio permit number 52646-1. Genomic DNA was extracted using the DNeasy® Blood & Tissue Kit (Qiagen), following the manufacturer’s protocol.
Geographic distribution of Tadarida brasiliensis sampling localities for COI and D-loop sequence data throughout the Americas. New sequences obtained in this study originated from Capão do Leão, Rio Grande do Sul, Brazil (solid circle in south Brazil). Guatemalan collection sites (El Rancho and El Progresso) are separated by 10 km..
COI and D-loop PCR and sequencing
A fragment of 687 bp for the cytochrome c oxidase subunit I gene (COI) and approximately 700 bp for the mitochondrial control region (D-loop) were analyzed. The COI fragment was amplified using the primers VF1D (5’-TTCTCAA CCAACCACAARGAYATYGG-3’) and VR1D (5’-TAGACTT CTGGGTGGCCRAARAAYCA-3’) (Ivanova et al. 2006). For D-loop, we used the primer pair F(mt) (5’-GTTGCTGGTTTCACGGAGGTAG-3’) and P(mt) (5’-TCCTACCATCAGCACCCAAAGC-3’) to amplify a variable region (Wilkinson & Chapman 1991). For both sequences, PCR reactions were carried out in 10µl volumes comprising approximately 100 ng of DNA, 0.4 µM of each primer, and 1x of HotStarTaq Master Mix (Qiagen). The amplification program followed 95°C for 5 minutes for initial denaturation; 30 cycles of 95°C for 1 minute for denaturation; 58°C for 50 seconds for COI primers annealing and 55°C for 50 seconds for D-loop primers annealing; 72°C for 1 minute for extension; with a final extension of 72°C for 5 minutes. The PCR products underwent 1% agarose gel electrophoresis and were subsequently purified using the PureLink Quick Gel extraction and PCR purification combo kit (ThermoFisher), following the manufacturer’s guidelines. Both forward and reverse strands of COI and D-loop sequences were sequenced at Macrogen (dna.macrogen.com/eng/). Consensus sequences were constructed by analyzing chromatograms using the Staden 4.11 package (Staden 1996). It was obtained 50 sequences for COI (PP806575 to PP806624) and 51 for D-loop (PP802915 to PP802965) (see Supplementary Table SI). The sequenced specimens were registered in the Brazilian National Management System for Genetic Heritage and Associated Traditional Knowledge (SISGEN) under the number AAEA29C.
In silico searches and sequence alignment
In silico searches were conducted to retrieve T. brasiliensis sequences for COI from the GenBank and the Barcode of Life Data System (BOLD) databases, and for the D-loop from the GenBank. Duplicate sequences were identified, with only one copy retained for analysis. Sequences containing stop codons were excluded, as these are indicative of nuclear mitochondrial DNA (Numt). The COI matrix has a total length of 687 bp and includes 128 sequences aligned using the Clustal algorithm. The D-loop matrix has a total length of 719 bp and comprises 271 sequences aligned using the Muscle algorithm. Both alignments were performed in MEGA11 (Tamura et al. 2021). For D-loop, three specimens from Capão do Leão have an exclusive 160 bp insertion (see Supplementary Table SI). The map in Figure 1 displays the collection sites of all sequences obtained from GenBank and from Capão do Leão. Supplementary Table SII lists the accession numbers and locations for all sequences retrieved from GenBank.
Bioinformatic analysis
DNA barcoding methods were applied to the COI dataset to assess whether distinct populations of T. brasiliensis form geographically determined clustering of sequences potentially representing separate (valid) species. These analyses were performed using Bayesian Inference analysis and genetic distances, applying the following approaches:
I. The identification of monophyletic clusters of sequences (Hebert et al. 2003, Herbert & Gregory 2005) was evaluated through COI sequences tree reconstructed by Bayesian analysis in MrBayes 3.2.7 (Ronquist et al. 2012). Sequence from Molossus molossus was used as outgroup (COI: JF455084). The HKY+I nucleotide substitution model was selected based on the Akaike Information Criterion test (AIC) (Akaike 1974) performed in jModelTest 2.1.8 (Darriba et al. 2012, Guindon & Gascuel 2003). The analysis involved two parallel runs with four MCMC chains, 1,000,000 generations, sampled every 1,000, with a burn-in of 25%. Convergence was assessed through the average standard deviation of split frequencies.
II. Species delimitation analysis was performed using coalescent-based analysis with the General Mixed Yule Coalescent (GMYC) model. Phylogenetic relationship for the COI gene was reconstructed through maximum likelihood (ML) analysis using IQTREE v.1.6.12 software (Nguyen et al. 2015), using the same nucleotide substitution as previous analysis, employing 10,000 ultra-fast bootstrap replicates. The resulting tree was used to construct ultrametric trees with non-zero length branches via the chronos function in the APE package (Paradis et al. 2004). The ultrametric tree was then used as input for GMYC analysis with the SPLITS package (Ezard et al. 2009) on R software v4.4.1 (R Core Team 2024).
III. The COI genetic distance analysis was conducted for sequences grouped by geographical location, both within and between each group. These groups were: Capão do Leão-Brazil (50 sequences); São Paulo-Brazil (20); Guatemala (06); Mexico (02); and USA (50). Genetic distances were calculated using the Kimura 2-parameter (K2P) nucleotide substitution model (Kimura 1980), as suggested by Hebert et al. (2003) and used by Clare et al. (2007) for bats. All ambiguous positions were removed for each sequence pair (pairwise deletion option). This analysis was performed in MEGA11 software.
IV. The differences in COI genetic divergences between intra and inter-geographical groups were analyzed to identify the presence of automatic grouping of COI sequences that reflect these geographical groups. In this analysis, sequences with lower barcode gap values (<1.5%) are expected to cluster together, indicating minimal differences among populations (Cai et al. 2010). This analysis was performed using the Automatic Barcode Gap Discovery (ABGD) method (Puillandre et al. 2011) within iTaxoTools (Vences et al. 2021, itaxotools.org).
COI and D-loop mitochondrial markers, with their distinct evolutionary dynamics, serve as complementary tools for reconstructing the population history of T. brasiliensis throughout the Neotropics. By combining their results, we achieve improved accuracy and resolution in population genetic studies, while maintaining a cost-effective and widely applicable approach for assessing demographic patterns across multiple timescales (Dias et al. 2017, Olímpio et al. 2025). Genetic diversity parameters of these markers, such as number of haplotypes (H), haplotype diversity (Hd), nucleotide diversity (pi), and polymorphic sites (S), were estimated for each genetic marker in DnaSP software (Rozas et al. 2017), excluding gaps/missing data and invariable sites.
Genetic differentiation among T. brasiliensis populations was quantified using an Analysis of Molecular Variance (AMOVA) based on Φ-statistics (analogous to Wright’s F-statistics; Excoffier et al. 1992). Population structure was evaluated at two hierarchical levels: (i) differentiation among geographical populations (Fst), and (ii) divergence between the combined Central/North American group (USA, Mexico, Guatemala) and the South American group (Brazil, Argentina, Chile) (Fct), to assess large-scale geographic isolation effects. Fst and Fct range from 0 (no differentiation) to 1 (complete differentiation), but statistically significant positive values equal or higher than 0.15 is consistent with high genetic differentiation (Holsinger & Wier 2009, Meirmans & Hedrick 2011, Weir & Goudet 2017). Separate AMOVA analysis were conducted for the COI and D-loop datasets, each with 10,000 permutations. The correlation between geographical distance and Fst data (YMatrix=Fst) was analyzed using Mantel’s test (Mantel & Greenhouse 1967) with 10,000 permutations. This analysis was performed to determine whether there is significant positive correlation between the population’s genetic distances and geographical distances, inferring the adequacy of an isolation by distance (IBD) model. The linearized version of Fst was used in Mantel’s calculations to consider the potential nonlinear relationship between Fst and geographic distances (Slatkin 1995). The geographical distance matrix was constructed using Geographic Distance Matrix Generator v1.2.3 (Ersts 2013) with geographical coordinates obtained from Google Earth Pro (Google Inc). The log-transformed values were used in Mantel’s analysis (Slatkin 1993), as performed in Excel®. For sequences from the northern distribution within the same country, obtained from multiple locations, the central geographical coordinate of the country was used for Mantel’s analysis (Supplementary Table SIII). AMOVA and Mantel tests were performed using Arlequin 3.5 (Excoffier & Lischer 2010). A Mismatch Distribution Analysis was conducted using DNAsp to test for signals of demographic or spatial population expansion events. To assess whether distinct patterns exist between the northern and southern populations, the analysis was performed separately for each group. The relationships between the resulting haplotypes were analyzed using a haplotype network inferred by the median-joining method, as implemented in the PopArt software (Leigh & Bryant 2015).
The neutrality tests Tajima’s D (Tajima 1989) and Fu & Li’s F (Fu 1997) were performed to assess demographic historical patterns deviations from the expectations of neutrality under a constant population size and random mating. Tajima’s D compares the average nucleotide differences between pairs of sequences (π) and the number of segregating sites (S); while Fu’s Fs neutrality test detects signs of selection, population expansion, or recent demographic events in a population comparing the observed number of haplotypes with the expected number of haplotypes under a neutral model of evolution. These tests were performed in DnaSP software (Rozas et al. 2017).
The Bayesian time-calibrated tree was constructed using separated COI and D-loop dataset in BEAST 2.6.3 (Bouckaert et al. 2019), applying HKY+I substitution models and a coalescent-constant size tree model for each analysis. Divergence times were estimated using a normal strict molecular clock with calibration (priors) set to 42 Mya, based on the dating of the oldest molossid fossil (Czaplewski et al. 2003). The analysis involved 200 million MCMC iterations, sampling every 10,000 steps. Convergence was ensured by checking if Effective Sample Size (ESS) was higher than 200 for all parameters using Tracer v.1.7.1 (Rambaut et al. 2018). Trees were summarized and annotated using the TreeAnnotator module in BEAST with a 10% burn-in. Molossus molossus sequences were used as outgroup (COI: JF455084; D-loop: KT721416).
Finally, the southern population size distribution of T. brasiliensis over time was estimated using the Extended Bayesian Skyline Plot (EBSP) analysis for both genetic markers. This analysis was performed in BEAST 2.6.3 (Bouckaert et al. 2019). The substitution models, clock models, and substitution rate priors were consistent with those used in the Bayesian time-calibrated analysis. The packages APE (Paradis et al. 2004), ggplot2 (Wickham et al. 2016) and tidyverse (Wickham et al. 2019) were used for graphical visualization of the results.
RESULTS
DNA barcode analysis
The Bayesian phylogenetic tree based on COI sequences reveals a geographical clustering of northern and southern populations, forming two main monophyletic groups, one consisting of all northern population sequences and the other comprising two distinct monophyletic clades from the southern populations (Figure 2; Supplementary Material - Figure S1). However, the GMYC coalescence analysis detected seven maximum likelihood (ML) clusters with high confidence interval of ML clusters (1 to 18, including the M. molossus outgroup sequence), negative threshold time (-0.314), and non-significant likelihood ratio between the null model and the GMYC model (1.07; p-value = 0.59).
Bayesian time-calibrated phylogenetic tree for the COI and D-loop mitochondrial data. Molossus molossus sequences were used as outgroup. Values above nodes refer to the mean divergence dating. Bars represent the 95% highest posterior density interval (HPD) for the divergence dating. Time in Mya. COI phylogeny: Clade 1 comprises sequences from São Paulo and Capão do Leão (Brazil); Clade 2 contains only Capão do Leão sequences. D-loop phylogeny has only sequences from Capão do Leão (Brazil). ** Clade comprising D-loop sequences with an exclusive insertion of 160 bp (KD18, KD24, and KD49). See Supplementary Figure S1 for sequence codes.
The mean intrapopulation values of the genetic distances for COI sequences ranged from 0.003 (to Brazil, Mexico, and USA populations) to 0.029 (Capão do Leão population), while the mean interpopulation values ranged from 0.004 (Mexico x Guatemala) to 0.062 (Capão do Leão x USA) (Table I). The interpopulation comparisons involving the studied population sequences showed the highest genetic distance values, while the comparisons among Guatemala, Mexico, and USA showed the lowest. Overall, the ABGD analysis identified groups corresponding to the southern and northern populations with no intermixing of sequences between the two regions (Supplementary Figure S2). With MinSlope = 1.5, only one partition was found. With MinSlope = 1.0, 12 clusters were identified, with priors for maximum intraspecific divergence ranging from 0.0017 to 0.0028.
Estimates of genetic distances for COI sequences grouped by geographical origin. Mean intrapopulation values (diagonal), mean interpopulation values (above diagonal), and minimum and maximum values (in parenthesis) are shown.
Genetic diversity
The mitochondrial markers COI and D-loop revealed a high genetic diversity within T. brasiliensis (Table II). For COI, 12 haplotypes (out of 128 sequences) were identified with an average haplotype diversity (Hd) of 0.782, a global nucleotide diversity (Pi) of 0.023, and a total of 10 polymorphic sites (S). Comparing the populations with 10 or more sequences sampled, Capão do Leão population exhibited the highest genetic diversity values. For D-loop, from 271 sequences, 216 haplotypes were identified with an average haplotype diversity (Hd) of 0.998, a nucleotide diversity (Pi) of 0.092, and 171 polymorphic sites (S). Comparing the populations with 10 or more sequences sampled, all of them showed high values of haplotype diversity.
Summary statistics of genetic diversity and neutrality tests performed for the COI and D-loop mitochondrial datasets. *: p<0.05; na: not applicable, at least four sequences are needed for the analysis. H: number of haplotypes. Hd: haplotype diversity: Pi: nucleotide diversity. S: polymorphic sites.
Population structure
Our analyses reveal contrasting evolutionary patterns between markers within and between populations. The COI haplotype network exhibits clear star-like pattern in both northern and southern populations, indicating recent population expansion within each region. In contrast, the D-loop network displays extensive mutational steps among haplotypes in both groups (consistent with its hypervariable nature), showing evenly distributed haplotypes overall but with localized star-like patterns restricted to specific southern haplotype clusters (Figure 3). These patterns suggest limited gene flow and an absence of major demographic expansions within population groups. The AMOVA results (Fst and Fct) reflect moderate to high genetic structure within and between northern and southern populations. Most Fst pairwise comparisons were greater than 0.150 (p<0.05), for both COI and D-loop data (Figure 4), with particularly strong genetic structuring among southern populations (Fst=0.089-0.415; p<0.05). The Fct indicate significantly genetic differentiation between northern and southern groups for for D-loop data (D-loop: Fct=0.598; p<0.05), while the COI showed high but non-significant differentiation (COI: Fct=0.566; p>0.05).
Matrix of pairwise COI and D-loop Fst (Wright’s F-statistics) values for Tadarida brasiliensis populations. USA: United States, Mex: Mexico, Gua: Guatemala, Chi: Chile, Arg: Argentina, Bra: Brazil, CL: Capão do Leão. *: p value < 0.05.
Haplotype network reconstructed from COI and D-loop sequences of Tadarida brasiliensis populations. Each circle represents a unique haplotype, with the size proportional to its frequency in the sample. Small lines in the lines connecting circles indicate mutational steps between haplotypes. The color of the circles denotes geographic origin.
A Mantel test was conducted to understand whether genetic variation is structured according to geographical proximity, assessing the correlation between genetic distances and geographical distances. For both COI and D-loop markers, the Mantel test revealed a strong positive correlation between geographical and genetic distances (r2: COI=0.67, p<0.05 and D-loop=0.63; p<0.05). However, when the test was performed for the southern and northern populations alone, the values were not statistically significant, despite the lower values for COI and D-loop for southern population and lower D-loop value for northern population (southern r2: COI=0.083, p=0.662 and D-loop=0.092, p=0.336; northern r2: COI=0.667, p=0.331 and D-loop=0.092,p=0.336).
Demographic background
The results of the Tajima’s D and Fu’s FS neutrality tests yielded negative values (Table II), indicating an excess of rare genetic variants within populations compared to what is expected under a neutral model of evolution for the genetic markers. For the mismatch distribution analysis, both COI and D-loop regions exhibited a bimodal distribution. While both northern and southern populations displayed this bimodal pattern, it was more pronounced in the southern populations (Supplementary Figure S3).
Based on the COI Bayesian time-calibrated tree (Figure 2), the southern T. brasiliensis population diverged from the northern population approximately 8.9 Mya. In contrast, the D-loop data suggest a more recent divergence time of 0.52 Mya between northern and southern populations. The discrepancy between divergence time estimates may reflect the distinct evolutionary dynamics of these markers, coupled with methodological considerations. Nonetheless, these markers provide complementary interpretations, while COI reveals deep phylogenetic splits, D-loop reflects recent population dynamics.
The COI data recovered two major clades within the southern populations, which diverged around 6.08 Mya: Clade 1, consisting of 29 sequences from Capão do Leão and all sequences from Brazil, and Clade 2, comprising 21 sequences, all from Capão do Leão (Supplementary Figure S1). For the D-loop marker, two clades were also recovered. However, one clade is performed by the three sequences from Capão do Leão that exhibit a 160 bp insertion (Figure 2; Supplementary Figure S4; alignment access: github.com/jlncdr/Dloop.git). This insertion is performed by a 80 bp duplication also found in all other T. brasiliensis D-loop sequences. But unlike the other specimens, the three specimens from Capão do Leão have five repeats, while all the others have only three repeats. Tandem repeats in D-loop represent a common and highly variable feature across Chiroptera, particularly within the Vespertilionidae family (Wilkinson & Chapmann 1991, Wilkinson et al. 1997, Gorobeyko et al. 2023). These repeats are consistently localized adjacent to the tRNA-Pro gene and exhibit both interspecific and intraspecific variation in length, sequence composition, and copy number. Although tandem repeats in the D-loop region have been well characterized in several bat families, their presence in Molossidae species has remained undocumented until now.
The effective population size analysis using the Extended Bayesian Skyline Plot (EBSP) for COI and D-loop sequences from southern populations indicated a trend of recent expansion. For COI, the expansion occurred within the last 0.5 Mya [mean 4.3113e-3; 95% HPD interval (4.2948e-3; 4.3284e-3)], while for the D-loop, the expansion occurred around 0.12 Mya [mean 0.4652; 95% HPD interval (0.4626; 0.4678)] (Figure 5).
HEffective population size (Ne) variation of Tadarida brasiliensis inferred for COI and D-loop sequences over time. The dashed black line represents the average variation in population size over time, with 95% confidence intervals (CPD) in green. Time in Mya.
DISCUSSION
In this study we provide a comprehensive assessment of the genetic structure of populations of the South American bat Tadarida brasiliensis, integrating population genetic approaches aimed at addressing gaps on its evolutionary history in south Brazil. The results from our DNA barcoding strategies (monophyly of sequences, coalescent-based species delimitation, and COI genetic distances analyzes) show that, despite the COI divergence between North and South American populations occurring approximately 8.9 Mya (Figure 2; Supplementary Figure S1), the COI genetic distance values are close to or below the estimated intraspecific threshold of 0.6% (Clare et al. 2007, 2011), ranging from 0.03% to 0.62% (Table I). This finding suggests that the various populations of T. brasiliensis should be regarded as a single evolving lineage rather than distinct subspecies (Morales et al. 2018). However, the highest values of intrapopulation COI genetic distances are observed in the comparisons involving Capão do Leão sequences. These results may reflect the split event of COI sequences in South American populations into two clades, which occurred during the geological and climatic events that occurred around 6 Mya (Late Miocene in South America) with the central Andes uplift, expansion of grasslands and reduction of rainforests took place, for example (Figure 2; Supplementary Figure S1). The D-loop data further corroborate the existence of North and South geographical genetic groups (Figure 2; Supplementary Figure S1). However, in this case, evolutionary events occurred approximately 0.5 Mya. The South American clade, composed exclusively by three D-loop sequences from Capão do Leão, contains the sequences that exhibit two additional tandem repeats compared to the repeat pattern already present on all T. brasiliensis sequences. The genetic-evolutionary events that led to this duplication occurred prior to the divergence events between the North and South populations for D-loop, which took place around 0.9 Mya (Figure 2; Supplementary Figure S1).
All populations studied exhibited high values for diversity indices for both the COI and D-loop markers (Table II; Figure 3). Migration strategies have been correlated with population structure and genetic diversity (Moussy et al. 2023), and the migration strategies for the northern populations of T. brasiliensis have been extensively studied (Villa & Cockrum 1962, Keeley & Keeley 2004, Russell & McCracken 2006, Morales et al. 2018). The northern populations of T. brasiliensis show a mix of resident and migratory behaviors, including long-distance and regional migrations (Russell & McCracken 2006). In contrast, information on the migration strategies of South American populations of T. brasiliensis remains limited. Evidence on density and circannual sex distribution in T. brasiliensis roosts in southern South America suggests that females disperse out of the roost, or migrate, in the winter season while males show a resident behavior (Marques & Fabián 1994, Fabián & Marques 1996, Botto-Nuñez et al. 2018, Boero et al. 2020). For the largest known colony of T. brasiliensis in South America, located in southern Argentina and estimated to comprise 12 million individuals, there is evidence of migration in both females and males (Díaz et al. 2021). Regarding the T. brasiliensis colony from Capão do Leão, although we recorded the sex of all tissue samples analyzed in this study, we currently lack accurate information about its circannual sex distribution. While initial observations suggest that this colony’s migration patterns may align with the annual cycle of migratory colonies in both hemispheres, this population roost is used year-round, including during both the reproductive and non-reproductive seasons.
Our data demonstrates pronounced genetic differentiation among populations, with the highest Fst values shown in pairwise comparisons among northern and southern populations (Figure 4). The Fct results indicate a clear genetic structure between northern and southern genetic groups, with minimal gene flow, at least for the D-loop marker. Thus, collectively, Fst and Fct results suggest restricted gene flow among populations, with isolation-by-distance (IBD; Slatkin 1993) as a plausible mechanism. Furthermore, the observed genetic differentiation in all pairwise comparisons, as indicated by Fst, appears strongly influenced by geographic distance, as confirmed by significant Mantel test results. Although this result could be improved by sampling additional populations across the Neotropical distribution of T. brasiliensis, particularly from northern South America. In contrast, the Mantel test was not significant when southern and northern populations were analyzed separately. This suggests that, at least for southern populations, genetic differentiation is primarily driven by limited gene flow, as reflected by the moderate to high Fst values, likely influenced by historical and demographic factors. For northern populations, the non-significant Fst values suggest high levels of gene flow, potentially indicative of panmictic population dynamics. This agrees with previous findings indicating that T. brasiliensis populations in North and Central America exhibit low genetic structure and high gene flow across large geographic distances, resembling a nearly panmictic population (McCracken et al. 1994, Russell et al. 2005a).
The negative values shown in the neutrality tests (Table II) might indicate positive selection or population expansion process, but the evidence is not strong enough to rule out the influence of genetic drift or other stochastic processes. However, the mismatch distribution analysis (Supplementary Figure S3) revealed a bimodal distribution in general for both datasets. The bimodal pattern indicates a demographic history of distinct expansion events or a population structure with distinct groups. Notably, when the analysis was performed separately for northern and southern populations, distinct trends emerged: the northern population exhibited signs of historical demographic stability without multiple expansion events, despite the recent population growth and colonization to new habitats (McCracken et al. 2018); while the southern population displayed evidence of more complex historical dynamics, involving more than one population expansion event. This pattern is corroborated by the EBSP results (Figure 5) which reveal a recent population expansion in southern populations with distinct time estimates between markers: 0.5 Mya for COI versus 0.12 Mya for D-loop. The temporal discrepancy parallels the earlier divergence estimates between northern and southern populations (8.9 Mya for COI versus 0.52 Mya for D-loop), providing insights into the species’ demographic history. The COI-inferred expansion (~0.5 Mya) coincides with the Mid-Pleistocene Warm Periods (Da et al. 2023), suggesting prolonged population growth following initial divergence, likely facilitated by range shifts during interglacial warming. In contrast, the D-loop expansion (~0.12 Mya) aligns with post-glacial recolonization during the Last Interglacial (Kukla et al. 2002), indicating rapid demographic recovery following Pleistocene glacial maxima. This suggests that the demographic histories of T. brasiliensis populations have been shaped by historical events, including recent population expansion with gene flow, large population sizes, and long-term stability. These processes, specific to each geographical region, have collectively contributed to the observed high genetic variability in this species.
The formation of the Isthmus of Panama (3–2.8 Mya, Late Pliocene) and the uplift of the Andes (10–5 Mya, Late Miocene to Pliocene) were two major geological events that shaped the geography, climate, and biodiversity of the Americas (O’Dea et al. 2016, Vázquez-López et al. 2024). For the widespread Neotropical bat Glossophaga soricina (Phyllostomidae), the major lineage splits occurred between 5–2.4 Mya (Late Mioceno to Late Plioceno), likely driven by these geological events (Dias et al. 2017). Similarly, the genus Molossus underwent rapid adaptive radiation, with most species emerging around 2.5 Mya (Pliocene–Pleistocene transition) (Olímpio et al. 2025). For T. brasiliensis, the higher environmental suitability and greater genetic diversity suggest that the species originated in central and eastern regions of South America, with subsequent dispersal to northern and southern areas (Zegarra et al. 2020, Amaral et al. 2023). Considering our results, the northern and southern lineages split occurred during the uplift of the Andes and the Isthmus of Panama formation, with southern population expansion during the Late Pleistocene.
Regarding the T. brasiliensis population from Capão do Leão the high genetic diversity observed may reflect mixed migratory influences. While this roost primarily harbors non-migratory individuals, it could simultaneously receive migratory bats from southernmost latitudes during winter or northern latitudes during summer, a seasonal pattern analogous to that documented in higher-latitude populations in North America (Morales et al. 2016). Indeed, studies conducted at similar southern latitudes suggest circannual variation in sex ratios within T. brasiliensis colonies, with activity persisting even during the cold season (Marques & Fabián 1994, Fabián & Marques 1996, Romano et al. 2015, Botto-Nuñez et al. 2018, Boero et al. 2020). Some of these colonies consist of reproductive females and their offspring, termed ‘maternal colonies’. Maternal colonies are typically large, and their female-biased seasonal dynamics suggest migratory behavior. However, activity ceases (or declines abruptly) during colder months (Romano et al. 2015, Boero et al. 2020, Llaven-Macías et al. 2021). The coexistence of these demographic patterns with high genetic diversity suggests a potential migratory corridor along the South America populations, where Capão do Leão may serve as a potential stopover (or coexistence zone) site between temperate and tropical populations.
Despite our study being based on new sequences from a single population in Southern Brazil, the genetic diversity detected in this population is comparable to that reported for North American populations, which included individuals sampled over a 2,800 km range (e.g., from California to South Carolina, USA) and over 100 individuals (McCracken et al. 1994, Russell & McCracken 2006, Morales et al. 2018), in contrast to the 50 individuals analyzed for the Capão do Leão population. Nonetheless, nuclear genetic markers are essential for elucidating and refining the evolutionary history of Tadarida brasiliensis in South America.
Acknowledgements
We thank Dr. Lizandra Robe and Dr. Jeferson Vizentin-Bugoni for valuable help during data analysis and to the reviewers for suggestions and improvement of the manuscript. This work was supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico, Brazil (CNPq, grant number 445498/2014-5).
References
- AKAIKE H. 1974. A new look at the statistical model identification. IEEE Trans Autom Control 19: 716-723.
- AMARAL IS, PEREIRA JB, VANCINE MH, MORALES AE, ALTHOFF SL, GREGORIN R, PEREIRA MJR, VALIATI VH & OLIVEIRA LR. 2023. Where do they live? Predictive geographic distribution of Tadarida brasiliensis brasiliensis (Chiroptera, Molossidae) in South America. Neotrop Biol Conserv 18 (3): 139-156.
- ARÉVALO‑CORTÉS J, TULCAN‑FLORES J, ZURC D, MONTENEGRO‑MUÑOZ SA, CALDERÓN‑LEYTÓN JJ & FERNÁNDEZ‑GÓMEZ RA. 2024. Description of the echolocation pulses of insectivorous bats with new records for Southwest Colombia. Mamm Res 69: 231-244.
-
BARQUEZ R, DIAZ M, GONZALEZ E, RODRIGUEZ A, INCHÁUSTEGUI S & ARROYO-CABRALES J. 2015. Tadarida brasiliensis The IUCN Red List of Threatened Species. Available at: https://www.iucnredlist.org/species/21314/22121621 (Accessed on 22 August, 2024).
» https://www.iucnredlist.org/species/21314/22121621 - BOERO L, POFFO D, DAMINO V, VILLALBA S, BARQUEZ RM, RODRÍGUEZ A, SUÁREZ M & BECCACECE HM. 2020. Monitoring and Characterizing Temporal Patterns of a Large Colony of Tadarida brasiliensis (Chiroptera: Molossidae) in Argentina Using Field Observations and the Weather Radar RMA1. J Remote Sens 12: 210-227.
- BOTTO-NUÑEZ G, GENTA M, DÍAZ M, RODALES AL & GONZÁLEZ EM. 2018. Circannual sex distribution of the Brazilian free-tailed bat, Tadarida brasiliensis (Chiroptera: Molosssidae), suggests migration in colonies from Uruguay. Mastozool Neotrop 25 (1): 213-219.
- BOUCKAERT R ET AL. 2019. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS Comput Biol 15: e1006650.
- CAI Y, YUE B, JIANG W, XIE S, LI J & ZHOU M. 2010. DNA barcoding on subsets of three families in Aves. Mitochondrial DNA 21(3-4): 132-137.
- CARRUTHERS T & SCOTLAND R. 2020. Uncertainty in divergence time estimation. Syst Biol 70(4): 855-861.
- CLARE EL, LIM BK, ENGSTROM MD, EGER JL & HEBERT PD. 2007. DNA barcoding of Neotropical bats: species identification and discovery within Guyana. Mol Ecol Notes 7(2): 184-190.
- CLARE EL, LIM BK, FENTON MB & HEBERT PD. 2011. Neotropical bats: estimating species diversity with DNA barcodes. PloS ONE 6(7): e22648.
- CLEVELAND CJ ET AL. 2006. Economic value of the pest control service provided by Brazilian free-tailed bats in south-central Texas. Front Ecol Environ 4: 238-243.
- CZAPLEWSKI NJ, MORGAN GS & NAEHER T. 2003. Molossid bats from the late Tertiary of Florida with a review of the Tertiary Molossidae of North America. Acta Chiropt 5(1): 61-74.
-
DARRIBA D, TABOADA GL, DOALLO R & POSADA D. 2012. jModelTest 2: more models, new heuristics and parallel computing. Nat Meth 9: 772. https://doi.org/10.1038/nmeth.2109.
» https://doi.org/10.1038/nmeth.2109 - DA S, ZHANG Z, LI Y, XU Q, FAN B, WANG S, DONG J, WANG Y & CHI Z. 2023. Pollen-based quantitative paleoclimatic record spanning the Mid-Brunhes Event in the Nihewan Basin, north China. Palaeogeogr Palaeoclimatol Palaeoecol 162: 111377.
- DIAS C, JÚNIOR J, PERINI F & SANTOS F. 2017. Biogeographic scenarios for the diversification of a widespread Neotropical species, Glossophaga soricina (Chiroptera: Phyllostomidae). Syst Biodiver 15: 440-450.
- DÍAZ M, MIOTTI D & ALURRALDE SG. 2021. Los murciélagos del dique Escaba y sus alrededores (Tucumán, Argentina): investigación, educación y conservación. In: TRACANNA BC (Ed), Escaba, un embalse en las Yungas del Noroeste Argentino. Serie Conservación de la Naturaleza 26. San Miguel de Tucumán: Fundación Miguel Lillo, Tucumán, ARG, p. 127-148.
-
ERSTS PJ. 2013. Geographic Distance Matrix Generator (version 1.2.3). American Museum of Natural History, Center for Biodiversity and Conservation. Available at: http://biodiversityinformatics.amnh.org/open_source/gdmg
» http://biodiversityinformatics.amnh.org/open_source/gdmg - ESCOBAR LE, JUAREZ C, MEDINA-VOGEL G & GONZALEZ CM. 2015. First Report on bat mortalities on wind farms in Chile. Gayana 79(1): 11-17.
- EXCOFFIER L & LISCHER HEL. 2010. Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Res 10(3): 564-567.
- EXCOFFIER L, SMOUSE PE & QUATTRO JM. 1992. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics 131: 479-491.
-
EZARD T, FUJISAWA T & BARRACLOUGH TG. 2009. SPLITS: SPecies’ Limits by Threshold Statistics. R package version 1.0-18/r45. Available at: http://R-Forge.R-project.org/projects/splits/ Accessed on: October, 2024.
» http://R-Forge.R-project.org/projects/splits/ - FABIÁN ME & GREGORIN R. 2007. Família Molossidae. In: REIS NR ET AL. (Eds), Morcegos do Brasil. Londrina, Universidade Estadual de Londrina, Paraná, p. 149-165.
- FABIÁN ME & MARQUES VR. 1996. Aspectos do comportamento de Tadarida brasiliensis brasiliensis (I. Geoffroy, 1824) (Chiroptera; Molossidae) em ambiente urbano. Biociências 4: 65-86.
- FU YX. 1997. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147: 915-925.
- GOROBEYKO UV, SHEREMETYEVA IN, KAZAKOV DV & GUSKOV VY. 2023. A new type of tandem repeats in Myotis petax (Chiroptera, Vespertilionidae) mitochondrial control region. Mol Biol Rep 50(6): 5137-5146.
- GUINDON S & GASCUEL O. 2003. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol 52(5): 696-704.
- HEBERT PDN, CYWINSKA A, BALL SL & DEWAARD J. 2003. Biological identifications through DNA barcodes. Proc R Soc Lond B Biol Sci 270: 313-321.
- HEBERT PDN & GREGORY TR. 2005. The promise of DNA Barcoding for Taxonomy. Syst Biol 54: 852-859.
- HOLSINGER KE & WIER BS. 2009. Genetics in geographically structured populations: defining, estimating and interpreting FST Nat Rev Genet 10: 639-650.
- HRISTOV NI, BETKE M, THERIAULT DEH BAGCHI A & KUNZ TH. 2010. Seasonal variation in colony size of Brazilian free-tailed bats at Carlsbad Cavern based on thermal imaging. J Mamm 91(1): 183-192.
- IVANOVA NV, DEWAARD JR & HEBERT PD. 2006. An inexpensive, automation-friendly protocol for recovering high-quality DNA. Mol Ecol Notes 6(4): 998-1002.
- KEELEY ATH & KEELEY BW. 2004. The mating system of Tadarida brasiliensis (Chiroptera: Molossidae) in a large highway bridge colony. J Mamm 85: 113-119.
- KIMURA M. 1980. A simple method for estimating evolutionary rate of base substitutions through comparative studies of nucleotide sequences. J Mol Evol 16: 111-120.
- KUKLA G ET AL. 2002. Last Interglacial Climates. Quat Res 58: 2-13.
- KUNZ TH, WHITAKER JRJO & WADANOLI MD. 1995. Dietary energetics of the insectivorous Mexican free-tailed bat (Tadarida brasiliensis) during pregnancy and lactation. Oecologia 101: 407-415
- LEIGH JW & BRYANT D. 2015. PopART: Full-feature software for haplotype network construction. Methods Ecol Evol 6(9): 1110-1116.
- LLAVEN-MACÍAS V, RUÍZ-MONTOYA L, LÓPEZ-GONZÁLEZ C, RICO Y & NARANJO EJ. 2021. Monthly fluctuation of colony size and composition of the free-tailed bat Tadarida brasiliensis in the southernmost roost of Mexico. Mamm Res 66: 339-348.
- MANTEL N & GREENHOUSE SW. 1967. Equivalence of maximum likelihood and the method of moments in probit analysis. Biometrics 23(1): 154-157.
- MARQUES RV & FABIÁN ME. 1994. Ciclo reprodutivo de Tadarida brasiliensis (I. Geoffroy, 1824) (Chiroptera, Molossidae) em Porto Alegre, Brasil. Iheringia, Série Zoologia 77: 45-56.
- MCCRACKEN G, BERNARD R, GAMBA-RIOS M, WOLFE R, KRAUEL J, JONES D, RUSSELL A & BROWN V. 2018. Rapid range expansion of the Brazilian free-tailed bat in the southeastern United States, 2008-2016. J Mammal 99: 312-320.
- MCCRACKEN GF. 2003. Estimates of population sizes in summer colonies in Brazilian free-tailed bats (Tadarida brasiliensis). Monitoring trends in bat populations of the United States and territories: problems and prospects. In: O’SHEA TJ & BOGAN MA (Eds), United States Geological Survey, Biological Resources Discipline, Information and Technology Report, USGS/BRD/ITR-2003-003, p. 21-30.
- MCCRACKEN GF, MCCRACKEN MK & VAWTER AT. 1994. Genetic structure in migratory populations of the bat Tadarida brasiliensis mexicana J Mammal 75(2): 500-514.
- MCCRACKEN GF, SAFI K, KUNZ TH, DECHMANN DK, SWARTZ SM & WIKELSKI M. 2016. Airplane tracking documents the fastest flight speeds recorded for bats. R Soc Open Sci 3(11): 160398.
- MEIRMANS PG & HEDRICK PW. 2011. Assessing population structure: FST and related measures. Mol Ecol Res 11(1): 5-18.
- MORALES AE, MORA DM & PIÑERO D. 2018. Spatial and environmental factors predict skull variation and genetic structure in the cosmopolitan bat Tadarida brasiliensis J Biogeogr 45: 1529-1540.
- MORALES AE, VILLALOBOS F, VELAZCO PM, SIMMONS NB & PIÑERO D. 2016. Environmental niche drives genetic and morphometric structure in a widespread bat. J Biogeogr 43(5): 1057-1068.
- MOUSSY C, HOSKEN DJ, MATHEWS F, SMITH GC, AEGERTER JN & BEARHOP S. 2013. Migration and dispersal patterns of bats and their influence on genetic structure. Mamm Rev 43(3): 183-195.
- NGUYEN LT, SCHMIDT HA, VON HAESELER A & MINH BQ. 2015. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol 32(1): 268-274.
-
O’DEA A ET AL. 2016. Formation of the Isthmus of Panama. Sci Adv 2(8): e1600883. DOI: 10.1126/sciadv.1600883.
» https://doi.org/10.1126/sciadv.1600883 - OLÍMPIO A, STEFANELLO F, DA NATIVIDADE B, BERNARDI I, LIMA AS, MENDES S, COSTA CS, FRAGA EC, BARROS M & SAMPAIO I. 2025. Cranial morphology reveals a lack of phylogenetic signal and rapid adaptive radiation in the bat genus Molossus (Chiroptera: Molossidae). PLOS One 20(4): e0320117.
- OWEN RD, CHESSER RK & CARTER DC. 1990. The systematic status of Tadarida brasiliensis cynocephala and Antillean members of the Tadarida brasiliensis group, with comments on the generic name Rhizomops Legendre. Occasional Papers Museum of Texas Tech University, n. 133. Texas: Texas Tech University Press, 18 p.
- PARADIS E, CLAUDE J & STRIMMER K. 2004. APE: analyses of phylogenetics and evolution in R language. Bioinformatics 20(2): 289-290.
- PUILLANDRE N, LAMBERT A, BROUILLET S & ACHAZ G. 2011. ABGD, Automatic Barcode Gap Discovery for primary species delimitation. Mol Ecol 21: 1864-1877.
- RAMBAUT A, DRUMMOND AJ, XIE D, BAELE G & SUCHARD MA. 2018. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst Biol 67(5): 901-904.
-
R CORE TEAM. 2024. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Version 4.4.1. Available at: https://www.R-project.org/
» https://www.R-project.org/ - RODRÍGUEZ-SAN PA & ALLENDES JL. 2016. Nuevos registros y extensión del rango geográfico latitudinal de Tadarida brasiliensis (Chiroptera: Molossidae) en Chile. Mastozool Neotrop 23(2): 567-569.
- ROMANO MC, MONTANI ME, CORDINI MC & AUIL S. 2015. First record of albinism in Tadarida brasiliensis (Chiroptera: Molossidae) in South America and new records of leucism in central Argentina. Chiropt Neotrop 21(1): 1312-1319.
- RONQUIST F, TESLENKO M, VAN DER MARK P, AYRES DL, DARLING A, HÖHNA S, LARGET B, LIU L, SUCHARD MA & HUELSENBECK JP. 2012. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol 61(3): 539-542.
- ROZAS J, FERRER-MATA A, SÁNCHEZ-DELBARRIO JC, GUIRAO-RICO S, LIBRADO P, RAMOS-ONSINS SE & SÁNCHEZ-GRACIA A. 2017. DnaSP 6: DNA Sequence Polymorphism Analysis of Large Datasets. Mol Biol Evol 34: 3299-3302.
- RUSSELL AL, COX MP, BROWN VA & MCCRACKEN GF. 2011. Population growth of Mexican free-tailed bats (Tadarida brasiliensis mexicana) predates human agricultural activity. BMC Evol Biol 11: 88-97.
- RUSSELL AL & MCCRACKEN GF. 2006. Population genetic structuring of very large populations: The Brazilian free-tailed bat Tadarida brasiliensis In: AKBAR Z ET AL. (Eds), Functional and Evolutionary Ecology of Bats, New York, Oxford University Press, New York, USA, p. 227-47.
- RUSSELL AL, MEDELLÍN RA & MCCRACKEN GF. 2005a. Genetic variation and migration in the Mexican free-tailed bat (Tadarida brasiliensis Mexicana Mol Ecol 14(7): 2207-2222.
- RUSSELL AL, TURMELLE AS, BROWN VA & MCCRACKEN GF. 2005b. Extremely variable di- and tetranucleotide microsatellite loci in Brazilian free-tailed bats (Tadarida brasiliensis). Mol Ecol Notes 5: 669-671.
- SCHWARTZ A. 1955. The status of the species of the brasiliensis group of the genus Tadarida J Mammal 36: 106-109.
- SIMMONS NB. 2005. Order Chiroptera. In: WILSON DE & REEDER DAM (Eds), Mammal species of the world: a taxonomic and geographic reference. 3rd ed., Baltimore: Johns Hopkins University Press, Maryland, USA, p. 312-529.
- SLATKIN M. 1993. Isolation by distance in equilibrium and non-equilibrium populations. Evolution 47: 264-279.
- SLATKIN M. 1995. A measure of population subdivision based on microsatellite allele frequencies. Genetics 139: 457-462.
- SPEER KA, PETRONIO BJ, SIMMONS NB, RICHEY R, MAGRINI K, SOTO-CENTENO JA & REED DL. 2017. Population structure of a widespread bat (Tadarida brasiliensis) in an island system. Ecol Evol 7(19): 7585-7598.
- STADEN R. 1996. The Staden sequence analysis package. Version 2.0. Mol Biotechnol 5: 233-241.
- TAJIMA F. 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123: 585-595.
- TAMURA K, STECHER G & KUMAR S. 2021. MEGA11: Molecular Evolutionary Genetics Analysis version 11. Mol Biol Evol 38(7): 3022-3027.
- VÁZQUEZ-LÓPEZ M, RAMÍREZ-BARRERA S, TERRONES-RAMÍREZ A, ROBLES-BELLO S, DE OCA A, RUEGG K & HERNÁNDEZ-BAÑOS B. 2024. Biogeographic factors contributing to the diversification of Euphoniinae (Aves, Passeriformes, Fringillidae): a phylogenetic and ancestral areas analysis. ZooKeys 1188: 169-195.
- VENCES M ET AL. 2021. iTaxoTools 0.1: Kickstarting a specimen-based software toolkit for taxonomists. Megataxa 6: 77-92.
- VILLA BR & COCKRUN EL. 1962. Migration in the guano bat Tadarida brasiliensis mexicana (Saussure). J Mammal 43: 43-64.
- WEIR BS & GOUDET J. 2017. A unified characterization of population structure and relatedness. Genetics 206: 2085-2103.
- WICKHAM H, CHANG W & WICKHAM MH. 2016. Package ‘ggplot2’. Create elegant data visualizations using the grammar of graphics. Springer-Verlag New York. ISBN 978-3-319-24277-4.
- WICKHAM H ET AL. 2019. Welcome to the Tidyverse. JOSS 4(43): 1686. Available at: https://tidyverse.tidyverse.org/articles/paper.html.
- WILKINS KT. 1989. Tadarida brasiliensis Mamm Species 331: 1-10.
- WILKINSON GS & CHAPMAN AM. 1991. Length and sequence variation in evening bat D-loop mtDNA. Genetics 128(3): 607-617.
- WILKINSON GS, MAYER F, KERTH G & PETRI B. 1997. Evolution of repeated sequence arrays in the D-loop region of bat mitochondrial DNA. Genetics 146(3): 1035-1048.
- YANG Z & RANNALA B. 2006. Bayesian estimation of species divergence times under a molecular clock using multiple fossil calibrations with soft bounds. Mol Biol Evol 23 (1): 212-226.
- ZEGARRA O, PACHECO J & PACHECO V. 2020. Distributional patterns of the Brazilian free-tailed bat Tadarida brasiliensis in the Peruvian territory. Therya 11(3): 495-507.
Edited by
-
Handling editor
Alexander Kellner
All data are available in the Supplementary Material and at https://github.com/jlncdr/Dloop.











