Picky City Bats: Variation in Urban Conditions Generates a Suitability Gradient that Structures the Bat Community.

The transformation of natural systems into urban cities represents a radical habitat modication and one that generates species-specic effects on bat communities. Cities present heterogeneous mosaics of urban conditions, which the bats can inhabit differentially by utilizing their intrinsic characteristics. We expected bat species richness, number of abundant species and feeding guilds to be positively impacted by an increased degree of suitability for bats in certain areas interspersed in the city. We also expected that the number of dominant species would follow an inverse pattern, with some species being fostered in less suitable areas for bats. We determined the levels of suitability for the presence of bats in urban-based basic geostatistical units (BGUs) in the Neotropical city of Morelia, in Michoacán state, Mexico. We identied three levels of suitability for bats, low, medium and high, using the percentage of arboreal vegetation and constructed and open areas in the BGUs. We monitored the bat assemblage monthly over an annual cycle using autonomous ultrasonic recorders to assess the abundance of each species at each suitability level. Across all three suitability levels throughout the year, we found a high dominance of three insectivorous bat species that forage in open areas. Diversity measures supported our predictions, agreeing with bat ight and foraging functional traits. These results highlight the importance of city features in driving bat diversity and how urban planning that fails to consider key city features, such as arboreal areas in the city, may reduce suitability for bats, leading changes in diversity.


Introduction
The transformation of natural habitats into cities constitutes one of the most radical changes to nature. These changes are so dramatic that the usual result is the generation of conditions with negative effects on biodiversity (Turner et al. 1994 These modi cations can dictate the availability of food resources for many animal species, as well as their activity patterns (Gonsalves and Law 2018;Salvarina et al. 2018). Hence, the presence and abundance of a given set of bat species in a particular site can be in uenced by seasonal changes over the course of the year. One example of this is the presence of the nectarivorous bat species Leptonycteris yerbabuenae, which can be in uenced by the seasonal occurrence of ower resources as it migrates south in the tropical areas of Michoacán state, in Mexico (Racey and Entwistle 2000). A positive relationship has also been found between the seasonal abundance of some insectivorous bat species and that of insects in tropical and temperate forests (Montiel et al. 2006;Novaes et al. 2017; Salvarina et al. 2018). Although seasonal change in food resource availability is considered the main driver in uencing the abundance of plantivorous and animalivorous bat species in non-urban areas (Cisneros et al. 2015), to our knowledge, there have been few evaluations conducted in cities, and those that have assessed the in uence of seasonality in cities have been performed for only a limited number of species at a time, and have not considered the bat assemblages as a whole (Rollinson et al. 2013;Novaes et al. 2017; Gonsalves and Law 2018).
Variation in urban conditions thus creates areas with different features that, in turn, generate different suitability levels for bat presence, which does not necessarily agree with the approach of distance from the center of the city. In this study, we aimed to evaluate the effect on the bat assemblage of a suitability gradient interspersed within the urban matrix of Morelia city, Mexico. We also aimed to evaluate the impact of seasonality on this bat assemblage along a suitability gradient. We hypothesized that intrinsic bat characteristics, such as ight and foraging functional traits, would enable differential use of the urban habitat and thus foster the presence of certain species in some areas, while excluding other species from less suitable areas. We have four predictions in two sets that highlight the overall agreement of the suitability level with the expected dimensions of diversity driven by the species-speci c differential effects of human impact on the city suitability levels for bats. The rst set of predictions considers the negative impact of urbanization on some bat species as a result of a reduction in the degree of suitability for bats, in which: a) the bats will present reduced species richness and number of relative abundant species, and b) feeding guilds such as frugivorous bats will be either reduced in species number or lost. The second set of predictions invokes the positive effects of habitat modi cation on some species, in which: c) there will be an increase in the number of dominant bat species, d) the aerial insectivorous guild will be widely distributed in the city given the openness of the urban matrix due to the presence of constructed and open areas. Furthermore, given the in uence of food resources on changes in bat abundance (Gomes et al. 2013;Cisneros et al. 2015), we also expect that abundance values will increase during the rainy season, particularly for the insectivorous bat species.

