Analyzing the adaptive potential of the long-nosed bat, Leptonycteris nivalis, in the face of climate change

Roberto-Emiliano Trejo-Salazar1,2 , Jaime Gasca-Pineda*1,3 , Rosalinda Tapia López3 , and Luis E. Eguiarte3 .

1Secretaría de Ciencias, Humanidades, Tecnología e Innovación (SECIHITI), Estancias Posdoctorales por México. E-mail: remilianotrejo@ciencias.unam.mx (RET-S)

2Departamento de Biología Evolutiva, Facultad de Ciencias, Universidad Nacional Autónoma de México (UNAM).

3Departamento de Ecología Evolutiva, Instituto de Ecología, Universidad Nacional Autónoma de México. Ciudad Universitaria, Ciudad de México, México. E-mail: rtapia@ecologia.unam.mx (RTL); fruns@unam.mx (LEE)

*Corresponding author: jaimegasca@yahoo.com

Tequila or magueyero bats (i.e., the genus Leptonycteris) play a very important role as pollinators in most Mexican and southern US ecosystems, which is why they have recently been given greater attention. Female members of Leptonycteris nivalis populations exhibit migratory behavior associated with nectar consumption and pollination of ecologically important Agave species in Mexico, and with their reproductive events. Such migratory behavior makes them susceptible to the consequences of climate change. In this study, we explored how the distribution and the adaptive potential of L. nivalis may be affected by possible global warming scenarios. We obtained samples from six locations across the distribution of L. nivalis. Using parallel massive-sequencing techniques, we obtained 13,112 filtered SNPs (single-nucleotide polymorphisms). Just two of these SNPs were associated with a climatic variable, but we detected associations between heterozygosity and climate change in two sites. Models of future distribution under climate change indicated a reduction in climatic stable areas favorable to this bat’s presence. Our results will support the implementation of protection strategies for bats and the ecosystems they inhabit, which are the most representative in Mexico. Also, it is very important to maintain connectivity between the localities and refuges occupied by the magueyero bat, as well as to maintain populations within areas of greatest thermodynamic stability, where the most favorable conditions for bat presence are predicted. Based on our analysis, we believe the species has the potential to withstand changes in temperature and precipitation patterns, which may support its conservation.

Keywords: Chihuahuan Desert, climate change, conservation genomics, Mexico, nectar-feeding bats, potential adaptive, potential distribution, population genomics, SNP

Los murciélagos del género Leptonycteris, a veces llamados tequileros o magueyeros, desempeñan un papel muy importante como polinizadores en la mayoría de los ecosistemas mexicanos y del sur de Estados Unidos, por lo que recientemente se les ha prestado mayor atención. La especie de murciélago Leptonycteris nivalis muestra un comportamiento migratorio en parte de sus hembras, que está asociado al consumo del néctar de especies ecológicamente importantes en México en particular del género Agave, y a los eventos reproductivos del murciélago. Estas características los hacen susceptibles a las consecuencias del cambio climático. En este estudio, exploramos cómo la distribución y el potencial adaptativo de Leptonycteris nivalis podrían verse afectados por posibles escenarios de calentamiento global. Para ello, obtuvimos muestras en seis localidades en la distribución de L. nivalis. Mediante técnicas de secuenciación masiva, obtuvimos 13,112 SNP (polimorfismos de un solo nucleótido) filtrados. En total sólo dos se asociaron con una variable climáticas, pero encontramos una relación estadística entre dos sitios y la susceptibilidad de la heterocigosis a cambios climáticos. Los modelos de distribución futura mostraron una reducción de las áreas climaticamente estables donde este murciélago podría sobrevivir. Pensamos que nuestros resultados respaldarán la implementación de estrategias de protección para los murciélagos y los ecosistemas que habitan. Es muy importante mantener la conectividad entre las localidades y los refugios ocupados por el murciélago magueyero. Es importante mantener espacios dentro del área de mayor estabilidad termodinámica, donde se prevén las condiciones más favorables para la presencia de murciélagos. Basándonos en nuestros análisis, creemos que la especie tiene el potencial de soportar cambios en los patrones de temperatura y precipitación, y estas circunstancias pueden favorecer la conservación de la especie.

Palabras clave: cambio climático, Desierto Chihuahuense, distribución potencial genómica de la conservación, genómica de poblaciones, murciélagos nectarívoros, potencial adaptativo, SNP

© 2026 Asociación Mexicana de Mastozoología, www.mastozoologiamexicana.org

Every species across all biomes is affected by global climate change, potentially harming the persistence and functioning of each ecosystem (Meek et al. 2023). Changes in temperature driven by global warming and the greenhouse effect directly affect organisms by challenging their physiological limits or altering their phenology or reproductivity cycles (Hayes et al. 2009). The consequences may include an expansion of the distribution range of pathogens, increased disease severity, and the introduction of invasive species into new ecosystems (Elad and Pertot 2014; Gona and More 2022; Mora et al. 2022).

The accelerated pace of climate change may prevent many populations from responding adequately, thus increasing the rate of decline and extinction (McLaughlin et al. 2002; Cahill et al. 2013; Bestion et al. 2015). Even if populations persist in their habitats and respond to these rapid changes, they still need to maintain levels of genetic diversity that allow the expression of alleles associated with adaptive characteristics (Waldvogel et al. 2020). The ability to respond to new environmental conditions is referred to as ‘evolvability’ or ‘adaptive potential’ (Houle 1992). In the context of climate change or global warming, there is a risk of reduced evolvability, as extinction could be accelerated by an inability to respond to selection pressures (Blows and Hoffmann 2005; Davis et al. 2005).

Species can avoid extinction by shifting their geographic distribution or moving to more favorable habitats, acclimating to stressful conditions through phenotypic plasticity (the ability of a genotype to express different phenotypes in different environments), or by adapting through genetic changes (Thomas 2010). However, understanding how and under what circumstances eco-evolutionary responses will occur, and differentiating among these responses to identify the adaptive potential in species with low adaptive capacity, is challenging (Urban et al. 2024).

One factor that promotes local adaptation is maintaining a species’ genetic diversity. Generally, greater genetic diversity reflects a greater adaptive potential of species to respond to extreme situations. This is why it is important to maintain genetic variants that can cope with different or changing ecological conditions (Whitlock 2015). Then, a relationship may also be established between climate change and genetic diversity and adaptive potential of natural populations (Pauls et al. 2013; Wanjala et al. 2023).

Nowadays, vulnerability to climate change is commonly assessed using forecasts generated by suitability models that project probable future scenarios (Levinsky et al. 2007; Rebelo et al. 2010; Deb et al. 2020; McGowan et al. 2021; Moura et al. 2023). Among the various genetic and epigenetic tools that enable the study of local adaptations, we used mass sequencing for SNP genotyping that allowed us to explore for potential local adaptations in the genome’s coding regions that could produce advantageous forms in response to changing environmental conditions, such as climate change, the greenhouse effect and global warming (Franks and Hoffmann 2012; Aguirre-Liguori et al. 2021; Hoffmann et al. 2021; Lancaster et al. 2022). This approach has recently become more common due to its importance in conservation biology, as it enables us to anticipate possible scenarios in the context of current climate change (Hayes et al. 2009; Kumar et al. 2016; Jordan et al. 2017; Luo et al. 2024; Klichowska et al. 2025).

Local adaptations are generally the result of divergent selection on one or more traits. These adaptations play a crucial role in determining populations’ ability to respond to environmental changes (Henriques et al. 2018; Cummins et al. 2019; Meek et al. 2023). However, this approach has not yet been applied to bats. Detecting candidate alleles under selection enables the use of statistical methods to associate their presence and frequency in natural populations with specific environmental variables (Beji et al. 2020; Guillaume et al. 2024), thus providing insights about the populations’ capabilities to face abrupt changes in environmental conditions (Hansen et al. 2012; Marková et al. 2023).

In this sense, maintaining the genetic diversity of populations is important because it provides the basis for adaptive responses to surrounding environments (Whitlock 2014; Kardos et al. 2021). Understanding the mechanisms behind the formation and maintenance of diversity and adaptive potential is essential for improving conservation and species management plans and programs.

Climate-related local extinctions have been recorded in hundreds of species (Wiens 2016). Wild mammals exhibit different degrees of vulnerability to climate change, some species have very specialized characteristics, which may cause them to disappear due to temperature changes that affect physiological processes, habitat loss, or resource loss (Boutin and Lane 2014; Hetem et al. 2014; McCain and King 2014; Chattopadhyay et al. 2019; Wells et al. 2022; de Castro et al. 2024).

