The bigger they are, the higher they go: Australian insectivorous bats confirm Bergmann’s 175-year-old prediction

Introduction

Insectivorous bats are living at the energetic edge of their physiological limits when considering energy requirements for flight, which is exacerbated in temperate areas and higher elevations. Although bats have a range of behavioural and physiological adaptations to cold conditions, they also show geographic variations that are expressed in a species’ morphometric traits and in distributions (Voigt et al. 2010; Stawski et al. 2014; Chenery et al. 2022). Given the vagaries of an insectivorous diet in areas where seasons influence food shortages, there is a possibility that insectivorous bat size distribution is complying with Bergmann’s 175-year-old prediction. Bergmann identified bats as a mammal taxon with very small forms that have a large surface area to body volume ratio, which is extreme when compared with other taxa (loosely translated from ‘… Ein Extrem von Ausbreitung der Körperoberfläche im Verhältniss zum Volumen bietet auf der andern Seite, eine Familie dar, welche meist kleine, zum Theil sehr kleine Formen zählt: die der Fledermäuse…’; Bergmann 1848, p. 44).

Bergmann made this prediction on the basis of what is in modern times known as Bergmann’s rule, which implies that larger animals of a species and larger species occur in colder environments as a result of improved energy conservation related to reduced surface to volume ratio in larger endotherm animals. The rule is often discussed in relation to its general applicability, although Bergmann identified many exceptions and variations associated with taxonomic differences (Bergmann 1848; Brown and Lee 1969; Meiri and Dayan 2003; Clauss et al. 2013; Alston et al. 2023; Baldwin et al. 2023; He et al. 2023). The controversy around the applicability of Bergmann’s rule is in part related to the use of a single English translation of the rule, which may have led to misinterpretation and misconceptions. When taking into account Bergmann’s extensive treatment of how the body surface area to volume ratio influences the temperature–size relationship (i.e. the mechanism behind the rule), many analyses with outcomes diverting from the rule become less controversial (Watt et al. 2010). So as to avoid misinterpretation, I provide the relevant translations and original for those passages from Bergmann that I believe are relevant for this paper.

Most studies investigating the relationship between many bat species and temperature have focused on mean body size and latitude associated mean temperatures (see e.g. Clauss et al. 2013; Jiang et al. 2019; Wang et al. 2020; Alston et al. 2023) and only one study investigated the body size of two bat species of two different genera with altitude (Castillo-Figueroa 2022). Bergmann also stated that closely related, sympatric species have morphological characteristics and distributions that show these size differences at species and intraspecies level. One such taxon is the genus Vespadulus (Vespertillionidae), where three tree-roosting species occur together in south-eastern Australia over a range of elevations. Vespadelus darlingtoni is the larger species, V. vulturnus the smaller, with V. regulus being the intermediate species. Studies have indicated that V. darlingtoni prefers elevations above 600 m asl, V. vulturnus prefers elevations below 1000 m, whereas there is in no indication for elevational preference for V. regulus (Law et al. 1998; Law et al. 2023). They differ in size and weight, with forearm length being a characteristic identifier and measure of size (Herr 1998; Churchill 2009; Meng et al. 2016; Law et al. 2022).

Previous work in the 1990s on the elevational distribution and size of these species was inconclusive, and this may have been related to the limited availability of highly intensive computational methods (see Herr 1998). An initial explorative linear regression analysis of the forearm length with elevation for this work showed small inconclusive trends for the three species and sex and also a high unexplained variance expressed in R² values (Fig. 1, first panel). However, forearm length counts indicated a trend for increasing forearm length in both adult sexes of the different species (Fig. 1, second panel).

Fig. 1.

Vespadelus forearm length with elevation. AF, adult female; AM, adult male.

