Heat over heritability: Increasing body size in response to global warming is not stabilized by genetic effects in Bechstein's bats

How well populations can cope with global warming will often depend on the evolutionary potential and plasticity of their temperature‐sensitive, fitness‐relevant traits. In Bechstein's bats (Myotis bechsteinii), body size has increased over the last decades in response to warmer summers. If this trend continues it may threaten populations as larger females exhibit higher mortality. To assess the evolutionary potential of body size, we applied a Bayesian ‘animal model’ to estimate additive genetic variance, heritability and evolvability of body size, based on a 25‐year pedigree of 332 wild females. Both heritability and additive genetic variance were reduced in hot summers compared to average and cold summers, while evolvability of body size was generally low. This suggests that the observed increase in body size was mostly driven by phenotypic plasticity. Thus, if warm summers continue to become more frequent, body size likely increases further and the resulting fitness loss could threaten populations.

higher adult mortality (Fleischer et al., 2017). It has been previously shown that during the last 25 years, female Bechstein's bats born in warmer summers grow to larger sizes (Mundinger et al., 2021).
During summer, female Bechstein's bats form maternity colonies to raise a single offspring per mother. In contrast to normothermic mammals that maintain a constant body temperature, bats of the temperate zone can lower their body temperature and enter daily torpor during adverse conditions (Willis & Brigham, 2005). This state of reduced physiological activity helps bats to save energy, but also inhibits growth in juveniles (Racey, 1973;Racey & Swift, 1981). The presumed attenuation of such daily torpor bouts in warmer summers likely explains why juvenile females grow larger in these years (Mundinger et al., 2021).
Taken together, the higher mortality of larger individuals and the increase in body size in warmer summers raise the concern that global warming may further threaten this species, which is already of high conservation concern (Petrov et al., 2018). However, energetic constraints imposed by thermoregulatory and flight costs might restrict morphological variation and evolvability of body size in bats (compare Cava et al., 2019;Rubalcaba et al., 2022).
Here, we study phenotypic plasticity, additive genetic variance (V A ), narrow-sense heritability (h 2 ), and the evolvability (I A ) of body size (measured as forearm length, FAL) in female Bechstein's bats, based on individualized morphological, demographic and genetic data from a 25-year monitoring of four wild Bechstein's bat colonies.
Estimating the heritability of body size is an important step to assess the future of body size changes and thus predict population persistence: if body size is highly heritable, selection could ultimately slow down the temperature-induced shift towards larger bats.
Using the framework of an 'animal model' (Kruuk, 2004;Wilson et al., 2010), we first address the question: how does additive genetic variance, narrow-sense heritability and evolvability of body size vary in relation to environmental conditions in Bechstein's bats? Specifically, we subdivided the data into quartiles based on the observed variation in ambient temperature during the critical time period for juvenile growth (Mundinger et al., 2021), to investigate whether V A , h 2 and I A of body size are reduced in more extreme environmental conditions. In general, morphological traits, such as body size, show higher heritability than other fitness-relevant traits (Hoffmann & Merilä, 1999;Teplitsky et al., 2009;Visscher et al., 2008), suggesting that developmental processes are robust to environmental challenges. However, during cold summers the frequent use of torpor, potentially coupled with lower insect prey availability (e.g. Welti et al., 2021), could prevent the genetic potential of individuals from being fully realized, thereby reducing the heritability of body size. In contrast, in warm summers, environmental conditions will cause juveniles to grow larger (Mundinger et al., 2021), again leading to lower heritability estimates, compared to average summers where juveniles show more average body sizes. Therefore, we expect h 2 and V A of body size to decrease in both the coldest and warmest years. In a second analysis, we estimated V A , h 2 and I A of body size in relation to the difference in temperature conditions experienced during the juvenile growth period of mothers and their daughters. Here, we expect heritability to decrease whenever the birth environments experienced by mothers and respective daughter differ strongly, as phenotypic differences should be more pronounced when growing up under very different temperature regimes. Due to the morphological constraints imposed by the capacity of flight, such as limits on the relationship of body mass to wing surface area (compare Cava et al., 2019;Rubalcaba et al., 2022), we additionally expect the evolvability of body size to be low.