Study site
The study was conducted in the city of Morelia, in Michoacán state, Mexico (19°52' to 19°26'N and 101°02' to 101°31'W). Morelia covers an area of 67.2 km 2 (IMPLAN 2018) and presents a minimum temperature of 10.4°C and maximum temperature of 27.5°C throughout the year, with an average annual rainfall of 796.5 mm (Servicio Meteorológico Nacional 2019). A total of 75% of the annual rainfall occurs from June to September, followed by a dry season that extends from October to May. Morelia is a medium-sized city with a population of 743,275 inhabitants in 2020 (INEGI 2021), and 400% growth towards the periphery recorded in the last three decades (López et al. 2001). Most of our understanding of the response of fauna in cities is derived from large urbanized areas such as Mexico City (Avila-Flores and Fenton 2005) and Sydney (Threlfall et al. 2011(Threlfall et al. , 2012a; however, the most common urban settings in which animals are found are medium-sized cities (

Selection of sampling units and estimation of suitability level gradient for bats in Morelia
For sampling, we used the units de ned by the Mexican Instituto Nacional de Estadística y Geografía (INEGI) that are used periodically for national human population censuses in Mexico (INEGI 1980(INEGI , 2010. These sampling units are known as Basic Geostatistical Units, or BGUs (also known as AGEBs, by their Spanish acronym) and are xed polygons, delimited by streets around a set of blocks, among which comparisons can then be made (INEGI 1980(INEGI , 2010. The use of these BGUs may enable us to delimit areas within the city that can be characterized by a set of parameters representing different levels of given urban characteristics. All BGUs in the city can then be associated along a gradient of suitability for bats, drawn from a set of parameters. Only urban BGUs within the limits of the city of Morelia were included in our sampling. To obtain a representation of the gradient of suitability for bats in Morelia city, we produced a land use cover map for the city. We used 2018 Landsat 8 satellite images to create this map. We classi ed the entire city into three cover types: a) areas with arboreal vegetation, b) open areas, presenting neither arboreal vegetation nor human constructions such as buildings and houses, and c) constructed areas that included houses and buildings. The map was obtained from training elds for all three cover types distributed throughout the city, using the Raster package (Hijmans 2019) for R 3.5.0 (R Core Team 2020). These training elds provided information at pixel level from satellite images. A supervised classi cation was subsequently performed using a classi cation and regression tree with the rpart package for R. The percentage of each cover type per BGU was obtained. The set of variables used to construct this map has been previously associated with the presence of bats (Li and Wilkins 2014).
The land use cover map was then used in a multicriteria analysis using ILWIS 3.3 to obtain the suitability gradient map for Morelia city. For this multicriteria analysis, we used the percentage of each cover type per BGU obtained from satellite images in raster format. All three variables were standardized to between 0-100, in which 100 represented the optimal suitability level for bats. Standardized variables were then weighted, based on expert opinion, in order to indicate how they could in uence the presence of bats as a whole, rather than on a species-by-species basis. Arboreal vegetation had the highest weighting value in the analysis (0.50). This analysis assigned BGUs in Morelia to one of three categories along a suitability gradient for bats. Low suitability level BGUs were represented to a greater extent to the north of the city, while the BGUs assigned to the medium suitability level were located mostly to the south of the city. High suitability level BGUs were located in the periphery and to the south of the city (Fig. 1). We selected four BGUs per suitability level, and acoustic monitoring was performed in each. Given that the selection of BGUs per suitability level was constrained to their location within the city (i.e. low suitability levels were mostly located to the north and medium and high suitability levels to the south of the city), as well as to the context in which they were located (i.e. we aimed to select BGUs that were surrounded by BGUs of a similar level category), their distribution per level was not homogenous across the city. The high suitability level category recorded the largest BGU average area (0.61 ± 0.22 km2, range = 0.38-0.81 km 2 ), compared to the low (0.45 ± 0.23, range = 0.17-0.68 km 2 ) or medium (0.54 ± 0.18; 0.42-0.81 km 2 ) suitability levels. We also obtained the percentage of land use cover type per suitability level. The BGUs that were classi ed as low suitability presented, on average, more open (average = 56%, range 48% -64%) and constructed (average = 42%, range 32% -52%) areas, and a lower percentage of arboreal vegetation (average = 2%, range 0% -5%), while those of high suitability recorded the highest percentage of arboreal vegetation (average = 52%, range 42% -63%) and a lower percentage of constructed areas (average = 15%, range 10% -19%) ( Table 1).

Experimental design for acoustic monitoring of bats
We monitored the bat assemblage at each of the four selected BGUs per suitability level in the city through passive acoustic monitoring, performed from November 2018 to October 2019 for four days each month, two days before and two days after the new moon in order to avoid the incidence of lunar phobia. Three BGUs were sampled simultaneously each day, one BGU per suitability level each night per month. Inclusion of each BGU in a given sampling day was random each month. Sampling sessions were conducted from 30 minutes before sunset and lasted for up to eight hours (Hourigan et al. 2006). No sampling was conducted on days with rain and wind. Sampling was performed using Song Meter SM4BAT FS autonomous ultrasound detectors in conjunction with an omnidirectional SMM-U1 ultrasonic microphone (Wildlife Acoustics 2016). The sites where these bat detectors were placed included the rooftops of houses and buildings, surrounded by more houses, which are the most common type of construction in the city, as well as buildings that are characteristic of the city. These sites were selected based on logistic accessibility, while also ensuring that the equipment would not be stolen. Microphones were placed in a site on average 4.99 ± 14.4 (range 4.95-5.45 meters, n = 12) m above ground level at a 45° angle facing open areas with no obstacles. Monitoring was performed with the ultrasonic recorders con gured such that an acoustic event was recorded when all of the following conditions were met: the sound reached at least 12 dB and had a duration of at least 1.5 ms. In addition, a bat pass event was recorded in independent les during surveys when the time elapsed between separate acoustic events was at least 3 s.
The microphone was set at 12 dB gain with no high-pass lter. Recordings were obtained in wav format at a sampling rate of 500 kHz.

Data sets for analyses
We used data from two sources: 1) a reference dataset with data of acoustic features of bat calls identi ed at species-level from a previous study (Zamora-Gutierrez et al. 2016) for training a Random Forest model, and 2) a eld dataset with acoustic parameters obtained from recordings from a monitoring survey of bats.

Data for training Random Forest model
We obtained a data set to train a Random Forest model. Random Forest (RF) is a tree-based modeling procedure for classi cation and regression problems to overcome the disadvantages of using a single Classi cation and Regression Tree (CART), including improvement of stability and reduction of over tting (Breiman 2001). RF produces a set of several independent trees using a random subset of sample data points along with independent variables at each split, and utilizes ensemble techniques to combine the output of the trees in order to improve predictive performance (Breiman 2001). RF has been previously used to successfully identify bat species using acoustic parameters from ultrasound bat recordings  Table 1). We then obtained a subset of records from the reference dataset (Zamora-Gutierrez et al. 2016) with the acoustic parameters for the list of bat species potentially present in the study area. The reference dataset included 27 acoustic parameters for search phase recordings/calls for 59 bat species occurring throughout Mexico, and is freely available from Dryad (https://datadryad.org/stash/dataset/doi:10.5061/dryad.760r8). However, only species with potential distribution in our study site were used.

Acoustic data from monitoring of bats for species identity assignation
We obtained a number of acoustic parameters from eld recordings to be used in the assignation of species identity. We automatically obtained seven acoustic parameters in recordings obtained from monitoring using the software Kaleidoscope Pro 3: i) constant frequency (Fc), ii) maximum frequency (Fmax), iii) minimum frequency (Fmin), iv) bandwith (Bndwdth), v) knee frequency (Fk), vi) duration (Dur), vii) initial slope (S1). Evaluation of the acoustic information provided by these parameters (Wildlife Acoustics 2016) indicated that they were equivalent to those used by the reference dataset. Additional acoustic variables obtained from recordings using Kaleidoscope included: average characteristic slope (Sc); average mean frequency (Fmean); average time between calls (TBC); average time to the knee (Tk); average time to the characteristic. However, these were not used in the species identity assignation since it was not possible to assign an equivalent to those used in the reference dataset. The variable of recording quality (Qual) was also obtained, but was used only during the depuration process and was excluded from the other analyses. Our database obtained from eld surveys was subject to a process of data depuration as described below.

