Tadarida brasiliensis, the Brazilian free-tailed bat, is an insectivorous molossid with a wide distribution in the Americas. It occurs in different ecosystems and uses varied shelters, from caves and crevices to human constructions, such as roofs and ceilings. Despite its wide distribution, there are several sampling gaps that make it difficult to identify the regions where the species occurs. This is a particular problem for the subspecies T. brasiliensis brasiliensis in South America, a region with few studies in comparison to North America. Considering these problems involved with identifying the distribution of T. b. brasiliensis in South America, we inferred its distribution based on 121 confirmed occurrences for the subspecies. We created a species distribution model (SDM) using the ensemble approach from the combination of BIOCLIM, SVM, GLM and MaxEnt algorithms. The resulting model suggested that the subspecies is unlikely to occur in the Amazon region and has a positive affinity with human population density, topography, a lower vegetation index, and the precipitation in the driest month. Our results show there is a large continuous area suitable for T. b. brasiliensis in central and eastern South America, with interruptions and narrow areas toward Central America. The population in this last area is separated from a smaller site in Chile by Andean deserts, snowy peaks, and high-altitude points. Our results demonstrated that along its distribution suitable habitat for T. b. brasiliensis is not continuous. The discontinuities in populations require further investigation to determine if there are phylogeographic consequences for the species.

Brazilian free-tailed bat, potential distribution, species distribution modeling, subspecies

The Brazilian free-tailed bat Tadarida brasiliensis (I. Geoffroy 1824) is an insectivorous molossid with wide distribution in North, Central, and South America (Reis et al. 2007; Barquez et al. 2015; Rodríguez-San and Allendes 2016). It occurs in different phytophysiognomies and varied shelters including caves and human structures, such as roofs, ceilings, and crevices (Schwartz 1955; Fabián and Gregorin 2007; Zegarra et al. 2020). The species can occur in colonies of hundreds to millions of individuals, with different colonies having different feeding, foraging, and migration habits (Schwartz 1955; Russell and McCracken 2006; Speer et al. 2017). Nine subspecies have been proposed for T. brasiliensis (Schwartz 1955): T. brasiliensis brasiliensis, that occurs widely in South America (Schwartz 1955; Russell 2003); T. b. cynocephala, that occurs in the southwestern and western United States, and throughout most of Mexico (Russell and McCracken 2006); T. b. mexicana, which occurs from Mexico to the western United States (Russell and McCracken 2006); T. b. intermedia, occurring in Central America, from Southern Mexico to Panama (Russell and McCracken 2006), and five subspecies occurring on 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).

Tadarida b. brasiliensis occurs from Peru (Zegarra et al. 2020) to mid Argentina (Patagonia), which is currently the southernmost occurrence known for the subspecies (Barquez et al. 2013). In the potential geographic distribution presented by Escobar et al. (2015), the most suitable regions identified for the occurrence of the subspecies were those related to the Andean ecosystems on both sides of the Andes Mountain Chain. Argentina and Chile were considered to contain the most suitable environments for the occurrence of T. b. brasiliensis (Escobar et al. 2015). However, suitable areas were also indicated in Peru close to the Andes. In a study focusing only on Peru, Zegarra et al. (2020) suggested there is a wide area suitable for T. b. brasiliensis on both east and west slopes of the Andes Mountain Chain and from sea level to 4105 m of altitude.

In eastern South America, the currently recognized geographic distribution was inferred by the interpolation of the most external points at which the species had been previously recorded (Barquez et al. 2015). However, in countries such as Brazil where large areas have not been sampled for bats, this methodology is biased towards the areas in which sampling was frequent (Delgado-Jaramillo et al. 2020). The available records of the subspecies are from areas geographically close to bat research groups or in urbanized and industrial regions that have been surveyed to create a fauna inventory as part of environmental licensing processes (Rondinini et al. 2006; Brito et al. 2009; Fernández and Nakamura 2015). In addition, occurrence records are often inaccurate due to misidentification of specimens, mainly with Geoffroy’s free-tailed bat, Nyctinomops laticaudatus (É. Geoffroy 1805) (Zortéa and Taddei 1995).

