Phyllostomid bats distribution and richness gradient in a subtropical Brazilian state

Phyllostomid bats (Chiroptera: Phyllostomidae) are key elements for the maintenance of New World forests, but little information on their distribution is available in some regions of Brazil. Here we use occurrence records and bioclimatic variables to model the distribution of phyllostomid bats in Santa Catarina, a subtropical Brazilian state. Estimates of geographic variation in species richness were then obtained by stacking the generated maps. Lastly, we tested how associated species richness is to ecoregions and Protected Areas. Our results suggest that the phyllostomid bats species richness is closely linked to the region’s climate gradient. Most species are restricted to the Serra do Mar ecoregion, where the temperature is high and varies less throughout the year. In contrast, the colder areas seem to house extremely simple communities, composed of a subset of the species present in the warmer areas. We found significant evidence that Protected Areas in Santa Catarina play an important role in the conservation of species, although there are still several places where species richness is high, but no Protected Area is available. The creation of new Protected Areas in these places can boost the species conservation, and, consequently, the ecological services provided by phyllostomid bats.


INTRODUCTION
Phyllostomid bats (Chiroptera: Phyllostomidae) constitute the most diverse group of Neotropical bats, with 62 genera and at least 223 species described (Simmons & Cirranello, 2020). All of this taxonomic diversity is reflected in a wide variety of morphological, physiological, behavioral and dietary adaptations, which makes these bats perform several crucial ecological services to New World (sub)tropical forests, including pollination, seed dispersal, and population control of insects and small vertebrates (Muscarella & Fleming, 2007;Kunz et al., 2011;Mikich et al., 2015).
When compared to other Neotropical bat families, phyllostomid bats have its biology considerably better known due to its greater detectability through mist nets, a method commonly used to capture bats in the region (Straube & Bianconi, 2002;Kunz & Parsons, 2009). Despite this facility, knowledge about bat fauna in some regions of Brazil is still insufficient, with large areas with no information on the occurrence of bats (Bernard et al., 2011a).
In Santa Catarina, a subtropical Brazilian state, the situation is similar Passos et al., 2010). In the last decades, some researchers Cimardi, 1996;Pacheco et al., 2007;Passos et al., 2010;Cherem et al., 2011;Cherem & Althoff, 2015;Bôlla et al., 2017;Cherem & Althoff, 2019) have inventoried the species and discussed ecological aspects of them, but this information does not cover the entire state and is highly concentrated in the eastern portion. Hence, most cities have never been the object of specific research (Fig. 1), which warns for the need for studies that fill these knowledge gaps, given the possibility that species or populations will become extinct in places where they were not even known (Costa et al., 2005;Bernard et al., 2011a).
Species Distribution Models (SDMs) represent a useful tool for addressing the lack of spatially comprehensive research (Guisan & Zimmermann, 2000;Guisan & Thuiller, 2005;. In general, SDMs use georeferenced occurrence data and environmental variables (usually climatic, topographic and/or land cover), to establish a statistical relationship between species and environment. Subsequently, this relationship can be used to make predictions, allowing researchers to make inferences in space (e.g., Rebelo & Jones, 2010, in search of new populations of rare species) and time (e.g., Aguiar et al., 2016, investigating the expected impacts of climate change). In recent years, driven by advances in computing and increased open data availability, SDMs have experienced a unique advance, thus becoming an important tool in ecology and conservation (Guisan & Thuiller, 2005;. In this work, we use georeferenced occurrence records and bioclimatic variables to generate Species Distribution Models for 20 species of phyllostomid bats occurring in Santa Catarina and then answer the following questions: (1) how are phyllostomid bats distributed throughout Santa Catarina? (2) Where are the species richness hotspots? (3) How effective are the Protected Areas of the Santa Catarina in protecting phyllostomid bats? Our first hypothesis (H1) is that the majority of phyllostomid bats in the state are present, exclusively or not, in the Serra do Mar ecoregion, in a longitudinal gradient similar to that observed in other regions of the Atlantic Forest (Batista et al., 2020). Besides, in order to investigate the effectiveness of Protected Areas, we also tested the hypothesis that Protected Areas are located in places with greater species richness than would be expected at random (H2).

Study Area
Santa Catarina state ( Fig. 1) is in the southern region of Brazil, covering a total area of 95,730.684 km², inserted fully below the Tropic of Capricorn. Elevation oscillates from sea level, on the coastal plains, to areas above 1,800 m a.s.l. (Maack, 2001). According to the Köppen classification, the climate ranges from Cfa (humid temperate climate with hot summer), in the lower areas, to Cfb (humid temperate climate with temperate summer), in the highest points of the state (Alvares et al., 2014). The whole area of Santa Catarina is inserted in the Atlantic Forest, with four main ecoregions (Olson et al., 2001): Restingas and Mangroves, restricted to coastal areas; Serra do Mar Coastal Forests (Serra do Mar), located in the east of the state, between 0 and 900 m a.s.l.; Araucaria Moist Forests (Araucaria), restricted to areas with elevation above 500 m a.s.l., in the central portion of the state; and Alto Paraná Atlantic Forest (Alto Paraná), located at intermediate elevations in the west of the state, mostly composed of semideciduous forests (Klein, 1978).

Occurrence Data
We extracted the georeferenced occurrence records from specialized literature (Appendix S1), online data-  Olson et al. (2001). Russi, C.H. et al.: Phyllostomid bats in subtropical Brazil Pap. Avulsos Zool., 2021;v.61: e20216160 2/12 bases, Zoological collection of the Regional University of Blumenau (CZFURB), and unpublished records made by F. Carvalho. Data from the specialized literature were found in the Web of Science (https://webofknowledge.com), Google Scholar (https://scholar.google.com) and SciELO (https://scielo.org), using the keywords "Phyllostomidae", "bats", "Chiroptera", "Atlantic forest", "subtropical Brazil" and "southern Brazil". Data from online databases came from Brazilian Biodiversity Portal (2021; https:// portaldabiodiversidade.icmbio.gov.br), SpeciesLink (2021; http://splink.cria.org), and Global Biodiversity Information Facility (GBIF, 2021; https://www.gbif.org), and were found using the keyword "Phyllostomidae". When only the name of the municipality of occurrence of a species was available, we used the centroid coordinates of this municipality, extracted through QGIS v.3.4.7 (QGIS Development Team, 2019). Although we recognize that this may introduce uncertainty into our analyses, the high spatial autocorrelation present in the environmental variables used tends to minimize this possible problem (Naimi et al., 2011). Besides, the algorithm used to generate our models seems to be one of the most robust to deal with this type of uncertainty (Graham et al., 2008).
As few records were available for Santa Catarina, and a minimal amount is required for analyses employed, records located within a 120 km buffer from the state were also considered. This distance was chosen based on an estimate of the maximum distances traveled by phyllostomid bats in the Atlantic Forest .
We checked the coordinates of all records obtained and excluded duplicates and records located at sea. Records in which the species nomenclature was incorrect (i.e., synonyms) were corrected following the most recent version of the nomenclature presented by the Brazilian Bat Research Society (Garbino et al., 2020; https://www. sbeq.net/lista-de-especies).
In total, we found 1,621 occurrence records (considering only one record of each species per location) in the entire buffer, related to 22 genera and 31 species (Table 1). When considering only the extent of Santa Catarina, the number of records was 983, related to 16 genera and 22 species. Among these, 20 species presented at least 8 occurrence records and were considered for further analyses. An electronic spreadsheet containing all the records found and their respective references is available in the GitHub repository (https://github.com/ chrussi/Phyllostomid_bats_2021).

Environmental Variables
In order to characterize the environmental conditions that phyllostomid bats are likely to experience in the study area, we used an initial set of 19 bioclimatic variables available on the WorldClim portal (version 2.0; Fick & Hijmans, 2017) at a resolution of 30 arc-seconds (≅ 0.86 km² at the Equator), which were cut out to the same extent as the occurrence records (Santa Catarina + 120 km buffer). The WorldClim bioclimatic variables represent the average of the 1970-2000 period for temperature and precipitation trends, such as average, seasonality, and extreme conditions, which are known to influence occurrence patterns of phyllostomid bats in the Atlantic Forest (Mello et al., 2009;Stevens, 2013;Batista et al., 2020). After visually inspecting all bioclimatic layers, we chose not to consider bio 8 and bio 9 due to anomalies observed at some points in the study area. Next, since many of the WorldClim variables are redundant and collinearity may decrease the predictive performance of the analyses (Dormann et al., 2013), we perform a Principal Component Analysis (PCA) using the "prcomp" function of the "stats" R package (R Development Core Team, 2021). We extracted the first four orthogonal axes of the PCA, which represented 97% of the variation of the bioclimatic variables (Appendix S2). These orthogonal axes were used for further analysis.

Species distribution and richness gradient
We used Maxent (Phillips et al., 2004), implemented in the "maxnet" R package , to generate Species Distribution Models (SDMs) of all species for which we found at least 8 occurrence records in the study area. Maxent compares the environmental conditions occupied by the species, extracted from the occurrence records, with the environmental conditions available in the study area, extracted from background points, to then estimate species responses to the environmental gradient, which are later used to make geographic predictions (Phillips et al., 2004;Elith et al., 2011).
Since the central and western regions of our study area have historically been subject to less sampling effort than the eastern region , we took a series of precautions to minimize this sample bias: (i) first, we customize the way that Maxent samples the background points. This is usually done by generating 10,000 random points across the study area. However, this generates a bias because not all areas were sampled equally , which tends to underestimate the occurrence of species in places with low sampling. Therefore, our 10,000 background points were spatially manipulated to represent the same sampling bias present in our presence records. The procedure used the "kde2d" function of the "MASS" R package (Venables & Ripley, 2002) to create a sampling density grid − based on the records of all phyllostomid bat species − and the "sample" function of the "base" R package, to choose proportionally more background points in places with higher sampling density. (ii) Next, we also remove records located very close (< 20 km) to each of all species with more than 15 records (Fourcade et al., 2014). We used this minimum distance because we noticed that it was enough to disaggregate clustered records without excessively reducing the number of records available; (iii) finally, we restrict the transformations that Maxent does on variables to only linear and quadratic. This implies simpler models, which give less importance to the noise caused by the sampling bias (Merow et al., 2014).
To test Maxent's ability to predict species presence/ absence in unsampled areas, we used a geographically structured cross-validation procedure with four iterations. For this, we use the "blockCV" R package (Valavi et al., 2019) to divide the presence records of each species and background points into four geographically distinct groups (Radosavljevic & Anderson, 2014). Thus, in each iteration, three groups were used to calibrate the models and one group was used to validate it. As measures of model performance, we use two distinct metrics: (i) Area Under the Curve (AUC) of the receiver operating characteristic curve (ROC; Fielding & Bell, 1997), a discrimination metric that measures the likelihood that a point of presence will receive a higher score on the models than a background point. Models with AUC values above 0.7 are considered good; and (ii) Continuous Boyce Index (CBI; Hirzel et al., 2006), a calibration metric that assesses how much model predictions differ from a random distribution of observed presences across prediction gradients. It varies from -1 to 1, with values close to 1 indicating good models. AUC and CBI were measured using the "SDMTools" (VanDerWal et al., 2014) and "ecospat" (Di Cola et al., 2017) R packages, respectively.
Subsequently, to obtain species richness gradient estimates, we cut out the continuous climate suitability maps generated by Maxent for the Santa Catarina extension and overlapped them, according to the approach recommended by Calabrese et al. (2014). The procedure was performed in R (R Development Core Team, 2021).

Hypothesis testing
We use Generalized Linear Models (GLM; Mccullagh & Nelder, 1989) with Poisson distribution to test our hypotheses. For the hypothesis that the Serra do Mar is the ecoregion where there are the highest values of species richness (H1), we generate 1,000 random points in each ecoregion, extract the values of the species richness predicted for each point, and execute the GLM, using the species richness as the response variable and ecoregions as a categorical explanatory variable. Since the "Restingas and Mangroves" ecoregion occupies a small part of the study area, we exclude it from this analysis. To test the hypothesis that the Protected Areas of the state are effective in protecting phyllostomid bats (H2), we generate 1,000 random points inside and 1,000 random points outside the Protected Areas, and then run the GLM using the species richness as the response variable and the origin of the point (inside or outside the Protected Area) as a categorical explanatory variable. We consider both Protected Areas with all levels of protection (Integral Protection and Sustainable Use) and Protected Areas with Integral Protection only. The analyses were based on the map of Protected Areas of Santa Catarina, made available by the Brazilian Ministry of the Environment (http://mapas.mma.gov.br/i3geo/datadownload.htm#). GLM was implemented in R (R Development Core Team, 2021). Both hypotheses (H1 and H2) were accepted only if the GLMs resulted in p-values below 0.05.

RESULTS
All SDMs had a satisfactory predictive performance (mean AUC = 0,784, sd = 0.06; mean CBI = 0.724, sd = 0.17), indicating a good fit to the data (for more details, see Appendix S3). The maps generated for all species are shown in Figs. 2 and 3. Except for C. auritus, D. rotundus, and S. lilium, most species were restricted to low elevation areas (< ≅ 900 m a.s.l.). We did not find any species restricted to high and cold areas, so the species abovementioned were also frequent in areas of lower elevation.  The stacking of species maps resulted in a strong gradient of species richness (Fig. 4). Estimates ranged from 3 species, in the highest and coldest areas of the Araucaria (dark gray regions on the map, Fig. 4), to 17 species (red regions on the map, Fig. 4), in the north of the Serra do Mar. Indeed, our GLMs supported the H1 that species richness is significantly greater in the Serra do Mar than in any other ecoregion (p-value < 0.001 when compared to Alto Paraná and Araucaria; Fig. 5). Besides, the H2 was also supported by GLMs (p-value < 0.001), so species richness is greater in Protected Areas than outside them, both for the Protected Areas with integral protection and sustainable use, and only those with integral protection.

DISCUSSION
Our analyses provided unprecedented high-resolution estimates for the distribution and geographic variation in the species richness of phyllostomid bats in Santa Catarina. Before, the only works that addressed the distri-bution of this family in the state were restricted to listing the cities Bôlla et al., 2017) or river basins (Pacheco et al., 2007) in which each species was registered. Besides, both hypotheses tested in this work were accepted. Therefore, our models predicted a greater species richness in the Serra do Mar than in other ecoregions, as well as a consequent greater species richness inside the Protected Areas than outside them, even when considering only Protected Areas with integral protection.
Despite being climatically heterogeneous and rich in other taxa, the subtropical Atlantic forest has lower levels of phyllostomid bats richness and endemism than rainforests located close to the equator, such as Amazonia (Bernard et al., 2011b), or even the central and northern portions of the Atlantic Forest (Stevens, 2013;Batista et al., 2020;Stevens et al., 2020). This has been attributed to the high elevations and latitudes observed and its consequent climatic conditions, which limit the occurrence of many species (Pereira & Palmeirim, 2013;Stevens, 2013;Batista et al., 2020). The same seems to be true in the Santa Catarina, where many species appear to occur only in the Serra do Mar, especially in the coastal plains of the northern part of the state (Fig. 4), where the temperature is higher and varies less throughout the year (Wrege et al., 2012).
The areas of Alto Paraná, although they have a hot and humid climate for several months of the year, have winters with temperatures low enough to limit several tropical plant species (Wrege et al., 2012;Gasper et al., 2015). These conditions also seem to limit the occurrence of some species of phyllostomid bats, thus boosting the longitudinal gradient observed both in Santa Catarina and in other regions of the subtropical Atlantic Forest (Batista et al., 2020). The situation is even more drastic in the Araucaria, where temperatures close to 0℃ are often observed in winter (Wrege et al., 2012) and the areas with the fewest species number have been identified.
Although we paid particular attention to the sampling bias (see Material and Methods), the high species richness in the Serra do Mar raises the question of how biased our results are, since it is precisely at these locations that the largest number of records were available (Fig. 1). However, we emphasize that a similar longitudinal gradient of species richness is also observed in bat communities from other regions of the Atlantic Forest (Batista et al., 2020;Delgado-Jaramillo et al., 2020), as well as in other typically tropical taxa, like angiosperms (Vibrans et al., 2020) and anuran amphibians (Vasconcelos et al., 2019). Thus, we believe that our results reflect with good accuracy the true gradient of species richness presented by the phyllostomid bat fauna of Santa Catarina.
C. auritus, D. rotundus and S. lilium seem to be the only phyllostomid bats to occupy most of the coldest areas of Santa Catarina (Fig. 2-3), which is a pattern similar to that of other regions of the Atlantic Forest (Cherem & Althoff, 2015;Carvalho et al., 2019;Miranda et al., 2019). Like Patterson et al. (1996), Carvalho et al. (2019), and Batista et al. (2020), our results suggest that the high and cold areas have a highly nested community of phyllostomid bats. In Santa Catarina, there are no species restricted to cold areas, and the three species abovementioned rep-  . Phyllostomid bats richness gradient predicted for Santa Catarina, Southern Brazil, and its association with Protected Areas. The map was obtained by stacking maps of species potential distribution generated through the Maxent algorithm. Russi, C.H. et al.: Phyllostomid bats in subtropical Brazil Pap. Avulsos Zool., 2021;v.61: e20216160 7/12 resent a subset of the species present in the lower and warmer areas. Hence, minimum extreme temperatures may act as a filter, selecting only a few phyllostomid bat species with physiological tolerance to these conditions (Stones & Wiebers, 1965;Stevens, 2013), or even restricting the occurrence of the organisms consumed by these bats, which would also make it impossible for the species analyzed to occur in the higher and colder areas.
Our results provide insights for species conservation by indicating the most species-rich areas -which should be the focus of actions aimed at the conservation of phyllostomid bats of Santa Catarina -, as well as demonstrating the importance of maintaining the Protected Areas existing in the region. In this sense, the Parque Estadual do Acaraí (Santa Catarina, 2005) seems to play a key role in the conservation of Santa Catarina's phyllostomid bat communities, as it is the main Protected Area with integral protection in the northeast of the state, in an area where at least 17 species seem to co-occur. The other Protected Areas in the northeast of Santa Catarina have a small extension or belong to the category of sustainable use (Fig. 4), which tends to be less restrictive and more susceptible to illegal deforestation. The situation is similar in the south of the state, where more than 12 species were predicted to occur in most locations and few Protected Areas with integral protection are available. The creation of new UCs in these regions seems to be important to protect the phyllostomid bats, especially when considering the intense degradation to which the Atlantic Forest has been subjected (Ribeiro et al., 2009).
The data collection resulted in only seven species for the Alto Paraná areas, although richness estimates indicate the presence of at least nine species in some places in the region (Fig. 4). Species such as Carollia perspicillata, Micronycteris megalotis, Diphylla ecaudata, and Platyrrhinus lineatus showed medium climate tolerance to these regions but were never registered there. Some species that depend on a less altered landscape may no longer occur in these areas due to the low number of well-preserved forest remnants available (see Appendix S4). In this context, the restoration of these environments may be important for these species to return to occupy the forests of the region.
Although our results demonstrate the effectiveness of Protected Areas, the species richness was relatively low in the three main Protected Areas with integral protection in Santa Catarina (Fig. 4): Parque Nacional da Serra do Itajaí (Brasil, 2006), Parque Estadual da Serra do Tabuleiro (Santa Catarina, 2009) and Parque Nacional de São Joaquim (Brasil, 1961). This is because these Protected Areas are located in mountainous regions, where the richness of phyllostomid bats tends to decline. However, these areas can be crucial for the conservation of other bat families typical of temperate regions, such as Vespertilionidae and Molossidae. Histiotus montanus (Vespertilionidae), for example, has been recorded exclusively in the coldest areas . Indeed, some studies have shown the existence of a greater richness of insectivorous bats in cold areas of South America  Miranda et al., 2019), but more fieldwork is still needed for accurate distribution estimates for these species in these areas.

CONCLUSION
Our results demonstrate that phyllostomid bats are distributed non-randomly in Santa Catarina and that bioclimatic variables, and, consequently, ecoregions, represent a good proxy for estimates of species distribution and richness. Based on these estimates, we show that Protected Areas are important for maintaining taxonomic diversity and ecological services provided by these species, although there are opportunities for improvement.

ACKNOWLEDGMENTS
The authors would like to thank Luís O.M. Giasson, Liu I. Orozco for the relevant suggestions made during the development of this work. We also thank Fábio L.V. Bones and Guilherme S. Grittz for reviewing the manuscript. This research did not receive any specific grant from funding agencies in the public, commercial, or notfor-profit sectors.