Random Forest model training and species identity assignation
We used the Random Forest (RF) analysis in the caret package (Classi cation and Regression Training) (Kuhn 2008) for R 3.5.0 (R Core Team 2020) to assign species identity to each bat pass event recorded in the bat acoustic monitoring. This assignation of species identity comprised two steps. In the rst step, we conducted a RF in which we rst used the subset of acoustic parameters from the reference dataset for previously identi ed bat species that potentially occur in the study area. We then split this subset into a training sub-dataset, which contained 80% records from the subset of the reference dataset, and a validation sub-dataset with the 20% remaining records. We trained the model using the training subdataset. The trained RF model was constructed with 1,000 trees. For model training, we used model parameters similar to those of the reference source. Three predictor variables were randomly sampled at each stage. The resampling method used was 5-fold cross-validation, using the ve groups previously de ned by the original authors from the reference dataset. We used the automatic grid search, with 1974 as the seed set. The variable to predict with the RF was bat species identity, using the 7 acoustic parameters as predictors. Finally, we validated predictions with the validation sub-dataset from the reference dataset. The predictive performance of the model for species identi cation was evaluated using accuracy. The resulting predictive performance of the trained model for species identi cation had an overall average accuracy of 65% in terms of correct classi cation. The worst performance in the training dataset occurred with Molossus sinaloae and the Corynorhinus phonotype (Supplementary material  Table 3). In the second step, we predicted species identity for each recorded bat pass, from eld data using the previously trained RF model (Supplementary material, Fig. 2).

