Species distribution modelling supports “nectar corridor” hypothesis for migratory nectarivorous bats and conservation of tropical dry forest

The Mexican long‐tongued bat (Choeronycteris mexicana), Mexican long‐nosed bat (Leptonycteris nivalis) and lesser long‐nosed bat (Leptonycteris yerbabuenae) (Phyllostomidae: Glossophaginae) undertake long‐distance migrations from south‐central Mexico to the south‐western United States. It is proposed that these bats migrate along a nectar corridor of columnar cacti and Agave species, but this has not been tested with independent data and the spatiotemporal nature of this relationship is poorly understood. Our goal was to test this nectar corridor hypothesis and determine the relative importance of food plant and abiotic variables to the distribution and seasonal movements of these migratory nectarivores.


| INTRODUC TI ON
Migration is the predictable, seasonal and synchronized movement between two regions in which different life-history events take place in order to maximize fitness (Dingle, 1996;Dingle & Drake, 2007;Greenberg & Marra, 2005;Lennox et al., 2016). Across a wide range of taxa, migration is favoured over residency when the distribution of resources changes seasonally (Dingle & Drake, 2007;Shaw & Couzin, 2012). Though migration takes place throughout the animal kingdom, most insights are from birds, fishes and insects (Dingle stable isotope analysis. Fleming et al. (1993) discovered that L. yerbabuenae feeds on plants of the C3 photosynthetic pathway in the southern portion of its range during the winter months and on those of the Crassulacean acid metabolism (CAM) photosynthetic pathway as they become available throughout the spring and summer months in the Sonoran and Chihuahuan deserts. Using mitochondrial DNA analysis, Wilkinson and Fleming (1996) theorized that L. yerbabuenae follows divergent migration routes along the Pacific coast of Mexico and along the foothills of the Sierra Madre Occidental. Rojas-Martínez et al. (1999) refined the theory using geographic analysis of plant phenology and capture records, suggesting that L. yerbabuenae populations occurring north of 30° migrate, while those around 21° are likely resident.
Though migratory routes ought to be spatially and temporally predictable , we still lack a spatially explicit framework with which to test the theory of generalized latitudinal migration. Hence, the overall influence of forage resource phenology on nectarivorous bat distribution and migration has not been substantiated. Simultaneously, recent extralimital occurrence records of all three species at the northern edge of their range suggest these bats may have a broader range than previously described (Bogan & Cryan, 2006;Ramsey & Whiteman, 2011;Bogan et al., 2017). A better understanding of the distribution and potential migratory pathways among seasonal habitats would inform conservation efforts for these bats.
Our objectives were to: (a) model the overall distribution of C. mexicana, L. nivalis and L. yerbabuenae to gain a better understanding of their geographic distribution and to test the relative importance of food plants vs. abiotic variables to their distribution; and (b) model seasonal habitats and migratory pathways between seasonal habitats to test the theory of latitudinal migration based on available food resources and to test the presence of divergent migration routes. We used presence-only data with species distribution and corridor modelling tools to test this theory. The models generated in this study shed light on the ecological relationships governing the distribution of migratory nectarivorous bats, enhance our ability to monitor populations throughout their range by providing a spatiotemporal framework with which to prioritize habitat conservation and allocate survey efforts, and provide a spatial model with which to test theories on migratory strategies.

| Study area
Our study area encompassed the entire potential migratory range of C. mexicana, L. nivalis and L. yerbabuenae including all of Mexico and the south-western United States. This region is typified by a transition in plant communities and increasing plant richness along elevation, temperature and humidity gradients (Sánchez-González & López-Mata, 2005). Elevations range from sea level to about 3,190 m (Silva-Flores, Pérez-Verdín, & Wehenkel, 2014). From the tropics of southern Mexico to the Sonoran and Chihuahuan deserts, biomes shift from tropical evergreen, semideciduous and deciduous forests to thorn forest and xeric shrubland with Neotropical savanna and warm temperate grasslands interspaced with sky islands of coniferous forests and Madrean montane grassland (Brown & Makings, 2014). Within these biomes, plant communities are highly variable based on aspect and soil characteristics, with high diversity of xerophytic plants along limestone escarpments (Gehlbach, 1967;Sánchez-González & López-Mata, 2005).

| General approach
We used maxent 3.3.3 to generate species distribution models (SDMs) (Phillips, Anderson, & Schapire, 2006). Species distribution modelling is widely used in ecology to better understand the relationship between organisms and environmental variables, shedding light on how this relationship influences the likelihood of presence across the landscape (Elith & Leathwick, 2009;Elith et al., 2010).
While there are many methods with which to generate SDMs (Guisan & Thuiller, 2005), maxent has become the most popular for both its ease of use and functionality (Morales, Fernández, & Baca-González, 2017). maxent is a machine learning algorithm that allows SDMs to be generated using presence-only data, making it an effective tool for predicting species distribution when obtaining presence-absence data is logistically impractical (Clemente et al., 2019;Kabir et al., 2017;Zhang, Yao, Meng, & Tao, 2018).
Modelling the distribution of species throughout multiple seasons at a continental or global scale can help reveal important patterns of seasonal habitat preferences at a scale more relevant to animal migration (Fink et al., 2010). Many studies assess the relationship between species occurrences and climatic variables to predict potential distribution across the landscape (Abolmaali, Tarkesh, & Bashari, 2018;Bosso et al., 2018;Thapa et al., 2018), but the specialized diets of nectarivorous bats may render food plants more biologically relevant to their distribution (Arita & Santos-del-Prado, 1999). Hence, we used the distribution of plants that have been documented in the diet of these bats as explanatory variables, in addition to more commonly used bioclimatic, topographic and land cover variables, to generate SDMs of C. mexicana, L. nivalis and L. yerbabuenae. We used the SDMs to test the relative importance of these variables to each bat and to predict migratory corridors by establishing resistance surfaces derived from the SDMs, which can adequately estimate wildlife corridors when more precise information, such as telemetry data, is unavailable (Zeller et al., 2018). We chose to use circuit theory over least-cost methods because although nectarivorous bats have strong spatial memory (Henry & Stoner, 2011), shifting phenology of food resources would prevent these bats from always following the most efficient route along the migratory corridor.

| Predictor layers
We used 51 food plant variables, along with food plant richness, bioclimatic variables, topographic and land cover variables to model the distribution of C. mexicana, L. nivalis and L. yerbabuenae. These included 19 bioclimatic variables (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005), distance to karst (Weary & Doctor, 2014), elevation (GTOPO30;Earth Resources Observation and Science Center 1997) and land cover (Colditz et al., 2014), amounting to 76 initial variables. Some studies have excluded bioclimatic variables that combine precipitation and temperature variables due to odd spatial patterns (Escobar, Lira-Noriega, Mdeina-Vogel, & Peterson, 2014). However, we chose to include these variables because we did not observe any odd spatial patterns and because of their combined influence on agave and cactus physiology (Nobel, 1988) and on bat foraging and migration (Frick et al., 2012;Pettit & O'Keefe, 2017). As distance to variables may improve model performance for central-place foragers such as bats (Rainho & Palmeirim, 2011) and nectarivorous bats rely on a matrix of foraging grounds and roosting habitat, we used karst as a surrogate variable for roosting habitat by creating a layer of distance to karst using the Euclidian distance tool in arcgis 10.4.1 (ESRI, 2016). We also used elevation and land cover as predictor variables due to their influence on habitat selection and the supposed segregation of these species based on elevational gradients (Arita, 1991;Baker & Cockrum, 1966).
We also generated food plant richness models from the SDMs to use as explanatory variables (Zhang et al., 2016). Species richness of food plants is positively associated with nectarivorous bat presence (Gómez-Ruiz & Lacher, 2016). We generated food plant richness models by applying a presence threshold to each food plant SDM using the maximum sum of sensitivity and sensitivity method to create binary presence-absence models (Liu, Newell, White 2016). We then used the raster calculator in arcgis 10.4.1 to aggregate binary models into an overall food richness model. We generated seasonal food richness models based on the seasonal availability of each food plant to nectarivorous bats to later use as covariates in our winter, spring, summer and fall seasonal bat distribution models (see Appendix S1).

Choeronycteris mexicana Leptonycteris nivalis
Leptonycteris yerbabuenae n, total n, rarefied n, total n, rarefied n, total n, rarefied TA B L E 1 Sample sizes of total occurrence records and spatially rarefied (50 km) occurrence records for overall species distribution models and seasonal distribution models

| Bat occurrence records
Because Leptonycteris spp. are frequently misidentified (Arita, 1991), we only used specimen records of this genus that were examined and verified by us or were included in the revision by Arita and Humphrey (1988). We differentiated L. nivalis and L. yerbabuenae using the length of the last digit of the third phalange (Hinman & Snow, 2003;Medellín, Arita, & Sanchez, 1997 (2011), which included measurements that allowed for reliable identification, as well as records in Bogan et al. (2017) and Bogan and Cryan (2006).

Choeronycteris mexicana is readily distinguished from Leptonycteris
based on the presence of a skirted uropatagium (Medellín et al., 1997 Cryan and Bogan (2003).
To maintain adequate sample sizes, we did not discard bat occurrence records with low precision because this has limited impact on model performance for volant animals with high dispersal potential (Hayes, Ozenberger, et al., 2015). However, many of our records were spatially clustered. To limit potential sampling bias, we spatially rarefied occurrence records for each species to 50 km to approximate daily dispersal distance (Horner et al., 1998;Medellín et al., 2018). This resulted in sample sizes which were sufficient for generating reliable maxent models (Van Proosdij, Sosef, Wieringa, & Raes, 2015) and allowed minimally for the implementation of hinge, linear and quadratic feature classes (Merow et al., 2013;Phillips & Dudík, 2008;Scheglovitova & Anderson, 2013; see Table 1).

| Generating bat SDMs
We defined the background extent for each SDM as all of Mexico, Texas, New Mexico, Arizona and California south of the latitude of the northern border of Arizona. We did not restrict this extent further because these highly mobile bats could theoretically reach any given area in the study region given their migratory capabilities of 1,200 km (Cockrum, 1991;Moreno-Valdez et al., 2000). Simultaneously, recent extralimital records of each species indicate poorly documented distributions and a broad background extent prevented us from inadvertently truncating potential ranges. We followed a series of steps to tune each model. First, we removed highly correlated variables (r > |0.90|) a priori based on biological relevance to limit redundancy (Fourcade, Engler, Rӧdder, & Secondi, 2014). We followed the methods of Warren et al. (2014)  After selecting the final variable set, we retuned the β parameter by generating 76 more models with β values ranging from 0 to 15 adjusted in increments of 0.2 and selected the optimal model using AIC c . Once final parameters were established, we bootstrapped a final model with the number of replicates equal to 20 when n > 20, or equal to n when n < 20 (Scheglovitova & Anderson, 2013).
We evaluated model performance using the area under the receiver operator characteristic curve (AUC), true skill statistic (TSS), sensitivity and specificity. To maintain sample sizes, test data were not separated from training data. To limit bias due to spatial clustering, we used a block-partitioned variation on k-fold cross-validation using the package "enmeVal" in program R to obtain AUC values (Muscarella et al., 2014).
We used TSS values from Landis and Koch (1977) to evaluate perfor-

| Seasonal models
To better understand shifts in habitat preferences throughout the year, we divided occurrence records by season to generate seasonally dynamic SDMs Smeraldo et al., 2018).
We divided records based on the month of the occurrence record. seasonal models for L. nivalis. However, our sample sizes for C. mexicana and L. yerbabuenae were adequate for each season of the year (see Table 1).
We generated seasonally dynamic models  for C. mexicana and L. yerbabuenae following the same model tuning and model evaluation methods as for the SDMs.
However, we only used predictor layers of food plants that would be available to a nectarivorous bat in the season being modelled.
We considered a plant available if it flowered in the given season.
We also considered cacti to be available if the species produced fruit in the season being modelled, given the importance of cactus fruit to migratory nectarivorous bats (Hernández & Herrera, 2016;Valiente-Banuet et al., 1996). We used herbarium records and the literature to determine seasonal availability (see Appendix S1). We also included a plant richness layer of seasonally available food plants, rather than overall plant richness. We evaluated seasonal models using block-partitioned AUC metrics, along with TSS, sensitivity and specificity.

| Predicting migratory corridors
We used circuitscape, a software which utilizes circuit theory Because we wanted to model connectivity between specific seasonal habitats rather than overall landscape connectivity, we averaged seasonal models for two seasons to represent the conductance of the landscape between those seasons. For example, if we wanted to estimate landscape conductance between spring and summer habitats, we averaged the spring and summer seasonal models. We then used spatially rarefied presence locations for the given seasons to use as focal nodes. Using resistance or conductance values of grid cells surrounding the focal nodes, circuitscape then produces a series of networks by which the species may travel between seasonal habitats based on its likelihood of crossing a given grid cell. We used the logistically transformed cumulative current map to denote potential migratory pathways for C. mexicana and L. yerbabuenae, with higher values representing more likely travel corridors.

| Bat SDMs
We generated SDMs with satisfactory performance based on their block-partitioned AUC evaluation metrics all >0.8 and TSS values of all moderate or substantial (see Table 2). Tuned β values were higher than maxent's default value of 1 (see Table 2). Variables in the final bat SDMs consisted of a mix of food, bioclimatic and topographic variables, with individual food plant variables and/or food plant richness variables having high contribution compared to other variables overall (see Table 3). The output of the SDM is a relative occurrence rate (ROR), which can be interpreted as the likelihood of occurrence or as environmental suitability. Given our use of potential food plant distribution as predictor variables, we interpret higher ROR values as reflecting a higher environmental suitability (Elith et al., 2010). In the SDMs for all three species of bats, there is high environmental suitability throughout south-central Mexico that varies with increasing latitude (see Figure 1).

| Seasonal models
Results of the seasonal models for C. mexicana and L. yerbabuenae were more nuanced than the overall SDMs (see Table 3

| Migration models
The

| Post hoc analysis
Given the importance of species richness of food plants in the SDMs for all three species of bats, we ran a post hoc analysis to test the relative importance of the components to this richness: agave richness, cactus richness and C3 plant richness. We did this by creating SDMs for each bat species using only these three richness variables (see Appendix S2). Each species exhibited varying reliance on the richness of these three plant groups, especially the agaves and cacti. The distribution of L. nivalis was only influenced by the richness of agaves. The distribution of Choeronycteris mexicana also was strongly influenced by the richness of agaves (63% contribution), but also moderately influenced by the richness of cacti (37% contribution). In contrast, the distribution of L. yerbabuenae was primarily influenced by the richness of cacti (85.8% contribution) with small influences by the richness of agave (7.7% contribution) and C3 plants (6.5% contribution). towards maternity roosts and eventual lactation (Fleming, Hooper, & Wilson, 1972;Hernández & Herrera, 2016). Similar to the findings of  with respect to continental scale bat migration, our models indicate that the distribution of migratory nectarivorous bats is the most uncertain and widespread during migration periods in spring and fall, especially for C. mexicana, highlighting the potential for increased vulnerability during these time periods, as well as the need for a better understanding of migratory corridors.
The migration models support the theory to preserve imperilled tropical dry forest ecosystems (Janzen, 1988;Quesada et al., 2009). Sustained plant richness in this ecosystem depends on the ability to preserve remaining tropical dry forest fragments via coordinated conservation efforts across the region (Banda et al., 2016). Important food plants for each bat species face threats from rapid deforestation of tropical dry forests (Janzen, 1988;Trejo & Dirzo, 2000), overutilization (Nabhan & Fleming, 1993), altered fire-return intervals from invasive species introductions (McDonald & McPherson, 2013) and climate change (Dávila, Téllez, & Lira, 2013); columnar cacti are especially susceptible to these threats. Given the strong importance of columnar cacti richness to L. yerbabuenae, mitigation of these threats may be especially important for the full recovery of this recently delisted species.
The SDMs and migration models from this study can be overlaid with protected lands databases to identify potential gaps in landscape connectivity, allowing researchers and land managers to prioritize landscape preservation and restoration efforts. The SDMs can also be used to allocate survey efforts to obtain more occurrence records for these species. Using the existing SDMs to target monitoring efforts could lead to new detections (Rebelo and Jones, 2010), allowing seasonal and migration models to be generated for L. nivalis and to be continually improved for C. mexicana and L. yerbabuenae.
The ability to target conservation efforts along migratory corridors for L. nivalis could be invaluable to the conservation and recovery of this species. Generating seasonal models for L. nivalis would allow us to more fully understand shifting habitat needs. A better understanding of migratory pathways for this species can help researchers identify stopover habitat, prioritize conservation of roosts and foraging grounds, inform restoration of food plants throughout migratory corridors and seasonal habitats, and identify flyways that could be impacted by human development. We recommend the use of these models to promote an overall transboundary approach to monitoring these species throughout their predicted ranges to enhance the conservation of these nectarivores and their habitats. Paso Biodiversity Collections for hosting us at museum collections during this project. We thank two anonymous referees for helpful comments that improved the paper.

DATA ACCE SS I B I LIT Y
Environmental layers used as predictor layers in this study are publicly available and sourced in the main text of this paper. All occurrence records used to generate the SDMs in this study (plant and bat SDMs) are publicly available via museum collections, publications or reports.
The authors can provide raster grids or kml files of SDMs or migration models generated in this study to land managers if requested, but due to the sensitive nature of these species, the models generated in this study have not been made publicly available.