Spatial and environmental variation in phyllostomid bat (Chiroptera, Phyllostomidae) distribution in Mexico

Spatial and environmental variation in phyllostomid bat (Chiroptera, Phyllostomidae) distribution in Mexico. Species’ spatial distribution patterns allow us to understand the establishment of different biotic components in different environmental conditions. This study analyzes the spatial distribution of the Phyllostomidae family in Mexico to identify groups of species that occur in similar sites, the environmental conditions associated with species distribution, and the percent of overlap with human–modified areas. The results suggest six groups of sites with particular species composition. The spatial variation in richness pattern was associated with species tolerance to environmental conditions, such as minimum temperature and tree cover. The convergence between species distribution and modified areas varied per species feeding guild. Insectivorous and nectarivorous bats were sensitive species because they occurred in narrow environmental conditions and their distributions overlapped with areas modified by human activities. The approach implemented here analyzes regional species distributions and estimates their environmental requirements, contributing to the development of optimal conservation strategies for susceptible bat species.


Introduction
Species diversity distribution is heterogeneous in space, so understanding and predicting this variation in different environments and taxa is a fundamental goal in ecology and biogeography (Lomolino et al., 2005). Species richness and community assembly variation may occur either gradually or abruptly, depending on the spatial scale analyzed and associated factors (e.g. availability of energy and water; Hawkins et al., 2003;Field et al., 2009). Species' spatial distribution patterns are useful to understand the establishment of different biotic components in different environmental conditions (Morrone, 2009). However, the persistence of biological diversity and the functionality of ecosystems are at risk due to the increase in anthropogenic activities such as agriculture, livestock, and urban activities (Klein-Goldewijk and Ramankutty, 2004). Changes in native vegetation, for example, have modified 12 % of Earth's land surface, affecting wildlife habitats and the diversity of species distribution (Ries et al., 2004). Therefore, the patterns of species composition vary depending on the environment, and high biodiversity areas that are threatened by human land cover changes must be evaluated to optimize conservation strategies, a key topic in conservation biology (Margules and Sarkar, 2009).
Species potential distribution models are a useful tool to identify environmental conditions related to species presence, high biodiversity areas and zones with similar species composition (Mateo et al., 2013). Despite being locally biased by the number of collections records, the type of sampling, the choice of predictors or the algorithm used (Peterson et al., 2011), the fact that species distribution models provide reliable results at a regional scale (Raxwhorty et al., 2007;Lee et al., 2012) makes this technique appropriate to study bat communities at larger scales. The superposition of information from a distribution model and land-use cartography provides an estimate of the degree of overlap between human-modified areas and can be a useful strategy to identify areas at risk from human activities (Wu et al., 2014). An integrative approach to distribution models allows helps to evaluate ecological suitability in biodiversity areas to identify species that are susceptible to environmental changes and important as conservation targets (Peterson et al., 2015).
The Phyllostomidae family is a Neotropical taxa that can be used to identify high biodiversity areas and zones at risk of human activities because it occurs in several different environments (Stevens, 2006) and has different specific feeding preferences (Giannini and Kalko, 2004). The regional distribution of phyllostomids has been associated with oscillations in temperature and humidity (McCain, 2007) and is locally determined by sensitivity to disturbance (Wordley et al., 2015). However, little is known about areas with a similar species composition, species environmental tolerances, and the degree of overlap between the distribution of bats and modified areas (López-González et al., 2011;Razgour et al., 2016). Phyllostomid bats provide ecosystem services as agents of seed dispersal, pollination, and pest regulation (Kunz et al., 2011).High diversity areas that may include species susceptible to constant human degradation must therefore be identified.
Mexico has a wide range of environmental conditions, produced by altitudinal variations, the influence of two oceans, and the convergence of two biogeographical regions (i.e., the Neotropical and the Nearctic; Morrone, 2005) that offer multiple environments for 60 phyllostomid bat species (Ramirez-Pulido et al., 2014). However, these bat species are at risk from human activities, mainly the conversion of tropical forests (Challenger and Dirzo, 2009). This study analyzes the distribution of Phyllostomidae family species in Mexico, as well their different trophic guilds, since this information is useful to understand the response of the community structure to human disturbance (Klingbeil and Willig, 2009). The objectives of this study were to identify the spatial richness patterns, define groups of sites with a similar species composition, estimate the environmental tolerance of each species, and calculate the convergence between species distribution and human-modified areas.