Process of data depuration
We performed a process of depuration on the les where a species was assigned in order to obtain recordings that were more suitable for analysis. We initially considered 38,625 bat pass events recorded in independent les during surveys. Prior to performing the analyses for species identity assignation, we rst excluded 16,240 les with less than or equal to three echolocation pulses from those bat passes recorded during surveys (Hourigan et al. 2010). We then used the variable Qual, obtained from Kaleidoscope (see above), as a proxy for the quality of recordings. We excluded 5,394 additional recordings for which the quality values were in the rst quartile, in order to ensure we were using the highest quality records in subsequent analyses. As a proxy for precision of species identi cation, we then obtained accuracy values for each species using the package RandomForest. Only species with accuracy values ≥ 0.4 were included in the analyses. Accuracy values ranged from 0-1.0, and 6,065 records with accuracy value < 0.4 were therefore assigned to an unknown species category and excluded from further analysis. We then obtained the relative abundance of each bat species by measuring its occurrence at one-minute intervals, rather than the number of times it was recorded during this one-minute period. The total number of one-minute intervals at which the species occurred for the complete sampling night per BGU was used as a proxy for the relative abundance of each species (Miller 2001). This procedure has been used to obtain abundance estimations in other animal groups, such as birds (Celis-Murillo et al. 2009, 2014), and has been correlated with visual estimations of abundance in bats (Gehrt and Chelsvig 2003). This procedure enabled us to not only reduce bias due to incorrect species assignation, but also to reduce overestimations (Miller 2001;Hourigan et al. 2010; Trejo-Ortiz 2011). Thus, we obtained a dataset of 10,188 records of the occurrence of bat species during the one-minute intervals.
Once a species identity was assigned to each bat pass, the acoustic records of the species were reviewed manually 90 min after initiating each recording session and again at 90 min after the end of each recording session. This was conducted in order to exclude urban sounds or those of non-bat species that are common in the urban soundscape during the sampling schedule. We then excluded 2,436 records from further analysis because these were not recordings of bat species. In addition, notable records, such as those for species with few records or belonging to the same genus, were also revised manually to ensure that the species were correctly assigned. These species included Eptesicus furinalis Artibeus-Sturnira phonotype, Molossus molossus, M. rufus and M. sinaloae. There are a number of species that have presented low performance in previous classi cation works (Pio et al. 2010; Zamora-Gutierrez et al. 2016), and some of these species identi ed in the species assignation procedure were therefore assigned to phonotypes. Records assigned to different species from the genus Artibeus and Sturnira were grouped into an Artibeus-Sturnira phonotype. Records assigned to Corynorhinus mexicanus and C. townsendii were indicated as a Corynorhinus phonotype. All records assigned to the different Myotis species potentially occurring in the study area were grouped into a Myotis phonotype. Finally, we excluded ve species that contained < 10 records from further analyses (Kloepper et al. 2016), although these species are considered in the taxonomic list. These excluded species comprised 18 records, which corresponded to Lasiurus blossevillii, Leptonycteris yerbabuenae, Molossus sinaloae and to the Corynorhinus and Myotis phonotypes. This process of quality veri cation yielded a nal database that comprised 7,726 records.