For species with a wide distribution it is often difficult to determine the real area of occurrence (Elith and Leathwick 2009; Rodrigues 2011). For T. b. brasiliensis this is particularly difficult due to taxonomic identification and sampling problems. Furthermore, the tendency of T. b. brasiliensis to form colonies of greatly varying size throughout South America, and its punctual presence in tropical regions, adds to the difficulty (Fabián and Gregorin 2007). Deficient sampling may lead to distribution polygons which are insensitive, or to identifying the species range in unsuitable areas or beyond its ecological limits (Elith and Leathwick 2009; Rodrigues 2011).

In this context, species distribution modeling (SDM) is a powerful tool to reduce both: false absences or presences in estimated species distributions (Rondinini et al. 2006; Bernard et al. 2014). SDM can indicate the potential distribution of a species based on occurrences and environmental variables, indicating suitable areas for the potential presence of a given species (Lobo et al. 2008; Phillips and Dudík 2008; Giannini et al. 2012). In response to the problems involved with identifying the distribution of T. b. brasiliensis in the South American continent, we here infer the predictive geographic distribution of subspecies throughout the South America continent and indicate the variables with the greatest influence in its potential distribution.

We used South America as our study area (Fig. 1) considering the previously known distribution and dispersal potential of T. b. brasiliensis (Schwartz 1955; Russell 2003). Based on Russell et al. (2005) and Morales et al. (2018) we chose to model only T. b. brasiliensis. This subspecies has a clear differentiation in relation to the subspecies of the northern hemisphere through microsatellite DNA and mtDNA analysis. The authors also suggest a taxonomic review of the South American genetic group in relation to other groups. In the SDM, we used bibliographic records of colonies, combined with the geographic location of sequences deposited in GenBank and data from specimens deposited in collections (GenBank 2004; RELCON 2010; Rodríguez-San and Allendes 2016; Nuñes et al. 2018; Lobato 2019; Zegarra et al. 2020). These data were combined with craniomandibular analysis of skull samples and molecular identification of tissues deposited in scientific collections in South America (D-loop gene, Amaral et al. unpublished data), totaling 136 locations. The coordinates at which each specimen was collected were compiled from the specimen identification tags. However, in the few cases where this information did not exist, the coordinates of the administrative headquarters of the municipalities where the specimen was collected were used.

Figure 1. 

Locations used in the species distribution model (SDM) of the Brazilian free-tailed bat, T. b. brasiliensis, in South America.

The geographic coordinates of the records were verified using QGIS version 3.16 (QGIS 2020). After that, we use a spatial filter the CoordinateCleaner R package (Zizka et al. 2019) to avoid sampling effort problems (Veloz 2009; Stevens 2013). To work around the error of spatial autocorrelation and sampling bias, we excluded the coordinates with a distance of less than 10 km using the spThin R package (Aiello-Lammens et al 2015). The final matrix of the dataset resulted in 121 coordinates of records of the subspecies T. b. brasiliensis (Fig. 1, Suppl. material 1).