Species potential geographic distribution
Potential geographic distribution models for phyllostomid bat species in Mexico were developed using the maximum entropy algorithm (MaxEnt ver 3.3; Phillips et al., 2006Phillips et al., , 2016, which has proven to be robust (Elith et al., 2006). The algorithm searches a combination of variables with maximum entropy and estimates the importance of each in species distribution with respect to sites where the species was recorded (Elith et al., 2011). Presence records were obtained from the database provided by the Instituto de Biología of the Universidad Nacional Autónoma de México. The database was complemented with records from the Global Biodiversity Information Facility (www.gbif.org). A total of 64,773 records were compiled, from which doubtful records or those with spatial redundancy were removed, resulting in 11,701 records.
The variables for predicting the species distribution included those that limit their presence at regional scales, such as temperature, precipitation, and elevation, and those variables related to niche requirements such as type of vegetation and a regionalization variable, at a spatial resolution of 0.0083º (~0.85 km 2 ; see table 1s in supplementary material). Climatic factors (such as temperature and precipitation) were used as direct variables that have physiological importance for bat species, but are not consumed; vegetation cover variables (such as trees or herbaceous plants) were related to resources used directly (López-González et al., 2011). Elevation was an indirect ecological variable that has no direct physiological relevance for species' persistence in response to future changes in environmental conditions (Guisan and Zimmermann, 2000), while the categorical variable of regionalization represents unique fauna and flora of the world's continents (Olson et al., 2001). The suite of predicted variables allows us to evaluate bat species' responses to the environmental gradient present in Mexico. The layers of climatic variables were modified to generate new layers (see table 2s in supplementary material) using the Image Calculator in IDRISI (edition Selva; Eastman, 2012). However, using a set of variables is biased to autocorrelation (Peterson et al., 2011) and so to reduce information redundancy, we deleted continuous variables with a high correlation (r 2 > 0.8) and performed a cluster analysis (1-Pearson r distance and Ward algorithm). The cluster analysis was performed with a sample of 500 random points, generated in ArcGIS (ver. 10.1; ESRI, 2010) using the Random Points Extension. The number of groups of correlated variables was defined from the amalgamation graph (threshold value = 1.1). From each group, we chose one variable and excluded the rest (see fig. 1s in supplementary material). Potential distribution models were performed with one categorical and 10 continuous low correlated variables to represent the environmental heterogeneity in the country (table 1).
Models were generated only for species with a minimum of 10 records (53 species; table 2), since MaxEnt is stable with this number of records (Wisz et al., 2008). Default MaxEnt settings were used, with 75 % of the records used to create the model and 25 % to test it. Due to this random component, 10 single models were generated for each species to compile a consensus map via the weighted average method (Marmion et al., 2009). This method provides a continuous map and to obtain the species presence/absence, we calculated the average occurrence threshold from the single models. The threshold was the fixed cumulative value for 10 % probability of occurrence. Species were considered present in values above the threshold and absent in values below the threshold. Thus, a binary consensus map was obtained, increasing accuracy and decreasing the uncertainty of single models.

Species richness variation
The spatial pattern of richness for phyllostomids in Mexico was obtained from the sum of the consensus presence/ absence (binary) models for each species in IDRISI (Eastman, 2012). Additionally, we generated a richness map for trophic guild classification (Giannini and Kalko, 2004;table 2). Richness maps show the variation in species richness and spatially to identify highly specific richness areas.

Sites with similar species composition
Binary maps were used to define sites with a similar species composition. We used 500 random points (sites) to extract the presence/absence from each model and to compile a binary matrix that contained all the possible species combinations. We defined groups of sites with similar species composition using a generalized k-means analysis, which finds the optimum 'partition' for dividing a number of objects into k clusters in function of categorical variables (presence-absence) in STATISTICA (StatSoft, 1984(StatSoft, -2013. The used k groups were defined from the amalgamation graph of a cluster analysis (Sorensen-Dice distance and Ward algorithm). The analysis searches for a combination of sites that maximize significant differences in local species composition between groups based on individual x 2 -tests for each species. The null hypothesis was that the frequency of sites per group was similar at a probability of 0.05. The result was the assignment of each site in a group as a function of local species composition. We present a map that shows the spatial tendency in the geographical distribution of groups of sites that differ in bat species composition.
In addition, we identified species that distinguish each group based on their proportional distribution affinity, obtained as the percentage of sites by group where each species occurs. The affinity values (from 0 to 100 %) indicate the species geographical distribution along the different groups of sites.

Species environmental segregation
Bat community-level response to environmental conditions was identified via the outlying mean index (OMI) to estimate the preference (marginality) and the amplitude (tolerance) of conditions used by each species (Dolédec et al., 2000). The analysis estimates the distance between the average conditions (centroid) in which each species was present with respect to the average conditions of the study area (center of gravity). The OMI generates canonical variables that account for the highest variability of sites, which were associated with sites where species were recorded by a canonical correspondence. As a result, the variability (inertia) associated with the species distribution is decomposed in three components (Dolédec et al., 2000): (1) marginality or OMI, which is the deviation of average conditions used for the species in the study area. High values for species were found in different conditions of the average evaluated and low values for species were found in the average; (2) tolerance, which is the number of sites with which species are associated and their location in an environmental gradient (i.e., niche breadth). Low values imply that a species occurs in a narrow range of conditions (i.e., specialist) and high values imply that a species occurs in a wide range (i.e., generalist); and (3) residual tolerance, which is the variation in species occurrence not explained by the variables.
The OMI analysis was performed in ADE-4 software (Thioulouse et al., 1997) using a faunal matrix with species presence/ absence and an environmental matrix with environmental characteristics; both were extracted from the 500 sites. The faunal matrix was used in the k-means analysis, while the environmental matrix included values of the 10 continuous variables used to generate the models (table 1). Additionally, we included a new column in the faunal matrix containing points where species absence was predicted. The column identifies environmental conditions that limit the species distribution. Statistical significance was estimated with a Monte Carlo permutation test (1,000 permutations; Metropolis and Ulam, 1949). We report eigenvalues, factor loadings, and a graph of species centroids.
To identify bat species that tend to co-occur in similar sites, we grouped the species according to the environmental conditions where their presence was predicted. The grouping was conducted with a Q agglomerative analysis for the faunal matrix via cluster analysis (Sorensen-Dice distance and Ward algorithm). The groups found are highlighted on the species marginality graph to determine whether a tendency exists for species to be homogeneously dispersed within the available environmental gradient. All cluster analyses were performed in STATISTICA (ver. 12;StatSoft, 1984StatSoft, -2013.

Species incidence with anthropic changes
The degree to which species occurrence overlapped with human activities was analyzed by estimating the percentage of area in which a species distribution converges with modified areas. Modified areas were obtained using a layer of land use/vegetation in Mexico (INEGI serie V, www.inegi.org.mx). The layer was reclassified into cropland, livestock, and urban zones. We obtained the percent of overlap between family and trophic guild richness using the classified land-use layer. A convergence percentage was calculated from the ratio of overlain pixels to the total pixels in the study area, multiplied by 100 (Venegas- Barrera and Manjarrez, 2011). Additionally, we present the convergence percentage by richness proportions for the family and guilds in groups of sites where species diversity was highest.

Spatial patterns in phyllostomid bat richness
The Phyllostomidae family is potentially distributed in 79.6 % of Mexico's land surface, and we found that the highest richness areas were in the southern tropical environments (> 40 spp., fig. 1A). Nectarivorous bats were the guild with the widest geographic distribution (73.9 %), whereas hematophagous species were found in 46.2 % of Mexico and frugivorous and insectivorous bats occurred in 51.5 % and 59 % of Mexico, respectively. Nectarivorous bats showed a potential highest richness on the Pacific coast, while the other guilds had a higher richness on the Atlantic coast ( fig. 1B-1E). Richness for all groups was lowest in the northern arid environments.

Groups of sites with similar phyllostomid bat composition
We found six groups of sites that contained a distinctive species (threshold value of 1.5 unities; fig. 2): (1) The Atlantic group, on the coast of the Gulf of Mexico and the Yucatan Peninsula, had the highest potential richness areas (10 to 49 species) with Mimon crenulatum, Lophostoma brasiliense, and Vampyrum spectrum as higher proportionality affinity species ( The second axis explained 10.8 % of the variation due to variations in the precipitation percentage in July and the herbaceous cover percentage ( fig. 3A). The first OMI axis separates the environments where the family occurs from those where it is absent. Phyllostomid bats were present in environments that have a yearly minimum mean temperature above 13.5 ºC (from 12 to 20 ºC; fig. 3B) and a mean tree-cover percentage higher than 24 % (from 20 to 80 %; fig. 3C). The second axis revealed an environmental gradient associated with variations in precipitation in which the species are distributed. Most species (e.g., Mimon crenulatum, Phylloderma stenops) occurred in environments where six to 12 % of the annual rainfall occurs in July ( fig. 3D) and the herbaceous cover percentage is from 20 to 40 % ( fig. 3E). Conversely, several species (e.g., Choeronycteris mexicana, Macrotus californicus) occurred in environments where 12 to 18 % of the annual precipitation occurs in July and the herbaceous cover constitutes 40 to 60 %. We found three groups of bat species that occurred in similar sites (threshold value of 2.5 units, fig. 3F). The groups were differentiated according to species distribution range: (1) Widespread distribution (across 17 to 51 % of Mexico's surface) was found throughout the entire study area, occurring in groups in all sites. This species group was composed of 24 species with mainly frugivorous or nectarivorous bats (e.g.,   Sturnira lilium, Leptonycteris yerbabuenae) and one hematophagous species (i.e., Desmodus rotundus).
(2) Regional distribution (in eight to 21%) was associated with the Neotropical region, with a limited presence in mountain systems, mainly present in the Atlantic group. This species group was composed of 18 species, mainly insectivorous bats (e.g., Overlap between phyllostomid bat distribution and human modified areas Native/secondary vegetation covers 72.2 % of Mexico, and the remaining 27.8 % has been human-modified for crops (17 %), livestock (9.9 %) and urban zones (0.9 %). We found 30.9 % of the entire Phyllostomidae family spatial distribution overlapped with human-modified areas (table 3). Cropland was the most broadly distributed cover type converging with the family and the four feeding guilds. In the Atlantic group, the high richness areas of the Phyllostomidae family and the guilds of hematophagous, insectivorous and frugivorous bats overlapped with livestock ( fig. 4). In the Pacific group, the high richness areas for nectarivorous bats converged similarly with all types of human changes.

Discussion
In this study, we provide a perspective about phyllostomid bat distribution in Mexico. We identified six site groups with distinctive species assembly, three species groups that co-occurred in similar sites, environmental tolerances of each species, and the degree of overlap with human activities. The distribution of phyllostomid bats showed a spatial richness pattern that related negatively to latitude, decreasing the number of species from tropical environments in southern to arid environments in northern Mexico ( fig. 1A; Stevens, 2006). The highest species richness areas were in the Tehuantepec Isthmus region, except for nectarivorous bats (figs. 1B-1E). The latitudinal pattern has been reported for groups    (Escalante et al., 2007). The   The species occurrences with common origin and differential environmental preferences in the same area are related to their evolutionary histories leading to their diversification and current distributions (Dumont et al., 2012).  The groups of sites identified as the Neotropic region, the Atlantic and Pacific groups, present distinctive environments (i.e., tropical and subtropical) and are composed of ancestral bat species, such as Micronycteris or Diphylla genus. The Atlantic group containe mostly insectivorous species from the subfamily Phyllostominae (e.g., Mimon crenulatum and Lophostoma brasiliense), with a high marginality (mean OMI = 6; table 2) that occur in specific environmental conditions, such as perennial forests. Insectivorous bats are narrowly distributed ( fig. 3) and have been associated with large forest fragments, where a high availability of resources is found (e.g., food and roost; Medellín et al., 2000). On the other hand, the Pacific group is composed principally of frugivorous and nectarivorous species with a narrow or widespread distribution (e.g., Musonycteris harrisoni and Chiroderma salvini), in deciduous forests. However, the species has a higher tolerance (mean Tol ~1) than species in the Atlantic group, allowing their persistence when resources are scarce through the variation in seasons found in deciduous environments (Chávez and Ceballos, 2001).
The groups of sites identified in the limits of the Neotropic region or in the Nearctic region present distinctive environments (i.e., temperate, arid) and are composed of the phylogenetically derived bat species, such as Artibeus and Leptonycteris genus ( fig. 2). In the Mountain and Plateau groups, an elevational gradient derived from the mountains promotes variable climatic conditions (i.e., temperature and humidity) that results in different vegetation types (e.g., temperate forests and grasslands; Ramanoorthy et al., 1998). The groups are composed of frugivorous and nectarivorous species (e.g., Artibeus aztecus and Leptonycteris nivalis) that have low marginality (mean OMI < 1.3, table 2). However, the species are distributed widely ( fig. 3) due to adaptations such as dietary changes and migratory movements that allow them to colonize temperate environments (Fleming et al., 2009). In Northern Mexico, the species most commonly present in the Arid and Desert groups are Macrotus californicus and Leptonycteris nivalis, species that have a widespread distribution and that are tolerant to average conditions in Mexico (mean Tol = 0.5). The species that occur in semi-arid environments, such as Macrotus californicus, tolerate desert conditions through behavioral and physiological adaptations (Bell et al., 1986), while Leptonycteris nivalis uses migration and mutualism with scrubland xeric plants (Fleming et al., 2009). The persistence of these species in arid zones is high. These zones are the best conserved ecosystems in Mexico because the low availability of water limits the establishment and growth of human populations (Challenger and Dirzo, 2009).
The potential distribution of the Phyllostomidae family in Mexico showed a 30 % overlap with areas modified by human activities (table 3). The percent of overlap with the modified areas varied among feeding guild species; for example, convergence of higher richness of hematophagous and nectivorous species with modified areas differed with insectivorous guild, so they may be affected differentially ( fig. 4). Frugivo-rous bats will be less vulnerable to human activities because they have a higher niche tolerance and a lower marginality (e.g., genus Sturnira and Artibeus, table 2). These species are generalist, and they exploit a broad amplitude of resources in both preserved and disturbed environments (Klingbeil and Willig, 2009). Common vampire bats (Desmodus rotundus) have a negative interaction with humans because they cause economic loss to the beef industry and play a role in the epidemiology of bovine paralytic rabies in rural areas (Dantas-Torres, 2008). However, future climate projections predict an increase in the distribution of D. rotudus (Lee et al., 2012). In contrast, other endangered hematophagous species (i.e., Diphylla ecaudata and Diaemus youngi) that do not interfere with human activities are under risk through population control of the common vampire bat (Vercauteren et al., 2012). Insectivorous bats (e.g., Micronycteris microtis and Mimon cozumelae) are vulnerable to agricultural practices because they show a low tolerance and a high marginality (Medellín et al., 2000) and because they converge with croplands where humans regulate pests with chemicals that are toxic to bats (table 3; Lawer and Darkoh, 2016). The nectarivorous bats, due to their widespread distribution and high richness areas, overlap with human-modified areas ( fig. 4) and their distribution can indirectly affect their vital biological processes, such as the loss of reproduction areas, mutualistic associations and modified migratory movements (Fleming et al., 2009).
The spatial-environmental distribution of the Phyllostomidae family in Mexico reveals a species response to niche requirements and human changes. For example, species in homogeneous environments (e.g., tropical forests) have smaller niche breadth and narrow distributions, being more sensitive to environmental changes (Brown, 2014). On the other hand, species in heterogeneous environments (e.g., deserts) have broad niche breadth and widespread distributions with tolerance to disturbance (Krebs, 2001). In the country, the highest native cover loss occurs in tropical forests of the Atlantic and Pacific coastal plains, because low slopes and precipitation provide optimal conditions for the development of agricultural activities (Challenger and Dirzo, 2009), while the human population is concentrated mostly in Meseta Central (Klein-Goldewijk and Ramankutty, 2004), areas where phyllostomid bat richness is also high ( fig. 1). Medium (2023) and long-trend (2033) scenarios predict that if conservation actions are not implemented, an unstable to very critical scenario will occur in these areas (Sánchez-Salazar et al., 2013); in this sense, Protected Areas play a fundamental role in conserving biodiversity (Myers et al., 2000).
Mexico has decreed 181 Protected Areas (www. conanp.gob.mx) to maintain the integrity of ecosystems and environmental services. The decrees are supported by laws that regulate anthropogenic activities. However, most of these areas are surrounded by modified zones, are in mountainous regions, do not include high diversity areas, or show some degree of deterioration (Fuller et al., 2006). In the case of phyllostomid bats, nectarivorous and insectivorous bats are susceptible to constant environmental degradation caused by humans. Thus, conservation strategies may complement the surface of protected areas through the delimitation of areas with suitable environmental conditions to promote the persistence and flow of biodiversity (Nori et al., 2016). The approach used in this work can be extrapolated to other regions and taxa, which in turn will provide a better understanding of regional patterns of richness and species composition (Peterson et al., 2015). Our results showing the areas with the highest bat richness may be useful to propose biological corridors and conservation priority areas using approaches such as systematic conservation planning (Margules and Sarkar, 2009