Data analysis
Data analysis was performed per BGU at two levels: 1) the complete annual sampling period, and 2) the dry and rainy seasons separately. The dry season comprised the period October to May, while the rainy season comprised June to September. We created rank relative abundance plots per suitability level. Estimations of diversity per BGU were obtained by calculating the rst three Hill numbers (Hill 1973 To evaluate the similarity in bat species composition among BGUs with different levels of suitability, we performed non-metric multidimensional scaling (NMDS). We used a Bray-Curtis similarity matrix, considering bat species richness and relative abundance. We used the stress value as a measure of the goodness of t in the data ordination, with a range of 0 to 1. Stress values closer to 0 indicate better ordination relative to the original data (Oksane et al. 2019). The NMDS graph represents the degrees of similarity among points in terms of bat richness and relative abundance. The points representing BGUs located further apart in the NMDS ordination space present greater differences than those of BGUs in greater proximity to each other. We plotted the rst two NMDS axes, and overlaid tted contours of the percentage of arboreal cover, constructed cover and open area cover per BGU onto the ordination space.
Fitted contours were obtained using the ordisurf function of the vegan package for R, and the values plotted in the ordination space indicated the degree of association between these BGUs and the variable.
Finally, to evaluate whether suitability level and season could help to explain variance in the diversity components of q0, q1, and q2, and relative abundance per bat species, we constructed generalized linear mixed models (GLMMs) independently. In all models, suitability (with three levels) and season (with two levels) were considered as xed effects. Since we were interested in the in uence of coarse-grain city features such as areas of suitability for bats, we concentrated our analysis on this level. The BGU and suitability level were considered as random factors in each model to account for the lack of independence. This enables us to take into account the systematic non-explained variation (Faraway 2016), due to the use of several sampling units per suitability level. For the q0 model, we used a negative binomial error distribution with the lme4 package (Bates et al. 2015). Both q1 and q2 models were performed using lognormal error distributions with the lme4 package. For the relative abundance model, we used a zeroin ated with negative binomial error distribution using the glmmTMB package (Brook et al. 2017). We used a signi cance level of p < 0.05, and all analyses were performed in R.

Bat richness and abundance in Morelia
We sampled for 48 nights in 12 BGUs in Morelia, for a total of 7,726 recording minutes. We recorded 14 species, belonging to 15 genera, ve families, and three phonotypes. We also recorded three feeding guilds, considering the classi cations proposed by Kalko et al. (1996): 1) aerial insectivorous bats; 2) gleaning insectivorous bats; and 3) frugivorous bats (Supplementary material, Table 2). The bat assemblage in Morelia was highly represented by insectivorous bat species, which accounted for almost 98% of the records.