| Study site, data collection on body sizes
We investigated a 25-year dataset (1996-2020) of individually RFIDtagged wild female Bechstein's bats living in four colonies (BS, GB2, HB and UA) within 15 km of one another that inhabit deciduous forests near the city of Würzburg, Bavaria, Germany (Kerth, 2022;Kerth et al., 2002). Although the composition of tree species and forest management is similar between sites, we added birth year as a random factor in all models and tested the effect of colony ID to account for any remaining spatio-temporal variance.
At least twice a year the bats were captured from their day roosts (bat boxes) to take biometric measurements, and a 3 mm wing tissue sample was taken upon first capture for genetic analyses (Kerth et al., 2000). During these capture events, unmarked females were marked with subcutaneously implanted RFID-tags (Trovan, Germany) (Kerth et al., 2011). FAL was measured to the nearest 0.1 mm using callipers. FAL is a reliable proxy of body size in bats as it highly correlates with total body length (R 2 = .93, Meng et al., 2016), but can be measured more precisely in the hand. To avoid including FAL measurements of individuals that were not fully grown, we only included measurements of females that recruited into the population in spring after their first hibernation. Males were not considered in this study, as they disperse from their natal colony (Kerth & Morf, 2004) and thus could not be sampled as adults.

| Genetic analysis and pedigree construction
After DNA extraction (Kerth et al., 2000), we genotyped each sample using 13 polymorphic microsatellite loci. Complete marker information can be found in Mundinger et al. (2022). Pedigrees were constructed as described in Mundinger et al. (2022) using the softwares coancestry (1.0.1.9) (Wang, 2011) and cervus (3.0.7;Kalinowski et al., 2007). For the heritability analysis, we only included motherdaughter pairs whose LOD scores (the natural logarithm of the likelihood-odds ratio; Kalinowski et al., 2007) were significant at the strict 95% confidence level and/or had 0 mismatches, yielding 570 mother-daughter pairs. Of these, adult body sizes for both individuals were available for 332 mother-daughter pairs, comprising 192 individual mothers across 25 cohorts. We applied the R package 'pedantics'  to search for errors (like missing mothers) in the pedigree and sort the birth chronology, so mothers were always listed before their offspring.

| Temperature conditions and changes in body size
Ambient temperature data were provided by the Bayerische ID 0570558). Based on previous analysis, we used the mean of minimum daily temperature for the period from 22 June to 16 July ( Figure 1a) as a proxy for environmental conditions, as this was estimated to be the sensitive time window during which body size is most susceptible to variation in minimum temperature (Mundinger et al., 2021). This sensitive time window is highly important, as warm temperatures during this period lead to larger adult body sizes, which has significant fitness-relevant consequences for reproduction and longevity in later life (Mundinger et al., 2022). We used this mean to split the 25 cohorts into four quartiles, ranging from cold to warm ('environment', Q1-Q4; Figure 1c) and tested for differences between all quartiles. These four environmental conditions were later used as a factor in the first part of the heritability analysis, where we assessed h 2 for the environments experienced by the juveniles during the critical growth phase.
To further evaluate whether an effect of summer temperature on maternal body size affects the heritability estimate for daughters born in differing environmental conditions, we calculated the size differences in FAL for each mother-daughter pair as well as the differences between the minimum temperatures in the respective birth years of mothers and their respective daughters. We then subdivided all mother-daughter pairs into three equally sized groups ('birth environment differentials'): (1) daughters born in colder conditions (on average at least −0.6°C colder) than the respective mothers, (2) daughters born in similar conditions (an average differences between > −0.6 and <0.9°C) and, finally, and (3) daughters born in warmer conditions (on average at least 0.9°C warmer) than their mothers ( Figure S1). As for some mothers the birth year was unknown, data for this part of the analysis were only available for 261 mother-daughter pairs. As in the first analysis, these three groups were then used as a factor to estimate h 2 .
F I G U R E 1 Boxplots in the upper row depict the variation over the years 1996-2020 of minimum temperature (°C) during the period of juvenile growth in summer (a) and forearm length (mm) in adult female Bechstein's bats (b). Red and blue lines depict the overall trend of temperature (a) and body size (b) as a smooth function. Figures are modified from Mundinger et al. (2021), with a slightly different data set. Colours represent the environmental conditions, ranging from cold years (Q1; light orange) to hot years (Q4; dark orange). Boxplots depict the average minimum summer temperature in the four environmental conditions Q1-Q4 (c) and the forearm length (mm) in adult female Bechstein's bats born in the respective environmental conditions (d). Count indicates the number of overlapping data points.

