Vulnerability of bat–plant pollination interactions due to environmental change

Plant–pollinator interactions are highly relevant to society as many crops important for humans are animal pollinated. However, changes in climate and land use may put such interacting patterns at risk by disrupting the occurrences between pollinators and the plants they pollinate. Here, we analyse how the co‐occurrence patterns between bat pollinators and 126 plant species they pollinate may be disrupted given changes in climate and land use, and we forecast relevant changes of the current bat–plant co‐occurrence distribution patterns for the near future. We predict under RCP8.5 21% of the territory will experience a loss of bat species richness, plants with C3 metabolism are predicted to reduce their area of distribution by 6.5%, CAM species are predicted to increase their potential area of distribution up to 1% and phanerophytes are predicted to have a 14% reduction in their distribution. The potential bat–plant interactions are predicted to decrease from an average of 47.1 co‐occurring bat–plant pairs in the present to 34.1 in the pessimistic scenario. The overall changes in suitable environmental conditions for bats and the plant species they pollinate may disrupt the current bat–plant co‐occurrence network and will likely put at risk the pollination services bat species provide.

decades with around one million species threatened with extinction (Díaz et al., 2019). Such biodiversity losses are thought to be due to changes in atmospheric and land use conditions, which pose a threat to the functioning of ecosystems and human well-being (Díaz et al., 2006;Travis, 2003). The tropics and subtropics have been strongly impacted by ongoing land use conversion, changes in precipitation patterns and strong extreme weather events (Malcolm et al., 2006;Myers et al., 2000;Sala et al., 2000). It has been shown that a changing environment may not only impact the distribution of species and reduce overall biodiversity levels (e.g. Aguirre-Gutiérrez et al., 2020) but can also change the composition of communities (Esquivel-Muelbert et al., 2019) and at the same time disrupt species interactions (e.g. plant-pollinator interactions; Diego et al., 2019;Gómez-Ruiz & Lacher, 2019). Importantly, the disruption of the interactions between species (e.g. plant-pollinator) in a community may also disrupt the functioning of the ecosystem and its Nature's Contributions to People (e.g. disruption of pollination services ;Díaz et al., 2015).
The effects that a changing environment has on ecosystems and on biodiversity per se are dependent on the functional trait composition of the species in the community and on how functionally diverse such communities are. Functional traits are those physiological, phenological and morphological characteristics of species that allow them to respond or have an effect on their environment (Violle et al., 2007). For example, the photosynthesis response to temperature, elevation, CO 2 concentration and water stress can be described in a plane curve having an optimum, where the photosynthesis is inhibited both at low and high values (Berry & Björkman, 1989), and the inherent ability for temperature acclimation of photosynthesis is different in C3, C4 and CAM plants (Kumar et al., 2017). The system C3 in plants is primitive and widely distributed in the geographic and taxonomic space, but tends to predominate in cooler habitats; C4 plants are thought to have originated in relatively arid regions in response to the low atmospheric CO 2 concentrations that arose sometime after the late Cretaceous, where high temperatures occur in combination with water stress; whereas desert CAM plants arose either in response to selection of increased water use efficiency or for increased carbon gain and now are adapted to drought in arid regions, where day and night temperatures can show a drastic shift, although the CAM system is also related to epiphytes in tropical environments (Ehleringer & Monson, 1993;Kumar et al., 2017).
Meanwhile, plant life forms are related to species responses to disturbances regimes and climatic factors (Cornelissen et al., 2003;Rowe & Speck, 2005). For example, the height and positioning of the foliage may be both adaptations and responses to grazing by different herbivores. Likewise, life form is associated with the relationship of the embryonic tissue to the ground surface that remains inactive during seasonal changes and then resumes growth with return of a favourable season. Since perennating tissue makes possible the plant's survival during an unfavourable season, the location of this tissue is an essential feature of plants adaptation to climate and disturbances (Cornelissen et al., 2003). In short, the different photosynthetic systems and the life forms of plants have a clear relationship with ecological interactions and ecosystem functionality.
Pollination is crucial for the functioning of ecosystems, it is responsible for the reproduction and evolution of angiosperms (flowering plants; Magnoliophyta), which represent ~90% of all worldwide land plants (Hernández-Hernández & Wiens, 2020), and it is of pivotal importance as a contribution to people's livelihoods (i.e. for crop production; Aizen et al., 2009). Yet, we are facing a pollination crisis as populations of insect and vertebrate pollinators are decreasing in abundance, occurrence and diversity across the world (Potts et al., 2010(Potts et al., , 2016Regan et al., 2015). The disruption of plant-pollinator interactions by a changing environment could thus have severe negative effects on ecosystems functioning, food security and human well-being at global scale (Potts et al., 2010(Potts et al., , 2016.
Bats are the most important mammal pollinators (Ratto et al., 2018), being responsible for the fertilization of at least 528 species of angiosperms distributed in 67 families; whereby contributing to ecosystem health, crop production and food security by pollinating ecologically and economically important plants through the tropics (Fleming et al., 2009;Fleming & Muchhala, 2008). There is no global assessment on the value of bat pollination services, but some studies have estimated that it generates revenues of US$117/ha for Durian in Indonesia (Sheherazade et al., 2019) and up to US$2,500/ha for the pitaya sector in Mexico (Tremlett et al., 2021). The tequila production, which is dependent on agave azul (Agave tequilana F.A.C. Weber), that is in fact bat pollinated, is responsible for profits of up to ~7000 million US dollar each year (Trejo-Salazar et al., 2016).
Moreover, bat pollination also plays a crucial role in maintaining the genetic diversity of plants, especially in fragmented and anthropized landscapes, as for instance for wild crop relatives such as banana (Musa spp.), agaves (Agave spp.) and dragon fruit (Hylocereus spp.; Williams-Guillén et al., 2016). Additionally, most wild plants having reproductive declines due to loss of pollinators are vertebrate pollinated (Ghazoul, 2005;Ratto et al., 2018), pointing to the possible effect of a decline or disruption on the distribution of vertebrate pollinators.
Bats are nowadays a matter of conservation concern as nearly 80% of bat species worldwide require either conservation or research attention (Frick et al., 2019). Environmental changes are projected to disproportionately affect tropical biodiversity, where bats are likely to experience considerable losses in their suitable habitat (Hughes et al., 2012;Zamora-Gutierrez et al., 2018). Jones et al. (2003 argue that geographic range sizes and wing morphology are among the most significant traits influencing bats extinction vulnerability. Species with small geographic ranges are less likely to escape threats, like environmental change (Pearson et al., 2014); while species with different wing morphology have different responses to changes in their habitat, in food availability, and even have different movement capabilities (Arita & Fenton, 1997;Norberg & Rayner, 1987). Other traits that could also influence bat responses to anthropogenic change are body size and reproductive rate, but their effect varies depending on several factors (Jones et al., 2003;Purvis et al., 2000).
Here, we investigate how the distribution of bat-pollinated plant species and bat pollinators will likely change as a result of changes in climate and land use conditions in the near future, and how the bat-plant co-occurrence patterns 'as a surrogate for possible interactions' may be disrupted under future environmental change. We use a large spatially explicit data set of observations of the presences of bats and the plants they pollinate for Mesoamerica, together with bat functional traits that are hypothesized relevant for their response to a changing environment. We then make predictions across Mexico as a main region of distribution of bat pollinators in Mesoamerica. We hypothesize that a changing climate and land use may cause a disruption in the co-occurrence patterns, and thus of the possible interactions, of bats and the plants they pollinate; that the strength of such disruption may depend on the plants metabolism and life form, and that bats responses to a changing environment, in climate and land use conditions, will be dependent on their intrinsic functional traits.

| Bats occurrence and trait data
We included all the species (N = 12) within the subfamily Glossophaginae distributed in Mexico as it includes all the bat pollinating species (Ramírez-Pulido et al., 2014). We used the occurrence database of Zamora-Gutierrez et al. (2018) which includes records from 1970 to 2014, but we updated their data set to include records up to 2020 following their same methods. We standardized scientific names following Ramírez-Pulido et al. (2014). After the cleaning process, we obtained 7458 bat occurrence records (Table S1). We obtained traits data from the PanTHERIA database (Jones et al., 2009; http://esapu bs.org/archi ve/ecol/E090/184/), published literature and unpublished data. We focused on traits that reflected extinction vulnerability, ecological plasticity, life history traits, population dynamics and dispersal/colonization potential (Jones et al., 2003;Pearson et al., 2014;Williams et al., 2008). The selected traits were litters per year, diet breadth, adult forearm length (mm), population group size, wing loading and wing aspect ratio (mm ; Table S2). We based the definition of these traits on Jones et al. (2009), except for wing loading and aspect ratio that followed the definitions used by Norberg and Rayner (1987; Data Table S1).

| Plant species selection, occurrences and trait data
We selected plant species following four criteria: (1) plant families presenting the highest richness of chiropterophilous species distributed in Mexico (Fleming et al., 2009); (2) within those families, we searched for all the species pollinated or visited by bats reported in the work of Fleming et al. (2009) and the database Bat-Eco Interactions (Geiselman & Sarah, 2020; http://www.batba se.org); (3) we complemented our search of literature in English and Spanish reporting plant species potentially pollinated by bats using the search engines of Google Scholar, Scielo and Web of Science using the search words "pollination", "polinización", "bats", "murciélagos", and "scientific name of each plant species"; and (4) from the plant species list generated, we confirmed that their flowers had characteristics of the chiropterophilous syndrome following the descriptions of Tschapka and Dressler (2002). After applying these selection criteria, we selected 126 species from the families Asparagaceae (N = 83), Bromeliaceae (N = 10), Cactaceae (N = 22) and Malvaceae (N = 11) to search for occurrence records (Data Table S2).
We searched for plant occurrence records through the Global Biodiversity Information Facility (GBIF; https://www.gbif.org), Portal  Table S2).
We searched for plant functional traits that could reflect environmental responses (hard traits) through the BIEN Botanical Information and Ecology Network and TRY Plant Trait databases (Kattge et al., 2020;Maitner et al., 2018). However, we did not find information for the plant species selected, so we had to focus on physiological and morphological traits (soft traits sensu Cornelissen et al., 2003). Soft traits are often good correlates of other hard traits, which may be more accurate indicators of plant functions responsible for responses or effects at broad ecological scales. We used two sets of soft traits: (1) plant life forms, which are associated with plant responses at environmental change, plant competitive strength and defence against herbivores and pathogens, and with effects on biogeochemical cycles and disturbance regimes (Cornelissen et al., 2003); (2) metabolism type, which is based on plant capacity to adjust their photosynthesis with climatic oscillations .  Table S2).

| Present and future environmental data
Climatic data from the present and future conditions were gathered from the Worldclim data set (Hijmans et al., 2005) with a spatial resolution of ~10 × 10 km at the equator. For the future conditions, we selected climatic data for the radiative concentration pathways (RCP) 2.6 and 8.5 to represent a sustainable and a worst-case scenario, respectively, for the year 2050. We included the RCP 2.6 as it assumes that CO 2 (carbon dioxide) levels start declining in the near future and get to zero by 2100, which may keep temperature from increasing above 2°C by 2100 (IPCC, 2014). As to investigate how a contrasting climate scenario may affect the distribution of species and their co-occurrence patterns, we also modelled the species distributions using the RCP 8.5 which assumes that CO 2 emission will continue rising through 2100 with a mean increase in temperature of 3.7°C. For the present and future conditions, we selected four main bioclimatic variables that have been shown to drive the distribution and persistence of plants MTemWarQua and MTemColQua depict the extreme of temperature across the year that may act as abiotic filters selecting for specific sets of species adapted to those temperature conditions.
MAnnPre and PreSeaso are key variables that depict the sum of water availability across the year but also how this varies temporally, which has been shown to strongly determine plant and ani- and obtained the proportion of each cell that is occupied by forest, shrubland, grassland and cropland area. We selected two land use projections to represent an optimistic and a pessimistic scenario: (1) a 'sustainable heaven' scenario based on a reduction on resources use, dependency on fossil fuels and deforestation within protected areas; and (2) a 'business-as-usual' scenario (for more details, see Zamora-Gutierrez et al., 2018).
We combined the two land use and the two climate change scenarios to obtain two combined extreme environmental change projections: (1) an optimistic combined scenario (RCP 2.6 and 'sustainable heaven' land use scenario) and (2) a pessimistic combined scenario (RCP 8.5 and 'business-as-usual' land use scenario). We used the same variables, spatial resolution and future scenarios for both bats and plants models. All analyses were carried out in the R statistical platform (R Development Core Team, 2013) using the raster package.

| Bats Joint Species Distribution Models
We modelled the bat species distributions using a Bayesian Joint Species Distribution Model (J-SDM) framework with the BORAL algorithm (Hui, 2016). We used J-SDMs as to be able to model the distribution of bat species also accounting for their co-occurrence patterns and functional traits. BORAL uses Bayesian Markov chain Monte Carlo methods and incorporates latent variables to model between species correlation (significant residual correlations). BORAL fits a generalized linear model using the species co-occurrence matrix, in this case of bat species, as response variable and the climatic and land use covariates as explanatory variables. We included the bat traits specified in the section above into the model to explain differences in how species respond to environmental conditions. This feature is of importance as we are interested in the way species composition changes across sites depending on the species traits available and the underlying environmental conditions. We modelled the bat species distributions using the same climatic and land use variables used to model the plant species. However, we here account for the co-occurrence patterns and functional traits of the bat species modelled as to improve their spatial predictions and gather insights about the trait-environment relationships which can at least partially determine where species occur. For the BORAL models we used 10,000 iterations to discard at the beginning of the MCMC sampler, we carried out 40,000 iterations and the thinning rate was set to 30. We used default uninformative priors for model fitting. We used 75% of the species presence data for model fitting and the remaining 25% for model testing and carried out 10 model repetitions for each algorithm. We transformed suitability scores into presence/ absence maps of each species using the 10% suitability threshold (Pearson et al., 2007). This threshold omits cells with habitat suitability lower than the suitability scores for the lowest 10% of occurrence records, and thus is a conservative manner of thresholding continuous suitability map predictions.
All analyses were carried out in the R statistical platform using the BORAL package.

| Plant species distribution modelling
As there is not widely available plant trait data for the species modelled, and building J-SDMs for all plant species together across the study area is not computationally feasible given the high number of species modelled (>100), we only modelled the species distributions based on land use and climatic conditions. We modelled the plant species distributions (SDM's) using the biomod2 package v. 3.1-64 (Thuiller et al., 2009). We fitted the SDMs using three modelling algorithms, generalized linear models (GLM; McCullagh & Nelder, 1989), Random Forest (RF; Breiman, 2001) and Maximum Entropy (MaxEnt; Phillips et al., 2006).
We used 75% of the species presence data for model fitting and the remaining 25% for model testing and carried out 10 model repetitions for each algorithm. The selected algorithms, except MaxEnt, make use of pseudo-absences for model fitting; therefore, we generated a random sample of pseudo-absences as five times the number of presences used. Then, we followed an ensemble modelling approach selecting only models with AUC (area under the curve; Hanley & McNeil, 1982) values above 0.5 from the three modelling algorithms to obtain more robust predictions for the potential distributions of the species (Araújo & New, 2007). We used the resulting 30 models from the three modelling algorithms as the basis for the final ensemble model for which we took the median prediction value as the ensemble rule. In this way, we avoid selecting a single model based only on their single AUC scores, which could render bias predictions (Aguirre- Gutierrez et al. 2013). The resulting ensemble model fitted under the present conditions was then projected to the future climatic and land use conditions (described above). We transformed suitability scores into presence/absence maps of each species using the 10% suitability threshold (Pearson et al., 2007).

| Potential bat-plant interaction and species richness data
The resulting presence/absence maps were aggregated to produce species richness maps for bats, plants and plant groupings by morphological and physiological characteristics. At each map cell, all occurring species of the given species group were counted. Similarly, we built potential interaction maps by summing all co-occurring batplant pairs in every single map cell. We assumed that the simultaneous presence of a bat and plant species from our restricted list is evidence that a pollination interaction could occur. The resulting maps were used to compare species richness and potential interactions across the different scenarios.
To estimate the large-scale territorial extent of an interaction between a bat-plant species pair, we tested different metrics. Firstly, we computed the area of distributional overlap, which provides a very rough estimate of the frequency of their interaction across all of Mexico. Secondly, we measured the proportion of the bat's distribution simultaneously covered by the plant's distribution. Similarly, we computed the proportion of the plant's distribution in which the bat is present. These proportional or relative measures of distribution intersections are especially suited for comparing the extent of plant or bat co-occurrence with respect to a particular bat or plant species. Finally, the ratio of the overlap area to the area covered by both species was tested as a measure of concordance of distributions (Nori et al., 2015). It is important to consider that our approach on interaction assessment does not reflect other types of interactions beyond those reflected by sympatric spatial distributions, and that the influence on species preferences is also not accounted for within this spatial analysis.

| Impacts of environmental changes on bats and plants distributions and richness
Overall the climatic variables related to temperature and precipitation appeared to be more important for determining the distribution of bat species than land use covariates ( Figure S1). The overall trend in climatic changes across Mexico is that temperatures are likely to increase for both the coldest and warmest quarters of the year across both scenarios (RCP 2.6 and RCP 8.5), precipitation is likely to decrease in scenario RCP 8.5 with more seasonal precipitations (a higher coefficient of variation; Figure S2a). Geographically, the biggest changes in temperatures are projected to be in the arid and semiarid regions, decreased precipitations are estimated along the west and east mountain ranges and more precipitation seasonality might be expected in the temperate and arid regions ( Figure S3).
Regards land use variables, the cover that had the biggest change was forest with a deep decrease in the business as usual scenario of change, followed with an increase of shrubland for the same scenario ( Figure S2b).
The bat species that lost a bigger proportion of their distributions under the pessimistic scenario were Choeroniscus godmani (−11%), Musonycteris harrisoni (−10%) and Lichonycteris obscura (−10%, Table 1). The highest current nectarivorous bats richness is predicted along the Mexican Pacific Coast, where also the highest losses of species richness are predicted for the pessimistic scenario ( Figure S4). Overall, 21% and 18% of the territory showed a loss of bat species richness, based on the pessimistic and optimistic scenarios, respectively (Figure 1), but most of the territory will maintain the same number of bat species per area (Figure 2). Some areas will gain above 100% of their current number of nectarivorous bat species based on the pessimistic scenario, like some portions of the Yucatán Peninsula ( Figure S4).  under the pessimistic scenario, while chamaephyte will likely show an increase in their distribution in both scenarios (Table 1).
For the analyses of all plants, the highest current species richness is predicted along the Sierra Madre Occidental, but this pattern of species richness varies depending on the plant group ( Figure S5). When all plants are analysed together, most of their territory (up to 71% for the pessimistic scenario) will show declines in the number of species (Figure 1), with the pessimistic scenario also showing the greatest shifts towards more area with less species ( Figure 2). The highest losses are likely to happen across the Mexican Plateau and north of the Sierra Madre Occidental ( Figure   S5). Based on plant metabolism, CAM plants will have the greatest losses on richness across their distribution (68%, Figure 1) with the pessimistic scenario showing the greatest shifts towards more area with less species (Figure 2). Those losses will be higher across most part of north Mexico ( Figure S5). For the case of plant life form, the rosulate chamaephyte group will lose species richness across 70% and 58% of their distributions based on the pessimistic and optimistic scenarios, and will lose all the species in 3% and 2% of their territory respectively (Figure 1). This group will also have the greatest shifts towards more areas with less species in the pessimistic scenario (Figure 2). Rosulate chamaephyte plants are likely to have severe losses in most of the territory except in the Yucatán Peninsula and in some parts of the Gulf coast and north of Sonora, which coincide with the regions with low current species richness ( Figure S5). Chamaephyte plants will have a total loss of species in 20% and 14% of their territory based on the pessimistic and optimistic scenarios, yet will also have the greatest gain in new territory mostly in south Mexico ( Figure S5). For the rest of the plant groups, most will lose richness across most of their distributions, rather than gaining species (Figures 1 and 2).
In both future scenarios, plant distributions showed a general elevation increase ( Figure S6a), although this effect is greater in the pessimistic scenario. Terrain at higher altitudes is usually constrained to narrow and ragged spaces thus restricting resident species to distribute along territories with these geometric characteristics. Consequently, the positive shift in mean altitude of some plant species distribution was generally accompanied by a positive change in the perimeter-area ratio, especially in the pessimistic scenario ( Figure S6b). The changes in area of distribution showed

F I G U R E 3
The plots show differences in the territorial distribution of diversity of bat-plant potential interactions in the present and future scenarios for each plant species group. Top panel: differences in territorial distribution for bat-plant potential interactions; Middle panels: differences in territorial distribution for Bat-C3 (left) and Bat-CAM plant potential interactions; Bottom panels: differences in territorial distribution for Bat-Chamaephyte, Bat-Phanerophyte, Bat-Rosulates, Bat-Succulents lianas and Bat-Succulent phanerophytes potential interactions a significant negative correlation with mean altitude, implying that distribution shrinkage is in part driven by species shifting towards higher altitudes ( Figure S6c).

| Present and future bat-plant potential interactions
Changes in bats and plants distributions had a substantial impact on their potential interactions (Figure 3). When comparing ranges for potential interaction of plant-bat pairs across the different scenarios, as measured by overlap of distribution area, predictions are equally dire. In both scenarios, 72.3% of plant species will see an overall decrease of bat interaction area (Table 1)

| DISCUSS ION
We analysed how a changing environment may impact the co- We showed that all plant groups analysed are shifting up in elevation tracking more suitable areas (see also Chen et al., 2011;Lenoir et al., 2008), which might be one of the causes of range contraction since mountain tops have a limited geographic extent. These shifts towards cooler environments could be partially explained by the changes observed in temperature and precipitation variables across the country, where temperate habitats mostly located in highland areas offer and are likely to continue to offer colder and moister environments compared to the lowlands Rzedowski, 2006). In fact, previous studies exploring the effects of  (Borges et al., 2016). Some of these plant groups find their greatest diversity in xerophilous scrub and seasonally dry tropical forest, where it is frequent to face water stress in the habitat and a selective advantage by nocturnal flowering. It is in these areas containing plants with high adaptations to be bat pollinated that will suffer the highest nectarivorous bats losses. Previous studies already highlighted these areas as potential hot spots of risks for biodiversity losses (Peterson et al., 2002;Zamora-Gutierrez et al., 2018). Although pollinating networks have some mechanisms to provide them with resilience to pollinators loses (Hegland et al., 2009), even the removal of one pollinating species from the community could result in changes on the whole network structure in terms of species composition, resources use and visitation rates with subsequent consequences on the biodiversity and the general functioning of the ecosystem (Brosi et al., 2017;González et al., 2009;Miller-Struttmann et al., 2015). Moreover, although nectarivorous bats could potentially pollinate a wide range of plants because of their generalized diet, species differ widely in their home ranges, thus in their abilities to connect different plant populations and to maintain gene flow (Rothenwöhrer et al., 2011).
Complete plant groups like CAM (agave and columnar cacti) and rosulate chamaephyte (agave), species adapted to arid and semiarid habitats, will lose local diversity in terms of species richness, in over 60% of the national territory, even though their aggregated distribu- These loses on plants adapted to dry environments are very likely to cause a drastic reduction in the diversity of food sources available for bats, especially for migratory species which mainly feed on these plant groups (Burke et al., 2019). This is very likely to result in a severe decrease of foraging sources for this species (Gómez-Ruiz & Lacher, 2019). Most of the bat species here analysed are cave-dwelling bats and thus, at different degrees, depend on fixed roosting sites (Gomez-Ruiz et al., 2015). This dependency is very likely to limit their ability to track the distribution changes of their food sources. Changes in plant distributions and richness could also cause a decrease in plant population densities with a subsequent decrease in food abundance.
Due to the particularly high metabolic rate of nectarivorous bats, their population dynamics and reproductive cycles are highly correlated with food availability (Peñalba et al., 2006;Stoner et al., 2003;Winter & von Helversen, 2003). Several studies across the Neotropics have recorded a high synchronization between peak abundances of nectarfeeding bats populations and their reproductive bouts with the peak abundance in flower resources availability (Heideman & Utzurrum, 2003;Stoner et al., 2003;Valiente-Banuet et al., 1996). Some bat communities have a seasonal structure driven by food abundance and the quality of the food sources available at different times (Tschapka, 2004). Changes in the structure of plant assemblages caused by environmental changes may have serious consequences on the competition between bat species thus altering the spatio-temporal dynamics of bat communities.
Another possible impact of environmental changes on nectarivorous bats are the alterations on the phenology of plants (Hegland et al., 2009;Memmott et al., 2007). Previous studies have revealed that even generalist pollinators could face a considerable loss of food sources due to plant phenological changes (Memmott et al., 2007). An observed trend in phenological changes is that flowers tend to bloom earlier with higher temperatures and animal-pollinated plants are more affected by this (Fitter & Fitter, 2002 Our study contributes to the existing knowledge gap on how species traits could potentially modulate species interaction, in this case pollination (Schleuning et al., 2020). We found that all the traits modelled here have a higher response over climatic than land use variables. Geographic range has been a predictor constantly highlighted as an important determinant of extinction risk in mammals, and particularly in bats (Jones et al., 2003;Pearson et al., 2014). Although we did not include this trait as part of the analysis as most of the modelled species are widely distributed across Mexico, the species that lost more distributions were those with smaller ranges. The thermal niche of species is also one of the most important physiological traits, which determines species responses to climatic changes (Gunderson & Stillman, 2015). Climatic models project a range of temperature increase of around 1. Wing morphology reflects an important ecological characteristic that can influence species vulnerability to anthropogenic changes like guild, habitat and dispersal capacity (see Jones et al., 2003). Low wing aspect ratios are linked to high extinction risks because they reflect higher energetic cost when flying; lower dispersal abilities and home ranges, and species with low wing aspect ratios also tend to have smaller colony sizes (Jones et al., 2003;Norberg & Rayner, 1987). In this study, species with high wing loading and aspect ratio are distributed towards the Nearctic regions of the country. Several of these species are migratory and can move long distances following the seasonality of nectar resources across the country. They also tend to have larger geographic ranges which further reduce their vulnerability to environmental changes (Jones et al., 2003). On the contrary, smaller pollinating bats with low wing loading and aspect ratios are distributed in the tropics (south Mexico). These species are mainly adapted to warm humid temperature with low seasonality.
Diet breath was greater in those species distributed in northern regions. Species can buffer the effects of climate change through behavioural changes and more generalist species have more ecological flexibility to adapt to those changes (Slatyer et al., 2013).
Species with more specialized diets are also more sensitive to plant phenological changes and changes in plant community compositions (Memmott et al., 2007). Therefore, Neotropical nectarivorous bats may be at higher risk of climatic changes due to their more specialized diets and in general due to their smaller sizes and wing morphology that do not allow them to follow changes in their food sources over extended geographical regions.
Ideally to incorporate the effects of a changing climate on the distribution of pollinating bats, a model accounting for the bat interactions with the plant species could be constructed. However, given the nature of the Bayesian modelling protocol followed which is highly resource demanding, and the number of plant species included in our study and all possible interactions with the bat species, we decided not to develop the models this way. We therefore use the modelled potential distribution of the bats and plants and their intersection as a proxy of possible interactions and pollination potential. This approach may actually overrepresent the possible interactions and thus underestimate pollination loss than expected if using actual plant-pollinator interaction networks. Therefore, our results can be seen as conservative in the effects that a changing climate and land use may have on the plant-bat pollination services because species that co-occur might not interact at all. As earlier suggested by Davis and collaborators (Davis et al., 1998), we find that predicting biotic responses to climate change alone and without including the effects of biotic interactions gives an incomplete picture of how species may respond to a changing environment. This is simply because species do not act alone and depend on more than just climatic conditions for their survival and persistence (Davis et al., 1998;de Souza Laurindo et al., 2017). Here, we tend to account for part of those biotic interactions by investigating the cooccurrence patterns between the different bats species modelled using the J-SDMs and by then analysing its co-occurrence with the plant species they pollinate.
In sum, we show that a changing environment may cause a disruption in the co-occurrence patterns, and thus of the possible in- Lim for providing information of bats traits and to Pedro Aguilar-Rodríguez for providing Bromeliaceae species occurrence data.

AUTH O R CO NTR I B UTI O N S
VZG and JAG conceived the study, designed and carried out the analysis. VZG, JAG, ANRV, ACC gathered species and traits data.
VZG, JAG and SMB carried out data analysis. VZG and JAG wrote the first draft of the paper and all co-authors commented on and approved the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The lists of bat and plant species used in this study can be found as supplementary information files. The climatic and land use data that support the findings of this study are openly available from their sources.