Bat diversity evaluation per suitability level
Three aerial insectivorous bat species were found to be the most abundant at all three suitability levels: Promops centralis was the most abundant (low suitability: 71% of all records per suitability level, medium suitability: 71% of all records per suitability level, and high suitability: 61% of all records per suitability level), followed by Eumops underwoodi (low suitability: 21% of all records per suitability level, medium suitability: 13% of all records per suitability level, and high suitability: 19% of all records per suitability level) and nally Molossus rufus (low suitability: 4% of all records per suitability level, medium suitability: 8% of all records per suitability level, and high suitability: 6% of all records per suitability level). The bat community per feeding guild and suitability level recorded speci c changes. The BGUs of low suitability level were characterized by presenting exclusively insectivorous bat species, while those of both medium and high suitability level presented frugivorous bat species, albeit in low numbers. Insectivorous gleaning species showed a positive increase in relative abundance with suitability level. The BGUs of high suitability level also presented an increase in the relative abundance of frugivorous bat species (Fig. 2a). Both speci c richness (q0) and the number of species with a similar abundance pattern (q1) increased concomitantly with increased suitability level (Fig. 3a, 3b). Average q0 and q1 values per BGU were lower at the low suitability level (average low: q0 = 7.25 species, range = 5-10; q1 = 2.22 species, range = 1.70-2.94) when compared to the high suitability level BGUs, which presented the highest average number of species (average high: q0 = 9.75 species, range 7-11; q1 = 3.30 species, range = 2.49-4.57, Fig. 3a and   3b). However, q2 showed an inverse pattern, with a reduction in the number of dominant species with increased suitability level for bats. The largest average q2 value per BGU was recorded at the low suitability level (q2 = 2.86 species, range = 1.91-4.84, Fig. 3c), compared to the high (q2 = 1.83 species, range = 1.45-2.38) suitability level (Fig. 3c).