In this work, I investigated whether bat body size (as inferred from forearm length) increases with elevation-related temperature decrease. It tests Bergmann’s prediction that temperature-related increase in body size applies to closely related bat taxa (Bergmann 1848). Mainly, I quantified the relationship between forearm length in V. darlingtoni, V. vulturnus, and V. regulus and the elevation (and, by implication, decreasing temperature) at which they live in a statistical sense, and also incorporated additional sources of variation for bat size and bat occurrence. The purpose of this was two-fold. First, I investigated whether there is a relationship between bat size and elevation by using more recent computational techniques of Bayesian multilevel modelling. Second, I provided an example of a modern Bayesian multilevel model (BMM). Advantages of BMMs are their ability to incorporate additional information such as measurement errors and correlations to better estimate the size–elevation relationship and to estimate the full range of uncertainties, thus providing uncertainty bounds (and precision estimates) for the magnitude of this relationship. The results of this analysis also have implications for climate change-related habitat changes. I briefly discuss the implications of increasing temperature for larger bats in the study area, by using an example of elevation-related potential habitat-suitability changes.

Materials and methods

For this study, I used the data from V. darlingtoni, V. regulus and V. vulturnus adult males and females. These bat measurements were from Herr (1998), who surveyed bats in the western slopes of the Australian Alps by using harp traps between the Murray River near Albury and the mountainous south-eastern, below-alpine areas towards the Alpine National Park. For this work, I used measurements from 775 bats that occured in eight areas of inland south-eastern Australia ranging in elevation from 100 m to 1200 m asl and covering nearly 7000 km². Further details on the study area, survey methods and timing are available from Herr (1998).

I focused on adult bats, because the growth development of subadults would introduce non-linearity, and, therefore, confound the analysis in size distribution of the three species. I used forearm length as an indicator for bat size in preference over weight and other length measurements. Unlike the weight of individual bats, which can change much over the course of their feeding as they forage large amounts of insects, forearm length provides the most reliable measure of body size because it does not vary in an adult bat (Kunz et al. 1995; Herr 1998; Gardner et al. 2009; Gonsalves et al. 2013; McGuire et al. 2018; Davy et al. 2022). Female Vespadelus species are generally larger than males, which is also reflected in their longer forearms. This aligns with the general female-biased sexual size dimorphism in insectivorous bats (e.g. Lentini et al. 2015; de Oliveira 2020).

Statistical model of bat forearm length at different elevations

I develop a BMM, which is an improvement on conventional linear regression models, because it can also account for variations in intercept, slope, correlations between grouping variables and measurement errors. I used the Bayesian regression model using Stan (BRMS) library associated with the freely available statistical software R. All Monte Carlo Markov-chain (MCMC) random-sampling includes four chains with 5000 iterations each to sample the posterior parameter values and the final model has a full MCMC convergence, with R^ values of <1.05 and an effective sample size ratio of at least 0.1 (Vehtari et al. 2017; R Core Team 2021; Bürkner et al. 2023a; Gabry et al. 2023). The high-performance computing enabled backbone for BRMS, cmdStanR, reduced performance time requirements of the MCMC algorithm from days to hours when using parallel within-chain processing.

In the analysis I related bat forearm length (response variable) to elevation (explanatory variable). Besides elevation, there are two variables that group the data into adult females and males (two groups) of the species V. darlingtoni, V. regulus, V. vulturnus (three groups), with sex grouped within species level. On the basis of Fig. 2, forearm length has an approximate normal distribution (N) with a linear predictor μ and a standard deviation σ.

Fig. 2.

Vespadelus forearm length. AF, adult female; AM, adult male.

The approach has an expectation that the response variable forearm length linearly changes with the explanatory variables plus an additional term ε that represents unexplained variance, such as errors in measuring bat forearm length. Specifically, the variation is with each elevation i, and the males and females of the adult bats from the three species (i.e. species) and sex nested within species (species|sex), where α represents the intersect, β_1−n the slopes for each of the explanatory variables. This leads to Eqn 1, as follows:

(1)Forearmi~f(α+β1elevationi+β3species|sex+εi)Forearmi~N(μi,σ)