The Chiroptera order is especially susceptible to ecological disturbances because most species have adapted to specific habitats and feeding habits, and most of them are susceptible to population decrease with an associated risk of extinction (Sherwin et al. 2012; Chattopadhyay et al. 2019; Gonçalves et al. 2021; McGowan et al. 2021; Festa et al. 2023). In Mexico, there are about 138 registered bat species, but only 21 are nectar-feeding (Medellín et al. 2008). Nectarivorous bats are responsible for pollinating many plant species that are considered ecologically important in the country’s most representative ecosystems, like the scrubland and low deciduous forest (Koleff et al. 2018). For instance, Leptonycteris nivalis, the long-nosed bat, murciélago magueyero mayor, or tequila bat, is one of the main pollinators of Agave plants (Hensley and Wilkins 1988). This chiropteran is distributed across central and northern Mexico, and in southern New Mexico, United States (Hensley and Wilkins 1988). It is the largest nectar-feeding bat in North America, and it is grouped into populations ranging from a few dozen to over 10,000 individuals (Hensley and Wilkins 1988). Besides many species of Agave, Leptonycteris nivalis pollinates several plant species of the Leguminosae and Cactaceae families, plants whose presence is found throughout Mexico and the southern United States in practically all ecosystems (Sánchez and Medellín 2007). In addition to its specialized nectar-consumption behavior, L. nivalis has another distinguishing characteristic: some of its females perform migratory movements from Central Mexico maternity shelters in the northern part of their distribution range, in Nuevo León, Mexico, and Texas, USA (Moreno-Valdez et al. 2000). The migration is directly related to their reproductive behavior during the spring and summer, coupling with the blooming season of various Agave species. This guarantees the availability of resources to feed during gestation and offspring rearing at the northern part of their distribution range (Bogan et al. 2017; Burke et al. 2019). However, these characteristics underscore the long-nosed bat’s vulne-rability to extreme temperature variations and habitat changes and degradation (Sánchez and Medellín 2007).

In Mexico, L. nivalis is classified as endangered in the official list of species, known as NOM-059-SEMARNAT-2010, which includes species under some form of wildlife protection, based on the Ley General de Vida Silvestre, Conservación y Aprovechamiento Sustentable (D.O.F. 2014). It has also been listed as endangered under the U.S. Endangered Species Act since 1988, and the International Union for Conservation of Nature (IUCN; Medellín 2016).

Leptonycteris nivalis evolutionary genetics and migration were recently studied using mitochondrial genetic markers (Trejo-Salazar et al. 2025). It was estimated that the species has undergone recent population growth, associated with an increase in global temperature following a decrease during the Pleistocene (Trejo-Salazar et al. 2023; Trejo-Salazar et al. 2025). This increase was supported using distribution models from the Last Interglacial, Last Maximum Glacial, Holocene, and Current epochs, in which the distribution extends to the north (Trejo-Salazar et al. 2025). Here, we describe a new set of 13,112 filtered SNPs (single-nucleotide polymorphisms) obtained using massive sequencing techniques from six locations throughout the distribution of L. nivalis. In addition, we generated models for the future potential distribution of L. nivalis under global warming. Thus, our main objective was to associate the genetic diversity of L. nivalis with changes in climatic conditions to support future management of conservation programs and maintain pollination services in Mexico’s most representative ecosystems. We expect that our results will enable the potential impacts of climate change on the tequila bat species, L. nivalis.

Materials and methods

Sampling. Forty-seven samples were taken from six localities in roosting caves and mist nests on field feeding-sites along L. nivalis distribution from 2014 to 2016 (Figure 1, Table 1), with sampling permit Secretaría del Medio Ambiente y Recursos Naturales (SEMARNAT) SGPA/DGVS/07161/15 (Supplemental File SD1), following the Animal Care and Use protocols of the American Society of Mammalogists (Sikes et al. 2016). Tissue samples were taken with a 3 mm2 biopsy wing punch in an area of the wing with no blood capillaries or nerve terminals. Wing biopsies were fixed in 90% ethanol at environmental temperature and then stored at -20°C until DNA extraction.

DNA Extraction and Genotyping. Total genomic DNA was extracted following Gasca-Pineda et al. (2015). Genomic DNA was genotyped by nextRAD sequencing libraries (SNPsaurus LLC, Eugene, OR, USA) as described by Russello et al. (2015). DNA was fragmented using Nextera reagent (Illumina, San Diego, CA, USA), which also ligates short adapter sequences to the ends of the fragments. The DNA was amplified for 27 cycles at 74°C, and libraries were sequenced on an Illumina HiSeq 4000 150 bp single‐end line (University of Oregon, Eugene, OR, USA).