| Heritability analysis using the 'animal model'
All analyses were performed in R (R Core Development Team, 2021).
We used the 'MCMCglmm' package (Hadfield, 2010) to fit univariate animal models with body size as the phenotypic response variable.
As body size was normally distributed (Shapiro-Wilk test p = .12), we used a Gaussian distribution in all models. Animal ID, mother ID and birth year of the daughter were included in all models as random effects.
A preliminary analysis revealed that adding the variables 'colony size' and 'colony ID' as additional fixed effects did not improve the model fit (Tables S1 and S2). Consequently, these two parameters were not included in the final analysis.
In the final analysis we ran two sets of models that had different fixed effect structures. The first set of models with 'environment' as fixed effects was designed to address environmental variation, by adding temperature quartile (Q1-Q4; samples sizes ranging from N = 66-101 mother-offspring pairs) as a fixed effect. The second set of models included as fixed factor the difference between mother's and daughter's birth environment ('birth environment differentials') categorized into three categories ('daughter warmer', 'similar' or 'daughter colder'; samples sizes for each level N = 87 motheroffspring pairs).
Our baseline model had the following structure: body size ~ Var1 + random(animal) + random(mother ID) + random(birth year), with body size modelled as Gaussian and Var1 = either environment From this baseline, we built and compared models with different structures (Table 1) using DIC. To quantify additive genetic variance and narrow-sense heritability across all environments and years in Bechstein's bats, we estimated V A , h 2 and I A of body size from model M1. To estimate how genetic variance and narrow-sense heritability vary among environments and among birth environment differentials we used model structure M2. Using DIC, we then compared models M3 and M5 to test for different V A between environments, and models M3 and M4 for the same question, but with the addition of a whole covariance matrix. Comparing models M4 and M5 helped us to assess genetic covariance between environments. While model M4 ran well for the birth environment differential with three categories, data were not sufficient for estimating the environment with four categories. Therefore, Model M4 was only estimated for the birth environment differential. All summary statistics are found in Notes S1-S9.
We built an inverse relatedness matrix with the 'inverseA' command from the MCMCglmm package to consider every relationship link in a pedigree. To deal with potential autocorrelation, we applied an expanded, uniform prior defined as an inverse-Gamma distribution, as is widely used in animal models (de Villemereuil, 2018), with V = 1, nu = 0.002, alpha.mu = 0 and alpha.V = 1000. Burn-in and thinning intervals were set to yield effective sample sizes >1000 and satisfy convergence criteria. This was achieved for the models without interaction at 700,000 iterations, with a burn-in period of 10,000 and a thinning interval of 500, and for the models with an environment interaction with the same burn-in period, 1,510,000 iterations and thinning of 1500.
Model convergence was further checked using the heidel.diag and gelman.diag functions.
As pointed out by Wilson (2008) andde Villemereuil et al. (2018), including fixed effects into heritability analyses will deprive the phenotypic variance of the portion of variance that would be explained by the fixed effects of the model. This can lead to an upward bias in h 2 , so we corrected for the variance of TA B L E 1 Overview of models compared with details on the fixed effect structure, random effect structure and residual variance structure. the fixed effect and always added V F . Analogous to h 2 , the other random effect sizes were calculated as their ratio to the V Phen , with m 2 as the ratio of the maternal effect, y 2 for the ratio of effect size of birth year of cohorts and e 2 as the ratio of the residual variance.

Narrow-sense heritability
In a final step, we calculated evolvability (I A ) by dividing V A by the mean square of the phenotype (i.e. FAL) multiplied by 100 (Houle, 1992;Postma, 2014).

| Body size and environment
Minimum temperatures during the sensitive time window in June-July ( Figure 1a) and mean FAL of the different female cohorts  (Table S2).
For the 332 female offspring assigned in the pedigree, adult FAL ranged from 38.9 to 45.9 mm with a mean of 42.7 mm (±1.2 SD).
Mean offspring FAL differed across the four quartiles defined by summer temperatures (Figure 1c,d) with successively larger mean offspring in Q3 and Q4 (Table S3), but notably SD did not differ (Table 1). Mean mother FAL were similar across all environmental quartiles (Table 1).