The BMM can incorporate the group-level variables sex and species in different ways, namely first, as components of the linear equation, and second, as group-level (random) parameters for intercept and slope of the equation. A ‘full model’ accounting for both ways provides the complete Eqn 2, where I expect a normal distribution for each intercept (α) and slope (β) parameter, with a standard deviation of τ. Incorporating random intercepts and slopes adds additional information into the modelling through using information from relationships of samples among and within the grouping variables (McElreath 2020).

In addition to the broad BRMS default priors, I used relatively uninformative priors for the MCMC sampling. The initial prior for the intercept α has a normal distribution, with a mean of 33 (the mean of all bat forearm length measurements) and a standard deviation of 1. Group-specific β intercept priors are of the same values. All other slopes have a normal prior, with a mean of 0 and a standard deviation of 1. Variations for the species|sex grouping intercept have a mean and standard deviation of 1, and the species-level group intercept variation has a mean of 2 and standard deviation of 1.25. This provides a sufficiently wider range of slopes in orders of magnitude higher than what one would expect from the data range. Variations around intersect and slopes are always held positive in the modelling. Additionally, the BMM accounts for group-specific correlations using Lewandowski-Kurowicka-Joe (LKJ) priors, Ω_k ~ LKJ(ζ), thus incorporating additional information into the estimation that would not be possible with a conventional linear regression model. Setting ζ prior to 1 assumes a uniform density for the correlation matrices (Lewandowski et al. 2009; Bürkner 2017; Nalborczyk et al. 2019).

I developed five model variations, allowing for different population effects and combinations in group-level variations with intercept and slope. I compared the effectiveness of these models with a leave-one-out analysis, using the theoretical expected log pointwise predictive density (elpd) to assess the predictive accuracy of both models (Vehtari et al. 2017).

Model 1 (Eqn 2) has all variables included as population-level effects and group-level (random) variation for intercept and elevation slopes, with the sex nested within species (species|sex); that is, I expect that the differences in sex are not the same among the species (i.e. male–female difference may result in different means and slopes for the different species).

(2)μi=(α+αspecies|sex)+(β1+βspecies|sex) elevationi+β2 species+β3 sexαspecies=N(0,τα, species)αspecies|sex=N(0,τα, species|sex)βspecies=N(0,τβ, species)βspecies|sex=N(0,τβ, species|sex)Priors:intersect α~N(33,1)variations τ for αspecies|sex~|N(1,1)|variations τ for αspecies~|N(2,1.25)|variations τ for βspecies|sex,βspecies~|N(0,1)|correlations Ω~LKJ(1)

A more limited Model 2 does not include the grouping variables as population effects. Instead, the model includes them just as group-level (random) parameters, which uses species by sex intercepts and species by sex variations on elevation slope (Eqn 3). Removing the elevation slope from the equation leads to Model 3 (Eqn 4).

(3)μi=(α+αspecies+αspecies|sex)+(β1+βspecies+βspecies|sex) elevationi

(4)μi=(α+αspecies+αspecies|sex)+β1 elevationi

For model comparison purposes, I included two more models, one that includes intercept and population effects for elevation (β₁) only (Eqn 5) and one model that includes only the intercept, that is, the α parameter (the grand mean intercept) is the only part of the right-hand side model equation. These last two models provide the benchmark for model comparison. Here an elpd difference (elpd diff) less than four with other models means that including these additional variables does not provide additional information, which would justify incorporating these variables.

(5)μi=α+β1 elevationi

Results