Raw data processing. Raw reads were trimmed by quality with trimmomatic v0.35 (Bolger et al. 2014), using the following parameters: LEADING:5 TRAILING:5 SLIDINGWINDOW:4:20 MINLEN:75. To remove Nextera Adapters, we used the bbduk module of the BBTools v37.87 software (https://sourceforge.net/projects/bbmap/files/). Filtered reads quality was evaluated using FASTQC v0.11.5 (Andrews 2010). Single Nucleotide Polymorphisms (SNPs) were obtained in ipyrad v0.979 (Eaton and Overcast 2020) using as a reference the previously published genome of Leptonycteris yerbabuenae (Gutierrez-Guerrero et al. 2020). We applied the following parameters: a minimum depth of six for statistical base calling, two maximum alleles per site in consensus sequences, 0.05 maximum uncalled bases in consensus, 0.05 maximum heterozygotes in consensus, four as the minimum number of samples per locus, 0.2 as the maximum SNPs per locus, and 0.5 as the maximum shared heterozygous sites per locus. After testing multiple filtering schemes, we implemented the following filtering parameters using VCFtools v 0.1.16 software (Danecek et al. 2011): a maximum missing data per locus of 20%, a Hardy-Weinberg test with p ≤ 0.01, and a maximum alleles per locus of two. Because outlier methods may be sensitive to linkage disequilibrium and rare variants (Frichot and François 2015; Luu et al. 2017), we implemented a more stringent filtering scheme, including a MAC filter of 3 to recover SNPs present in at least 2 individuals and a thinning of 450 pb. This filtering recovered 13,112 SNPs.

Genetic Diversity and Population Genetic Structure. We estimated the basic summary statistics, He, Ho, and FIS, using the adegenet 2.1.11 (Jombart 2008; Jombart and Ahmed 2011) and hierfstat 0.5-11 (Goudet and Jombart 2022) R 4.5.1 packages (R Core Team 2025). To detect possible genetic structuring, we implemented Admixture v1.3.0 (Alexander et al. 2009) to assess the number of genetic groups (K) in the collected samples. We tested for 1-20 K groups, and the optimal K value was determined by the lowest 50-fold CV error (Alexander and Lange 2011). As the genetic structure of this species is unknown, we opted to use a superior K value that accounts for the possible genetic groups represented in the sample. We also performed a principal component analysis (PCA) using the R package adegenet to summarize genetic similarities among individuals.

To detect genetic structuring between sampling locations, we performed a paired FST analysis. The test was performed using a 10,000 bootstrap test; if 95% of the intervals did not reach zero, we considered the value significant. As there was only one individual in the Puebla and Zacatecas populations, these sites were removed from the analyses to obtain summary statistics and pairwise FST. In turn, we used Nei’s distances to represent the genetic differences between sampling locations using the Ward clustering method. An isolation-by-distance (IBD) analysis was performed to assess whether geographic distance between sampling locations influences genetic differentiation in Leptonycteris nivalis. This analysis was based on the correlation between genetic and geographic distances between pairs of sites, using only the set of previously identified neutral markers (SNPs) to avoid bias from potential natural selection. Nei’s distance, calculated from the allele frequencies at each location, was used as a measure of genetic distance. To strengthen the analysis at the individual level, we also used Euclidean distances between individuals derived from the first two principal components (PCs) of a principal component analysis (PCA). Geographic distance was calculated as the straight-line Euclidean distance (in kilometers) between the sampling sites’ geographic coordinates. The correlation between the genetic and geographic distance matrices was evaluated using a Mantel test with 10,000 permutations to determine the p-value. A p-value of less than 0.05 was considered to be evidence of a correlation, indicating that genetic differentiation increases with geographic distance (see below).

Additionally, we used the kinship analysis in VCFtools to evaluate kinship relationships within the study sample. Specifically, the ‘--relatedness2’ command was used to calculate a matrix of kinship estimates for all pairs of individuals from the genotype data contained in the source VCF file. This estimator is based on the correlation of allele states across multiple genetic markers between individuals, providing a robust measure that is relatively insensitive to population structure. VCFtools generates a tabulated file (.relatedness2) as output for this analysis, containing the estimated kinship values (commonly denoted as RELATEDNESS_PHI) for each pair of individuals. These values quantify the proportion of genetic material two individuals are expected to share by descent. This allows family relationships to be identified, such as duplicates/sampling errors, full siblings, and first- and second-degree relatives.

Suitability models and Future Scenarios. We obtained the environmental data from the WorldClim database (Fick and Hijmans 2017). To avoid redundancy among variables, we applied a Variance Inflation Factor (VIF) test and retained variables with a VIF below 10 using the AlleleShift R package (Kindt 2020). The final data set included: Bio.2 (Mean Diurnal Temperature Range), Bio.3 (Isothermality), Bio.8 (Mean Temperature of Wettest Quarter), Bio.9 (Mean Temperature of Driest Quarter), Bio.14 (Precipitation of Driest Month), Bio.15 (Precipitation Seasonality), Bio.16 (Precipitation of Wettest Quarter), Bio.18 (Precipitation of Warmest Quarter) and Bio.19 (Precipitation of Coldest Quarter). Then, we used the R package FactoMineR (Le et al. 2008) to implement a PCA analysis to evaluate the relative contribution of each variable.

Distribution models were constructed using a set of nine bioclimatic variables. For future projections, we used MIROC6 from the Coupled Model Intercomparison Project Phase 6 (CMIP6) for the intervals 2041-2060, 2061-2080, and 2081-2100, available from the WorldClim database (www.worldclim.org). To evaluate different climate change scenarios, we included the Shared Socio-economic Pathways models ssp126 and 585. The first is a more conservative model that predicts a lower increase in global temperatures, while the second predicts the highest increase (Fick and Hijmans 2017). Models and future projections were generated using the R library biomod2 v 4.2-6-2 (Thuiller et al. 2025). The ensemble models were created using the committee-averaging criteria (Araújo and New 2007; Qiao et al. 2015). Models were computed using the generalized boosted model (GBM; Friedman 2001) and random forest (RF; Breiman 2001). For each computed algorithm, five independent pseudo-absence sets of 1000 points were generated with two replicates, using 70% of the records to train the model and 30% to evaluate performance. We assessed model performance using the area under the receiver operating characteristic curve (AUC; Swets 1988) and true skill statistic (TSS; Allouche et al. 2006). The ensemble computation was performed using the committee average criterion restricted to models with AUC > 0.9 and TSS > 0.8. Ensembles were transformed to presence/absence values using the TSS metric (in our runs, ROC and TSS performed almost equally; we selected TSS binary because it yielded more conservative models).

Outlier detection and GenomeEnvironment Association analysis. We implemented three different approaches to evaluate loci under selection. First, we used the R package pcadapt 4.4.1(Luu et al. 2017; Privé et al. 2020). This method identifies outlier loci using principal component analysis (PCA), computing K vectors of z-scores to infer the statistical association of each SNP with the K principal components. To select the optimal value of K, we used the “screeplot” option, followed by inspection of the paired plots of the recovered components to ensure they reflected sample genetic grouping. Outlier loci were selected using the Bonferroni correction for multiple comparisons, considering an alpha threshold of ≤ 0.1.

Second, we implemented a GenomeEnvironment Association (GEA) method using Redundancy Analysis (RDA) to identify outlier SNPs associated with environmental variables, following Capblancq and Forester (2021). Initially, we used the same variable data set as in the suitability models; however, to avoid overfitting, we excluded variables identified as collinear by the RDA analysis. The RDA was performed using the R package vegan (Oksanen et al. 2025) using SNPs as response variables and the three selected environmental variables as explanatory variables (Bio.3, Bio.9, Bio.14). RDA requires a genotype matrix without missing data; therefore, missing values were imputed using the frequency of the most common allele (the mode). Outlier SNPs were identified on each of the first three constrained axes using a threshold of 2 standard deviations. Then, we evaluated the association between each outlier and the environmental variables, estimating the Pearson correlation coefficient. We retained SNPs with correlation values of 0.5 or greater.

Third, we used a latent factor mixed model (LFMM; Frichot et al. 2013) implemented in the LEA v3.20.0 R package (Frichot and François 2015; Gain and François 2021). This method performs statistical tests to identify loci associated with environmental gradients (Frichot et al. 2013). The analysis was run using the environmental variables previously selected for each analysis. Outliers were selected using the Bonferroni correction with a threshold of ≤ 0.1.

We opted for these methods because they do not require an a priori definition of genetic groups, unlike FST-based methods. Considering all loci detected by the three methods, we generated a neutral data set for the genetic structure analysis. To annotate outlier loci, we considered only loci shared by at least two methods. Then, loci were mapped to the reference genome and extracted 500 bp in both directions (a total of 1,000 bp). We also performed a blast search (Altschul et al. 1997) against the UniProtKB database (https://www.uniprot.org/) and only hits with an e-value < 1x10-5 were considered.

We calculated the risk of non-adaptation (RONA) using pyRONA v0.4.3 (Rellstab et al. 2016; Pina-Martins et al. 2019). This method is based on the association between current allele frequencies and current climatic variables, and the estimation of the differences between current and future expected allele frequencies. For this analysis, we used the outlier loci detected by the three methods, all based on the same data set of environmental variables used in the LEA and RDA outlier detection methods. A high RONA value suggests a greater likelihood of maladaptation to climate change. For this analysis, we implemented the MIROC6 Coupled Model Intercomparison Project Phase 6 ssp585 for the interval 2081-2100.

Results

Genetic Diversity and Population Structure. We obtained an overall of 2,776,591 raw SNPs. After filtering, we identified 13,112 SNPs distributed across 2,065 scaffolds in the Leptonycteris yerbabuenae reference genome. The observed heterozygosity for the overall sample was Ho = 0.1866 (s.d. = 0.259), whereas the average expected heterozygosity was He = 0.2112 (s.d. = 0.1815) and the FIS was 0.1097 (s.d. = 0.3101). The genetic-distance dendrogram revealed a genetic grouping irrespective of the geographical distance between sampling locations (Figure 2). Furthermore, the analysis of FST between the different locations showed similar patterns (Figure 3), with low but > 0.05 values (except for Michoacan vs. Hidalgo; Table SD2, supplementary materials). The principal component analysis revealed a homogeneous population, consistent with the FST-based analyses (Figure 4).

Relatedness analysis. Our analysis of relatedness among all individuals revealed minimal kin relationships among two pairs of individuals in the same locality (Figure 5). Only individuals ln70-ln73, and individuals ln64-ln69, who are all from Nuevo León, and the degree of relationship in the two pairs of individuals reflects a close kinship similar to that between uncles and nephews. The first pair exhibits a higher degree of relatedness between its members (Figure 5).

To complete these results, the Admixture analysis showed no genetic structuring (K = 1, Figure SD2 supplemental materials). Moreover, worth noting that the plots with higher K values showed no admixture between individuals (SD3 supplemental material). The distance-based isolation analysis for the samples (SD4 supplemental material) did not show a correlation between geographic distance and genetic distance (p > 0.3).

Potential Distribution Model in Future Scenarios. Our models projecting the potential distribution of the migratory bat, L. nivalis, in future climate change scenarios show contraction in areas of greater climate suitability. In the most optimistic scenario (Figure 6a, SSP1-2.6), this contraction is minimal and concentrated in the northern part of L. nivalis’ distribution, where temperatures would not rise by more than the global average increase of approximately 1.8 °C by the year 2100 compared to pre-industrial levels.

By contrast, in the most pessimistic scenario (SSB585; Figure 6b), which projects a temperature increase of 4.4 °C (with a probable range of 3.3 to 5.7 °C) above pre-industrial levels by 2100, the contraction of climatic stable areas mainly occurs in the north of the distribution range, particularly affecting areas where maternity caves are located, primarily in Nuevo León in Mexico and Texas in the United States. The area that showed the greatest climatic suitability across the three time periods is in the center and southern parts of the current distribution. This means that even in the scenario of the greatest global warming, the center of distribution and part of the south would remain suitable for the species.

Outlier detection and SNP annotation. The three methods yielded an overall of 582 outlier SNPs. The PCAdapt analysis recovered 123 outlier SNPs, RDA 254 SNPs, and LEA 213; from them, only two SNPs were shared by at least two methods. (Figure 7) Of the outlier loci, only one had a match with uncharacterized proteins from species of the Vespertilionidae family.

Risk of non-adaptedness. RONA analysis shows an association between heterozygosity and changes in three climate variables in the context of projected global change (Figure 8). In other words, current heterozygosity maintains adaptive potential under present conditions; however, changing some of these variables may compromise this potential in response to environmental conditions. The most significant risk relationship in the context of climate change is that between heterozygosity and the Bio.9 variable, ‘Mean Temperature of Driest Quarter’. However, this risk is more pronounced in the northern part of the species’ distribution, specifically in Texas (0.93415) in the United States and, to a lesser extent, in Nuevo León, Mexico (0.5362). The Bio.3 variable, ‘Isothermality’, and the Bio.14 variable, ‘Precipitation of Driest Month’, showed average RONA values at sites in the center and south of the long-nosed bat’s distribution range. For Bio.3, the three sites at greatest risk were Texas (0.45894), Zacatecas (0.44286) and Puebla (0.3934). Meanwhile, the highest values for Bio.14 were found in Puebla (0.55603) in the south and in Texas (0.38008) in the north (see all values in Supplemental Data SD5).

Discussion

The genetic diversity, expressed as expected heterozygosity (Ho), is considered a parameter that can help us predict the adaptive potential of species in the face of rapidly emerging adverse conditions (Schmidt and Russello 2025), as is currently being observed during the Anthropocene. While this parameter does not precisely identify the genes and alleles under local selection and adaptation, it provides an overview of the populations’ adaptive potential. We initially assumed that L. nivalis, if it had low levels of heterozygosity, would be vulnerable to the abrupt changes expected in the near future, in particular if their distribution in maternity caves found in the northern part of the species’ range-- primarily in Nuevo León in Mexico and Texas in the United States during spring-summer seasons-- would be affected by the climatic changes. For this reason, it is important to explore the optimal conditions for the species’ migratory and reproductive dynamics.

Maintaining heterozygosity is therefore extremely important in L. nivalis, given that its genetic diversity is low (Ho = 0.1866) compared to that of other bat species considered for national and/or international protection, such as the sister species, the lesser long-nosed bat, L. yerabuenae (Ho= 0.292; Peralta et al. 2025), even with a larger sample analyzed in our study. Heterozygosity of L. nivalis is also lower than in Myotis lucifugus (Ho= 0.247-0.263; Lilley et al. 2020) Myotis grisescens (Ho= 0.183; Nagel et al. 2023). Myotis sodalis (Ho= 0.083) (Nagel et al. 2023), as well as Miniopterus schreibersii (Dufresnes et al. 2023), and Myotis septentrionalis (Ho = 0.000086 -0.000115; Grimshaw et al. 2024).

Given the geographic distribution of the species, we consider our sampling to be adequate of a genomic study (were the large number of genetic makers can compensate the smaller samples of individuals, see simulations of Aguirre-Liguori et al. 2020), although we recognize that there is a bias towards northern locations, Texas and Nuevo León, as well as shortage in one locality at the center and one locality at southern of L. nivalis distribution. However, bias and a small sample size may have prevented the detection of genetic differences between populations or genetic structure.

Based on previous work, we expected to detect population genetic structure, given that evidence of distinct genetic groups had been reported in a recent study using two mitochondrial markers and one associated-chromosome Y gene data (Trejo-Salazar et al. 2025). This is important for associating outlier alleles with local adaptations and for proposing conservation strategies that prioritize sites where individuals carry these alleles. However, our results did not detect any genetic structure. Genetic structure analyses in this study and previous studies (Pourshoushtari and Ammerman 2021; Trejo-Salazar et al. 2025) suggest that genetic flow remains strong and constant, preventing genetic differentiation between refugees and locations inhabited by the species. Therefore, it appears that females maintain this genetic flow when they return to the south-central part of the species’ distribution during the autumn-winter season, keeping the sites or populations as a single unit.

Despite the absence of genetic population structure and the relatively large sample size, we detected only two candidate alleles under selection that are statistically correlated with the climatic variables used for distribution model projections. But they are not clearly representatives of maladaptations or of any potential adaptation related to climate change. In other words, we identified alleles related to climatic conditions, temperature, and precipitation, we did not find local adaptations that would allow us to prioritize the protection of one area over another. We also detected a positive correlation between genetic diversity and four climatic variables. This result is relevant, as it enables us to predict the future of tequila bat populations and prioritize connectivity between sites, migratory routes, and the conservation of heterozygosity levels for the species.

Thus, environmental pressures or resource scarcity can affect the entire population of a species, pushing it into a complex situation. Based on our results, Nuevo León is among the most susceptible locations to climate change. According to RONA’s analysis, there is a closer relationship between diversity and the precipitation-related climate variable (Bio9) in this location. In Nuevo Léon, at the El Infierno cave, both sexes gather to form shelters for mating and raising offspring. The long-nosed bat is present there alongside at least five species of the genus Agave during their flowering stage. Furthermore, the presence of L. nivalis in this area is positively related to ambient temperature (Moreno-Valdez et al. 2004), so any variation in temperature could negatively affect the species’ ecological dynamics and, consequently, its primary food source.

On the other hand, we think that the lack of population structure helps to explain the distant kinship relationships between individuals in our sample. In other words, the absence of population structure means that it is less likely that close relatives will be found in small population samples (Weir and Goudet 2017). A larger sample size would likely yield more detailed information on both genetic population structure and kinship relationships, and apparently, genetic flow between refuges has remained constant. However, gene flow must be maintained to ensure that the adaptive potential and the outlier loci remain part of the species’ gene pool (Chhina et al. 2024). Nevertheless, intense selective pressures, such as climate change or global warming, could cause a bottleneck (Jangjoo et al. 2016) and affect the parameters addressed in this study, thereby impacting genetic diversity and the adaptive potential stored in this genomic variation.

We consider our analysis for detecting outlier loci under selection to be representative, as, in general terms, this sampling reflects the species’ distribution. Nevertheless, among the 13,112 SNPs used for outlier detection, fewer than 5% were identified as outliers by any of the three methods, and only two were shared across methods. We considered that our results are similar to those found in previous studies in mammals. For instance, for the sister species, L. yerbabuenae, Peralta et al. (2025) reported 17,732 SNPs after filtering, of which 7,172 were useful for demographic analyses with 32 individuals, compared with our sample of 47 individuals. In 288 individual samples of the American pika (Ochotona princeps), Henry and Russello (2013) detected 68 outlier loci under selection, and 15 of these were statistically associated with at least one climatic variable. In the arboreal Australian mammal Petauroides volans, Knipler et al. (2023) detected 66 loci associated with climatic variables using an 85-individual sample.

As for studies focused on the Chiroptera group, despite many species being considered vulnerable due to ecological characteristics, morphological features, and physiological specializations, there is insufficient information to predict likely scenarios using genomic tools and to establish preventive or corrective measures in response to climate change. One of the few studies to address the adaptive potential of the European bat, Plecotus austriacus, was Razgour et al. (2018), which examined a sample of 83 individuals and detected 24 single-nucleotide polymorphisms (SNPs) associated with extreme temperature and summer precipitation variables. In contrast, we detected statistically significant associations with just two SNPs and climatic variables: temperature during the driest month and precipitation of the warmest quarter. However, these were merely statistical correlations, and we did not find a biological association.

Regarding the potential distribution of the species under global warming scenarios, the most pessimistic scenario predicts a global temperature increase of around 4.4°C relative to pre-industrial levels. This would lead to a contraction of L. nivalis current distribution, mainly affecting migratory areas during the breeding season (Moreno-Valdez et al. 2000). Nevertheless, this species has already experienced similar dynamics in the past (Trejo-Salazar et al. 2025), during the climate changes from the Last Interglacial to the Holocene and Current time. Thus, we consider that L. nivalis populations can adjust to seasonal changes in climatic conditions through migratory movements. However, these changes have occurred in the past over longer periods than those we are expecting in the Anthropocene, so we cannot be sure that species will be able to withstand such predicted radical changes in less than 100 years.

Despite their obvious importance, relatively few studies have incorporated local adaptation into predictions of how species’ distributions and abundances will be affected by climate change (Smith et al. 2019) or have taken local adaptation into account when designing conservation and recovery plans (Peterson et al. 2019). Managed (assisted) gene flow and the incorporation of local adaptation into conservation plans can provide options for populations that have become increasingly maladapted to their environment, due to reduced genetic diversity or rapidly changing environments, or for restoring resilient populations in places where they have become locally extinct.

Conclusions

The future distribution of the long-nosed bat may be compromised in areas that are crucial for its reproductive behavior, particularly in places where it establishes maternity roosts during the spring and summer seasons. By inferring and predicting the distribution patterns and consequences for the adaptive potential of this nectarivorous bat, we can suggest possible scenarios of particular relevance, considering, for instance, that some plant species depend on pollination by these mammals. In other words, this way we can predict possible consequences at the ecosystem level due to climate change in Mexico.

It is possible that L. nivalis will adapt to these future climate scenarios. However, it must still be considered that current global warming is very rapid, and although bats of the genus Leptonycteris have survived and adapted to similar -but slower- climatic changes of the planet, it is not possible to ensure that in such a short time the current populations will manage to maintain a large enough size.

Regarding climate change scenarios, two predictions were made. The one we consider most optimistic forecasts an increase of approximately 1.8°C by the end of this century. In that case, the changes in climatic conditions favorable to bat presence are minor, and therefore, we believe the consequences may be barely noticeable. However, in the most pessimistic scenario, with a warming of around 4.4°C above pre-industrial global temperatures, the situation is more catastrophic, both for bats and, as will surely be the case, for the plants that depend on pollination by these bats.

Acknowledgments

This work is dedicated to Dra. Livia León Paniagua in recognition of her contributions to Mexican mammalogy. She played a fundamental role in training human resources and studying Mexican bats and rodents. We would like to express our gratitude for her hard work curating one of the country’s largest zoological collections for over 40 years. On a personal level, Dra. Livia supported Roberto Trejo’s academic, professional and personal development, always offering words of encouragement and demonstrating her humanity. Without her presence and support, we would not have been able to finish this work or other previous publications. We admire her as a person and an academic because she has always supported and collaborated with those in need, promoting an atmosphere of respect and professionalism with her positive and often maternal attitude. Author thanks Dra. Erika Aguirre Planter for technical support in laboratory, Dra. Edén Rodríguez and Mc S. Dulce Hernández for supporting DNA extraction. Gasca-Pineda J. received a post-doctoral fellowship from “Estancias Posdoctorales por México para la Formación y Consolidación de las y los Investigadores por México” (SECIHTI , formerly CONAHCYT, CVU: 172925); Trejo-Salazar RE received a post-doctoral fellowship from “Estancias Posdoctorales por México para la Formación y Consolidación de las y los Investigadores por México” (SECIHTI , formerly CONAHCYT, CVU: 249584). Proyecto Conacyt Investigación en Fronteras de la Ciencia 2015, número 177, Genómica de la Diversidad de Vertebrados Mexicanos 1: Leptonycteris y la evolución de la nectarivoría en murciélagos y aves.

Declaration of Artificial Intelligence use

The authors declare that we have not used any artificial intelligence tools at any stage of the preparation of this manuscript.

Author contributions

Roberto-Emiliano Trejo-Salazar designed the project, collect samples, database construction, analyses and drafting the manuscript; Jaime Gasca-Pineda participated in the database construction, computational analyses and drafting manuscript; Rosalinda Tapia López contributed with DNA extraction and laboratory support and drafting manuscript; Luis E. Eguiarte helped with the logistics and drafted and corrected the manuscript.

Supplementary data

Supplementary materials are available at Therya online and https://doi.org/10.5281/zenodo.17613255

SD1. Permit to collect biological samples.

SD2. Pairwise FST values between the sampled locations and Bootstrap values for paired comparisons of FST between sampled locations.

SD3. CV values for different K values in the admixture analysis.

SD4. Distance isolation analysis.

SD5. RONA analysis values for each location.

Literature cited

Alexander DH, Novembre J, and Lange, K. 2009. Fast model-based estimation of ancestry in unrelated individuals. Genome Research 19:1655-1664. https://doi.org/10.1101/gr.094052.109

Alexander DH, and Lange K. 2011. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinformatics 12:246. https://doi.org/10.1186/1471-2105-12-246

Aguirre-Liguori JA, Luna-Sánchez JA, Gasca-Pineda J, and Eguiarte LE. 2020. Evaluation of the minimum sampling design for population genomic and microsatellite studies: An analysis based on wild maize. Frontiers in Genetics 11:870. https://doi.org/10.3389/fgene.2020.00870

Aguirre-Liguori JA, Ramírez-Barahona S, and Gaut BS. 2021. The evolutionary genomics of species’ responses to climate change. Nature Ecology & Evolution 5:1350-1360. https://doi.org/10.1038/s41559-021-01526-9

Allouche O, Tsoar A, and Kadmon R. 2006. Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). Journal of Applied Ecology 43:1223–1232. https://doi.org/10.1111/j.1365-2664.2006.01214.x

Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, et al. 1997. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Resources 25:3389–3402. https://doi.org/10.1093/nar/25.17.3389

Andrews S. 2010. FastQC: a Quality Control Tool for High Throughput Sequence Data. [Accessed October 10th, 2024]. http://www.bioinformatics.babraham.ac.uk/projects/fastqc

Araújo M, and New M. 2007. Ensemble forecasting of species distributions. Trends in Ecology & Evolution 22:42–47. https://doi.org/10.1016/j.tree.2006.09.010

Beji S, Fontaine V, Devaux R, Thomas, M, Negro, SS, Bahrman, N, et al. 2020. Genome-wide association study identifies favorable SNP alleles and candidate genes for frost tolerance in pea. BMC Genomics 21:536. https://doi.org/10.1186/s12864-020-06928-w

Bestion E, Teyssier A, Richard M, Clobert J, and Cote J. 2015. Live fast, die young: experimental evidence of population extinction risk due to climate change. PLoS Biology 13:e1002281. https://doi.org/10.1371/journal.pbio.1002281

Blows MW, and Hoffmann AA. 2005. A reassessment of genetic limits to evolutionary change. Ecology 86:1371–84. https://doi.org/10.1890/04-1209

Bolger AM, Lohse M, and Usadel B. 2014. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30: 2114–2120. https://doi.org/10.1093/bioinformatics/btu170

Boutin S, and Lane JE. 2014. Climate change and mammals: evolutionary versus plastic responses. Evolutionary Applications 7:29–41. https://doi.org/10.1111/eva.12121

Bogan MA, Cryan PM, Weise CD, and Valdez EW. 2017. Landscape movements by two species of migratory nectar-feeding bats (Leptonycteris) in a northern area of seasonal sympatry. Western North American Naturalist 77:317–330. https://doi.org/10.3398/064.077.0305

Burke RA, Frey JK, Ganguli A, and Stoner KE. 2019. Species distribution modelling supports “nectar corridor” hypothesis for migratory nectarivorous bats and conser-vation of tropical dry forest. Diversity and Distributions 25:1399–1415. https://doi.org/10.1111/ddi.12950

Breiman L. 2001. Random forests. Machine Learning 45:5–32. https://doi.org/10.1023/A:1010933404324

Cahill AE, Aiello-Lammens ME, Fisher-Reid MC, Hua X, Karanewsky CJ, Yeong R, et al. 2013. How does climate change cause extinction? Proceedings of the Royal Society B: Biological Sciences 280(1750): 20121890. https://doi.org/10.1098/rspb.2012.1890

Capblancq T, and Forester BR. 2021. Redundancy analysis: A Swiss Army knife for landscape genomics. Methods in Ecology and Evolution 12:2298–2309. https://doi.org/10.1111/2041-210X.13722

Chattopadhyay B, Garg KM, Ray R, and Rheindt FE. 2019. Fluctuating fortunes: genomes and habitat reconstructions reveal global climate-mediated changes in bats’ genetic diversity. Proceedings of the Royal Society B 286(1911):20190304. https://doi.org/10.1098/rspb.2019.0304

Chhina AK, Abhari N, Mooers A, and Lewthwaite JMM. 2024. Linking the spatial and genomic structure of adaptive potential for conservation management: a review. Genome 67:403–423. https://doi.org/10.1139/gen-2024-0036

Cummins D, Kennington WJ, Rudin‐Bitterli T, and Mitchell NJ. 2019. A genome‐wide search for local adaptation in a terrestrial‐breeding frog reveals vulnerability to climate change. Global Change Biology 25:3151–3162. https://doi.org/10.1111/gcb.14703

Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. 2011. The variant call format and VCFtools. Bioinformatics 27:2156–2158.

Davis MB, Shaw RG, and Etterson JR. 2005. Evolutionary responses to changing climate. Ecology 86:1704–1714. https://doi.org/10.1890/03-0788

Deb JC, Forbes G, and MacLean DA. 2020. Modelling the spatial distribution of selected North American woodland mammals under future climate scenarios. Mammal Review 50:440–452. https://doi.org/10.1111/mam.12210

de Castro Evaldt BH, Leite YLR, and Loss AC. 2024. Climate change impact on small mammals from two Neotropical hotspots. Biological Journal of the Linnean Society 143(3): blae014. https://doi.org/10.1093/biolinnean/blae014

Diario Oficial de la Federación. 05-03-2014. Norma Oficial Mexicana – 059 (NOM-059). Acuerdo por el que se da a conocer la lista de especies y poblaciones prioritarias para la conservación. Estados Unidos Mexicanos. Mexico City (MEX): Secretaría de Medio Ambiente y Recursos Naturales.

Dufresnes C, Dutoit L, Brelsford A, Goldstein-Witsenburg F, Clément L, López-Baucells, et al. 2023. Inferring genetic structure when there is little: population genetics versus genomics of the threatened bat Miniopterus schreibersii across Europe. Scientific Reports 13:1523. https://doi.org/10.1038/s41598-023-27988-4

Eaton DA, and Overcast I. 2020. ipyrad: Interactive assembly and analysis of RADseq datasets. Bioinformatics 36:2592–2594. https://doi.org/10.1093/bioinformatics/btz966

Elad Y, and Pertot I. 2014. Climate change impacts on plant pathogens and plant diseases. Journal of Crop Improvement 28:99–139. https://doi.org/10.1080/15427528.2014.865412

Festa F, Ancillotto L, Santini L, Pacifici M, Rocha R, Toshkova N, et al. 2023. Bat responses to climate change: a systematic review. Biological Reviews 98:19–33. https://doi.org/10.1111/brv.12893

Fick SE, and Hijmans RJ. 2017. WorldClim 2: new 1‐km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37:4302–4315. https://doi.org/10.1002/joc.5086

Franks SJ, and Hoffmann AA. 2012. Genetics of climate change adaptation. Annual Review of Genetics 46:185–208. https://doi.org/10.1146/annurev-genet-110711-155511

Friedman JH. 2001. Greedy function approximation: a gradient boosting machine. Annals of Statistics 29:1189–1232. https://www.jstor.org/stable/2699986

Frichot E, Schoville SD, Bouchard G, and François O. 2013. Testing for associations between loci and environmental gradients using latent factor mixed models. Molecular Biology and Evolution 30:1687–1699. https://doi.org/10.1093/molbev/mst063

Frichot E, and François O. 2015. LEA: An R package for landscape and ecological association studies. Methods in Ecology and Evolution 6:925–929. https://doi.org/10.1111/2041-210X.12382

Gain C, and François O. 2021. LEA 3: Factor models in population genetics and ecological genomics with R. Molecular Ecology Resources 21:2738–2748. https://doi.org/10.1111/1755-0998.13366

Gasca-Pineda J. 2015. Genética de la conservación del bisonte de la pradera (Bison bison bison) y del borrego cimarrón (Ovis canadensis) en México. [PhD thesis]. [Mexico City (MEX)]: Universidad Nacional Autónoma de México.

Gona PN, and More AF. 2022. Bacterial pathogens and climate change. The Lancet 400:2161–2163.

Gonçalves F, Sales LP, Galetti M, and Pires MM. 2021. Combined impacts of climate and land use change and the future restructuring of Neotropical bat biodiversity. Perspectives in Ecology and Conservation 19:454–463. https://doi.org/10.1016/j.pecon.2021.07.005

Goudet J, and Jombart T. 2022. hierfstat: Estimation and Tests of Hierarchical F-Statistics. <ttps://doi.org/10.32614/CRAN.package.hierfstat>, R package version 0.5-11, https://CRAN.R-project.org/package=hierfstat

Grimshaw JR, Donner D, Perry R, Ford WM, Silvis A, Garcia CJ, et al. 2024. Disentangling genetic diversity of Myotis septentrionalis: population structure, demographic history, and effective population size. Journal of Mammalogy 105:854–864. https://doi.org/10.1093/jmammal/gyae056

Guillaume AS, Leempoel K, Rogivue A, Gugerli F, Parisod C, and Joost S. 2024. Integrating very high resolution environmental proxies in genotype–environment association studies. Evolutionary Applications 17:e13737. https://doi.org/10.1111/eva.13737

Gutierrez-Guerrero YT, Ibarra-Laclette E, Martínez del Río C, Barrera-Redondo J, Rebollar EA, Ortega J, et al. 2020. Genomic consequences of dietary diversification and parallel evolution due to nectarivory in leaf-nosed bats. Gigascience 9:giaa059. https:// doi.org/10.1093/gigascience/giaa059

Hansen MM, Olivieri I, Waller DM, Nielsen EE, and GeM Working Group. 2012. Monitoring adaptive genetic responses to environmental change. Molecular Ecology 21:1311–1329. https://doi.org/10.1111/j.1365-294X.2011.05463.x

Hayes BJ, Bowman PJ, Chamberlain AJ, Savin K, van Tassell CP, Sonstegard TS, et al. 2009. A Validated Genome Wide Association Study to Breed Cattle Adapted to an Environment Altered by Climate Change. PLoS One 4:e6676. https://doi:10.1371/journal.pone.0006676

Henriques D, Wallberg A, Chávez-Galarza J, Johnston JS, Webster MT, and Pinto MA. 2018. Whole genome SNP-associated signatures of local adaptation in honeybees of the Iberian Peninsula. Scientific Reports 8:11145. https://doi.org/10.1038/s41598-018-29469-5

Henry P, and Russello MA. 2013. Adaptive divergence along environmental gradients in a climate‐change‐sensitive mammal. Ecology and Evolution 3:3906–3917. https://doi.org/10.1002/ece3.776

Hensley AP, and Wilkins KT. 1988. Leptonycteris nivalis. Ma-mmalian Species 307:1–4. https://doi.org/10.1644/797.1

Hetem RS, Fuller A, Maloney SK, and Mitchell D. 2014. Responses of large mammals to climate change. Temperature 1:115–127. https://doi.org/10.4161/temp.29651

Hoffmann AA, Weeks AR, and Sgrò CM. 2021. Opportunities and challenges in assessing climate change vulnerability through genomics. Cell 184:1420–1425. https://doi.org/10.1016/j.cell.2021.02.006

Jangjoo M, Matter SF, Roland J, and Keyghobadi N. 2016. Connectivity rescues genetic diversity after a demographic bottleneck in a butterfly population network. Proceedings of the National Academy of Sciences 113:10914–10919. https://doi.org/10.1073/pnas.160086511

Houle D. 1992. Comparing evolvability and variability of quantitative traits. Genetics 130:195-204. https://doi.org/10.1093/genetics/130.1.195

Jombart T. 2008. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24: 1403–1405. https://doi:10.1093/bioinformatics/btn129

Jombart T, and Ahmed I. 2011. adegenet 1.3-1: new tools for the analysis of genome-wide SNP data. Bioinformatics 27: 3070–3071. https://doi:10.1093/bioinformatics/btr521

Jordan R, Hoffmann AA, Dillon SK, and Prober SM. 2017. Evidence of genomic adaptation to climate in Eucalyptus microcarpa: Implications for adaptive potential to projected climate change. Molecular Ecology 26:6002–6020. https://doi.org/10.1111/mec.14341

Kardos M, Armstrong EE, Fitzpatrick SW, Hauser S, Hedrick PW, Miller JM, et al. 2021. The crucial role of genome-wide genetic variation in conservation. Proceedings of the National Academy of Sciences 118:e2104642118. https://doi.org/10.1073/pnas.2104642118

Kindt R. 2020. AlleleShift: an R package to predict and visualize population-level changes in allele frequencies in response to climate change. PeerJ 9:e11534. https://10.7717/peerj.11534

Knipler ML, Gracanin A, and Mikac KM. 2023. Conservation genomics of an endangered arboreal mammal following the 2019–2020 Australian megafire. Scientific Reports 13:480. https://doi.org/10.1038/s41598-023-27587-3

Klichowska E, Wróbel A, Nowak A, and Nobis M. 2025. Eco‐evolutionary genomics reveal mountain range‐specific adaptation and intraspecific variation in vulnerability to climate change of alpine endemics. Molecular Ecology 34:e70113. https://doi.org/10.1111/mec.70113

Koleff P, Urquiza-Haas T, Ruiz-González SP, Hernández-Robles DR, Mastretta-Yanes A, Quintero E, et al. 2018. Biodiversity in Mexico: State of knowledge. In: Pullaiah T, editor. Global Biodiversity Volume 4: Selected Countries in the Americas and Australia 1st Edition. New York (EEUU): Apple Academic Press; p. 285–337.

Kumar S, Banks TW, and Cloutier S. 2016. SNP discovery through next-generation sequencing and its applications. In: Kumar S. editor. Crop Breeding Bioinformatics and Preparing for Climate Change. New York (EEUU): Apple Academic Press; p. 209–244.

Lancaster LT, Fuller ZL, Berger D, Barbour MA, Jentoft S, and Wellenreuther M. 2022. Understanding climate change response in the age of genomics. Journal of Animal Ecology 91:1056–1063. https://doi.org/10.1111/1365-2656.13711

Le S, Josse J, and Husson F. 2008. FactoMineR: An R Package for Multivariate Analysis. Journal of Statistical Software 25:1-18. https://doi.org/10.18637/jss.v025.i01

Levinsky I, Skov F, Svenning JC, and Rahbek C. 2007. Potential impacts of climate change on the distributions and diversity patterns of European mammals. Biodiversity and Conservation 16:3803–3816. https://doi.org/10.1007/s10531-007-9181-7

Lilley TM, Wilson IW, Field KA, Reeder DAM, Vodzak ME, Turner GG, et al. 2020. Genome-Wide Changes in Genetic Diversity in a Population of Myotis lucifugus Affected by White-Nose Syndrome. G3 Genes|Genomes|Genetics 10:2007–2020.

Luo Y, Qin W, Yan Y, Yin K, Zang R, and Du FK. 2024. Climate change vulnerability and conservation strategies for tertiary relict tree species: Insights from landscape genomics of Taxus cuspidata. Evolutionary Applications 17(9):e13686. https://doi.org/10.1111/eva.13686

Luu K, Bazin E, and Blum MG. 2017. pcadapt: an R package to perform genome scans for selection based on principal component analysis. Molecular Ecology Resources 17:67–77. https://doi.org/10.1111/1755-0998.12592

Marková S, Lanier HC, Escalante MA, da Cruz MO, Horníková M, Konczal M, et al. 2023. Local adaptation and future climate vulnerability in a wild rodent. Nature Communications 14:7840. https://doi.org/10.1111/1755-0998.12592

McCain CM, and King SR. 2014. Body size and activity times mediate mammalian responses to climate change. Global Change Biology 20:1760–1769. https://doi.org/10.1111/gcb.12499

McGowan NE, Roche N, Aughney T, Flanagan J, Nolan P, Marnell F, et al. 2021. Testing consistency of modelled predictions of the impact of climate change on bats. Climate Change Ecology 2:100011. https://doi.org/10.1016/j.ecochg.2021.100011

McLaughlin JF, Hellmann JJ, Boggs CL, and Ehrlich PR. 2002. Climate change hastens population extinctions. Proceedings of the National Academy of Sciences 99:6070–6074. https://doi.org/10.1073/pnas.052131199

Medellín RA, Arita HT, and Sánchez O. 2008. Identificación de los murciélagos de México: Clave de campo (2ª ed.). Mexico City (MEX): Asociación Mexicana de Mastozoología, A.C., Instituto de Ecología, Universidad Nacional Autónoma de México.

Medellín R. 2016. Leptonycteris nivalis. The IUCN Red List of Threatened Species ٢٠١٦: e.T١١٦٩٧A٢٢١٢٦١٧٢. [Accessed on 30 June 2025]. https://dx.doi.org/10.2305/IUCN.UK.2016-1.RLTS.T11697A22126172.en

Meek MH, Beever EA, Barbosa S, Fitzpatrick SW, Fletcher NK, Mittan-Moreau, et al. 2023. Understanding local adaptation to prepare populations for climate change. Bioscience 73:36–47. https://doi.org/10.1093/biosci/biac101

Mora C, McKenzie T, Gaw IM, Dean JM, von Hammerstein H, Knudson TA, et al. 2022. Over half of known human pathogenic diseases can be aggravated by climate change. Nature Climate Change 12:869–875. https://doi.org/10.1038/s41558-022-01426-1

Moura MR, Oliveira GA, Paglia AP, Pires MM, and Santos BA. 2023. Climate change should drive mammal defaunation in tropical dry forests. Global Change Biology 29:6931–6944. https://doi.org/10.1111/gcb.16979

Moreno-Valdez A, Grant WE, and Honeycutt RL. 2000. A simulation model of Mexican long-nosed bat (Leptonycteris nivalis) migration. Ecological Modelling 134:117–127. https://doi.org/10.1016/S0304-3800(00)00253-2

Moreno-Valdez A, Honeycutt RL, and Grant WE. 2004. Colony dynamics of Leptonycteris nivalis (Mexican Long-Nosed Bat) related to flowering Agave in Northern Mexico. Journal of Mammalogy 85:453–459.

Nagel JJ, Nelson DM, and Gugger PF. 2023. Population genetic structure and effective size of two endangered cave bat species. Acta Chiropterologica 25:203–211. https://doi.org/10.3161/15081109ACC2023.25.2.001

Oksanen J, Simpson G, Blanchet F, Kindt R, Legendre P, Minchin P, et al. 2025. Package ‘vegan’. Community Ecology Package, version, 2:1–295

Pauls SU, Nowak C, Bálint M, and Pfenninger M. 2013. The impact of global climate change on genetic diversity within populations and species. Molecular Ecology 22:925–946. https://doi.org/10.1111/mec.12152

Peralta DM, Jaramillo-Correa JP, Hernández-Rosales HS, Túnez JI, Gasca-Pineda J, Medellín RA, et al. 2025. To Migrate or not to migrate? Exploring the genomic basis of partial migratory behavior in bats Genome Biology and Evolution 17: evaf203.

Peterson ML, Doak DF, and Morris WF. 2019. Incorporating local adaptation into forecasts of species’ distribution and abundance under climate change. Global Change Biology 25:775–793. https://doi.org/10.1111/gcb.14562

Pina-Martins F, Baptista J, Pappas GJr, and Paulo OS. 2019. New insights into adaptation and population structure of cork oak using genotyping by sequencing. Global Change Biology 25:337–350.

Pourshoushtari RD, and Ammerman LK. 2021. Genetic variability and connectivity of the Mexican long-nosed bat between two distant roosts. Journal of Mammalogy 102:204–219. https://doi.org/10.1093/jmammal/gyaa138

Privé F, Luu K, Vilhjálmsson BJ, and Blum MG. 2020. Performing highly efficient genome scans for local adaptation with R package pcadapt version 4. Molecular Biology and Evolution 37:2153–2154. https://doi.org/10.1093/molbev/msaa053

Qiao H, Soberón J, and Peterson AT. 2015. No silver bullets in correlative ecological niche modelling: insights from testing among many potential algorithms for niche
estimation. Methods in Ecology and Evolution 6:1126–1136. https://doi.org/10.1111/2041-210X.12397

R Core Team (2025). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/

Razgour O, Taggart JB, Manel S, Juste J, Ibañez C, Rebelo H, et al. 2018. An integrated framework to identify wildlife populations under threat from climate change. Molecular Ecology Resources 18:18–31. https://doi.org/10.1111/1755-0998.12694

Rebelo H, Tarroso P, and Jones G. 2010. Predicted impact of climate change on European bats in relation to their biogeographic patterns. Global Change Biology 16:561–576. https://doi.org/10.1111/j.1365-2486.2009.02021.x

Rellstab C, Zoller S, Walthert L, Lesur I, Pluess AR, Graf R, et al. 2016. Signatures of local adaptation in candidate genes of oaks (Quercus spp.) with respect to present and future climatic conditions. Molecular Ecology 25:5907–5924. https://doi.org/10.1111/mec.13889

Russello MA, Waterhouse MD, Etter PD, and Johnson EA. 2015. From promise to practice: pairing non-invasive sampling with genomics in conservation. PeerJ 3:e1106. https://doi: 10.7717/peerj.1106

Sánchez R, and Medellín RA. 2007. Food habits of the threatened bat Leptonycteris nivalis (Chiroptera: Phyllostomidae) in a mating roost in Mexico. Journal of Natural History 41:1753–1764. https://doi.org/10.1080/00222930701483398

Schmidt DA, and Russello MA. 2025. Genomic Vulnerability of a Sentinel Mammal Under Climate Change. Molecular Ecology 34:e17688. https://doi.org/10.1111/mec.17688

Sherwin HA, Montgomery WI, and Lundy MG. 2012. The impact and implications of climate change for bats. Mammal Review 43:171–182. https://doi: 10.1111/j.1365-2907.2012.00214.x

Swets JA. 1988. Measuring the accuracy of diagnostic systems. Science 240:1285–1293. https://doi.org/10.1126/science.3287615

Sikes RS. the Animal Care and Use Committee of the American Society of Mammalogists. 2016. Guidelines of the American Society of Mammalogists for the use of wild mammals in research and education. Journal of Mammalogy 97:663–688. https://doi.org/10.1093/jmammal/gyw078

Smith AB, Beever EA, Kessler AE, Johnston AN, Ray C, Epps CW, et al. 2019. Alternatives to genetic affinity as a context for within-species response to climate. Nature Climate Change 9:787–794. https://doi.org/10.1038/s41558-019-0584-8

Thomas CD. 2010. Climate, climate change and range boundaries. Diversity and Distributions 16:488–495. https://doi.org/10.1111/j.1472-4642.2010.00642.x

Thuiller W, Georges D, Gueguen M, Engler R, Breiner F, Lafourcade B. et al. 2025. biomod2: Ensemble Platform for Species Distribution Modeling. https://doi:10.32614/CRAN.package.biomod2

Trejo-Salazar RE, Gámez N, Escalona‐Prado E, Scheinvar E, Medellín RA, Moreno‐Letelier A, et al. 2023. Historical, temporal, and geographic dynamism of the interaction between Agave and Leptonycteris nectar‐feeding bats. American Journal of Botany 110:e16222. https://doi.org/10.1002/ajb2.16222

Trejo-Salazar RE, Gasca-Pineda J, Hernández-Bolaños K, Hernández-Rosales DC, Tapia-López R, Aguirre-Planter E, et al. 2025. High genetic variation, low differentiation, and Pleistocene expansions of the migratory and endangered long-nosed tequila bat, Leptonycteris nivalis, inferred using both maternal and paternal genetic markers. PLoS One 20:e0316530. https://doi.org/10.1371/journal.pone.0316530

Urban MC, Swaegers J, Stoks R, Snook RR, Otto SP, Noble DW, et al. 2024. When and how can we predict adaptive responses to climate change? Evolution Letters 8:172–187. https://doi.org/10.1093/evlett/qrad038

Waldvogel AM, Feldmeyer B, Rolshausen G, Exposito-Alonso M, Rellstab C, Kofler R, et al. 2020. Evolutionary genomics can improve prediction of species’ responses to climate change. Evolution Letters 4: 4–18. https://doi.org/10.1002/evl3.154

Wanjala G, Astuti PK, Bagi Z, Kichamu N, Strausz P, and Kusza S. 2023. A review on the potential effects of environmental and economic factors on sheep genetic diversity: Consequences of climate change. Saudi Journal of Biological Sciences 30:103505. https://doi.org/10.1016/j.sjbs.2022.103505

Wells CP, Barbier R, Nelson S, Kanaziz R, and Aubry LM. 2022. Life history consequences of climate change in hibernating mammals: a review. Ecography 2022: e06056. https://doi.org/10.1111/ecog.06056

Weir BS, and Goudet J. 2017. A Unified Characterization of Population Structure and Relatedness. Genetics 206:2085–2103. https://doi.org/10.1534/genetics.116.198424

Wiens JJ. 2016. Climate-related local extinctions are already widespread among plant and animal species. PLoS Biology 14:e2001104. https://doi.org/10.1371/journal.pbio.2001104

Whitlock R. 2014. Relationships between adaptive and neutral genetic diversity and ecological structure and functioning: a meta‐analysis. Journal of Ecology 102:857–872. https://doi.org/10.1111/1365-2745.12240

Whitlock MC. 2015. Modern approaches to local adaptation. American Naturalist 186:S1–S4. https://doi.org/10.1086/682933

Associated editors: Giovani Hernández Canchola and Pablo Colunga Salas

Submitted: November 18, 2025; Reviewed: December 17, 2025

Accepted: May 8, 2026; Published online: May 29, 2026

THERYA, 2026, Vol. 17(2):259-276

DOI: 10.12933/therya.2026.6252 ISSN 2007-3364

Figure 1. Map of sampled locations; the number associated with each blue dot corresponds to the location number in Table 1.

Table 1. Sample sizes by location are shown, along with genetic diversity, observed heterozygosity (Ho), and expected heterozygosity (He) and Inbreeding coeficient Fis.

Locality

n

Ho

He

FIS

1

Nuevo León (nl)

14

0.1913

0.2123

0.0757

2

Michoacán (mich)

7

0.1867

0.2113

0.0793

3

Puebla (pue)

1

-

-

-

4

Zacatecas (zac)

1

-

-

-

5

Texas (tex)

19

0.2041

0.2157

0.0457

6

Hidalgo (hid)

5

0.1900

0.2122

0.0502

Overall

47

0.1866

0.2112

0.1097

Figure 2. Dendrogram of Euclidean distances illustrating genetic relationships among collection sites. The dendrogram was constructed using Euclidean distances calculated from genome-wide SNP data. Each terminal branch represents a site. The tree structure depicts the hierarchical clustering of sites based on genetic similarity, with shorter branches indicating greater genetic proximity. This analysis provides insights into population structure and the genetic differentiation among sampled localities.

Figure 3. Heatmap of pairwise FST values among sampled localities. FST estimates are shown for comparisons involving Nuevo León, Michoacán, Texas, and Hidalgo. Color intensity reflects the degree of genetic differentiation, with darker shades indicating higher FST values. Localities represented by a single individual (Zacatecas and Puebla) were excluded from the analysis to ensure reliable estimates. The heatmap provides a visual summary of population structure and genetic divergence across sampling sites.

Figure 4. Principal Component Analysis (PCA) of 47 individuals based on 13,112 SNPs. Each point represents an individual, colored by sampling locality: Hidalgo, Michoacán, Nuevo León, Puebla, Texas, and Zacatecas. The first two principal components (Axis1 and Axis5) are shown, explaining 2.758% and 2.698% of the total genetic variance, respectively. The analysis reveals population structure and genetic relationships among localities, with individuals from the same geographic region tending to cluster together.

Figure 5. Heatmap of pairwise relatedness estimates among 187 individuals. Relatedness coefficients were calculated using VCFtools based on genome-wide SNP data. The color gradient represents the degree of genetic relatedness, with warmer colors (e.g., red) indicating higher relatedness and cooler colors (e.g., white) indicating lower relatedness or unrelated individuals. Individuals are ordered along both axes by sample ID, and the matrix is symmetric, reflecting bidirectional pairwise comparisons. This analysis provides insights into familial relationships and population structure within the sampled cohort.

Figure 6. Maps show the suitability models for the L. nivalis species for the present (in green) and three future time periods, with major climate stability, until the end of the century. The colors blue, brown, and yellow represent the time periods that overlap.

Figure 7. A) Venn diagram showing the number of SNPs identified by three different filtering methods. RDA, PCAdapt, and LEA yielded 254 (43%), 121 (21%), and 211 (36%) SNPs, respectively. The overlapping regions indicate the number of SNPs shared between methods. Only two SNPs (<1%) were common to all three approaches. B) Adjusted p-values (-log10) are shown for each method, with values of 0 indicating strong statistical support after correction. Percentages represent the proportion of SNPs relative to the total number of unique SNPs identified across all methods.

Figure 8. Risk of non-adaptedness (RONA) analysis across sampling localities for three bioclimatic variables (Bio.9, Bio.3, Bio.14). A) Bar plots showing RONA values per locality, representing the statistical association between heterozygosity and the risk of maladaptation under future climate change scenarios for each bioclimatic variable. B) Geographic distribution of RONA values across the sampled range. Warmer colors (closer to red) indicate localities with a higher risk of non-adaptedness if the associated climatic variable changes. Values for selected representative localities (Bio.9, Bio.3, Bio.14) are shown to illustrate spatial variation in climate-associated genetic vulnerability.