| Change in body size between mothers and daughters
When examining the size difference between daughters and their respective mothers in relation to the mother's size (Figure 2a) Figure S1). On average, daughters born in summers colder than their mothers' birth environment reached a smaller size than their respective mothers ( Figure 2b), whereas daughters born during warmer summers than their mothers, grew respectively larger (linear model for FAL offspring y = 0.07 + 0.28 × ΔTº between mother and offspring, p < .001, R 2 = .15). When birth environments of mothers and daughters had a similar average temperature, no directional change in FAL was detected, but variance was high.

F I G U R E 2
Changes in body size between mother and offspring as a function of body size of the mother (a) and mean summer temperature (b) Left: linear regressions for the different environmental conditions (orange hues).
Genetic and nongenetic maternal effects were overall low, accounting for only 1%-2% of the total variance. Birth year of the cohorts explained 10% (Crl: 0.02-0.19) of the variance of h 2 after the summer temperature was controlled for as a fixed factor. The remaining residual variance accounted for 26% of the observed variance.
Across the four environmental quartiles, we found that h 2 estimates for the three coldest quartiles (Q1, Q2 and Q3) were similar (Figure 3a), ranging from an average of 0.48 to 0.56 (Table 2A), while h 2 during the hottest summers (Q4) was markedly lower (0.25). This reduction was caused by both a reduction in additive genetic variance (V A ), and an increase in residual variance (e 2 ). The distribution of h 2 estimates (Figure 3a) was substantially broader in both the coldest and hottest quartiles, compared to the two middle quartiles.
Comparing h 2 and V A among birth environment differentials (Figure 3b and Table 2B), revealed a similar pattern: h 2 was substantially lower when daughters grew up in comparatively warmer summers with a narrow distribution of posterior estimates. In contrast, when daughters experienced colder environments during the critical growth phase, h 2 was above average. Under similar birth conditions, h 2 estimates were broadly distributed. TA B L E 2 Means, variances and estimated random effect sizes for animal models of adult daughter and mother female forearm length FAL (mm) for (A) the four quartiles based on the summer temperature profiles of each year (Env = Environment) and (B) the difference in summer temperature experienced between mother and respective daughter (Birth Diff = Birth environment differentials). the 95% posterior credible interval. The overall value was estimated using the baseline models with model structure M1, the differences for either the four-level factor environment (A) or the three-level factor of birth environment differentials (B) were estimated using model structure M2.

F I G U R E 3
Boxplots show the differences in heritability estimates (a) during four different environmental conditions (Q1-Q4, Q1 as light orange and Q4 as dark orange), and (b) in birth environment comparison between mothers and respective daughters, as proportion of the phenotypic variance (blue = daughter born in colder conditions; white = daughter born in similar conditions; red = daughter born in hotter conditions). Violin plots depict the distribution of h 2 estimates as density curves. The dark grey area depicts the median of the overall heritability (solid black line) with the respective SDs (dotted lines), the light grey area the 95% posterior intervals.

Evolvability of FAL was low in all environmental conditions
(I A = 0.04 ± 0.02), suggesting that the evolutionary potential of body size in the study population is low.
Comparing the models M3 and M5 we found different V A , both for the four environments as well as in the birth environment differentials (Table 3 and Notes S3 and S4; S7 and S9). Comparing models M4 and M5 in the birth environment differentials showed that covariance improved model fit as indicated by the lowest DIC.