The results of the analysis showed that the bats comply with Bergmann’s rule, that is, bat size increases with elevation, and, by implication, decreasing temperature. All models provided similar mean elevation estimates but only Model 3 had an acceptable estimate error. The errors for Models 1 and 2 were of the same order of magnitude as the estimate (rounded to the nearest 1/10). Model 3 had an elevation error that was one order of magnitude smaller than the estimate (Table 1). When comparing models, Models 1–3 had a much larger elpd diff value than the benchmark models, indicating that the inclusion of elevation and nested grouping of species with sex provide improved model information (Table 2). The absolute elpd diff differences among Models 1–3 were small and <3, which suggests that these models have similar predictive performance (Sivula et al. 2022). Only Model 3 provided satisfactory estimates from a statistical perspective, and it used the least amount of parameter estimations. The additional information gained from including the grouping variables species and species|sex as population effects was limited and they resulted in an overly large standard error. This is fortuitous because I am mainly interested in the relationship of forearm length with elevation, and not in species specific or sex-related variations. Hence, I subsequently discuss the details of Model 3 only, where I further implement measurement errors.

The priors for the measurement error model are broad and provide a wide range of potential forearm length estimates, as the prior densities in Fig. 3 show. Although having low-informative priors, the model provides fairly tight posterior estimates of forearm length as the posterior densities in Fig. 3 attest. With measurement errors for forearm length and elevation included, the residual standard deviation estimate (σ) is 0.82 mm, leading to a residual variance of 0.67. The model estimates a forearm length increase of 0.11 mm for every 100 m elevation increase. The standard error for all estimates is at least one order of magnitude smaller for the population-level effects and less than half the size for residual standard deviations (Table 3), with the latter reflecting the overlap in forearm length between species and sex. This is also apparent in the overlapping posterior densities in Fig. 4 and the overlapping measured forearm length values (‘+’ in Fig. 5) and their posterior predictions (confidence intervals in Fig. 5). Despite this overlap, and even with the inclusion of measurement errors, the model attributed the elevation-related forearm length increase of 1.1 mm over a 1000 m elevation range (Fig. 5) fairly confidently, as the low standard error indicates (Table 3). In the context of Table 3 and Fig. 5, the reviewers struggled with the discrepancies of the reported intercept estimate in the table and location of the visual line in the figure. This difference is a result of summarising aspects of the model; the single grand mean intercept in the table is calculated as the mean of all posterior intercept estimates (i.e. including all intercept estimates at all elevations), whereas the line in the figure is a result of estimating the 5% quantile interval around the median. Although this is an interesting side effect of using mean posteriors predictions in model summaries, the intercept is not of much other interest for the final model. What is important to consider are the slopes (β) in Eqn 4 to ascertain how much a unit of elevation increase influences forearm length. This slope is a result of the relationship of elevation–forearm length and the variation around the grand mean intercept, which also includes the species and species by sex specific variations incorporated in the posteriors.

Fig. 3.

Comparison of prior and posterior densities of forearm length for three Vespadelus species and their sex. The dashed vertical line indicates the average bat forearm length. AF, adult female; AM, adult male.

Fig. 4.

Posterior densities and predictions of Vespadelus forearm length with elevation. AF, adult female; AM, adult male.

Fig. 5.

Posterior densities and predictions of Vespadelus forearm length with elevation. CI, credible interval as grey band; grey line, median of posterior prediction.

However, the intercept standard deviations of the group-level variables (species, α_species and sex within species, α_species|sex, see Eqn 4) show a somewhat larger standard error (Table 3). This indicates that these variables do not influence the relationship as clearly as does the elevation. This is reasonable, given the large range of forearm length among species and within sex of adult bats. Here, the model estimates an upper and lower range of 3.49 mm and 1.83 mm (estimate ± s.e) for species-related deviation from the grand mean intercept respectively. For sex within species, this deviation ranges from 1.62 mm (upper limit) to 0.6 mm (lower limit).

Discussion

This analysis showed a clear increasing relationship between bat size and elevation as a proxy for decreasing temperature, in line with Bergmann’s prediction and based on his original assertion that the surface to volume ratio of body size in animals influences their distribution. There is a positive statistical relationship of bat size (represented by their forearm length) with elevation for a group of three taxonomically closely related, sympatric species. The BMM provided estimates of the relationship with high confidence, as the low standard errors in Table 3 show. While the 1.1 mm increase per 1000 m is smaller than species- and sex-related changes, and 1/8 of the residual standard deviation, the BMM readily estimated residuals and the contribution of each explanatory variable to the prediction. The advantage of the BMM approach over conventional linear regression is that it can readily incorporate measurement errors for both the predictor and the explanatory variables. This leads to more robust estimates and provides support for an increasing relationship of bat size with elevation.