The extension of all variables used was adjusted to South American territorial boundary using the raster R package (Hijmans and Etten 2012; R Core Team 2021), through the crop, mask and aggregate functions (Hijmans et al. 2023). To build the models, 24 variables were obtained initially with 5’ arc-minutes of spatial resolution (~10 km²) and, which are related to the occurrence of the subspecies (Schwartz 1955; Fabián and Gregorin 2007; Barquez et al. 2015; Lobato 2019; Zegarra et al. 2020). Nineteen bioclimatic variable layers from the dataset are available in WorldClim v. 2.1 (Fick and Hijmans 2017; https://www.worldclim.org/) and five other variables: land cover, human population density, vegetation index, land surface temperature at night and topography were obtained from the Neonasa system (https://neo.gsfc.nasa.gov/) (Suppl. material 2).

We select the least correlated variables for the SDM using the Variance Inflation Factor (VIF) with a value of < 0.7 (Dormann et al. 2013) using the “vifcor” function of the usdm R package (Naimi et al. 2014; Naimi and Araújo 2016; R Core Team 2021). The selection indicated 10 less correlated variables, which were used in the SDM: mean diurnal range (mean of monthly (max temp - min temp)) (var02), isothermality ((BIO2/BIO7) (x100)) (var03), precipitation of driest month (var14), precipitation seasonality (coefficient of variation) (var15), precipitation of warmest quarter (var18), precipitation of coldest quarter (var19), human population density (var21), vegetation index (var22), land surface temperature at night (var23), and topography (var24).

In the present study, in addition to presence data, pseudo-absence data were used (VanDerWal et al. 2009; Ribeiro-Souza et al. 2022). Pseudo-absence points were randomly generated equal to number to presence points throughout South America, except for the pixels that contained presence data, through the argument “method = gRandom” in sdm R package (Naimi and Araújo 2016; R Core Team 2021). For modeling the SDM we use the ensemble of forecasts approach (Araújo and New 2007). The predictive map based on the species distribution model (SDM) was generated using four algorithms: climate-envelope-model - BIOCLIM (Nix 1986), Support Vector Machines - SVM (Vapnik et al. 1995), Maximum Entropy - MaxEnt (Phillips et al. 2006), and Generalized Linear Models - GLM (McCullagh and Nelder 1989). BIOCLIM is based only on presence records (Pereira et al. 2020), the SVM method is based on data including points at which T. b. brasiliensis is known to be present, combined with pseudo-absence, weather and environmental information to create a new hyperspace, in which complicated patterns can be represented in a less complex way (Drake et al. 2006). The GLM was fitted using presence and pseudo-absences (Segurado and Araújo 2004). MaxEnt uses the principal of maximum entropy to predict species distributions taking into account the fundamental niche of the species, and correlating presence and background data with environmental variables (Martínez and Di Cola 2011). The algorithms were applied with different techniques using the dismo (Hijmans et al. 2011), raster (Hijmans and Etten 2012) and sdm (Naimi and Araújo 2016) R packages (R Core Team 2021).

To calibrate and evaluate the models, we used the dplyr (Wickham et al. 2023), raster (Hijmans and Etten 2012) and dismo (Hijmans et al. 2011) R packages (R Core Team 2021). Occurrence locations were divided into two subsets (training and test) that contained 70% and 30% of occurrence locations, respectively (Ribeiro-Souza et al. 2022). This process was repeated 20 times for each algorithms combination, totaling 80 different predictions (20 randomizations × 4 algorithms) to estimate the potential distribution of T. b. brasiliensis.

The performance of each model was evaluated by the estimated values of Area Under the Curve (AUC), which calculates the probability of correct classification (Deleo 1993; Fielding and Bell 1997; Shabani et al. 2018) and, True Skill Statistic (TSS) for threshold that maximizes the sum of sensitivity and specificity (Max SSS) (Liu et al 2013), which is the number of correct classifications minus those assigned randomly (Allouche et al. 2006; Shabani et al. 2018). The value of AUC varies between 0 and 1, where values close to 1 indicate that the data were well distinguished, there is no overlap between presence and pseudo-absence data. If the value is 0.5 or < 0.5, it indicates overlap between the data (Fielding and Bell 1997). TSS values are calculated using sensitivity and specificity and, vary between -1 to + 1. Values close to 1 indicate excellent models and values zero or < 0 indicate a prediction performance random (Allouche et al. 2006). Only the predictions that had AUC and TSS values greater than the average of all the adjusted models generated by the four algorithms were used and combined in the ensemble. Lastly, we created ensembles with “weighted averages” by TSS values (Ribeiro-Souza et al. 2022) previously computed using the sdm R package (Naimi and Araújo 2016; R Core Team 2021).

Finally, this method resulted in a map being generated with suitable areas for the distribution of T. b. brasiliensis, where values close to 1 mean suitable areas for the occurrence of the subspecies. Finally, we used the threshold that maximizes the sum of sensitivity and specificity (Max SSS) (Liu et al. 2013). This resulted in a binary map (presence = 1 and absence = 0), in which the red color represents the places of presence.

The distribution model generated for T. b. brasiliensis in South America presented a high value of AUC (0.90 ± 0.06) and TSS (0.79 ± 0.09). The most important variables in terms of contribution to the suitability model (Suppl. material 3) were: isothermality (var03 - 35%), human population density (var21 - 13.1%), topography (var24 - 7.4%), precipitation in the driest month (var14 - 5.7%), and vegetation index (var22 - 5.1%).

Fig. 2 shows the predicted distribution maps generated, where the orange to red colors (≥ 0.6) in map A indicated the regions with the highest suitability and, red color in map B indicate the regions with presence for T. b. brasiliensis in South America. The response curves (Suppl. material 4) showed a greater affinity with regions where the isothermality (var03) variation is negative. In turn, according to the ensemble model, the species has a positive affinity with human population density (var21), topography (var24), lower vegetation index (var22), and with the precipitation in the driest month (var14). It was also possible to observe a negative relationship with precipitation in the warmest quarter (var18), with greater affinity for neutral values for this variable.

Figure 2. 

Suitable areas for the occurrence of T. b. brasiliensis in South America A “warmer” colors (yellow to red) indicate higher suitability of subspecies occurrence B binary map - red color indicates presence locations of subspecies. Black dots are localities of occurrence.

In Brazil, the SDM showed high environmental suitability for T. b. brasiliensis in the south and southeast regions. This presents a wide and continuous area suitable for the occurrence of this subspecies that extends to neighboring countries, reaching the territory of Peru in its northern portion. In Uruguay, a wide distribution was observed throughout the country. However, in Argentina T. b. brasiliensis is more likely to occur in the north, northeast and west, with no suitable regions for the occurrence of the subspecies in the south of the country. According to the model, in Chile, the subspecies has a wide area of habitat suitable for its occurrence, with regions separated to the east and southwest by the barrier formed by the Andes Mountain chain. Furthermore, according to the model, the subspecies is likely to occur along the east face of the Andes, forming a corridor with strong narrowing and interruptions in its northern portion along the borders of Ecuador, Colombia and Venezuela. It was also evident that the subspecies is less likely to occur in the Amazon region.

The model generated for the potential distribution of T. b. brasiliensis showed that the subspecies has the potential to occupy a wide area in most regions of South America, being less likely in the Amazon and Cerrado biomes. Areas suitable for the subspecies include the south and southeast of Brazil, all of Uruguay, and the north, west and northeast of Argentina. In addition to these regions, a region on the east side of the Andes was evidenced, which seems to be a corridor towards Central America, which includes part of Bolivia, Paraguay and Peru. A suitable region for the subspecies was also evident on the western side of the Andes, which includes Chile.

The influence of environmental factors on the presence and activity of bats has been reported in several studies, especially air temperature and relative humidity (Lacki 1984; Maier 1992; Adam et al. 1994; Hayes 1997; Russ et al. 2003). Our results point to these as the most significant factors defining the potential distribution of T. b. brasiliensis in South America. In previous studies on T. brasiliensis in the northern hemisphere, these environmental characteristics were found to be determinants for the activity of bats, but they were also related to the availability of food resources (insects) (Akasaka et al. 2009; Wang et al. 2010).

In addition, in the present study, a negative relationship with precipitation was found in the hottest period of the year (austral summer). This result could be related to the reduction of bat activity on rainy nights, probably due to the difficulty imposed by the reduced efficiency of echolocation and the effects of rain on thermoregulation (Baskaran et al. 2016; Russo and Voigt 2016). The reduction in bat activity in the summer could cause an energy deficit precisely in the period of greatest energy demand for the species, during pregnancy, lactation, and the first flights of the young (Wilkins 1989; Marques and Fabián 1994). In Molossus molossus (also an insectivore species) it was observed that adult females, out of the reproductive period, and after 48 hours of fasting,tend to significantly decrease the levels of glucose in the organism in comparison to male specimens (Goulart 2008; Freitas et al. 2010). In addition, it was verified in females of another insectivorous species (Eptesicus fuscus) that, after an event of energetic deficit during pregnancy and lactation, negative effects occur in the female’s organism, reflecting directly on the offspring, with weight loss and potential less survival capacity (Hamilton and Barclay 1994).

The SDM also suggests that the variable human population density is related to the occurrence of suitable places for the subspecies to roost, thus corroborating the habit of the subspecies to use human structures as daytime shelters (Marques and Fabián 1994). Another variable positively related to habitat suitability that reflected the behavior of the subspecies was the areas of undergrowth or mosaic of forest fields, a characteristic explained by the morphology and behavior of the Molossidae family, which has adaptations for foraging in open environments (Schnitzler and Kalko 2001; Jung et al. 2014).

Our results do not corroborate the predictive distribution proposed by Escobar et al. (2015) for T. b. brasiliensis. The model generated by the authors showed a wide region of inadequacy for the occurrence of the subspecies in eastern South America, mainly in Brazil and Uruguay, countries with confirmed occurrence of the subspecies. In the SDM, information on the occurrence of all subspecies of T. brasiliensis were used, without considering the differences in the variation in the niches of the genetic groups. According to Morales et al. (2018) there is a relationship between the environmental niche variation and the genetic groups of T. brasiliensis that are influenced by the climatic differences between the Nearctic, Neotropical and Antilles regions. Furthermore, no spatial filter was used to avoid autocorrelation or sample bias in the presence of the species. In a previous study, it was shown that the spatial autocorrelation of presence records can lead to a lack of independence between data sets, probably causing bias in model predictions (Veloz 2009).

In our study, among the regions in eastern South America identified as being unsuitable for the subspecies, is the state of Bahia in Brazil. There were records of T. b. brasiliensis for this state (Barquez et al. 2015), but we found no specimens for the region in any collection analyzed in the present study. After reviewing the external dental and morphological characteristics of some individuals (MNRJ6556, MNRJ6557, MNRJ6558, MNRJ6562, MNRJ6564, MNRJ6565, and MNRJ6566) identified as T. b. brasiliensis, it was found that the records for Bahia were incorrect, because they were all specimens of N. laticaudatus. In this context, we emphasize the need to review the real northeastern limit of the subspecies in Brazil, mainly through the search for new colonies and shelters. This is necessary since our records and SDM point to the state of Minas Gerais as the northeastern limit for the subspecies in country.

The species distribution model generated in the present study suggests that the main discontinuities in the distribution of the T. b. brasiliensis, besides the Amazon and Cerrado biomes, are in the equatorial region, and the high altitude areas of the Andes Mountain Chain. Barquez et al. (2015) also identified the Atacama Desert as an unsuitable habitat within the subspecies distribution polygon, probably due to its extreme environmental conditions. Although there are some records of subspecies in the surrounding arid areas, this evidences that the T. b. brasiliensis occurs in specific environments in the biomes. Our results reaffirm this, highlighting the need to consider information about the heterogeneity of environments and formations to describe the regions (microecosystems) where the subspecies occurs in future studies on the distribution of T. b. brasiliensis.

Moreover, our results also did not indicate the environmental suitability of the extreme south of Argentina for the distribution of the subspecies. Possible conditions resulting from climate change may expand the environmental suitability of this region where so far there are a few records of solitary individuals (Zapata et al. 2015; Hill et al. 2016; Teixeira 2019). Some studies have indicated that the area of habitat suitability for many species tends to increase towards the poles as a response to future scenarios where a reduction in environmental stress caused by cold may occur due to climate change. In contrast, the environmental suitability of equatorial regions may be reduced (Hill et al. 2016; Qin et al. 2019).

The barrier formed by the Andes Mountain Chain in some areas creates a substantial barrier to the east/west flow of many species, including bats (Koopman 1978; Patterson et al. 2012). Although T. brasiliensis is a species capable of overcoming major geographical barriers through flight (Ditchfield 2000; Russell and McCracken 2006), its main impediment to gene flow is still environmental discontinuities (Wiens and Donoghue 2004). These barriers reinforce the need for phylogenetic studies that make it possible to understand the influence of these discontinuities on the genetic diversity of the species and its taxonomic implications.

The results of the present study support the need to integrate taxonomic and population studies with studies on species distribution. Although T. brasiliensis is a molossid widely distributed across the American continent (Reis et al. 2007; Barquez et al. 2015; Rodríguez-San and Allendes 2016), and is considered a common bat in much of its distribution (Barquez et al. 2015), there are still many questions related to the nine subspecies, mainly taxonomic (Russell 2003; Morales et al. 2018). T. b. brasiliensis is currently the only one with taxonomic support (Morales et al. 2018), however it is one of the least studied subspecies. Our study showed that it is poorly sampled in some regions, especially those that are not close to study groups with mammals, and confused with N. laticaudatus, making the veracity of information on the presence of the subspecies even more difficult. We suggest caution with the use of occurrence databases available on the internet, as many of the individuals reviewed during this study were mistakenly identified as T. brasiliensis and appeared on different platforms. We recommend that future studies on the distribution of the species first consider genetic groups and use reliable presence information, such as colony records and specimens with genetic information since environmental differences between genetic groups must be considered prior to SDM.

We are grateful to the curators and technicians who sent us tissue samples or welcomed us into their collections during the Sars-Cov-19 pandemic period: Dr. Nilton Cáceres from the Mammal Collection of the Federal University of Santa Maria (CMUFSM), Dr. Marcelo Weksler and Dr. João Oliveira from the Mammals Collection of the National Museum of Rio de Janeiro (CMMNRJ), Dra. Márcia Jardim from the Mammals Collection of the Zoobotany Foundation of Rio Grande do Sul (CMFZB), Dr. Adriano Peracchi, Dr. Marcelo Nogueira and Biol. Karen Toledo from the Adriano Lúcio Peracchi Collection (UFRRJ), Biol. Fábio da Costa from the Natural Sciences Laboratory of the University of Caxias do Sul (MUCS), Me. Flávia da Silva e ao Biol. Élinton Rezende of the Augusto Ruschi Zoobotanical Museum (MUZAR) and Dr. Guillermo D’Elía of the Universidad Austral de Chile. In particular, we would like to thank Dr. Rubén Barquez and Dr. Monica Díaz, responsible for the Colección Mamíferos de Lillo (CML) of the Universidad Nacional de Tucumán and Dr. Pablo Teta and Dr. Sergio Lucero from the collection of the Museo Argentino de Ciencias Naturales -Bernardino Rivadavia. We also want to thank the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) by the Programa de Apoio à Pós-Graduação em Instituições de Ensino Superior Comunitárias (PROSUC) that provided the grant to IA, the Brazilian National Council for Scientific and Technological Development (Conselho Nacional de Desenvolvimento Científico e Tecnologico - CNPq) from Brazil, that provided the Produtivity grant to L.R.O. and financial resources for this research (CNPq N°. 303813/2011-3, 308650/2014-0, 310621/2017-8 and 315361/2021-2) and R.G. thanks The (CNPq N°. 304907/2019-7; 311204/2020-1), and thanks and Fundação de Amparo à Pesquisa do Estado de Minas Gerais (process PPM 00203-18) M.H.V, thanks The São Paulo Research Foundation - FAPESP (grant number 2022/01899-6). We also thank the anonymous reviewers who contributed with comments that were reflected in improvements to the final version of this article.