Bat diversity evaluation per season
Overall, all three suitability levels, in both the dry and rainy seasons, were dominated by the three aerial insectivorous bats (P. centralis, E. underwoodi and M. rufus). However, the relative abundance of P. centralis was higher in the dry season than in the rainy season. The bat assemblage showed seasonal changes in guild level per suitability level. The low suitability level presented only insectivorous bat species in both seasons, while a higher number of gleaning bat species was recorded during the dry season, compared to the rainy season (Fig. 2b). Although both the medium and high suitability levels presented frugivorous bat species during both seasons, the high suitability level BGUs presented a larger relative abundance of frugivorous bat species during the rainy season (Fig. 2c, 2d). Moreover, a set of insectivorous bats, such as Mormoops megalophylla and M. molossus, presented increased relative abundance during the rainy season in the high suitability level BGUs (Fig. 2d). The average species richness (q0) and the average number of dominant species (q2) were both lower in the rainy season compared to the dry season (Fig. 3d, 3f). In contrast, the number of species with similar abundance pattern values (q1) showed an inverse pattern, with higher average values presented in the rainy season (Fig. 3e).
Results from the GLMMs indicated that the suitability level and seasonality explained the variation in our data, but in different ways. This analysis indicated that the low suitability levels and the high suitability level in the rainy season signi cantly explained variation for q0, q1, and q2 as a measure of bat diversity in Morelia for the annual sampling period. These results also indicated that variance in the diversity components of q0, q1, and q2 could be explained by the dry season. The largest effect was caused by the dry season for speci c richness (q0). The largest effect for the number of species with a similar abundance pattern (q1) was due to a high suitability level in the rainy season, and both dry season and low suitability level equally have the lowest effect. The dry season had the greatest effect in q2, while relative abundance was explained positively solely by the high suitability level in the rainy season (Table 2). Marginal R 2 was low in all models.
Similarity in bat biodiversity across the city Results from the non-metric multidimensional scaling indicate that species composition and relative abundance varied among suitability levels. For the annual evaluation, the BGUs were segregated in the ordination space according to suitability level. Medium and high suitability level BGUs were more similar in composition and relative abundance in comparison to those from the low suitability level (Stress = 0.093, Fig. 4a, 4b, and 4c). This pattern, in which high and medium suitability level BGUs overlapped and segregated from the low suitability level, was repeated in the rainy season (Stress = 0.092, Supplementary material Fig. 1a, 1c, and 1e). However, during the dry season, the NMDS indicated that the BGUs per suitability level were more similar in species composition and relative abundance among each other, since all three suitability levels overlapped (Stress = 0.128, Supplementary material Fig. 1b, 1d, and 1f). The BGUs distributed in the ordination space according to suitability level also agreed with the overlaid tted contours of the percentage of open area and constructed area cover types per BGU onto the ordination space for the complete sampling period (open areas: Fig. 4a; constructed areas: Fig. 4b), and for the rainy species that increased using the autonomous recorders. Although we consider that the use of this sampling method reduces underestimation of bat diversity in species richness and abundance, the use of autonomous recorders is not exempt from bias. We acknowledge that this sampling approach could have led to an underestimation of species richness. For instance, given that areas with closed arboreal vegetation are preferentially avoided when these bat recorders have been used (Wildlife Acoustics 2016), any species that do use these areas would not be recorded as occurring in the city. In addition, bat species from the Phyllostomidae family can produce low intensity calls (Kalko 2002), which could complicate recording these individuals with bat detectors, and we acknowledge that the abundance estimations for members of this family in our study may have been underestimated as a result. This could potentially lead to sampling bias in relative abundances and species richness due to systematic error across the BGUs, and the numbers of records for this family would not necessarily indicate whether these species are absent or in low numbers in our sampled BGUs. Nevertheless, these recorders do present lower bias than traditional sampling methods such as mist nets. In addition, security considerations prevented us from sampling in arboreal areas in the city. Morelia recorded the dominance of three bat species across all three suitability levels. Such a high dominance of few species is a pattern that is consistently found in modi ed habitats and has been observed in cities not only in bats but also other animal groups, such as birds (Ortega-Álvarez and MacGregor-Fors 2009; Dixon 2012; Moretto and Francis 2017). This is more evident at the low suitability level, since this presented the highest number of dominant species when compared to the other suitability levels. The increase in the number of dominant species is normally the result of changes in the habitat, creating more favorable local conditions for these bat species and then fostering their increase (Threlfall et al. 2012b). Low suitability BGUs, characterized by presenting the largest percentage of their area with no vegetation, could therefore favor the foraging strategy of the aerial insectivorous species that were recorded in these BGUs. Aerial insectivorous bat species foraging in open areas could then more e ciently exploit locally modi ed conditions. One example of this is the dominance of the aerial insectivorous species P. centralis, E. underwoodi and M. rufus. These three bat species belong to the Molossidae family, which is considered less susceptible to more open urban environments since most of the species in this family present higher body mass as well as a larger wing load and wing aspect (Jung and Kalko 2011).
These features permit them to perform sustained ights at greater heights and higher speeds, enabling them to move over larger distances in one night (Norberg and Rayner 1987 Threlfall 2016, 2018). This guild could therefore be the most resilient to changes in conditions of suitability for bats in the city and be more widely distributed throughout the city.
Conversely, we found an increase in the richness and abundance of frugivorous bat species with an increased suitability level. This suggests that the low percentage of arboreal vegetation favoring aerial insectivorous bat species could also act as a limiting factor for other bat species in the city, such as the frugivorous bat species. Frugivorous bat species have been captured in Morelia city using mist nets (Martínez-Mijares 2016). Although there is no previous sampling in the high suitability level BGUs that could con rm the prior presence of frugivorous bats in this suitability level, sampling in these previous studies using mist nets did show frugivorous bats occurring in the city. Considering the suitability level that we assigned to each BGU in the city, sampling sites for these studies where frugivorous bats were recorded could be assigned to medium suitability level BGUs, given their location. One of the main drivers increasing the abundance of frugivorous bat species in urban parks is the amount and quality of fruits they use, as occurs with A. jamaicensis (Jara-Servín et al. 2017). The increase in richness and abundances of frugivorous and gleaning bat species in areas of medium and high suitability levels highlights the importance of the presence of areas of arboreal vegetation to the interior of cities. Of particular importance in fostering the presence of frugivorous bat species is the medium suitability level, despite its similarity to the low suitability level in terms of proportion of vegetation, given the large proportion of this suitability level BGUs in Morelia city.
During the rainy season, we recorded an increase in the number of individuals of insectivorous species such as P. hesperus, M. molossus, E. brasiliensis and M. megalophylla. The most common explanation for these seasonal changes is increased food availability due to the rainfall during this season, particularly in terms of the availability of prey for the insectivorous bat species. There are no studies in the city of Morelia evaluating seasonal uctuations in the availability of ying insects on which bats could forage; however, the abundance of ying insects in areas near to Morelia characterized by oak forest, indicate that they increase in numbers during the months of highest rainfall (Jurado-Vargas 1990; Báez-SantaCruz 2011). An alternative explanation for the increase in abundance during the rainy season could be related to the synchronization of the timing of the independence of pups from their mothers with that of the highest availability of food resources during this season. This synchrony has been recorded in both Pipistrellus nathusii and Nictalus noctule, where an increase in abundance co-occurs with pup independence not only in forested areas in the Americas, but also in farmlands in Europe (Racey and Entwistle

City management implications
The assignment of BGUs to a suitability level for bats showed geographic strati cation according to suitability level. The BGUs to the north of the city were mostly of low suitability for bats, whereas those to the south of the city were mostly medium with a few of high suitability level. This is re ected in the distribution of urban parks in the city, since Morelia has nine urban parks, of sizes ranging from 10 to 118 ha, with an important presence of arboreal vegetation. Seven of these nine parks are in the southern part of the city, whereas only two are in the northern part. This indicates not only the lower presence of vegetation in the north of the city, but also that the north hosts fewer urban parks, where there is a dominance of low suitability level BGUs, and suggests that land use cover could be a key factor driving the suitability of areas in the urban matrix for bats and creating a north-south strati cation of areas within the city.
The most marked difference between the low/medium and the high suitability levels is the percentage of arboreal vegetation areas in each BGU. Low and medium suitability areas were more similar in terms of the percentage of open and constructed areas, and thus over 90% of the city presents these features. This suggests that the sampled low and medium suitability level areas should be more similar in terms of both richness and abundance, and it could be expected that these segregate from the high suitability areas in the ordination space. Despite this, it was the areas of medium and high suitability levels that were more similar in diversity and were segregated from the low suitability areas, as indicated by NMDS plots. This could suggest that areas of medium and high suitability levels could be functionally similar as a result of the high proportion of medium and high suitability level areas close one to each other. Alternatively, the suitability level may in uence not only species richness and abundance but also species composition, since there are species, such as the frugivores, that are restricted to medium and high suitability levels.
Thus, while the medium and high suitability levels located in the south of the city could be similar in terms of diversity dimensions, they would not be functionally interchangeable due to differences in species composition.
The low suitability level represents one extreme of the suitability gradient, with gleaning and aerial insectivorous bat species occurring at this level of suitability. On the other side of the spectrum is the high suitability level, where we recorded an increase in relative abundance of frugivorous and gleaning insectivorous bat species. This suggests that some bat species will apparently be favored by the conditions for foraging provided by the presence of forested areas. Some other species, such as frugivorous bat species in areas of high suitability level, would require particular conditions to occur in this highly modi ed area. This acts to create a gradient in patterns of richness and abundance, but more markedly generates a potential spatial segregation to the interior of the city in bat species composition as the result of their different requirements and variability in the level of susceptibility to habitat modi cation based on the intrinsic characteristics of bat species.
Our study therefore indicates that suitability levels consistently in uence the structuring of the bat assemblage in the city of Morelia. We found that species richness and feeding guilds, as well as the number of abundant species, increased in line with suitability, while the number of dominant species showed an inverse pattern. Morelia is considered a medium-sized city that presents accelerated growth (López et al. 2001). If this trend continues, the negative impact of this growth could be more marked in the long run. This is particularly true for those species using high suitability areas where there is a higher occurrence of individuals from the frugivorous feeding guild. This would mean that Morelia could keep losing feeding guilds not only from the high suitability areas, but also from the medium suitability areas to the south. This would also mean an increase in the abundance of some bat species that have been favored by the current modi cation status, as is the case of aerial insectivorous bat species. Several bat species require the presence of both open areas and those with arboreal vegetation, so it is important to not only maintain the existing arboreal vegetation areas in urban paci cation plans for the city, but also to promote the creation of more urban parks and arboreal forested areas, particularly in the northern part of the city.  Table 2. Results from the generalized linear mixed models (GLMMs) for q0, q1, q2, and relative bat abundance for suitability level and the dry and rainy seasons. High suitability level and rainy season are included in the intercept. 95% con dence intervals are shown in parenthesis.   Rank abundance plots for bat species: a) Complete sampling period by suitability level: b) low bat suitability level per season, c) medium bat suitability level per season, d) high bat suitability level per season. Name of bat species according to Supplementary Material, Table 2.