Evolutionary and environmental influences on cranial shape
and forearm length variation in Corynorhinus mexicanus
1Red de Biología y Conservación de Vertebrados, Instituto de Ecología A.C., Xalapa-Enríquez, Veracruz, Mexico. E-mail: isachar26@hotmail.com (ILL-C)
2Departamento de Zoología, Instituto Politécnico Nacional, Escuela Nacional de Ciencias Biológicas, Mexico City, Mexico. E-mail: artibeus2@aol.com (JO)
3Centro de Investigación en Biodiversidad y Conservación, Universidad Autónoma del Estado de Morelos, Cuernavaca, Morelos, Mexico. E-mail: ospinagarcess@gmail.com (SMO-G)
4Centro de Investigaciones Tropicales, Universidad Veracruzana, Xalapa-Enríquez, Veracruz, Mexico.
*Corresponding author: cmacswiney@uv.mx
Phenotypic variation in bats reflects the combined influence of environmental and evolutionary factors. Body size often responds to thermal gradients (Bergmann’s rule), whereas cranial and mandibular morphology are often associated with evolutionary history and diet, which are further shaped by resource seasonality. Corynorhinus mexicanus, an endemic bat of Mexico’s temperate forests, was analyzed to assess if environmental factors, evolutionary history, and their interaction influence forearm size and cranial and mandibular shape. Linear mixed models were fitted using forearm measurements and principal components derived from environmental variables, incorporating previously documented intraspecific lineages. Models were compared under annual and monthly temperature variation, accounting for reproductive timing and sex-specific energetic investment. Cranial and mandibular variation were examined using a two-dimensional geometric morphometric protocol to assess the influence of lineage and environmental seasonality as a proxy for dietary composition. Forearm length did not differ between lineages in females but did in males, with those from the Sierra Madre Occidental lineage being smaller. In both sexes, body size followed Bergmann’s rule. Cranial differences between lineages were restricted to the lateral view of the braincase and the ventral view of the rostrum. Overall, cranial variation appeared to be primarily driven by evolutionary history, whereas forearm size reflected a combined influence of evolutionary and environmental factors. Maximum temperatures during the reproductive period influenced forearm length in males, while minimum temperatures affected females. Finally, we highlighted the importance of considering temporal extent and environmental resolution in eco-evolutionary studies.
Key words: Bergmann’s rule, body size, diet, heat conservation, heat dissipation, Vespertilionidae
La variación fenotípica en murciélagos refleja la influencia conjunta de factores ambientales y evolutivos. El tamaño corporal suele responder a gradientes térmicos (regla de Bergmann), mientras que la morfología craneal y mandibular suele relacionarse con la historia evolutiva y la dieta, modulada por la estacionalidad. Corynorhinus mexicanus, un murciélago endémico de los bosques templados de México, fue analizado para evaluar si factores ambientales, la historia evolutiva y su interacción influyen en el tamaño del antebrazo y en la forma y tamaño del cráneo y la mandíbula. Se ajustaron modelos lineales mixtos con medidas de antebrazo y componentes principales de variables ambientales, incorporando los linajes intraespecíficos conocidos. Los modelos se compararon considerando la variación anual y mensual de la temperatura, en relación con los periodos reproductivos y la asincronía energética entre sexos. La variación craneal y mandibular se evaluó mediante un protocolo de morfometría geométrica en dos dimensiones, analizando el efecto del linaje y la estacionalidad como un proxy de la dieta. La longitud del antebrazo no difirió entre linajes en hembras, pero sí en machos, siendo los del linaje de la Sierra Madre Occidental más pequeños. El tamaño corporal de ambos sexos siguió el patrón de la regla de Bergmann. Las diferencias craneales entre linajes se restringieron a la vista lateral de la caja craneal y a la ventral del rostro. En conjunto, los resultados indican que la variación craneal está determinada principalmente por la historia evolutiva, mientras que el tamaño del antebrazo responde de forma combinada a factores evolutivos y ambientales. Las temperaturas máximas durante el periodo reproductivo influyeron principalmente en la longitud del antebrazo en machos, mientras que las mínimas lo hicieron en las hembras. Finalmente, se destaca la importancia de considerar la escala temporal y la resolución ambiental en los estudios eco-evolutivos.
Palabras clave: Conservación del calor, dieta, disipación del calor, regla de Bergmann, tamaño corporal, Vespertilionidae
© 2026 Asociación Mexicana de Mastozoología, www.mastozoologiamexicana.org
Variation in phenotypic traits within species is a common phenomenon with important ecological consequences (e.g., community structure and demographic patterns; Bolnick et al. 2003; 2011). Such variation can arise from several factors, including adaptation to environmental conditions (e.g., Stevens et al. 2016; Maluleke et al. 2017), exploitation of specific resources (e.g., Encarnação et al. 2005; Hedrick 2021), interspecific competition (Salinas-Ramos et al. 2020), sex (e.g., Myers 1978; Camargo and de Oliveira 2012; Ulpian and Rossi 2017), genetic characteristics of individuals (e.g., Majerus and Mundy 2003), and the evolutionary history of populations (e.g., Uvizl and Benda 2021; Gorobeyko et al. 2025). Exploring how phenotypic traits relate to intrinsic and extrinsic factors can yield valuable insights into species’ responses to environmental change, such as global climate shifts (Paltrinieri et al. 2025), and can support biodiversity conservation by serving as potential indicators of metapopulation persistence (Gibert 2016).
Intraspecific phenotypic variation in bats has been widely documented, particularly in body size (e.g., Myers 1978; Williams and Findley 1979; Ulpian and Rossi 2017; Russo et al. 2024) and cranial and mandibular morphology (e.g., Marchan-Rivadeneira et al. 2012; Morales et al. 2018; Hedrick 2021), both of which exhibit diverse responses associated with evolutionary history and environmental conditions.
Body size variation within species has traditionally been attributed to adaptive responses to environmental clines (Ashton 2002; Clauset and Erwin 2008; Meiri 2008; Chown and Gaston 2010). In bats, the association between body size and environmental factors has been explained through several mechanisms, often linked to physiological constraints (Safi et al. 2013; Castillo-Figueroa 2022; Alston et al. 2023). In a physiological context, variation in bats’ body size has been negatively correlated with the environmental temperature where the organism lives (Bergmann 1848; Safi et al. 2013). This pattern referred to as Bergmann’s rule, although originally proposed to describe interspecific patterns, has also been applied at the intraspecific level, predicting larger individuals in colder regions as a result of enhanced heat conservation (Ashton et al. 2000; Watt et al. 2010). Despite the main reason behind body size variation according to the Bergmann’s rule is the heat conservation (Bergmann 1848), other reasons have been proposed to produce the same pattern. One particularly is relating to heat stress, where larger individuals have lower critical thermal maxima, favoring smaller body sizes in warmer environments (Smith et al. 1995; Peralta-Maraver and Rezende 2021). Both mechanisms (heat conservation and heat dissipation) may produce the same pattern of body size variation, but the environmental constraints are different, and they obey opposite boundaries of thermal tolerance of individuals. Another phenomenon behind body size variation is related to trophic resource availability, suggesting that individuals tend to be larger in more productive environments (McNab 2010), where higher primary productivity, temperature, and precipitation provide greater energy and resource inputs (Chu et al. 2016). As noted, the pattern exhibited here is opposite to Bergmann’s rule.
Regarding intraspecific variation in cranial and mandibular morphology, both structures are often influenced by several interacting factors, including evo-lutionary history and dietary adaptations (e.g., Marchan-Rivadeneira et al. 2012; Morales et al. 2018; Hedrick 2021), with the former reflecting evolutionary history effects and the latter shaped by diet influence.
In insectivorous bats, diet composition may depend, among other factors, on prey availability, which is strongly influenced by environmental conditions, as insect diversity and abundance fluctuate with seasonal changes in temperature and precipitation, as well as primary productivity (e.g., Borer et al. 2012; Moosman et al. 2012; Lind et al. 2017; Zhu et al. 2024). These fluctuations are most evident along latitudinal and altitudinal gradients (Alston et al. 2023), where seasonal variations in primary productivity may alter insect abundance and diversity, potentially forcing bats to broaden their prey spectrum (e.g., Moosman et al. 2012) and driving local cranial and mandibular morphological adaptations (e.g., Mendes et al. 2024).
Corynorhinus mexicanus G.M. Allen, 1916, is an insectivorous bat endemic to temperate forests located on the slopes of the Sierra Madre Occidental (SMOc) and the Trans-Mexican Volcanic Belt (TMVB) in Mexico (López-Cuamatzi et al. 2024). Previous study suggested this species comprises two mitochondrial lineages, each restricted to each mountain range, showing low genetic divergence (López-Cuamatzi et al. 2024). Variations in body size (measured as forearm length) and cranial morphology have been documented and attributed to evolutionary history (mitochondrial lineages) (López-Cuamatzi et al. 2024), although the influence of environmental factors has not been formally tested.
Recognizing environmental conditions and evolutionary history as primary drivers of intraspecific phenotypic variation, this study aims to elucidate the relationship between forearm length and cranial and mandibular shape variation in C. mexicanus with thermal and precipitation variables, mitochondrial linages, as well as the interaction between these components.
Specifically, for body size (measured as forearm length), we tested the following non-mutually exclusive hypotheses and predictions (Figure 1):
However, a methodological limitation in studies testing the heat conservation and heat-stress mechanisms of Bergmann’s rule is the assumption that environmental pressures –commonly represented by monthly, quarterly, or annual temperature extremes (e.g., Castillo-Figueroa 2022; Alston et al. 2023; Paltrinieri et al. 2025)– affect all individuals uniformly and continuously, regardless of their phenological stage. This simplification could bias interpretations of temperature–body size relationships. In temperate zones, maximum and minimum values of temperatures and precipitation occur in specific seasons that may not coincide with critical periods such as reproduction, when the energy demands of Holarctic bats are highest (Kurta and Kunz 1987; McLean and Speakman 1999). This complexity is further amplified by sex-specific phenological patterns, as males and females of some Nearctic species often reach peak reproductive energy demands at different times of the year (Gustafson 1979; Oxberry 1979). Considering this, we test a modification of GED hypothesis but considering the phenological demands of individuals:
Materials and methods
Collection of forearm data. We compiled forearm measurements of individuals from various localities representing the two genetic groups identified in C. mexicanus (Figure S1 of Supplementary Data SD1). These data were obtained from the specimen database of taxidermy-preserved individuals reviewed in López-Cuamatzi et al. (2024). Individuals catalogued as juvenile were excluded from the analyses.
Collection of cranial and mandible shape data. To characterize the cranial and mandibular shape variation, we photographed the mandible and dorsal, ventral and lateral views of the cranium using a Nikon D3000 reflex camera (Nikon Corporation, Tokyo, Japan) and a Nikkor 2.8F 60 mm macro lens (Nikon Corporation, Tokyo, Japan). Two-dimensional landmark and semi-landmark configurations were digitized on the cranial and mandibular images using the software tpsUtil and tpsDig2 (Rohlf 2017, 2019). The number of landmarks and semi-landmarks varied between views (see López-Cuamatzi et al. 2024 for details).
Landmark configurations were aligned using Generalized Procrustes Analysis to remove non-shape variation by standardizing origin, scale, and rotation (Rohlf 1990). Bone contours were digitized as semi-landmarks and aligned by minimizing Procrustes distances to the consensus curve (Pérez et al. 2006; Gunz and Mitteroecker 2013). Because skull views are bilaterally symmetrical, shapes were averaged across sides using bilateral symmetry analysis (Klingenberg 2015). Shape variables were represented as Procrustes coordinates, and size was estimated as centroid size (CS), calculated as the square root of the summed distances from landmarks to the centroid (Bookstein 1997). Analyses were performed in the "R-package" geomorph v.4.0.1 (Adams et al. 2025) on R v.4.4.2 (R Core Team 2024).
To reduce sample size requirements, each cranial view was divided into two developmentally and functionally recognized modules (rostral and basicranial), commonly described in mammals, including bats (Marroig et al. 2009; Porto et al. 2009). Although cranial modularity has previously been reported in C. mexicanus (López-Cuamatzi et al. 2024), modularity was re-evaluated because the current analysis excludes observations from Sierra Madre Oriental (formally C. leonpaniaguae) included in the previous study. Modularity was evaluated using the covariance ratio (CR) with 1,000 permutations, implemented in the geomorph R-package v.4.0.1 (Adams et al. 2025). When the CR indicated a lack of modularity, the entire structure was analyzed as a single unit. To compare shape independently of size, we fitted linear models between shape and centroid size for each view or module and used the residuals as size-free shape coordinates for subsequent analysis, using Procrutes ANOVA models.
Collection of environmental data. We used the geographic coordinates of each voucher specimen to extract environmental data from climatic raster files. To test the General Environmental–Driven Hypothesis (GED), we compiled a dataset with values from WordClim’s bioclimatic variables (Fick and Hijmans 2017): Bio5 (Max Temperature of Warmest Month) and Bio10 (Mean Temperature of Warmest Quarter) as proxies for heat dissipation and Bio6 (Min Temperature of Coldest Month) and Bio11 (Mean Temperature of Coldest Quarter) as proxies for heat conservation. We used minimum and maximum temperature data instead of mean temperature because our objective was to assess the thermal stress experienced by individuals, and mean temperature does not capture the minimum and maximum thermal exposures that may generate such stress. Additionally, thermal variability also influences physiological tolerance, affecting individual fitness (see Bozinovic et al. 2011) and spatial distribution (i.e., Zimmermann et al. 2009).
Environmental data were extracted at three spatial resolutions, 30-arcseconds (~1 km), 2.5-arcminute (~4.5 km), and 5-arcminute (~9.25 km), to account for resolution-related uncertainty. The 5-arcminute resolution approximates the dispersal capacity reported for Corynorhinus species (Fellers and Pierson 2002) but may introduce environmental error due to its coarse resolution. In contrast, 30-arcsecond raster offer finer environmental detail but may not adequately capture broader gradients relevant to Corynorhinus movements and habitat use.
Following (Alston et al. 2023), to evaluate the Resource Availability–Driven Hypothesis (RAD), we used the MODIS monthly Net Primary Productivity (NPP) dataset (NASA Earth Observatory 2025) at a 0.1° spatial resolution (~10 km) as a proxy for insect abundance. Primary productivity is strongly associated with insect biomass (Lind et al. 2017) and has previously shown a moderate positive correlation (r = 0.66; Borer et al. 2012). For each specimen locality, we extracted 24 years of data (2001–2024) and calculated the median NPP value. Additionally, we used Bio18 (Precipitation of the Warmest Quarter) at 30-arcsecond and 5-arcminute resolutions as a proxy for trophic resource availability for bats, reflecting insect abundance during summer (see Frick et al. 2010).
To evaluate the Reproductive-Cost by Thermal Regulation (RCTR) hypothesis, we compiled historical monthly climate data from WorldClim (Fick and Hijmans 2017) for each individual. This dataset consists of raster files at 5-arcminute resolution, including average minimum and maximum temperature (°C) from 1950 to 2024. We extracted monthly data for each of the 75 years available and calculated the median for each month, resulting in 75-year monthly medians. Median values were used to minimize the influence of outliers, which are common in environmental datasets. Two sex-specific datasets were generated to reflect periods of peak energetic demand: August–November for males, when C. mexicanus invests in spermatogenesis and sperm maturation (León-Galván et al. 2005), and March–July for females, corresponding to gestation and lactation (López-Wilchis 1989).
To test the Dietary–Driven Hypothesis (DD), we generated a dataset including WorldClim’s variables Bio4 (Temperature Seasonality), Bio15 (Precipitation Seasonality) (Fick and Hijmans 2017), together with Bio 18 and NPP, as proxies for insect richness and abundance. We used Bio4 and Bio15 because insect diversity can be influenced by thermal and humidity seasonality along altitudinal gradients, such as in mountainous environments (i.e., Wilson et al. 2007; Montañez-Reyna et al. 2022). As with the GDE, RAD and RCTR datasets, environmental data for DD were extracted at resolutions of 30-arcseconds and 5-arcminute.
Treatment of environmental variables. For the GED dataset, we compared environmental variable values (Bio5, Bio6, Bio10, and Bio11) across three spatial resolutions (30-arcseconds, 2.5-arcminute, and 5-arcminute) using Mantel tests based on Spearman’s rank correlation to reduce sensitivity to outliers and non-linear relationships. Euclidean distances were used to construct distance matrices for each environmental variable, and Mantel tests were performed with 9,999 permutations using the “vegan” R-package in R (Oksanen et al. 2025). This analysis evaluated whether datasets of differing resolutions produced consistent results and whether resolution influenced model outcomes. We then performed a Principal Component Analysis (PCA) on each dataset of bioclimatic variables to reduce dimensionality and summarize environmental information. The number of components to retain was determined using Horn’s parallel analysis (Dinno 2009), as implemented in the fa.parallel function of the R-package psych (Revelle and Revelle 2015). Retained components were subsequently used in further analyses. For the RCTR and DD datasets, PCAs were conducted as described above. In the RCTR dataset, PCAs were performed separately by sex due to differences in the number of monthly climatic variables considered (Figure 2). For the RAD hypothesis, Bio18 and NPP were used directly without dimensionality reduction through PCA; however, both variables were standardized to a mean of zero and unit variance. Additionally, we assessed their correlation using a Spearman correlation test to avoid collinearity in subsequent modeling.
Sexual dimorphism in forearm, and cranial and mandible shape. Because the RCTR hypothesis implies sex-specific associations with climatic variables, we tested for sexual dimorphism. Forearm length was analyzed using Welch’s t-test, appropriate for unbalanced designs, after confirming normality (Shapiro-test, P > 0.05) and homoscedasticity (Levene’s test, P > 0.05). For each cranial and mandibular view or module, we evaluated the effect of sex on size-free shape using Procrustes ANOVA in Geomorph, and analyzed centroid size using linear models with random permutations in the “RRPP” R-package (Collyer and Adams 2018). Both analyses were conducted with 1,000 permutations.
Modeling unit. For the forearm analysis, multiple specimens were collected from some localities (i.e., identical geographic coordinates). Therefore, each individual was treated as the unit of analysis, while acknowledging that some individuals were spatially clustered. To account for this structure, we fitted linear mixed models (LMMs) (Korner-Nievergelt et al. 2015), including locality as a random effect, with individuals grouped according to shared geographic coordinates.
Because several localities were represented by a single individual, the dataset contained a substantial proportion of singleton groups (i.e., clusters with only one observation). Such imbalance can lead to random-effect variance estimates approaching zero. Nevertheless, simulation studies indicate that LMMs are generally robust to high proportions of singleton clusters, even under moderate sample sizes (Bruyndonckx et al. 2018). Although some approaches recommend removing random effects when variance estimates are negligible, we retained the hierarchical structure given the inherent spatial grouping of the data. We therefore report the LMM results while acknowledging that, in some cases, the estimated random-effect variance was close to zero.
In contrast, for the morphometric analyses we used the mean cranial and mandibular shape per locality as the modeling unit. Although individual-level variation provides valuable insight into intraspecific morphological patterns, incorporating individuals directly into the models could bias inference regarding environmental or lineage effects. By aggregating shape data at the locality level, we reduced pseudo-replication and ensured that the modeling unit matched the scale at which environmental variables were measured.
Statistical analyses: forearm. We fitted linear mixed-effects models (LMMs) using the lmer function from the “lme4” R-package (Bates et al. 2015) to test the proposed hypotheses. We used forearm length as the response variable, locality as a random effect, and the LD, GED, RAD, and RCTR datasets as predictors. Because the LD hypothesis is not mutually exclusive with the rest of the hypotheses, we also fitted models including both additive (e.g., LD + GED) and interaction terms (e.g., LD * GED) to evaluate their combined and interactive effects. Model assumptions, including dispersion, residual uniformity, and outlier detection, were evaluated using the simulateResiduals function from the “performance” R-package (Lüdecke et al. 2021).
Model comparison was conducted using AIC, AICc, and ΔAICc, with selection based on the lowest AIC and AICc values (Korner-Nievergelt et al. 2015). The significance of each term of the model was assessed considering (P < 0.05). An intercept-only null mixed-effects model was fitted as a baseline to assess whether the observed patterns differed from those expected under random variation (Gotelli and Graves 1996). When multiple models showed comparable support (ΔAICc < 2) and were not nested, all top-ranked models were considered. However, when competing models were nested, we used likelihood ratio tests (through the anova function in the “stats” R-package [R Core Team 2024]) to evaluate whether the more complex model provided a significantly better fit. If the null hypothesis was not rejected, the simpler model was retained.
The intraclass correlation coefficient (ICC) and marginal and conditional R² values of models selected were computed using the icc and r2 functions in the “performance” R-package (Lüdecke et al. 2021). In instances where variance components of the random effects could not be reliably estimated, conditional R² was not computed; instead, marginal R² and variance components are presented. Additionally, when interaction terms in the model selected were statistically supported, simple slopes of continuous predictors within each factor level were estimated using the emtrends function from the “emmeans” R-package (Lenth and Piaskowski 2025).
Statistical analyses: cranial and mandible shape. The LD and DD hypotheses for cranial and mandibular shape were tested using Procrustes ANOVA in the Geomorph. Linear models were fitted on size-free shape coordinates derived from the residuals of the allometric regression between centroid size and shape. Predictors variables included the DD dataset and mitochondrial linage (LD hypothesis). As the two hypotheses are not mutually exclusive, we also tested models with interaction terms between lineage and environmental variables.
Models were compared using ANOVA, and the significance of each term was evaluated (P < 0.05). Model selection followed a nested structure, for example, LD, LD+DD, and LD*DD formed one set, while DD, LD+DD, and LD*DD formed another. However, the simple LD and DD models are not nested within each other. When the best-fitting models were simple, both were discussed. If a more complex model was supported in one set but not the other, the simpler nested model was preferred. A null model was also included for baseline comparison. Finally, centroid size was analyzed with linear models implemented in “RRPP” R-package (Collyer and Adams 2018), following the same logic as using 1,000 permutations.
After identifying significant predictors, we conducted individual tests to assess the association with size-free shape and centroid size. When lineage had a significant effect on size-free shape, we performed paired comparisons between the group means, and their significance was tested by permutation testing (10,000 permutations) and Bonferroni’s correction, comparing the observed Procrustes distance (procD) with that obtained from the random assignment of observation to groups in the R-package “Morpho” v. 2.7 (Schlager 2017). Differences in mean shape configurations among lineages were assessed using a Discriminant Function Analysis (DFA) based on the first ten principal components (PCs) obtained from a prior Principal Components Analysis. For views and modules exhibiting significant differences between lineages, Mahalanobis distances and corresponding P-values were calculated using permutation tests of the original data matrix with 1,000 replicates. Classification accuracy was evaluated using cross-validation, where each specimen was sequentially excluded from the training dataset and classified using discriminant functions built from the remaining individuals. For environmental variables we applied two-block Partial Least Squares (PLS) analyses with 999 iterations. For centroid size, significant lineage effects were evaluated using pairwise-distance tests with 999 iterations, while environmental variables were tested using Pearson or Spearman correlations depending on the normality of centroid size residuals.
Results
Environmental data summary. Regarding the environmental variable datasets used for GED hypothesis, Mantel tests revealed a strong correlation between the 2.5-arcminute and 5-arcminute environmental datasets (r = 0.95, P < 0.05), whereas the correlation between the 30-arcsecond and 5-arcminute datasets was notably lower (r = 0.67, p < 0.05). Based on these results, subsequent analyses were conducted using the 30-arcsecond and 5-arcminute datasets. Monthly climatic data for testing the RCTR hypothesis were analyzed at 5-arcminute resolution, while data for the DD hypothesis were analyzed at both resolutions to ensure comparability across analyses. Under the RAD hypothesis, Bio18 and NPP were only weakly correlated (NPP-Bio18 at 30-arcseconds, rho = -0.3, P < 0.001; NPP-Bio18 at 5-arcminute, rho = -0.29, P = 0.001), allowing both variables to be included as predictors in subsequent analyses.
Sexual dimorphism and modularity in cranial shape. The forearm dataset used to test sexual dimorphism consisted of 131 observations (63 females and 68 males). Welch’s t-test revealed significant sexual dimorphism, with females exhibiting greater forearm length than males (t125.35 = -7.059, P < 0.05; Figure S2 of Supplementary Data SD1). Consequently, all subsequent models were performed separately for each sex.
The number of observations for cranial and mandibular shape analyses varied among datasets (Table S1, Supplementary Data SD1). The modularity test supported the division into rostral and basicranial modules in the lateral (CR = 0.7578, P = 0.001) and ventral (CR = 0.7665, P = 0.001) views, whereas the dorsal view lacks statistically significant modularity (CR = 0.964, P = 0.44). Procrustes ANOVA detected sexual dimorphism in size-free shape only in the dorsal cranial view and lateral basicranial view (Table S1, Supplementary Data SD1). Post hoc tests indicated that in the dorsal view, females and males differed (procD = 0.003, P = 0.049), primarily in nasal bone width and maxilla length, while in the lateral basicranial view (procD = 0.008, P = 0.013), differences were concentrated at the base of the braincase (Figure S3, Supplementary Data SD1). Centroid size also exhibited sexual dimorphism in lateral rostral, ventral basicranial, and mandibular views (Table S2, Supplementary Data SD1). The trend of dimorphism varied across structures: females were larger than males in the lateral rostrum and mandible, whereas the opposite pattern occurred in the ventral braincase (Figure S4, Supplementary Data SD1).
PCA environmental variables. For the GED hypothesis, the first two principal components (PCs) explained 99.1% (PC1 = 59.87%; PC2 = 39.23%) of the total variance at the 30-arcsecond resolution and 98.83% (PC1 = 64.67%; PC2 = 34.16) at the 5-arcminute resolution. Horn’s parallel analysis suggested retaining two components in each dataset.
At the 30-arcsecond resolution, PC1 represented a general temperature gradient, with all bioclimatic variables loading positively, and Bio10 and Bio11 showed the strongest contributions to this component. PC2 differentiated variables associated with lower temperature contrast (lower maximum temperatures; positive loadings) from those linked to higher temperature contrast (higher maximum temperatures; negative loadings) between maximum and minimum temperature conditions, with Bio5 and Bio6 contributing most strongly to this axis. At the 5-arcminute resolution, we observed a similar arrangement of components to that obtained from the 30-arcsecond resolution PCA. The only difference was the opposite sign of the PC1 loadings (negative loadings), which reflects the arbitrary orientation of principal component axes and does not indicate structural differences between analyses.
For the RCTR hypothesis, the first two principal components explained 96.42% of the total variance in males (PC1 = 79.8%; PC2 = 16.5%) and 94.39% in females (PC1 = 65.0%; PC2 = 29.2%). Although Horn’s parallel analysis indicated that one component should be retained in males and two in females, the first two components were retained in both sexes to maintain comparability. In males, PC1 represented a general thermal gradient across months, with all maximum and minimum temperature variables loading positively, indicating that this axis captured overall temperature conditions during the August–November period. Maximum temperatures across months, together with minimum temperatures of August and September, showed strong contributions to this component. PC2 distinguished maximum temperatures (positive loadings) from minimum temperatures (negative loadings), with November and October minimum temperatures exhibiting the strongest negative contributions. This pattern reflects a contrast between daytime and nighttime thermal conditions, particularly emphasizing the influence of late-season minimum temperatures, and therefore represents a gradient of increasing dominance of nocturnal cooling toward the end of the period.
In females, PC1 represented a general thermal gradient across months, with all maximum and minimum temperature variables loading positively. The strongest contributions to this component were observed for minimum temperatures, particularly May and June, together with maximum temperatures from early spring months (April and March), indicating that this axis captured overall thermal conditions during the March–July period. PC2 contrasted high summer maximum temperatures (June and July, which showed the strongest positive loadings, followed by May) with lower minimum temperatures of early spring months (especially March and April, which exhibited the strongest negative loadings). This pattern reflects a seasonal contrast between peak summer heat and cooler spring nights, representing the degree of thermal intensification from early spring to midsummer. In simple terms, PC2 indicates how midsummer and early spring temperatures were similar.
For the DD hypothesis, the first two principal components (PCs) at 30-arcsecond resolution explained 83.9% of the total variance (PC1 = 68.3%; PC2 = 15.6%), whereas the first three PCs at 5-arcminute resolution accounted for 84.7% (PC1 = 70%; PC2 = 14.7%) of the variance. At 30-arcsecond resolution, PC1 was primarily associated positively with Bio4 and negatively with NPP, while PC2 was negatively associated with Bio15 and Bio18. This pattern reflects a dominant environmental gradient in which increasing seasonality is associated with decreasing productivity. Similar patterns of variable contribution were observed at the 5-arcminute resolution.
Testing hypotheses: forearm. The forearm dataset used for modeling included 131 individuals grouped into 65 localities, which were treated as a random effect in the LMMs. When separated by sex, the female dataset comprised 63 individuals (31 from TMVB and 32 from SMOc), whereas the male dataset included 68 individuals (48 from TMVB and 20 from SMOc).
Based on AICc and ΔAICc criteria, the best-supported model for males corresponded to the LD hypothesis in the 30-arcsecond dataset and to the GED hypothesis in the 5-arcminute dataset. However, in both cases, alternative models also showed comparable support (ΔAICc < 2), specifically the model representing the RCTR hypothesis in the 30-arcsecond dataset, and the LD and RCTR hypotheses in the 5-arcminute dataset (Table 1). Because the models retained under the ΔAICc criterion were not nested, all of them were considered in subsequent interpretations.
For females, the best-supported model included the GED effect in the 30-arcseconds dataset and the LD × RCTR interaction in the 5-arcminute dataset. However, in the 30-arcseconds dataset, alternative models also showed comparable support (ΔAICc < 2), particularly LD + GED, LD * GED, and LD * RCTR. Because GED, LD + GED, and LD * GED are nested models, likelihood ratio tests were performed. These tests indicated that the more complex models did not provide a significantly better fit than the simpler ones (GED vs. LD + GED: X21 = 0.48, P = 0.48; GED vs. LD * GED: X23 = 5.74, P = 0.12). Therefore, only the GED and LD * RCTR models were retained for subsequent interpretation. All fitted models for both sexes and datasets, including those not selected, met the assumptions of dispersion, uniformity, and absence of outliers.
In males, the LD model revealed a significant effect of lineage, with SMOc individuals showing lower values of forearm than TMVB (β = -0.71, SE = 0.25, CI95% = -1.23 – -0.22, t14.71 = -2.84, P = 0.013) (Figure 3A). For this model, random-effect variance was negligible (σ² = 0.001, SD = 0.04), with ICC = 0 and marginal R² = 0.11. In the GED model (5-arcminute dataset) for males, PC2 had a significant positive effect (β = 0.25, SE = 0.11, CI95% = 0.03 – 0.47, t30.25 = 2.19, P = 0.035) (Figure 3B), whereas PC1 was not significant (β = 0.14, SE = 0.08, CI95% = -0.02 – 0.3, t24.19 = 1.74, P = 0.09). Since positive values of PC2 indicate localities with greater influence of cold conditions and less dominance of extreme heat, these results indicates that males tend to be larger in colder locations and smaller in locations dominated by extreme heat. For this model, random-effect variance was low (σ² = 0.06, SD = 0.26; ICC = 0.036), with marginal and conditional R² values of 0.15 and 0.18, respectively.
Under the RCTR model, PC1 had a significant negative effect on forearm length (β = -0.13, SE = 0.048, CI95% = -0.23 – -0.03, t17.5 = -2.79, P = 0.012) (Figure 3C), whereas PC2 was not significant (β = -0.13, SE = 0.1, CI95% = -0.33 – 0.06, t33.05 = -1.29, P = 0.2). This result indicate that males inhabiting warmer environments, particularly during August–September, tend to exhibit shorter forearm length. For this model, random-effect variance was low (σ² = 0.06, SD = 0.24; ICC = 0.014), with marginal and conditional R² of 0.13 and 0.14, respectively.
In females, the GED model showed a significant negative effect of PC1 on forearm length (β = -0.30, SE = 0.09, CI95% = -0.48 – -0.11, t60 = -3.24, P = 0.001) (Figure 3D), while PC2 was not significant (β = 0.04, SE = 0.1, CI95% = -0.15 – 0.25, t60 = 0.48, P = 0.63). These findings suggest that females inhabiting colder environments tend to exhibit greater forearm length, while individuals from warmer environments show reduced forearm size. For this model, the random-effect variance was zero, resulting in ICC = 0 and a marginal R² of 0.15.
The LD × RCTR model for females revealed significant effects of lineage, PC1, and PC2, as well as significant interactions between lineage and both principal components (Figure 3E-F). The SMOc lineage exhibited lower forearm values than the reference lineage (TMVB) (β = -2.39, SE = 1.02, CI95% = -4.25 – -0.50, t36.22 = -2.34, P = 0.025). PC1 was negatively associated with the forearm length (β = -0.76, SE = 0.22, CI95% = -1.16 – -0.35, t26.76 = -3.40, P = 0.002), whereas PC2 showed a positive association (β = 1.28, SE = 0.50, CI95% = 0.36 – 2.16, t28.9 = 2.58, P = 0.015). Significant interactions were detected between lineage and PC1 (β = 0.76, SE = 0.24, CI95% = 0.32 – 1.18, t27.2 = 3.18, P = 0.004) and between lineage and PC2 (β = -1.33, SE = 0.57, CI95% = -2.32 – -0.26, t28.49 = -2.35, P = 0.026), indicating lineage-specific relationships between principal components and the forearm length (Figure 3G). Random-effect variance was small (σ² = 0.11, SD = 0.34), with ICC = 0.012 and marginal and conditional R² values of 0.24 and 0.25, respectively.
Simple slope analyses showed that PC1 was negatively associated with the forearm length in TMVB (β = -0.75, CI95% = -1.21 – -0.30) but not in SMOc (β = -0.005, CI95% = -0.18 – 0.16) (Figure 3E). Similarly, PC2 had a positive effect in TMVB ( = 1.25, CI95% = 0.25 – 2.26) (Figure 3F) but no significant association in SMOc (β = -0.018, CI95% = -0.57 – 0.54). Because PC1 represents a general thermal gradient during March–July, the negative association indicates that females inhabiting generally warmer localities exhibited shorter forearms. In contrast, PC2, which captures seasonal changes from early spring to midsummer, suggests that females from localities with a more pronounced seasonal transition between these periods had longer forearms.
Testing hypotheses: cranial and mandible shape and size. For most cranial views, fitted shape models using either the 30-arcsecond or 5-arcminute datasets did not outperform the null models. Only the lateral-view modules and ventral rostrum view showed better fit than the null models. None of the centroid size models outperformed the null models.
In the lateral braincase module, for females, the best-fitting model in both datasets (30-arcseconds and 5-arcminute) was the LD*DD interaction. None of the individual terms were significant, but PC1 (F1, 18 = 1.85, P = 0.053) and PC2 (F1, 18 = 1.91, P = 0.058) approached significance. In males, no model outperformed the null model in the 30-arcsecond resolution, but at 5-arcminute resolution, the best-fitting model was also LD*DD, with lineage (F1, 26 = 2.46, P = 0.017) and lineage*PC2 (F1, 26 = 2.35, P = 0.024). Pairwise analyses showed differences between TMVB and SMOc males in the frontal and basal regions of the braincase (procD = 0.013, P = 0.017) (Figure 4A). Discriminant Function Analysis (DFA) based on the first 10 PCs supported these differences (Mahalanobis D = 2.06, p < 0.001), with cross-validation correctly classifying 84% of TMVB and 50% of SMOc individuals (Figure 4B). PLS analyses revealed a positive, though non-significant, trend between shape and PC2 within lineages (Table 2).
In the lateral rostrum module, models using the 30-arcsecond dataset did not outperform the null model. In contrast, for the 5-arcminute dataset, the best-fitting model included LD*DD interactions, with significant lineage*PC1 (F1, 41 = 2.96, P = 0.023) and lineage*PC2 (F1, 41 = 4.13, P = 0.007) effects. PLS analyses within lineages indicated positive but non-significant relationships for both interactions (Table 2).
In ventral rostrum view, the best fitting model included only the lineage term (F1, 43 = 2.36, P = 0.022), which was significant. Pairwise analyses showed that TMVB and SMOc individuals differed in nasal bone size and length (procD = 0.012, P = 0.012) (Figure 4C). DFA based on the first 10 PCs confirmed these differences (Mahalanobis D = 1.26, P = 0.003), with cross-validation correctly classifying 75% of TMVB and 37.5% of SMOc individuals (Table 2, Figure 4D).
Discussion
In our study, models associated with the LD and RTCR hy-potheses provided the most consistent explanation for fore-arm variation in C. mexicanus. Although the GD hypothesis was similarly supported in both sexes, its performance varied depending on the spatial resolution of the environmental dataset. Given that hypothesis support contingent upon predictor resolution may reflect scale-dependent effects rather than clear biological mechanisms, we limit further discussion of the GED hypothesis. Nonetheless, we acknowledge that when GD hypothesis was supported, the relationship between forearm length and environmental variables followed the pattern predicted by Bergmann’s rule.
By the other hand, despite the comparative suitability of forearm length and body mass as proxies of intraspecific body size in bats is beyond the scope of this study, we provide the follow rationale to clarify our decision to use forearm length as a structural measure of body size in C. mexicanus in the present discussion.
Since size refers to the spatial extent occupied by an object in a three-dimensional space, size is fundamentally the volume, which depends on linear dimensions (height, length, and width). Although mass is often correlated with volume, it represents the amount of matter rather than the space occupied; therefore, mass and size are related but conceptually distinct properties. In bats, both body mass and forearm length have been widely used as proxies of body size at the interspecific level, as they generally scale in the same direction (Safi et al. 2013). However, at the intraspecific level, forearm length has been questioned as a proxy because it does not necessarily scale proportionally with body mass (McGuire et al. 2018; Mellado et al. 2024). We argue that dismissing forearm length on the basis that body mass represents a “truer” measure of size reflects a conceptual conflation between mass and volume. While body mass is appropriate for interspecific comparisons, its use at the intraspecific level often reflects variation in body condition or nutritional state rather than structural size.
In this context, we highlight two considerations of our study. First, the absence of support for the RAD hypothesis may be related to the use of forearm length rather than body mass as the response variable, given that resource availability is more directly associated with body condition and nutritional state. Second, because Bergmann’s rule is grounded in surface-to-volume dynamics (Bergmann 1848), we consider forearm length an appropriate structural proxy of size, even at the intraspecific level, as changes in linear dimensions necessarily affect volume. In other words, affect the amount of space that individuals occupied, which by definition is their size.
Variation of forearm. Our results indicate a sex-specific pattern of forearm length variation: in females, forearm seems to be more strongly constrained by environmental selective pressures than by the species’ evolutionary history, whereas in males, although environmental factors also play a role, evolutionary history also has an influence on forearm length. Previously, López-Cuamatzi et al. (2024) reported forearm differentiation between the SMOc and TMVB lineages in males but not in females, suggesting that the “Big Mother Hypothesis” (Stevens et al. 2013) could explain the absence of differentiation in females. According to this hypothesis, sexual dimorphism in forearm length in vespertilionid bats, characterized by a bias toward larger sizes in females, is the result of a positive selection favoring females due to the cost-benefit balance of energy expenditure during reproduction and lactation (Stevens et al. 2013). Our findings support this interpretation, indicating that environmental conditions during the reproductive period may select for specific body sizes (measured here as forearm length), leading to phenotypic convergence among females regardless of lineage.
Nevertheless, the lack of differences in forearm size among females’ lineage does not necessarily imply an absence of influence from evolutionary history. Our results for females for both resolution datasets reveal a significant interaction between lineage and environment during reproductive period: in the TMVB lineage, female forearm size is influenced by environmental variables, whereas in the SMOc lineage it is not. This could be interpreted as a reduced or negligible effect of temperature in SMOc, either due to an unknown biological factor or because temperatures in that region may be more homogeneous and span a narrower range. However, our data do not allow us to discern which of these hypotheses is more plausible, indicating that the lineage-dependent nature of this response requires further investigation.
Regarding specific environmental effects, in TMVB females, PC1 was negatively related to forearm size, exhibiting the classic Bergmann’s rule pattern. This component was mainly associated with minimum temperatures during the reproductive months, but also with maximum temperatures in March and April, the latter corresponding to the highest values of the maximum temperature set (Figure 5A). This pattern suggests that forearm size variation in females may be mainly shaped by selection for heat conservation, underscoring the role of efficient thermoregulation during gestation (March-April) (López-Wilchis 1989), when energy resources are partitioned between maternal maintenance and fetal development. But also, the high temperatures of March and April may also indicate selective pressure for heat dissipation during the early stages of gestation. Elevated thermoregulatory demands during this period –whether conserving or dissipating heat– could impose fitness costs on C. mexicanus females at TMVB.
Heat dissipation and conservation patterns in TMVB females are also reflected in the relationship between PC2 and forearm length. In localities where early spring (early gestation) and midsummer (lactation) temperatures are similar, forearm length tends to be shorter, likely due to relatively homogeneous and warmer conditions where avoiding overheating is important. Conversely, where early spring is cooler and the seasonal contrast is stronger, females exhibit longer forearms, possibly reflecting greater selective pressure for heat conservation during early gestation (Wolf et al. 2025).
In males, differentiation between lineages, also reported by López-Cuamatzi et al. (2024), underscores the relevance of evolutionary factors in forearm variation. However, our analyses indicate that environmental conditions during the reproduction (RCTR hypothesis) also play an important role. It is important to note that the hypotheses tested in this study (lineage and environment) are not mutually exclusive. Therefore, interpretation should not focus on determining which hypothesis better explains forearm size variation in males, but rather on recognizing that such variation reflects environmental effects acting in both lineages within their respective morphological ranges.
The PC1 of RCTR model was negatively related to forearm length and exhibited the classic Bergmann pattern. This component was associated with maximum temperatures during the reproductive period and with minimum temperatures in August and September, the warmest months in the analyzed series (Figure 5B). This suggests that, unlike in females, forearm variation in males may be influenced by thermal stress due to overheating under higher temperatures. Komar et al. (2022) have demonstrated that male bats tend to reduce the maturation speed and production of sperm at high temperatures, reducing the chances of mating. Since August and September coincide with key reproductive processes in males of C. mexicanus –epididymal development, recrudescence of testes, and accessory sexual glands (Léon-Galván et al. 2005)– additional energy expenditure for heat dissipation could affect sperm quality and, consequently, reproductive success.
Paltrinieri et al. (2025) previously demonstrated a sex-specific response of forearm length to environmental temperature, with females of European bats exhibiting adaptations toward heat conservation and males showing evidence of both heat conservation and dissipation. However, unlike their conclusion that both mechanisms are equally plausible for males, our results suggest a stronger selective pressure for heat dissipation.
Variation in the shape and size of the cranium and man-dible. Our results indicate that the skull and mandible are relatively conserved structures in C. mexicanus, with morphological differences so subtle that, in most analyzed modules, they did not deviate from a random pattern. López-Cuamatzi et al. (2024) previously attributed this morphological conservatism to low genetic differentiation and the recent divergence between lineages, a pattern also observed in Myotis velifer peninsularis (Nájera-Cortazar et al. 2015). Our data support this conclusion and further show that environmental variables –considered as proxies for insect abundance and diversity– play a minimal or negligible role in the morphological variation of these structures.
The only significant differences observed in some skull modules were primarily associated with mitochondrial lineage. Although López-Cuamatzi et al. (2024) previously reported differentiation in the lateral view of the cranial vault in females, our analysis unexpectedly detected this differentiation in males, while it was absent in females. This discrepancy may be explained by the analytical approach: López-Cuamatzi et al. (2024) used the individual as the sampling unit, whereas in our study we used the locality as the unit, averaging the cranial shape of multiple individuals from the same population. This procedure led to a loss of variability that our analysis could not recover. Nonetheless, we considered a locality-based approach more suitable for assessing ecological patterns at a spatial scale, as individual-level variation in our dataset could obscure the effects of environmental factors.
Similar to the lateral cranial view, the ventral view of the rostrum showed differences associated with mitochondrial lineage. Although this differentiation had not been previously reported, its subtlety and cross-validation values suggest an emerging signal that should be re-evaluated in the future with a larger sample size. While lineage appears as a factor associated with variation in both modules (cranial vault and rostrum), this does not necessarily imply causality. Therefore, if this variation reflects the species’ evolutionary history, its interpretation must be approached cautiously.
Regarding the lateral view of the rostrum, López-Cuamatzi et al. (2024) previously reported a trend toward differentiation between the SMOc and TMVB lineages, although it was not statistically significant. Our results suggest that this variation may be associated with lineage or its interaction with environmental variables; however, regression analyses conducted separately for each lineage did not recover statistical significance. This is consistent with statistical reasoning, as a term may be significant in a global model but not in separate analyses, especially when effects are weak (Pike 2011). Furthermore, this result was obtained only using environmental data at 5-arcminute resolution and was not consistent across both datasets (30-arcseconds and 5-arcminute), indicating that the finding should be interpreted cautiously, as it may be an artifact of the data.
Consideration of temporal extent and spatial resolution of the database. The performance and predictive capacity of a geospatial model are affected by three types of uncertainty across spatial and temporal dimensions: uncertainty in estimating the response variable, uncertainty in estimating predictor variables, and uncertainty arising from model formulation and structure (Svensson et al. 2013). At both dimensions, two factors are particularly relevant for predictor variable uncertainty: resolution and extent (Svensson et al. 2013).
Although high-resolution predictor variables often improve predictive performance of geospatial models (e.g., Franklin et al. 2013; Chauvier et al. 2022; Cohen and Jetz 2025; but see Lowen et al. 2016), they may not be the most biologically appropriate. Because bats exhibit considerable interspecific variation in dispersal abilities and home-range sizes (Wood et al. 2024), using environmental data at a spatial resolution finer than a species’ dispersal capacity may underestimate the influence of environmental features across its effective area of movement. Conversely, coarser resolutions, although potentially more representative of the environmental conditions across an individual’s movement area, may overestimate environmental characteristics by incorporating suboptimal habitats into the averaged values. Thus, selecting the appropriate resolution in spatial analyses constitutes a trade-off that should fundamentally depend on species-specific biological knowledge.
In our study, information provided by the 30-arcsecond resolution was particularly suitable, as the distribution of C. mexicanus in montane ecosystems exhibits climatic variation determined by altitudinal gradients (see Oke and Thompson 2015; Scherrer et al. 2019). Fine resolution allows these altitudinal gradients to be captured more accurately. However, although useful in the case of C. mexicanus, it should not be considered ideal for all species.
Results on lateral rostrum shape indicate that model term significance, particularly for lineage, depended on the resolution of environmental data used. This demonstrates that methodological decisions can influence conclusions and lead to interpretations that may or may not accurately reflect species biology. Therefore, we recommend that studies analyzing ecological phenomena using spatial data include analyses at multiple resolutions to minimize methodological bias in biological interpretations. This recommendation is especially relevant for research on highly mobile species, such as bats.
Regarding the temporal extent of environmental variables, studies testing Bergmann’s rule have commonly relied on monthly, quarterly, or annual temperature ranges to evaluate the relationship between temperature and body size (e.g., Barros et al. 2014; Castillo-Figueroa 2022; Alston et al. 2023; but see Paltrinieri et al. 2025). This approach implicitly assumes that the influence of thermal conditions on body size operates equally throughout the year. Under this assumption, thermoregulatory and energetic demands are considered temporally uniform, thereby justifying the use of annual or temporally aggregated variables. However, thermoregulatory and energetic demands are not constant over time and often peak during specific periods depending on sex and life-history stage (Gustafson 1979; Oxberry 1979; Kurta and Kunz 1987; McLean and Speakman 1999). This is especially evident during the reproductive season, when energetic requirements increase (Kurta and Kunz 1987; McLean and Speakman 1999). Males allocate energy to spermatogenesis and gonadal maturation, whereas females invest in fetal development and lactation (Kurta and Kunz 1987, León-Galván et al. 2005). If individuals must simultaneously allocate additional energy to thermoregulation during these periods, this energy is diverted from reproduction. Because thermoregulatory efficiency is influenced by body size, thermal conditions experienced during reproduction are likely to impose stronger selective pressures on body size than conditions encountered during periods of lower energetic demand.
For example, two hypothetical females of different body sizes inhabiting the same colony would experience similar overall thermal conditions. However, if temperatures during lactation are relatively cool, the larger female would conserve heat more efficiently due to a lower surface-area-to-volume ratio. In contrast, the smaller female would lose heat more rapidly and would need to allocate additional energy to thermoregulation, thereby reducing the energy available for milk production. A similar rationale applies to males: under warm autumn conditions, larger males may need to invest more energy in heat dissipation, potentially decreasing the energy available for gonadal development and sperm maturation. Because these energetic trade-offs directly affect reproductive performance, they may have fitness consequences and thus be subject to natural selection.
Therefore, we argue that the influence of environmental conditions on body size is likely strongest during reproductive periods. Assessing temperature effects at an annual scale using mean yearly values may obscure these relationships, as such measures incorporate thermal conditions from periods of lower energetic demand. Consequently, this approach may introduce thermal noise and mask biologically meaningful patterns. Our results indicate that models based on environmental variables restricted to the reproductive period had greater explanatory power than those using annual variables, supporting the idea that the strength of environmental influence on body size is not uniform throughout the year. In this context, our study provides a framework for future research exploring this pattern in bats more broadly, particularly among Holarctic species that exhibit asynchronous reproductive energy investment.
Conclusions
Our study examined phenotypic variation in forearm length as a proxy of body size, as well as cranial and mandibular morphology, in C. mexicanus, exploring their associations with evolutionary history and environmental factors. Phenotypic variation in C. mexicanus displays sex-biased patterns, with distinct factors driving variation in each sex. Forearm length aligns with the LD and RCTR hypotheses in males and the LD*RCTR hypothesis in females. However, the temperatures associated with forearm variation in the RCTR hypothesis suggest different underlying mechanisms such as conservation in females and heat dissipation in males. In contrast to forearm length, cranial and mandibular shape and size exhibit morphological conservatism, with only slight differences in certain views among lineages, consistent with the LD hypothesis.
Our findings also contribute to ongoing discussions about the appropriate temporal scale of temperature data used to test Bergmann’s rule and geographic patterns of body size variation, particularly at the intraspecific level. We propose that, within species, environmental thermal conditions exert their strongest influence on body size during the reproductive season, when energetic investment is highest. During this period, physiological pressures related to heat conservation or dissipation may intensify, allowing environmental conditions to play a critical role due to their direct effects on individual fitness.
Nevertheless, these conclusions require further evaluation through additional studies to determine whether this pattern is universal or broadly applicable. Moreover, the spatial resolution of environmental data should be carefully considered, as uncertainty in predictor estimates may influence model performance and, consequently, the interpretation of results.
Acknowledgments
All the authors wish to honor, through this text, the tireless work of Dr. Livia León Paniagua. Each of us has had the privilege of collaborating closely with her on projects, shared advisorships, and conjoined publications, and on every occasion, we have witnessed her exceptional qualities as both an academic and a human being, as well as her ability to inspire students and young researchers close to her.
We are grateful to the following curators and museums for providing access to their collections: Dr. Celia López González from the Colección de Mamíferos, CIIDIR-Durango, Instituto Politécnico Nacional; Dr. Livia León Paniagua from the Colección de Mamíferos, Museo de Zoología “Alfonso L. Herrera”, FC-UNAM; Dr. Fernando A. Cervantes and Yolanda Hortelano-Moncada from the Colección Nacional de Mamíferos, IB-UNAM; Dr. Cynthia Elizalde Arellano from the Colección Mastozoológica, ENCB-IPN; Dr. Rogelio Rosas Valdez from the Laboratorio de Colecciones Biológicas y Sistemática Molecular de la UACB, UAZ; Dr. Santiago Espinosa Andrade from the Colección de Vertebrados, Facultad de Ciencias-UASLP; Dr. Christian A. Delfín and MSc. Nallely Verónica Rodríguez Santiago from the Colección de Mamíferos, IIB-UV; Dr. José Ramírez Pulido and Noé González Ruiz from the Colección de Mamíferos de la Sierra Volcánica Transversal de México, UAM-Iztapalapa, as well as the Colección de Mamíferos de UAEH. We are also grateful for the assistance of Dr. Pedro A. Aguilar-Rodríguez, MSc. Melany Aguilar López and Dr. Melina del Real Monroy during specimens’ revisions. The first author gratefully acknowledges the Secretaría de Ciencia, Humanidades, Tecnología e Innovación (SECIHTI) for the Master’s and SNI-III Research Assistant fellowships (No. 1018336), and the Centro de Investigaciones Tropicales, UV, for support through a fieldwork grant.
Declaration of Artificial Intelligence use
The authors declare that AI was used to assist with style suggestions, textual clarity, and revision of grammar and syntax in English.
Author contributions
Issachar L. López-Cuamatzi: Conceptualization, Formal analysis, Investigation, Methodology, Visualization, Writing – Original Draft Preparation. Jorge Ortega: Conceptualization, Funding Acquisition, Methodology, Project administration, Resources, Writing- Reviewing and Editing. Sandra M. Ospina-Garcés: Conceptualization, Methodology, Writing- Reviewing and Editing. M. Cristina MacSwiney G.: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing- Reviewing and Editing.
Supplementary data
SD1. Results of sexual dimorphism analysis for forearm length and cranial and mandibular shapes.
Adams D, Collyer M, Kaliontzopoulou A, and Baken E. 2025. “Geomorph: Software for geometric morphometric analyses. R-package version 4.0.10.” https://cran.r–project.org/R-package=geomorph.
Alston JM, Keinath DA, Willis CK, Lausen CL, O’Keefe JM, Tyburec JD, et al. 2023. Environmental drivers of body size in North American bats. Functional Ecology ٣٧:١٠٢٠–١٠٣٢. https://doi.org/10.1111/1365-2435.14287
Ashton KG. 2002. Patterns of within‐species body size variation of birds: strong evidence for Bergmann’s rule. Global Ecology and Biogeography 11:505–523. https://doi.org/10.1046/j.1466-822X.2002.00313.x
Ashton KG, Tracy MC, and Queiroz AD. 2000. Is Bergmann’s rule valid for mammals? The American Naturalist ١٥٦:٣٩٠–٤١٥. https://doi.org/10.1086/303400
Barros LAVD, Fortes RDR, and Lorini ML. 2014. The Application of Bergmann´ s Rule to Carollia perspicillata (Mammalia, Chiroptera). Chiroptera Neotropical 20:1243–1251.
Bates D, Mächler M, Bolker B, and Walker S. 2015. Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software 67:1–48. https://doi.org/10.48550/arXiv.1406.5823
Bergmann C. 1848. Über die Verhältnisse der Wärmeöko-nomie der Thiere zu ihrer Grösse. Göttingen (GER): Vandenhoeck und Ruprecht.
Bolnick DI, Amarasekare P, Araújo MS, Bürger R, Levine JM, Novak M, et al. 2011. Why intraspecific trait variation matters in community ecology. Trends in Ecology and Evolution 26:183–192. https://doi.org/10.1016/j.tree.2011.01.009
Bolnick DI, Svanbäck R, Fordyce JA, Yang LH, Davis JM, Hulsey CD, et al. 2003. The ecology of individuals: incidence and implications of individual specialization. The American Naturalist ١٦١:١–٢٨. https://doi.org/10.1086/343878
Bookstein FL. 1997. Morphometric tools for landmark data: Geometry and biology. Cambridge (UK): Cambridge University Press.
Borer ET, Seabloom EW, and Tilman D. 2012. Plant diversity controls arthropod biomass and temporal stability. Ecology Letters 15:1457–1464. https://doi.org/10.1111/ele.12006
Bozinovic F, Bastías DA, Boher F, Clavijo-Baquet S, Estay SA, and Angilletta Jr. MJ. 2011. The mean and variance of environmental temperature interact to determine physiological tolerance and fitness. Physiological and Biochemical Zoology 846:543–552. https://doi.org/10.1086/662551
Bruyndonckx R, Hens N, and Aerts M. 2018. Simulation‐based evaluation of the linear‐mixed model in the presence of an increasing proportion of singletons. Biometrical Journal 601:49–65. https://doi.org/10.1002/bimj.201700025
Camargo NFD, and de Oliveira, HF. 2012. Sexual dimorphism in Sturnira lilium (Chiroptera, Phyllostomidae): can pregnancy and pup carrying be responsible for differences in wing shape? PLoS One ٧:e٤٩٧٣٤. https://doi.org/١٠.١٣٧١/journal.pone.٠٠٤٩٧٣٤
Castillo–Figueroa D. ٢٠٢٢. Does Bergmann’s rule apply in bats? Evidence from two neotropical species. Neotropical Biodiversity ٨:٢٠٠–٢٢١. https://doi.org/10.1080/23766808.2022.2075530
Chauvier Y, Descombes P, Guéguen M, Boulangeat L, Thuiller W, and Zimmermann NE. 2022. Resolution in species distribution models shapes spatial patterns of plant multifaceted diversity. Ecography 2022:e05973. https://doi.org/10.1111/ecog.05973
Chown SL, and Gaston KJ. 2010. Body size variation in insects: a macroecological perspective. Biological Reviews 85:139–169. https://doi.org/10.1111/j.1469-185X.2009.00097.x
Chu C, Bartlett M, Wang Y, He F, Weiner J, Chave J, and Sack L. 2016. Does climate directly influence NPP globally? Global Change Biology ٢٢:١٢–٢٤. https://doi.org/10.1111/gcb.13079
Clauset A, and Erwin DH. 2008. The evolution and distribution of species body size. Science ٣٢١:٣٩٩–٤٠١. https://doi.org/10.1126/science.1157534
Cohen JM, and Jetz W. 2025. Fine‐Grain Predictions Are Key to Accurately Represent Continental‐Scale Biodiversity Patterns. Global Ecology and Biogeography 34:e13934. https://doi.org/10.1111/geb.13934
Collyer ML, and Adams DC. 2018. RRPP: An package for fitting linear models to high‐dimensional data using residual randomization. Methods in Ecology and Evolution 9:1772–1779. https://doi.org/10.1111/2041-210X.13029
Dinno A. 2009. Exploring the sensitivity of Horn’s parallel analysis to the distributional form of random data. Multivariate Behavioral Research 443: 362–388. https://doi.org/10.1080/00273170902938969
Encarnação JA, Kierdorf U, Holweg D, Jasnoch U, and Wolters V. 2005. Sex‐related differences in roost‐site selection by Daubenton’s bats Myotis daubentonii during the nursery period. Mammal Review 35:285–294. https://doi.org/10.1111/j.1365-2907.2005.00066.x
Fellers GM, and Pierson ED. 2002. Habitat use and foraging behavior of Townsend’s big–eared bat (Corynorhinus townsendii) in coastal California. Journal of Mammalogy 83:167–177. https://doi.org/10.1644/1545-1542(2002)083%3C0167:HUAFBO%3E2.0.CO;2
Fick SE, and Hijmans RJ. 2017. WorldClim 2: new 1km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37:4302–4315. https://doi.org/10.1002/joc.5086
Franklin J, Davis FW, Ikegami M, Syphard AD, Flint LE, Flint AL, et al. 2013. Modeling plant species distributions under future climates: how fine scale do climate projections need to be? Global Change Biology 19:473–483. https://doi.org/10.1111/gcb.12051
Frick WF, Reynolds DS, and Kunz TH. 2010. Influence of climate and reproductive timing on demography of little brown myotis Myotis lucifugus. Journal of Animal Ecology 79:128–136. https://doi.org/10.1111/j.1365-2656.2009.01615.x
Gibert JP. 2016. The effect of phenotypic variation on metapopulation persistence. Population Ecology ٥٨:٣٤٥–٣٥٥. https://doi.org/10.1007/s10144-016-0548-z
Gorobeyko UV, Kazakov DV, Kadetova AA, Sheremetyeva IN, Guskov VY, Kartavtseva IV, et al. 2025. Intraspecific structure of Myotis petax Hollister, 1912 (Chiroptera: Vespertilionidae) based on mitochondrial DNA and morphological data. Vertebrate Zoology 75:87–106. https://doi.org/10.3897/vz.75.e134683
Gotelli NJ, and Graves GR. 1996. Null models in ecology. Washington D.C. (EEUU): Smithsonian Institution Press.
Gunz P, and Mitteroecker P. 2013. Semi landmarks: a method for quantifying curves and surfaces. Hystrix 24:103–109. https://doi.org/10.4404/hystrix–24.1–6292
Gustafson AW. 1979. Male reproductive patterns in hibernating bats. Reproduction ٥٦:٣١٧–٣٣١. https://doi.org/10.1530/jrf.0.0560317
Hedrick BP. 2021. Inter–and intraspecific variation in the Artibeus species complex demonstrates size and shape partitioning among species. PeerJ 9:e11777. https://doi.org/10.7717/peerj.11777
Klingenberg CP. 2015. Analyzing fluctuating asymmetry with geometric morphometrics: concepts, methods, and applications. Symmetry 7:843–934. https://doi.org/10.3390/sym7020843
Komar E, Fasel NJ, Szafrańska PA, Dechmann DK, Zegarek M, and Ruczyński I. 2022. Energy allocation shifts from sperm production to self–maintenance at low temperatures in male bats. Scientific Reports ١٢:١–11. https://doi.org/10.1038/s41598-022-05896-3
Korner-Nievergelt F, Roth T, Von Felten S, Guélat J, Almasi B, and Korner-Nievergelt P. 2015. Bayesian data analysis in ecology using linear models with R BUGS and Stan. Massachusetts (EEUU): Academic Press. https://doi.org/10.1016/C2013-0-23227-X
Kurta A, and Kunz TH. 1987. Size of bats at birth and maternal investment during pregnancy. Symposia of the Zoological Society of London 57:79–106.
Lenth R, and Piaskowski J. 2025. emmeans: Estimated Marginal Means aka Least-Squares Means R package version 201. https://githubcom/rvlenth/emmeans
León-Galván MA, López-Wilchis R, Hernández-Pérez O, Arenas-Ríos E, and Rosado A. 2005. Male reproductive cycle of Mexican big–eared bats, Corynorhinus mexicanus (Chiroptera: Vespertilionidae). The Southwestern Naturalist 50:453–460.
Lind EM, Pierre KJL, Seabloom EW, Alberti J, Iribarne O, Firn J, et al. 2017. Increased grassland arthropod production with mammalian herbivory and eutrophication: A test of mediation pathways. Ecology 98:3022–3033. https://doi.org/10.1002/ecy.2029
López-Cuamatzi IL, Ortega J, Ospina-Garcés SM, Zúñiga G, and MacSwiney G MC. 2024. Molecular and morphological data suggest a new species of big-eared bat (Vespertilionidae: Corynorhinus) endemic to northeastern Mexico. PLoS One 19:e0296275. https://doi.org/10.1371/journal.pone.0296275
López–Wilchis R. 1989. Biología de Plecotus mexicanus (Chiroptera: Vespertilionidae) en el estado de Tlaxcala, México. [PhD thesis]. [Ciudad de México (MEX)]: Universidad Nacional Autónoma de México.
Lowen JB, McKindsey CW, Therriault TW, and DiBacco C. 2016. Effects of spatial resolution on predicting the distribution of aquatic invasive species in nearshore marine environments. Marine Ecology Progress Series ٥٥٦:١٧–٣٠. https://doi.org/10.3354/meps11765
Lüdecke D, Ben-Shachar MS, Patil I, Waggoner P, and Makowski D. 2021. performance: An R package for assessment comparison and testing of statistical models. Journal of Open Source Software 660: 1–8. https://doi.org/10.21105/joss.03139
Majerus ME, and Mundy NI. 2003. Mammalian melanism: natural selection in black and white. Trends in Genetics ١٩:٥٨٥–٥٨٨. https://doi.org/10.1016/j.tig.2003.09.003
Maluleke T, Jacobs DS, and Winker H. 2017. Environmental correlates of geographic divergence in a phenotypic trait: A case study using bat echolocation. Ecology and Evolution ٧:٧٣٤٧–٧٣٦١. https://doi.org/10.1002/ece3.3251
Marchan-Rivadeneira MR, Larsen PA, Phillips CJ, Strauss RE, and Baker RJ 2012. On the association between environmental gradients and skull size variation in the great fruit–eating bat, Artibeus lituratus (Chiroptera: Phyllostomidae). Biological Journal of the Linnean Society ١٠٥:٦٢٣–٦٣٤. https://doi.org/10.1111/j.1095-8312.2011.01804.x
Marroig G, Shirai LT, Porto A, de Oliveira FB, and De Conto V. 2009. The evolution of modularity in the mammalian skull II: evolutionary consequences. Evolutionary Biology 36:136–148. https://doi.org/10.1007/s11692-009-9051-1
McGuire LP, Kelly LA, Baloun DE, Boyle WA, Cheng TL, Clerc J, et al. 2018. Common condition indices are no more effective than body mass for estimating fat stores in insectivorous bats. Journal of Mammalogy 995:1065–1071. https://doi.org/10.1093/jmammal/gyy103
McLean JA, and Speakman JR. 1999. Energy budgets of lactating and non‐reproductive brown long‐eared bats (Plecotus auritus) suggest females use compensation in lactation. Functional Ecology 13:360–372. https://doi.org/10.1046/j.1365-2435.1999.00321.x
McNab BK. 2010. Geographic and temporal correlations of mammalian size reconsidered: a resource rule. Oecologia ١٦٤:١٣–٢٣. https://doi.org/10.1007/s00442-010-1621-5
Meiri S. 2008. Evolution and ecology of lizard body sizes. Global Ecology and Biogeography 17:724–734. https://doi.org/10.1111/j.1466-8238.2008.00414.x
Mellado B, de Oliveira Carneiro L, Nogueira MR, and Monteiro LR. 2024. Not all size measures are created equal: Different body size proxies are not equivalent fitness predictors in the bat Carollia perspicillata. Journal of Mammalian Evolution 311:1–12. https://doi.org/10.1007/s10914-024-09702-x
Mendes SB, Stefanello F, Costa CLDS, Lima ACDS, Olímpio APM, Pires W MDM, et al. 2024. Morphological and molecular data combined reveal inter–and intraspecific cranial shape variations in bats of Artibeus Leach, 1821 (Chiroptera: Phyllostomidae). Biological Journal of the Linnean Society 143:blae031. https://doi.org/10.1093/biolinnean/blae031
Montañez-Reyna M, León-Cortés JL, Infante F, Naranjo EJ, and Gómez-Velasco A. 2022. Diversity and climatic distribution of moths in the tribe Arctiini Lepidoptera: Erebidae: Arctiinae in Mexico. Annals of the Entomological Society of America 1153:253–266. https://doi.org/10.1093/aesa/saac002
Moosman Jr PR, Thomas HH, and Veilleux JP. 2012. Diet of the widespread insectivorous bats Eptesicus fuscus and Myotis lucifugus relative to climate and richness of bat communities. Journal of Mammalogy ٩٣:٤٩١–٤٩٦. https://doi.org/10.1007/s13364-023-00678-2
Morales AE, De la Mora M, and Piñero D. 2018. Spatial and environmental factors predict skull variation and genetic structure in the cosmopolitan bat Tadarida brasiliensis. Journal of Biogeography ٤٥:١٥٢٩–١٥٤٠. https://doi.org/10.1111/jbi.13243
Myers P. 1978. Sexual dimorphism in size of vespertilionid bats. The American Naturalist 112:701–711.
Nájera-Cortazar LA, Álvarez-Castañeda ST, and De Luna E. 2015. An analysis of Myotis peninsularis (Vespertilionidae) blending morphometric and genetic datasets. Acta Chiropterologica 17:37–47. https://doi.org/10.3161/15081109ACC2015.17.1.003
NASA Earth Observatory. 2025. Net primary productivity (1 month – Terra/MODIS). Washington D.C. (EEUU): National Aeronautics and Space Administration; July, 2025. [Accessed July 21th, 2025]. https://neo.gsfc.nasa.gov/view.php?datasetId=MOD17A3H_Y_NPP
Oke OA, and Thompson KA. 2015. Distribution models for mountain plant species: the value of elevation. Ecological Modelling ٣٠١:٧٢–٧٧. https://doi.org/10.1016/j.ecolmodel.2015.01.019
Oksanen J, Simpson G, Blanchet F, Kindt R, Legendre P, Minchin P, et al. 2025. vegan: Community Ecology R-package. R-package version 2.8-0, https://vegandevs.github.io/vegan/.
Oxberry BA. 1979. Female reproductive patterns in hibernating bats. Reproduction ٥٦:٣٥٩–٣٦٧. https://doi.org/١٠.١٥٣٠/jrf.٠.٠٥٦٠٣٥٩
Paltrinieri L, Razgour O, Santini L, Russo D, Aihartza J, Aizpurua O, et al. ٢٠٢٥. The effects of climate on bat morphology across space and time. Ecography 2025:e07663. https://doi.org/10.1002/ecog.07663
Peralta–Maraver I, and Rezende EL. 2021. Heat tolerance in ectotherms scales predictably with body size. Nature Climate Change ١١:٥٨–٦٣. https://doi.org/١٠.١٠٣٨/s٤١٥٥٨-٠٢٠-٠٠٩٣٨-y
Pérez SI, Bernal V, and González PN. ٢٠٠٦. Differences between sliding semi–landmark methods in geometric morphometrics with an application to human craniofacial and dental variation. Journal of Anatomy 208:769–784. https://doi.org/10.1111/j.1469–7580.2006.00576.x
Pike N. 2011. Using false discovery rates for multiple comparisons in ecology and evolution. Methods in Ecology and Evolution 2:278–282. https://doi.org/10.1111/j.2041-210X.2010.00061.x
Porto A, de Oliveira FB, Shirai LT, De Conto V, and Marroig G. 2009. The evolution of modularity in the mammalian skull I: morphological integration patterns and magnitudes. Evolutionary Biology 36:118–135. https://doi.org/10.1007/s11692-008-9038-3
R Core Team. 2024. R: A language and environment for statistical computing R Foundation for Statistical Computing Vienna Austria https://wwwR-projectorg/
Revelle W, and Revelle MW. 2015. Package ‘psych’. The comprehensive R archive network. https://personality-project.org/r/psych/
Rohlf FJ. 1990. Rotational fit Procrustes methods. In: Rohlf FJ, Bookstein FL, editors. Proceedings of the Michigan morphometrics workshop. Ann Arbor (EEUU): University of Michigan Museum of Zoology Special Publication 2; p. 227–236.
Rohlf FJ. 2017. TpsDig2 v.2.3.6. Department of Ecology and Evolution, State University of New York. New York.
Rohlf FJ. 2019. TpsUtil ver. 1.7.8. Department of Ecology and Evolution, State University of New York. New York.
Russo D, Jones G, Martinoli A, Preatoni DG, Spada M, Pereswiet‐Soltan A, and Cistrone L. 2024. Climate is changing, are European bats too? A multispecies analysis of trends in body size. Ecology and Evolution 14:e10872. https://doi.org/10.1002/ece3.10872
Safi K, Meiri S, Jones KE, Smith FA, and Lyons K. 2013. Evolution of body size in bats. In: Smith FA, and Lyons SK, editors. Animal body size: Linking pattern and process across space, time, and taxonomic group. Chicago (EEUU): The University of Chicago Press; p. 95–151.
Salinas‐Ramos VB, Ancillotto L, Bosso L, Sánchez‐Cordero V, and Russo, D. 2020. Interspecific competition in bats: state of knowledge and research challenges. Mammal Review ٥٠:٦٨–٨١. https://doi.org/10.1111/mam.12180
Scherrer D, Christe P, and Guisan A. 2019. Modelling bat distributions and diversity in a mountain landscape using focal predictors in ensemble of small models. Diversity and Distributions ٢٥:٧٧٠–٧٨٢. https://doi.org/10.1111/ddi.12893
Schlager S. 2017. Morpho and Rvcg–Shape Analysis in R: R-packages for geometric morphometrics shape analysis and surface manipulations. In: Zheng G, Li S, Szekely G, editors. Statistical shape and deformation analysis: methods, implementation and applications. London (UK): Academic Press p. 217–256. https://doi.org/10.1016/B978-0-12-810493-4.00011-0
Smith FA, Betancourt JL, and Brown JH. 1995. Evolution of body size in the woodrat over the past 25,000 years of climate change. Science ٢٧٠:٢٠١٢–٢٠١٤. https://doi.org/10.1126/science.270.5244.2012
Stevens RD, Johnson ME, and McCulloch ES. 2016. Geographic variation of wing morphology of great fruit–eating bats (Artibeus lituratus): environmental, genetic and spatial correlates of phenotypic differences. Biological Journal of the Linnean Society 118:734–744. https://doi.org/10.1111/bij.12787
Svensson JR, Jonsson L, and Lindegarth M. 2013. Excessive spatial resolution decreases performance of quantitative models, contrary to expectations from error analyses. Marine Ecology Progress Series ٤٨٥:٥٧–٧٣. https://doi.org/10.3354/meps10307
Ulpian CM, and Rossi MN. 2017. Intraspecific variation in body size and sexual size dimorphism, and a test of Rensch’s rule in bats. Acta Zoologica, 98:377–386. https://doi.org/10.1111/azo.12183
Uvizl M, and Benda P. 2021. Intraspecific variation of Myotis emarginatus (Chiroptera: Vespertilionidae) inferred from mitochondrial and nuclear genetic markers. Acta Chiropterologica 23:285–300. https://doi.org/10.3161/15081109ACC2021.23.2.002
Watt C, Mitchell S, and Salewski V. 2010. Bergmann’s rule; a concept cluster? Oikos ١١٩:٨٩–١٠٠. https://doi.org/10.1111/j.1600-0706.2009.17959.x
Williams DF, and Findley JS. 1979. Sexual size dimorphism in vespertilionid bats. American Midland Naturalist 113–126. https://doi.org/10.2307/2425072
Wilson RJ, Gutierrez D, Gutierrez J, and Monserrat VJ. 2007. An elevational shift in butterfly species richness and composition accompanying recent climate change. Global Change Biology 139:1873–1887. https://doi.org/10.1111/j.1365-2486.2007.01418.x
Wolf JM, Lehmann P, and Kerth G. 2025. Field respirometry in a wild maternity colony of Bechstein’s bats Myotis bechsteinii indicates high metabolic costs above but not below the thermoneutral zone. Journal of Experimental Biology. 228:JEB249975. https://doi.org/10.1242/jeb.249975
Wood MR, de Vries JL, Monadjem A, and Markotter W. 2024. Review and meta-analysis of correlates of home range size in bats. Journal of Mammalogy ١٠٥:١٠٤٤–١٠٥٦. https://doi.org/10.1093/jmammal/gyae036
Zhu D, Liu Y, Gong L, Si M, Wang Q, Feng J, and Jiang T. 2024. The consumption and diversity variation responses of agricultural pests and their dietary niche differentiation in insectivorous bats. Animals 14:1–22. https://doi.org/10.3390/ani14050815
Zimmermann NE, Yoccoz NG, Edwards Jr. TC, Meier ES, Thuiller W, Guisan A, and Pearman PB. 2009. Climatic extremes improve predictions of spatial patterns of tree species. Proceedings of the National Academy of Sciences 106:19723–19728. https://doi.org/10.1073/pnas.0901643106
Associated editors: Giovani Hernández Canchola and Pablo Colunga Salas
Submitted: October 30, 2025; Reviewed: January 8, 2026
Accepted: March 9, 2026; Published online: May 29, 2026
THERYA, 2026, Vol. 17(2):143-162
Issachar L. López-Cuamatzi1,4 , Jorge Ortega2 , Sandra M. Ospina-Garcés3 , and M. Cristina MacSwiney G.4* .
DOI: 10.12933/therya.2026.6245 ISSN 2007-3364
Figure 1. Schematic representation of the hypotheses tested for forearm, cranial, and mandibular shape variation.
Figure 2. Workflow for processing environmental data. The diagram depicts the steps from initial input data (coordinates and environmental layers) to the final datasets used for linear model fitting. During preprocessing for the RCTR hypothesis, PCA was performed separately by sex. Final datasets (Outputs) were selected based on Mantel test results.
Figure 3. (A) Differences in forearm length between lineages in male C. mexicanus. (B) Positive association between forearm length and PC2 of the GED (5-arc-minute) dataset in males; positive PC2 values indicate lower temperature contrast (lower maximum temperatures), whereas negative values indicate higher contrast (higher maximum temperatures). (C) Negative association between forearm length and PC1 of the RCTR dataset in males; positive PC1 values correspond to higher temperatures and negative values to lower temperatures. (D) In females, forearm length was negatively associated with PC1 of the GED (30-arc-second) dataset, where positive values represent higher temperatures and negative values represent lower temperatures. (E–F) In the LD+RCTR dataset, forearm length in females showed significant interactions with lineage for both PC1 and PC2. For PC1, positive values indicate higher temperatures and negative values lower temperatures; for PC2, positive values reflect greater temperature contrast between early spring and midsummer, whereas negative values indicate more similar seasonal temperatures. SMOc individuals showed no detectable relationship between temperature and forearm length (G), displaying a homogeneous response to climatic PCs (uniform color), whereas TMVB individuals exhibited forearm variation associated with both PC1 and PC2 (color gradient).
Figure 4. (A) Lateral braincase shape variation in males of C. mexicanus, amplified using “mag = 3” for clearer visualization. (B) Linear Discriminant Function (LDF) scores for males from TMVB and SMOc based on lateral braincase shape. (C) Ventral rostrum shape variation, amplified using “mag = 3”. (D) LDF scores for individuals from TMVB and SMOc based on ventral rostrum shape.
Table 2. Summary of the best Procrustes ANOVA models fitted for cranial and mandibular shape. Sex and resolution categories are indicated below each view, along with sample size (n).
|
View |
n |
Factor |
d.f. |
SS |
MS |
R2 |
F |
Z |
P |
|
Lateral braincase |
24 |
Lineage |
1 |
0.0005849 |
0.00058489 |
0.04811 |
1.3176 |
0.69787 |
0.249 |
|
Females |
PC1 |
1 |
0.0008224 |
0.00082245 |
0.06765 |
1.8528 |
1.62565 |
0.053 |
|
|
30 sec |
PC2 |
1 |
0.0008496 |
0.00084957 |
0.06988 |
1.9138 |
1.56436 |
0.058 |
|
|
Lineage*PC1 |
1 |
0.0007242 |
0.00072415 |
0.05956 |
1.6313 |
1.29356 |
0.111 |
||
|
Lineage*PC2 |
1 |
0.0005303 |
0.00053025 |
0.04362 |
1.1945 |
0.55129 |
0.295 |
||
|
Residuals |
18 |
0.0079903 |
0.0004439 |
0.67523 |
|||||
|
Total |
23 |
0.0121575 |
|||||||
|
Lateral braincase |
24 |
Lineage |
1 |
0.0008136 |
0.00081358 |
0.06692 |
1.7084 |
1.18269 |
0.119 |
|
Females |
PC1 |
1 |
0.0005085 |
0.00050846 |
0.04182 |
1.0677 |
0.36244 |
0.362 |
|
|
5 min |
PC2 |
1 |
0.0005309 |
0.00053088 |
0.04367 |
1.1148 |
0.46211 |
0.32 |
|
|
Lineage*PC1 |
1 |
0.0004954 |
0.00049538 |
0.04075 |
1.0402 |
0.27462 |
0.395 |
||
|
Lineage*PC2 |
1 |
0.0009379 |
0.00093788 |
0.07714 |
1.9694 |
1.52011 |
0.06 |
||
|
Residuals |
18 |
0.008572 |
0.00047622 |
0.70508 |
|||||
|
Total |
23 |
0.0121575 |
|||||||
|
Lateral braincase |
32 |
Lineage |
1 |
0.0010894 |
0.00108942 |
0.07129 |
2.4677 |
2.12624 |
0.017 |
|
Males |
PC1 |
1 |
0.0004505 |
0.00045047 |
0.02948 |
1.0204 |
0.12911 |
0.452 |
|
|
5 min |
PC2 |
1 |
0.000559 |
0.00055896 |
0.03658 |
1.2661 |
0.5972 |
0.28 |
|
|
Lineage*PC1 |
1 |
0.0007205 |
0.00072055 |
0.04715 |
1.6322 |
1.25127 |
0.113 |
||
|
Lineage*PC2 |
1 |
0.0010416 |
0.00104155 |
0.06816 |
2.3593 |
2.04082 |
0.024 |
||
|
Residuals |
26 |
0.0114781 |
0.00044147 |
0.75115 |
|||||
|
Total |
31 |
0.0152807 |
|||||||
|
Lateral rostrum |
47 |
Lineage |
1 |
0.001616 |
0.0016159 |
0.03872 |
1.9171 |
1.29774 |
0.102 |
|
Both sexes |
PC1 |
1 |
0.001298 |
0.0012977 |
0.0311 |
1.5395 |
0.98793 |
0.168 |
|
|
5 min |
PC2 |
1 |
0.000408 |
0.0004083 |
0.00978 |
0.4844 |
-0.78256 |
0.764 |
|
|
Lineage*PC1 |
1 |
0.002497 |
0.0024966 |
0.05983 |
2.9619 |
2.0349 |
0.023 |
||
|
Lineage*PC2 |
1 |
0.003489 |
0.0034887 |
0.0836 |
4.1389 |
2.43528 |
0.007 |
||
|
Residuals |
41 |
0.034559 |
0.0008429 |
0.82818 |
|||||
|
Total |
46 |
0.041729 |
|||||||
|
Ventral rostrum |
45 |
Lineage |
1 |
0.0014626 |
0.00146264 |
0.05208 |
2.3627 |
1.9463 |
0.022 |
|
Both sexes |
Residuals |
43 |
0.0266191 |
0.00061905 |
0.94792 |
||||
|
30 sec |
Total |
44 |
0.0280817 |
Figure 5. Trends of maximum and minimum temperatures by sex and lineage. The reproductive period (colored rectangle) used to test the RCTR hypothesis is shown for females (A) and males (B), with mean values (circles) and standard deviations (intervals) calculated from historical monthly records spanning 1950-2024.