The findings of increasing forearm length (as a reliable size indicator in bats) with elevation (and, by implication, decreasing temperature) conforms with the generality of Bergmann’s rule, where larger animals can tolerate colder temperatures as a result of reduced surface to volume ratio. The result also supports Bergmann’s prediction that the smaller species of the same genera, where the main difference is their size, require a warmer climate and this requirement is directly related to the size difference (‘… Gäbe es nun Genera, deren Species sich so weit als möglich (eine Beschränkung dieser Möglichkeit ist im Anhange angedeutet) nur durch die Grösse unterschieden, so würden die kleinen Arten derselben durchweg ein wärmeres Klima fordern und zwar nach einem aus der Grössendifferenz genau bestimmten Maasse…’; Bergmann 1848, p. 46).

Accepting this conformance for this study relies on the following two causal relationships: (1) the relationship between body size and temperature that Bergmann’s rule implies, and (2) the linear relationship of average temperature increase with elevation. The causal relationship between the need for increased body size in declining temperature is based on the unequal surface:volume-ratio change (i.e. where volume increases with the power of 3 and surface with the power of 2). As a consequence, larger warm-blooded animals can regulate their core temperature more easily than can smaller animals. This is true for animals of closely related taxa, where the smaller trait differences are less likely to compensate temperature-regulation deficiencies (Brown and Lee 1969; Watt et al. 2010; Clauss et al. 2013; Baldwin et al. 2023). The more unrelated the taxa are, the broader are the organisational modifications that can over-ride the principle (loosely translated from ‘… Die Breite dieser Modificationen, welche sämmtlich auf der Organisation beruhen, muss natürlich um so mehr abnehmen je ähnlicher die Thiere einander in ihrer Organisation sind. So differente Grössen, wie sie in den Extremen in den verschiedenen Klimaten, besonders im tropischen gegeben sind, setzen grosse Verschiedenheiten der Organisation voraus…’) (Bergmann 1848, p. 46). This study used three closely related Vespadelus taxa for analysing the body size to elevation relationship, with the assumption that elevation is a proxy for temperature.

The second assumption is based on the (region-specific, very generalised) linear relationship between elevation and temperature known as the lapse rate. This is as a rule of thumb, which, on average, is 0.75°C decrease for every 100 m of elevation (Galloway 1988; Rolland 2003). However, the lapse rate implies a general relationship that varies geographically with temperature, precipitation and with air pressure, and this variation is exacerbated in mountainous regions (Nunez and Colhoun 1986; Minder et al. 2010; Pan et al. 2021; Rosas-Chavoya et al. 2021). In the Australian Alps region, for example, this rate is applicable over the winter period, but can be 0.91°C degrees per 100 m over summer and varies with microclimate locally (Green 2014). Hence, the model in this paper incorporates elevation variation, which is a way of addressing this variability in temperature, in addition to allowing for bat flight distances increasing the elevation range.

The size–elevation relationship also has implications for insectivorous bats in the study area under climate change. If insectivorous bats align with Bergmann’s rule, anthropogenic-related changes in climate may influence their survival ability at different elevations through a warming effect (see also Law et al. 2018). Climate change is likely to influence the bode size, and there is a risk that animals may not have sufficient time for adaption if the current environmental change is more rapid than their adaptaptability through phenotypic plasticity (Charmantier et al. 2008; Teplitsky and Millien 2013; Merilä and Hendry 2014; Thackeray et al. 2016; Radchuk et al. 2019). Current evidence around this risk for V. darlingtoni is inconclusive because previous studies did not find a consistent trend for a change over recent decades in density, survival or condition from one montane study area (Law et al. 2018; Law et al. 2023). Nevertheless, the risk still applies to bats in the study area, and, in addition, they are exposed to climate change-induced suitable habitat reduction. In areas where larger bats occur at higher altitudes, this may push them further into higher areas where more extreme temperature changes are prevalent, thus limiting their suitable-habitat availability (Moritz et al. 2008; Chen et al. 2011; Spence and Tingley 2020; Spence et al. 2022).