| DISCUSS ION
We find that heritability of body size overall was moderate, but clearly affected by environmental conditions in a 25-year dataset of wild Bechstein's bats. Heritability (h 2 ) and additive genetic variance (V A ) were substantially lower in the hottest summers (Q4), and when daughters were born in a warmer environment than their mother ('daughter warmer'). In hot summers, daughters tended to grow disproportionally larger, likely as a consequence of several factors: changes in bat behaviour (i.e. reduced use of torpor), environmental effects (e.g. high food availability), and lower physiological costs (i.e. ambient temperatures are closer to the bats' thermoneutral zone).
At the other end of the spectrum, in cold summers (Q1), average h 2 and additive genetic variance (V A ) did not differ with respect to their mean from the values in more average environments (Q2 and Q3), and h 2 estimates were highest when daughters were born in colder environments than their mothers. This contradicts our expectation of lower body size heritability in cold conditions when the use of torpor should have increased phenotypic variation. However, flexible social thermoregulation is a major driver of sociality in this and other temperate-zone bat species, and might alleviate the effect of cold spells (e.g. Kerth & van Schaik, 2012). Coupled with a preference of lactating females for boxes with a warm microclimate (Kerth et al., 2001), this suggests that Bechstein's bats are well adapted to cope with colder conditions during the critical phase of juvenile growth. Nevertheless, the density distribution of the h 2 estimates in Q1 was far broader, with a considerable proportion of estimates with low values. This suggests that for some of the individuals born in particularly cold years, frequent torpor use (Willis et al., 2006) or low food availability (Welti et al., 2021) substantially curtails juvenile growth.
In summary, while our expectation of low heritability and additive genetic variance in response to the hottest temperatures was met, our predictions for cold summers were refuted. In this context, similarly to our results, body size in house mice was also found to show very little plasticity in response to cold environments (Ballinger & Nachman, 2022). Regarding hot environments, Chapuis et al. (2021) Median h 2 ± SD Mean y 2 ± SD Mean m 2 ± SD Mean e 2 ± SD Mean I A ± SD reported that in desert locusts, additive genetic variance in traits closely related to fitness was largely unaffected by thermal stress. In addition, a study of the effects of thermal stress on genetic variance in six species found that thermal stress generally did not affect genetic variance, but in the case that it did, the effect was highly species specific (Fischer et al., 2021). Similarly, also Rowiński and Rogell (2017) did not find a clear direction in how heritability changes across benign or stressful conditions in their meta-analysis. Since thermal stress can lead to fundamentally different responses, predicting responses to temperature change is nearly impossible and must be evaluated in a species-specific context. Furthermore, our results indicate that the observed increase in body size in Bechstein's bats during the last decades is primarily a consequence of a high phenotypic plasticity and a particularly low heritability in warm summers.
Heritability in the heritability estimate, as these are known to reduce the comparability of heritability estimates across studies Kruuk et al., 2008;Wilson et al., 2010). Thus, we conclude that the increased plasticity of offspring growth during hot summers, is the factor driving the low heritability of body size in Bechstein's bats. Indeed, when daughters are born in warmer environments than their mothers, h 2 is particularly low. This notion is further supported by the strong influence of summer temperature on FAL found in Mundinger et al. (2021).  (Kerth & König, 1999). Group size has a strong effect on social thermoregulation, and thus can be used to buffer changes in ambient temperatures (Pretzlaff et al., 2010).

Evolvability of body size in
Moreover, roost switching allows bats to select warmer roosts during cold temperature periods, or vice versa, thus enabling bats to buffer the effect of ambient temperatures (Alston et al., 2022;Bohnenstengel, 2012;Kerth et al., 2001). Our results suggest that this behavioural flexibility in grouping and roost selection seems to be particularly well developed as an adaptation to cold temperatures, as h 2 remained comparable high in colder and average years.
In contrast, similar behavioural adaptations to warm temperatures are either not present in this population, or ineffective. Finally, prey availability, which also plays a substantial role in shaping juvenile development, and hence adult body size, is most likely affected by additional environmental conditions and not only by the ambient temperature in the time window used in our models.
With global warming increasing the frequency of warm summers (IPCC, 2021), we expect a further increase in bat body size. This is reinforced by that fact that behavioural adaptations of bats are mostly targeted to mitigate the negative effect of cold temperatures, but are inconsequential in hot environments. While this finding raises concerns in light of larger females exhibiting higher adult mortality, increasing FAL at the same time also leads to a faster pace of life in Bechstein's bat (Mundinger et al., 2022). Under current conditions, the faster rate of reproduction of larger bats can somewhat compensate for the observed mortality increase in our study populations.
However, there is a discernible trend towards reduced lifetime reproductive success in the largest individuals (Mundinger et al., 2022).
Our new findings of very low h 2 under warmer conditions suggest that the observed trend towards larger females will be more pronounced in future years. This is likely to endanger Bechstein's bat populations as increased stochasticity of extreme years, or low habitat quality, will most likely predominantly affect large individuals with faster life histories, as their lifetime reproductive success will be more strongly affected by foregoing reproduction in a poor year.
Moreover, with maximally one offspring per year, female Bechstein's bats are strongly limited in how much they can accelerate their pace of life. Under the current climate change scenarios, this gives rise to concern about future population persistence, as increased mortality in large bats will at some point likely not be compensated for by positive fitness effects associated with larger body sizes.

ACK N O WLE D G E M ENTS
We thank numerous helpers in the field and the 'Bayerische Landesanstalt Open Access funding enabled and organized by Projekt DEAL.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no potential conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available