In the context of the study area, if global warming results in a competitive advantage for smaller bats, one could expect the following two overlapping outcomes for the larger species: (1) the average forearm length (by implication, size) of the population would reduce, thus exposing bats to more competition with the smaller species; and (2) the competition from the smaller species would push the larger species into higher (by implication, cooler) areas. The three bat species have overlapping echolocation call-frequency ranges and they will feed on the full size ranges of available insects, with V. vulturnus having a potential upper size limit of prey it can ingest when compared with the two larger species. However, even though a smaller bat may not be able to ingest as large a prey owing to its physical size, it has the advantage of accessing smaller tree hollows for roosting, as well as being able to use the hollows larger bats rely on (see e.g. Herr 1998). Being pushed into higher habitats may have negative consequences. These higher habitats are already limited; within the study area, V. darlingtoni has a stronghold in the higher elevations and mountain peaks, so there is nowhere higher to go. Further afield, their distribution is mostly limited to below the subalpine areas, to an altitude of approximately 1500 m, where large trees are in low numbers (Slatyer 2010). This limit is related to their preference for larger roost trees (above 80 cm in diameter at breast height) with hollow development that provides appropriate shelter and microclimate for maternity groups (Taylor and Savva 1988; Herr 1998; Hoye et al. 2008).

A very coarse back of the envelope (BOTE) calculation for the study area can provide an estimate of suitable habitat reduction. When assuming that current elevation provides the current suitable temperature range for the bats, a 1.5°C and 3°C temperature increase would change suitability to a 1.5°C (or 3°C) higher elevation equivalent. Under climate change, tree-line upward movement may occur, but the availability of larger roosting trees would require significant time to occur, so I ignore upslope habitat improvement from an upward-shifting subalpine boundary. Further, I assume that (1) competition from the smaller species makes habitat for V. darlingtoni less suitable. Although competition among bats is very difficult to prove directly, ecological and evolutionary theory dictates its validity, and studies have generally accepted it (Roeleke et al. 2018; Morningstar et al. 2019; Salinas-Ramos et al. 2019; Novella-Fernandez et al. 2020). I also assume that (2) the average temperature lapse rate for the study area does not change under future climate (i.e. ignoring lapse-rate change at lower altitudes, see e.g. Mokhov and Akperov (2006), Kröner et al. (2016)).

The result is an elevation change of 200 m for 1.5°C and 400 m for 3°C, that is, 1300 m current elevations are equivalent to 1500 m elevations under 1.5°C increase and 1100 m current elevations are equivalent to 1500 m under 3°C increase. Under such circumstances, there is a reduction of suitable habitats for V. darlingtoni in those areas that in the future would be the 1500 m or above equivalent of the current elevations. This area is readily calculated for the study area on the basis of the Australian 9-s DEM (see e.g. Hutchinson 2008), to a 0.6 km² (0.9%) suitable-habitat loss for a 1.5°C average temperature increase and a 1.7 km² (3%) suitable-habitat loss for a 3°C average temperature increase. Climate change also may result in other changes such as habitat loss from, for example, increased fires and changes in vegetation as a response to reduced precipitation or a change in the tree line through climate tracking (Körner 2012; Zylstra 2018; Young and Littleboy 2019; Naccarella et al. 2020; Ji et al. 2022; Verrall et al. 2023). Although this BOTE does not consider such changes, they are still useful in estimating potential climate-change impacts on larger bats. Here, the causal link is the mechanism behind Bergmann’s rule, and its consequences that allow smaller bats to access higher areas under a warmer climate.