Introduction

Seasonal food scarcity challenges animal energy budgets. A positive or less negative energy balance across seasons can be achieved by a more constant net energy intake than predicted solely by food abundance (Sol, 2009), or by investing less in costly organs (Heldstab et al., 2018). A link between both strategies is the brain: A relatively large brain can improve cognitive ability and behavioral flexibility (Benson-Amram et al., 2016; Lefebvre et al., 2004; Reader and Laland, 2002; Sol et al., 2005), enabling animals to effectively locate diverse and dispersed food sources to buffer environmental fluctuations in seasonal habitats (‘cognitive buffer hypothesis’) (Allman et al., 1993; Sol, 2009; van Woerden et al., 2012). However, brain tissue has high metabolic costs (Aiello and Wheeler, 1995; Lukas and Campbell, 2000; Mink et al., 1981) that may not be temporarily reducible (Mink et al., 1981), constraining brain size evolution under periodic food scarcity (‘expensive brain hypothesis’) (Isler and van Schaik, 2009).

Food scarcity can also affect physiological responses. For example, a longer digestive tract may permit more efficient resource uptake during a short active period and be favored by selection in species with prolonged hibernation (Sibly, 1981). The evolution of the digestive tract could thus parallel that of the brain. By contrast, there could be an evolutionary trade-off between the two organs (‘expensive tissue hypothesis’; Aiello and Wheeler 1995; also see Isler and van Schaik (2006) for a more general ‘energy trade-off hypothesis’). Furthermore, physiological buffering often involves a seasonal reduction in metabolic rate or activity (e.g., hibernation), with energy drawn from stored fat reserves (Heldstab et al., 2016). Buffering lean periods by fat stores and reduced activity contrasts with cognitive abilities to cope with food scarcity, suggesting differential investment in adipose and brain tissue between species in highly seasonal environments (‘fat−brain trade-off hypothesis’; Heldstab et al., 2016; Navarrete et al., 2011). Finally, allocating fat stores to sustained brain function throughout prolonged hibernation might compete with investments in other tissues (Isler, 2011; Isler and van Schaik, 2009).

Species that breed soon after hibernation, such as some mammals (Place et al., 2002; Psenner, 1957) and many amphibians (Fei and Ye, 2001; Wells, 1977), face further unique challenges, as their reproductive tissue must regrow before emergence when stored resources are most limited (Isler, 2011; Isler and van Schaik, 2009). Testes, however, may be subject to intense selection by sperm competition resulting from female multiple mating (Lüpold et al., 2020), a widespread phenomenon where males cannot monopolize their mates (Lüpold et al., 2014). Indeed, across anurans (frogs and toads) that often breed in small water bodies, males invest relatively more in their testes and less in their forelimbs (used in pre-mating competition) as population density increases (Buzatto et al., 2015; Lüpold et al., 2017). If a shorter active period leads to more synchronized breeding and thus increases the risk of sperm competition or sperm depletion (Vahed and Parker, 2012), selection for relatively larger testes would be stronger where fat stores need to last longer, affecting resource demands and allocation during hibernation. Including reproductive investments and breeding patterns in studies of allocation trade-offs in response to hibernation and environmental seasonality would thus seem critical but remains to be done, particularly in the context of brain evolution.

The opposing selection pressures on brain or gonad size (i.e., cognitive or fitness benefits versus metabolic costs), varying degrees of seasonality and diverse strategies for buffering periodic food scarcity (e.g., cognitive versus physiological) between species render environmental fluctuations an ideal context for studying brain size evolution. The different hypotheses for the coevolution of brain size with other organs were developed separately for different mammalian or avian taxa, the two vertebrate classes with the largest brains relative to body size (Jerison, 1973). These hypotheses have yet to be directly tested against one another in a single taxon and ideally in the immediate context of seasonal activity, considering the extent, rather than the mere presence/absence, of hibernation. This last point is important because ‘hibernation’ can range between brief inactive bouts and long dormancy, with different energetic and life-history constraints. Furthermore, understanding the generality of the patterns reported in mammals or birds needs to be validated in other taxa, ideally with smaller brains and different energy demands. Such a generalization would help contextualize brain evolution in the two largest-brained taxa in relation to selection on encephalization and metabolic constraints.

A particularly suitable system is presented in ectothermic (‘cold-blooded’) species whose metabolism and activity are tightly linked to ambient temperature (Wells, 2007), possibly leading to stricter physical boundaries to behavioral flexibility than in the endothermic (‘warm-blooded’) mammals or birds. Compared to the divergent groups of mammals, for which most hypotheses on brain evolution have been proposed, anurans are also more homogeneous in body size and shape, diet and locomotion (Kardong, 2019), while still showing strong environmental effects on brain, reproductive and other investments (Liao et al., 2022; Lüpold et al., 2017; Wells, 2007). Being less influenced by Bauplan and lifestyle, evolutionary trade-offs with brain size could thus be easier to isolate in anurans than in mammals.

Here, we studied across 116 anuran species how varying ‘hibernation’ periods and their environmental correlates affect male investments in brain, testes, body fat, limb muscles, and the main visceral organs (see Fig. 1).Unlike mammals, anurans do not actively depress their metabolism but rather stop activity when ambient temperatures drop below the activity range (referred to here as ‘brumation’ for distinction; Pinder et al. 1992; Wells 2007). However, in both taxa, hibernators use fixed energy stores across sequential investments, albeit at a lower metabolic rate (Staples, 2016), whilst non-hibernators can partly replenish resources but need more food even when it is scarce (Heldstab et al., 2016). These contrasting strategies may be associated with a varying cost-benefit balance between organs, including the brain.

Figure 1.

Although anurans halt general activity during brumation, this does not mean complete brain inactivity. Several species move in their burrows or underwater hibernacula in response to changes in soil temperature or oxygen concentration (Holenweg and Reyer, 2000; Stinner et al., 1994; van Gelder et al., 1986), or when disturbed (Niu et al., 2022; Tattersall and Ultsch, 2008). Moreover, brumating frogs can renew more brain cells than active frogs, possibly to prevent brain damage (Cerri et al., 2009), and this could increase with brain size. Larger brains could also be less tolerant to low oxygen conditions (Pinder et al., 1992; Tattersall and Ultsch, 2008; Wells, 2007) if the findings from other ectotherms extend to anurans (Sukhum et al., 2016). Overall, maintaining a relatively larger brain while brumating may entail higher costs that could constrain brain evolution compared with species with short or no seasonal inactivity.

Instead of the typical indirect proxies or scores (Heldstab et al., 2018, 2016; Luo et al., 2017; Navarrete et al., 2011), we quantified seasonal changes in tissue size and estimated brumation duration as the time that temperatures were continuously below the species-specific activity threshold. These estimates aligned with periods of inactivity recorded during field surveys, and were robust across multiple sensitivity analyses (e.g., considering temperature buffering by underground burrows). We then tested the different hypotheses on brain evolution and further integrated variation in breeding context and reproductive evolution in this resource allocation framework. By examining ectothermic species with diverse energy demands and brain sizes, our findings provide a crucial context for understanding brain evolution and resource allocation.

Results

Determinants of brumation duration

For a mean (±SD) of 3.41 ± 0.95 males in each of 116 anuran species, we combined experimentally determined thermal activity thresholds with multi-year temperature fluctuations at their collection sites to estimate the species-specific brumation periods (details on methodology and various sensitivity analyses for validation in Material and Methods). These periods averaged between 0.6 ± 0.5 and 250.5 ± 16.7 days across the five years examined, with high repeatability within species (R = 0.95 [95%CI: 0.93, 0.96]; Fig. S1; Table S1). These brumation periods largely fell into three groups with ≤9 (N = 22 species), 18−27 (N = 3), and ≥47 days of expected brumation (N = 91), respectively. As the ground microclimate may buffer some of the fluctuations in air temperature and frogs can endure short cold spells without dormancy, we conservatively considered those 25 species with ≤27 days below their experimental temperature threshold unlikely to show any sustained brumation.

We employed phylogenetic generalized least-squares (PGLS) models (Freckleton et al., 2002; Ho and Ané, 2014) to dissect the factors driving brumation duration. Our analyses indicated that the brumation period increased with the latitude and elevation of study sites, along with the extent of annual temperature fluctuations (r ≥ 0.42, t₁₁₄ ≥ 5.00, P < 0.001; Table S2). In contrast, brumation duration was inversely correlated with the annual mean temperature and precipitation, and dry season length (r < −0.29, t₁₁₄ < −3.24, P < 0.002), but not significantly covarying with longitude or annual fluctuations in precipitation (r < 0.09, t₁₁₄ < 0.95, P > 0.34; Table S2). Except for the duration of the dry season, these results persisted for those 91 species with a high probability of sustained brumation (see above; Table S3). Among these species, those from cooler and more seasonal climates exhibited lower entry and emergence temperatures from their inactive state (Table S4), suggesting heightened cold tolerance to maximize their active period.

Impact of brumation on tissue investments

To investigate the repercussions of brumation on individual tissue investments, we used the same males from our 116 species. In separate sets of PGLS models, neither snout-vent length (SVL) nor body mass covaried with brumation duration or any other environmental variable (|r| ≤ 0.15, |t₁₁₄| ≤ 4.03, P ≥ 0.11; Table S5), except for a weak, non-statistically significant trend toward reduced body mass at higher elevations (r =−0.18, t₁₁₄ = −1.91, P = 0.06, phylogenetic scaling parameter λ = 0.94 [95%CI: 0.85, 1.00]). The different tissues increased with body size, albeit with different allometric slopes. Body fat and limb muscles showed steeper scaling (all β ≥ 1.10 [1.03, 1.18]), brain size exhibited a shallower slope (β = 0.49 [0.44, 0.54]; Table S6; Fig. 1A), and the remaining tissues did not deviate from proportionate scaling (i.e., 95%CI including 1.00; Table S6; Fig. 1A). Hence, the evolution of brain size appears to be more constrained than that of other organs when selection favors larger body size.

Given these distinct allocation patterns between brain and other tissues, we next tested whether extended brumation constrains brain size evolution, akin to suggestions for mammals based on presence/absence of hibernation (Heldstab et al., 2018). In PGLS models, absolute brain size was independent of brumation duration (r = −0.10, t₁₁₅ = −1.08, P = 0.28, λ = 0.89 [0.73, 0.97]), but accounting for SVL as a body size proxy, species with prolonged brumation had relatively smaller brains, quantified during breeding (partial r, r_p = −0.31, t₁₁₃ = −3.50, P < 0.001, λ = 0.35 [0.00, 0.61]; Fig. 2A, Table S7) or pre-hibernation (N = 50 species means based on 2.64 ± 0.94 males each: r_p = −0.51, t₄₇ = −4.03, P < 0.001, λ = 0.00 [0.00, 0.56]; Table S7). This trend was distinct from body size changes in response to brumation, as SVL remained independent of brumation duration (Table S5). On average, brumating species tended to possess relatively smaller brains compared to those less likely to experience prolonged brumation (Table S8), aligning with mammalian studies utilizing hibernation presence/absence (Heldstab et al., 2018). This trend extended to those 91 species categorized as brumating for some period (Table S9), reinforcing the link between brumation and brain evolution beyond coarse binary classification. This pattern further persisted when recalculating brumation periods using conservative thresholds of 2°C or 4°C below their experimentally derived thresholds, simulating shelter buffering (details and validation in Material and Methods; Table S10).

Figure 2.

Species displaying extended brumation periods further exhibited relatively more body fat (r_p ≥ 0.25, t₁₁₃ ≥ 2.72, P ≤ 0.008; Fig. 2B) and, particularly during breeding, had relatively larger testes (r_p = 0.36, t₁₁₃ = 4.06, P < 0.001, λ = 0.77 [0.40, 0.90]; Fig. 2C), along with relatively smaller hindleg muscles (r_p = −0.22, t₁₁₃ = −2.37, P = 0.02, λ = 0.22 [0.00, 0.51]; Tables S7). These patterns generally held true when using the presence/absence of brumation (Table S8), buffered temperature fluctuations (Tables S10), or excluding the 25 species unlikely to brumate (except for the non-significant effect on body fat; Table S9). Other tissue sizes remained independent of brumation (Tables S7 to S10) and showed not significant changes between sampling periods (Fig. 2).

Comparing pre-and post-brumation males, averaged across all species, we found a roughly 50% reduction in fat tissue and a 100% increase in testis size, indicating resource depletion and testicular regrowth during brumation, respectively (Fig. 1B). Only brain size deviated from zero among the remaining tissues, but this change was minimal compared to the body fat and testes, and within the range of many unchanged tissues. Thus, the biological significance of this putative increase in brain size during brumation remains uncertain, possibly reflecting general differences between sampled individuals. We will thus refrain from further interpretation. Among the two tissues with considerable change, the extent of fat depletion increased significantly with the brumation period (r = −0.36, t₄₈ = −2.68, P = 0.01, λ = 0.00 [0.00, 0.55]), as did that of testis regrowth (r = 0.39, t₄₈ = 2.89, P = 0.006, λ = 0.82 [0.16, 0.98]; Fig. 1C).

To examine whether the increase in relative testis size with prolonged brumation might be mediated by a shorter, more synchronized mating season (Wells, 2007), we tested for links between brumation duration and different breeding parameters. Prolonged brumation notably shortened the breeding season (r = −0.57, t₄₁ = −4.47, P < 0.001, λ = 0.00 [0.00, 0.38]; Fig. S2A), exerting a stronger effect than climatic variables (Table S11), particularly when considered together (Table S12). Hence, the effect of these climatic variables may be mediated by brumation. Furthermore, a shorter breeding season increased the probability of dense breeding aggregations (phylogenetic logistic regression: N = 42, z = −3.03, P = 0.002, α = 0.02; Fig. S2B). That brumation might mediate climatic impacts on breeding aggregations through the duration of the breeding season was also supported by a phylogenetic confirmatory path analysis (Gonzalez-Voyer and von Hardenberg, 2014; von Hardenberg and Gonzalez-Voyer, 2013; Fig. S3; Table S13). Finally, when combining these data with published data on the density of breeding populations (Lüpold et al., 2017; N = 8 species overlapping), a trend emerged toward higher mean population densities in species with shorter breeding seasons (r = −0.69, t₆ = −2.37, P = 0.06, λ = 0.00 [0.00, 1.00]), albeit with a small sample size (Fig. S2C).

To explore possible causal links between breeding parameters and relative testis size, we employed directional tests of trait evolution (Pagel, 1994; Revell, 2012), assessing whether changes in two binary traits are unilaterally or mutually dependent, or independent (Pagel, 1994). To this end, we converted our continuous to binary variables (see Material and Methods). Relating small/large testes to short/long breeding seasons, the best-supported scenario based on the Akaike Information Criterion (AIC) was independent evolution. Yet, the model with changes in relative testis size dependent on those in the breeding season found similar support (ΔAIC = 0.83, w_AIC = 0.33 compared to independent model with w_AIC = 0.50; Fig. S4A), differing from the remaining two models with clearly higher AIC scores (ΔAIC ≥ 3.40, w_AIC ≤ 0.09). In another directional test, between relative testis size and aggregation formation, the independent model was the best-supported (w_AIC = 0.66), followed by the scenario of increased relative testis size in response to aggregation formation (Fig. S4B), albeit above a ΔAIC cut-off of 2 (ΔAIC = 2.38, w_AIC = 0.20), with the remaining models being clearly less supported (ΔAIC ≥ 3.64, w_AIC ≤ 0.11). Hence, it is at least possible that breeding conditions mediate the positive relationship between brumation duration and relative testis size in our relatively small sample of species, whilst a response of breeding conditions to variation in relative testis size is clearly rejected by our analyses.

Covariation between tissues

As all tissues rely on the same finite resources, their responses to brumation are likely interconnected. Pairwise partial correlations, controlling for SVL and phylogeny, revealed positive covariation or non-significant associations among all tissue masses (Fig. S5). As such, our data do not support the expensive tissue (Aiello and Wheeler, 1995) or the more general energy trade-off hypotheses (Isler and van Schaik, 2006), which predict brain size trade-offs with the digestive tract or other costly organs, respectively. However, since brain size differed from fat, hindlimb muscles and testes in allometric relationships and responses to brumation, pairwise correlations may not capture more complex allocation. To examine the relative investments in these four most informative tissues simultaneously, we represented total body mass as proportions of these tissues and the remaining mass, generating a five-variable compositional dataset (van den Boogaart and Tolosana-Delgado, 2013). The combined mass of all four focal tissues scaled proportionately with body size (allometric β = 1.02 [0.96, 1.07], λ = 0.01 [0.00, 0.11]), confirming size-independent proportion of the total resources allocated to the four focal tissues combined. However, the resource distribution among these tissues varied considerably between species. Pairwise correlations between the four focal tissues, transformed to centered log ratios (van den Boogaart and Tolosana-Delgado, 2008) and controlling for phylogeny, showed brain mass to covary negatively with fat and testis mass, whilst testis mass covaried negatively with hindlimb muscle mass but not with body fat (Fig. 3A-D; Table S14).

Figure 3.

To further examine the effect of brumation duration on all five variables simultaneously, we conducted a phylogenetic multivariate regression analysis (Clavel et al., 2015) on the same compositional data, but now transformed to isometric log ratios as recommended for multivariate models (van den Boogaart and Tolosana-Delgado, 2013, 2008). Brumation duration significantly affected anuran body composition (Pillai’s trace = 0.33, effect size ξ² = 0.30, P = 0.001). The back-transformed coefficients of body fat (0.23) and testis mass (0.26) exceeded the expected coefficient of 0.20 if brumation had no effect, whilst brain mass, hindlimb muscles and the rest of the body had lower coefficients (0.16, 0.17, and 0.18, respectively; also see Fig. 3).

Finally, a phylogenetically informed principal component analysis (Revell, 2012) confirmed the negative associations of brain size with body fat and testis mass, and that of testis mass with hindlimb muscles (Fig. 3E,F, Table S15). Here, the first three principal components (PC1 to PC3) explained 84.7%, 8.0% and 5.1% of total variance, respectively. Although PC2 and PC3 explained a relatively small proportion of the total variance, they separated the different tissues. PC2 loaded mainly by brain size (0.30) and testis mass (−0.44), and PC3 by brain size (0.26) and body fat (−0.31). Biplots indicated negative associations between the vectors of the same traits as above within the multivariate trait space (Fig. 3F). Further, brumation duration covaried negatively with PC2 (r = −0.53, t₁₁₄ = −6.59, P < 0.0001, λ = 0.74 [0.49, 0.92]), consistent with reduced brain size and increased in testis mass towards longer brumation, but it was not associated with PC3 (r = −0.12, t₁₁₄ = −1.30, P = 0.20, λ = 0.41 [0.00, 0.70]).

Direct and indirect effects revealed by path analysis

To untangle the evolutionary links between tissue sizes, and to test competing brain evolution hypotheses in the context of brumation, we conducted a phylogenetic confirmatory path analysis (Gonzalez-Voyer and von Hardenberg, 2014; von Hardenberg and Gonzalez-Voyer, 2013) with 28 predetermined path models (Fig. S6; Table S16). The averaged model (Fig. 4) confirmed the negative effect of prolonged brumation on relative brain size (β = −0.15 [−0.22, −0.07]), along with direct (β = 0.16 [0.07, 0.26]) and indirect positive effects on relative testis size (Fig. S7). These effects on testis size were mediated by the relative amount of adipose tissue, which increased with brumation duration (β = 0.13 [0.06, 0.20]) and digestive tract size (β = 0.51 [0.37, 0.65]), and in turn positively influenced relative testis size (β = 0.39 [0.18, 0.59]).

Figure 4.

Discussion

Our study on anurans, with validated brumation periods and direct measures of expensive tissues, provides novel insights into brain and reproductive evolution in ‘cold-blooded’ organisms exposed to environmental seasonality. We found that species with longer brumation exhibited relatively smaller brains and allocated greater fat reserves primarily in reproduction, possibly due to the shorter breeding season with its socio-ecological implications.

Environmental and resource considerations in brain evolution

We demonstrated that species in cooler and more seasonal climates were more cold tolerant, thereby likely optimizing their active period. Yet, low temperatures impair foraging and digestion in ectotherms (e.g., Fontaine et al., 2018; Riddle, 1909), such that high seasonality may lead to longer brumation and smaller brains. These results confirm that, unlike birds (Sol, 2009) and some mammals (van Woerden et al., 2010), challenging and unpredictable environmental conditions select for physiological rather than cognitive buffering in anurans (Luo et al., 2017). Supporting a large brain may not be sustainable without continued resource intake, or larger brains could be less tolerant to hypoxic conditions during brumation (Sukhum et al., 2016). However, selection for relatively larger brains may also simply be stronger in species with longer active (and short brumation) periods owing to extended cognitive benefits such as predator evasion (Kotrschal et al., 2015; Liao et al., 2022) or exploitation of better and more diverse food sources (Jiang et al., 2022; Lefebvre et al., 1997).

In pairwise comparisons, the relative sizes of the tissues examined here, including the brain, were generally positively correlated. These results reject both the expensive tissue and energy trade-off hypotheses (Aiello and Wheeler, 1995; Isler and van Schaik, 2006), which predict trade-offs of brain size with the size of the digestive tract or other costly organs, respectively. This lack of support in anurans aligns with a previous report in mammals (Navarrete et al., 2011) despite their smaller brains and vastly different ecology and physiology, including a lower metabolic rate and largely lacking physiological thermoregulation. When focusing jointly on the four tissues (brain, body fat, testes, hindlimb muscles) that covaried with brumation duration, however, relative brain size covaried negatively with the relative mass of both fat tissue and testes, supporting the fat−brain trade-off (Navarrete et al., 2011) or expensive sexual tissue hypotheses (Pitnick et al., 2006), respectively.

Brumation and fat tissue

Species with longer brumation periods exhibited relatively more total body fat and a higher degree of its depletion, supporting the hypothesis that anurans buffer lean periods by metabolizing stored fat (Huang et al., 2020; Luo et al., 2017). Although adipose tissue may not itself be metabolically expensive, transporting it adds costs to locomotion, particularly when jumping away from predators (Moreno-Rueda et al., 2020) or climbing trees compared to moving horizontally on land or in water (Alexander, 2003; Hanna et al., 2008). Consistent with this notion, arboreal species tended to be leaner compared to (semi)aquatic or terrestrial species (Table S17), controlling for brumation duration and relative brain size, both of which we had shown to covary with body fat (Figs. 2 and 3).

Brumation and testis evolution

Species with prolonged brumation also had relatively smaller hindleg muscles and larger testes. The negative relationship between hindleg muscle mass and brumation duration may be linked to more movement during a longer active period, including predator evasion (Liao et al., 2022; Marchisin and Anderson, 1978). Larger testes may result from a shorter breeding season, leading to denser and more synchronous mating activity (Wells, 2007), as suggested by our path analysis. Breeding aggregations increase male-male competition over fertilization and thus enhanced investments in sperm production (Liao et al., 2018; Lüpold et al., 2020). Our results thus reveal how brumation patterns, influenced by environmental variation and physical constraints, affect the socio-ecological context of breeding, the mode and degree of sexual selection, and ultimately the evolution of mating systems, broadening Emlen & Oring’s (1977) general predictions.

In addition to the average size of the testes, their seasonal change also varied with the brumation period. Seasonally breeding anurans regress and regrow their testes between mating seasons (Ogielska and Bartmańska, 2009). Non-brumating species can use energy uptake to compensate for testicular recrudescence, whilst those with a short breeding season after a prolonged inactive period depend on the stored fat to regrow their testes before or immediately after emergence from their hibernaculum. Hence, resources are diverted away from the brain and other organs, especially in species such as Brachytarsophrys spp., in which the fully developed testes combined weigh 12−14 times more than the brain (Data S1).

Brain-testis trade-off

A phylogenetic path analysis confirmed the negative effect of brumation duration on relative brain size and revealed its direct and indirect positive effects on relative testis size, mediated by the amount of adipose tissue, which responded to variation in the inactive period (energetic demand) and the size of the digestive tract (energy uptake). That body fat did not contribute to brain size evolution in this more comprehensive analysis compared to pairwise correlations suggests that the fat−brain trade-off may not be direct. Rather, longer brumation, and thus a short active period, may enhance selection on fat storage for testicular investments in addition to starvation avoidance, while reducing selection for larger brains due to a shifted balance between cognitive benefits and energetic costs (Fig. 3).

A brain-testis trade-off has been reported for bats (Pitnick et al., 2006), but not replicated later in the same (Dechmann and Safi, 2009) or other mammalian taxa (Lemaître et al., 2009). In anurans, the apparent trade-off may result indirectly from opposing selection on brain and testis sizes via environmental seasonality and relative durations of the active and inactive periods. The testes may evolve in response to increased sperm competition and depletion during the shorter and more synchronized breeding season. The brain, while also responding to sexual selection (Mai et al., 2020), is central to various activities, including feeding (Lefebvre et al., 1997) or predator avoidance (Kotrschal et al., 2015; Liao et al., 2022) that are themselves subject to climatic conditions and may independently influence brain evolution. Additionally, whereas testes can regress to save energy when inactive (Ogielska and Bartmańska, 2009), brain metabolism may be less reducible (Mink et al., 1981), resulting in different cost-benefit balances between these organs in relation to seasonality.

Conclusions

In conclusion, our analyses show how brumation in anurans, influenced by high environmental seasonality, may affect resource allocation between costly tissues, directly or through its environmental correlates. The non-independent selective processes promoting diversification in different traits emphasize the need to study the evolutionary trajectory of a given trait such as brain size in the immediate context of both simultaneous investments to other tissues and the species-specific ecology. Our findings also have important implications in the context of sustained environmental change, exacerbated by climate change with its effects on temperature and precipitation patterns around the globe (IPCC, 2022). Such climatic shifts may disrupt the timing of brumation and breeding in many anurans, change their exposure to pathogens or predators through modified activity patterns, affect their ability to find food or suitable breeding sites, or change their population dynamics through resource availability, intra-and interspecific competition, dispersal capabilities, and gene flow (Alves-Ferreira et al., 2022; Blaustein et al., 2010, 2001; Carey and Alexander, 2003). In species with temperature-dependent sexual differentiation, thermal shifts may further change operational sex ratios in breeding populations and thus likely mating dynamics and sexual selection (Eggert, 2004; Lüpold et al., 2017; Ruiz-García et al., 2021). All these factors impose intense environmental pressure on resource acquisition and allocation patterns, and it remains to be seen to what extent variation in the adaptability, and thus resilience, between species exposed to environmental change is attributable to such competing needs between investments and species-specific constraints.

Materials and Methods

Sample collection and preparation

Between 2010 and 2020 and as part of concurrent studies, we collected a total of 396 sexually mature males from 116 anuran species (3.41 ± 0.95 males each) in post-brumation breeding condition and an additional 132 adult males from 50 of these species (2.64 ± 0.94 males each) shortly before entering their hibernacula (Data S1 and S2). For each species, we sampled all males at a single location in southern and western China with known longitude, latitude, and elevation (Data S3). Upon transfer to the laboratory, we sacrificed the individuals by single-pithing, measured their snout-vent length (SVL) to the nearest 0.01 mm with calipers and then preserved them in 4% phosphate-buffered formalin for tissue fixation.

After two months of preservation, we weighed each complete specimen to the nearest 0.1 mg using an electronic balance to obtain body mass before dissecting them following a strict protocol. We separately extracted the brain, heart, liver, lungs, kidneys, spleen, digestive tract, testes, limb muscles, and fat stores, cleaned these tissues and immediately weighed them to the nearest 0.1 mg with an electronic balance. We additionally measured the length of the digestive tract to the nearest 0.01 mm using calipers. We excluded emaciated individuals or those exhibiting visible organ pathologies from our analyses.

Environmental seasonality

For each collection site, we retrieved from the 30-year climate history of https://www.meteoblue.com the monthly mean temperature (in °C) and total precipitation (in mm) (Data S3) and used these values to calculate location-specific annual means and coefficients of variation. We also determined the duration of the dry season, P2T, as the number of months, for which the total precipitation was less than twice the mean temperature (Walter, 1971).

Brumation period

One way that anurans can physiologically respond to seasonality is by adjusting their thermal sensitivity and thus brumation period (Wells, 2007), which in turn could directly or indirectly affect the evolution of brain size (Heldstab et al., 2018). Hence, we estimated the brumation period for all 116 species. To this end, we visited the field sites for 30 of our species daily around the expected start and end times of brumation (based on prior experience). For each species, we recorded the dates and temperatures (using a Kobold HND-T105 high-precision thermometer to the nearest 0.1°C) when the last frogs of a given species were seen at the end of their active period (with no further activity detected for at least seven days) and when the first individuals were detected in the spring. For the same 30 species (and using the same individuals as for morphological measurements), we then experimentally simulated brumation using a Q18 temperature-controlled refrigerator in Shenzhen Pioneer (SAST). We gradually lowered and raised the temperature at a rate of 0.5°C/hour and recorded the temperature at which test subjects entered and left brumation. When entering brumation, anurans drop their heart rate, become sluggish, draw their nictitating membranes across the eyes for protection, spread their legs for stability, and change their physiology to avoid freezing, start breathing through their skin, or switch to anaerobic metabolism (Fei and Ye, 2001; Pinder et al., 1992; Tattersall and Ultsch, 2008; Wells, 2007). To minimize disturbance in our experiment, we used the motion-less four-point stance with the nictitating membranes drawn across the eyes as our proxy of brumation. Our experimental threshold temperatures were tightly associated with the corresponding field measurements both for the start (r = 0.97, t₂₈ = 22.26, P < 0.0001, λ = 0.04 [0.00, 0.47]) and end of the inactive state (r = 0.98, t₂₈ = 28.05, P < 0.0001, λ = 0.04 [0.00, 0.43]). Hence, we assessed the corresponding temperatures for all remaining species in the laboratory and estimated the brumation period based on the daily mean temperatures at the corresponding collection sites as retrieved from Chinese Meteorological Stations (http://www.lishi.tianqi.com) between 2012 and 2016.

We defined the brumation period as the number of consecutive days in each year that remained below this threshold. For simplicity we determined the active rather than brumation period, starting with the first day that the mean daily temperature rose above the activity threshold and remained there for at least five consecutive days, and ending with the last day before the temperature dropped below the activity threshold and remained there until the end of the calendar year. The brumation period then represented the difference between the activity period and the total number of days in each calendar year. Across these five years, the measured temperature thresholds yielded highly repeatable species-specific estimates of the number of days below the activity range (R = 0.95 [95%CI: 0.93−0.96]), as determined by the rpt function in the rptR package (Stoffel et al., 2017) across all 116 species (Fig. S1; Table S1). Further, across the 30 species that were examined both in the lab and the field (see above), these predicted brumation periods were also correlated with the observed brumation periods in the field (r = 0.96, t₂₈ = 18.03, P < 0.0001, λ = 0.05 [0.00, 0.48]; Fig. S8A), which themselves were highly repeatable between years within species (R = 0.98 [0.96−0.99]; Table S1).

Based on this data validation, we used for each species the mean brumation period predicted from our experimentally simulated temperature thresholds. However, to test for potential buffering effects of burrowing in the soil relative to the air temperatures reported by the meteorological stations, we also repeated these estimates by using more conservative thermal thresholds. Here, we restricted the putative brumation days to those with a reported air temperature of either 2°C or 4°C below the experimentally derived inactivity thresholds, simulating prolonged activity by seeking shelter in burrows. The 2°C threshold was based on a pilot study comparing direct measurements of air and burrow temperatures for four different burrows in each of five of our study species (burrow depths: 32.0 ± 3.2 to 121.0 ± 17.8 cm; Fig. S9). Across these species, the burrow-to-air temperature difference reached 1.03 ± 0.35°C to 2.45 ± 0.60°C in measurements around the peak of the brumation period (i.e., early January; Fig. S9). However, since these temporal snapshots were based on sites at relatively low elevation (≤320 m a.s.l.) due to accessibility of burrows during winter, we also used a second, more conservative buffer (4°C below activity range) for comparison. These temperature buffers shortened the predicted brumation periods to a varying degree between species (Fig. S8B); yet the predicted periods covaried strongly between the different temperature thresholds (all r > 0.90, t₁₁₄ > 21.96, P < 0.0001, all λ < 0.01).

Phylogeny reconstruction

To reconstruct the phylogeny, we obtained the sequences of three nuclear and six mitochondrial genes from GenBank (for accession numbers and sequence coverage see Data S4). The three nuclear genes included the recombination-activating gene 1 (RAG1), rhodopsin (RHOD) and tyrosinase (TYR). The six mitochondrial genes were cytochrome b (CYTB), cytochrome oxidase subunit I (COI), NADH dehydrogenase subunits 2 and 4 (ND2 and ND4), and the large and small subunits of the mitochondrial ribosome genes (12S/16S; omitting the adjacent tRNAs as they were difficult to align and represented only a small amount of data). We aligned the sequences by multi-sequence alignment (MUSCLE) in MEGA v.10.2.2 (Tamura et al., 2013) before comparing possible nucleotide substitution models. The best substitution model, as determined by the function modelTest() in the R (R Core Team, 2022) package phangorn (Schliep, 2011) based on the corrected Akaike Information Criterion, AICc, was GTR+Γ+I for all genes except RHOD, for which HKY+Γ had stronger support.

Using BEAUTi and BEAST v.1.10.4 (Suchard et al., 2018), we then constructed the phylogeny with unlinked substitution models, a relaxed uncorrelated log-normal clock, a Yule speciation process, and the best-supported nucleotide substitution models. We omitted time calibration due to a lack of fossil dates. We ran the Markov Chain Monte Carlo (MCMC) simulation for 55 million generations while sampling every 5,000^th tree with a 10% burn-in. Most effective sample size (ESS) values by far exceeded 375 (i.e. all well above the recommended threshold of 200) for all but two tree statistics in the program Tracer v.1.7.2 (Rambaut et al., 2018), thus indicating satisfying convergence of the Bayesian chain and adequate model mixing. Finally, we generated a maximum clade credibility tree with mean node heights and a 10% burn-in using TreeAnnotator v.1.10.4 (Suchard et al., 2018), presented in Fig. S10.

Breeding conditions

To test if a prolonged brumation period reduces the time available for reproduction, thereby changing the level of competition over mates and fertilizations (Lüpold et al., 2017), we extracted the start and end dates of the breeding season from our field notes of concurrent studies on species-specific life histories. These data were available for 43 of our species (Data S3). We used dates when the first and last clutches were observed in focal ponds as a proxy of mating activity, given that males release their sperm during oviposition in these external fertilizers. For each species, dates from at least two years were combined and averaged to obtain the mean duration of the breeding season.

We further recorded whether dense mating aggregations are typically observed in these species. We have previously shown that larger mating clusters, with multiple males clasping the same females, have a significant effect on the evolution of testis size due to the resulting competition among sperm for fertilization (Lüpold et al., 2017). Here, we had no detailed data on the sizes of aggregations and so were only able to code the typical presence or absence of aggregations as a binary variable (Data S3).

Finally, we used our direct estimates of species-specific population densities from our previous study (Lüpold et al., 2017) to test whether a shorter breeding season results in denser breeding populations. Although population density is a more direct measure than the occurrence of aggregations, such data were available for only eight of our species, each based on multiple populations per species (Lüpold et al., 2017). All these data were not necessarily derived from the same years or populations of our main dataset, but given the within-species repeatability in breeding populations (Lüpold et al., 2017) and in the duration of the breeding season (R = 0.88 [0.79−0.93]; Table S1), these differences should be relatively small compared to the interspecific variation and mostly introduce random noise.

Data analyses

General methods

We conducted all statistical analyses in R v.4.2.0 (R Core Team, 2022), using log-transformed data for all phenotypic traits, and for the CV in temperature among the ecological variables. To account for non-independence of data due to common ancestry (Freckleton et al., 2002; Pagel, 1999), we conducted phylogenetic generalized least-squares (PGLS) or phylogenetic logistic regressions (e.g., for occurrence of breeding aggregations), using the R package phylolm (Ho and Ané, 2014) and our reconstructed phylogeny. To account for variation around the species means, we bootstrapped for each model (at 100 fitted replicates) the standardized regression coefficients along with the phylogenetic scaling parameter λ and calculated their corresponding 95% confidence intervals. The λ values indicate phylogenetic independence near zero and strong phylogenetic dependence near one (Freckleton et al., 2002).

Unless stated otherwise, all PGLS models focusing on the relative mass of tissues as the response included snout-vent length (SVL) as a covariate in addition to the focal predictor variable(s). We chose SVL instead of body mass because it is the commonly used measure of body size in anurans and independent of seasonal fluctuations in tissues such as body fat, testes, or limb muscles. One exception, however, was the analysis of phylogenetically informed allometric relationships, for which we cubed SVL such that a slope of 1 equaled unity (isometry). For these allometric relationships we calculated ordinary (generalized) least-squares rather than reduced major-axis regressions, because their greater sensitivity to changes in the steepness, but lower sensitivity to changes in scatter, capture allometric slopes more adequately (Kilmer and Rodríguez, 2017).

Pairwise correlations between tissues

To examine the covariation between different tissues across species, we first calculated pairwise partial correlations controlling for SVL and phylogeny. To this end, we calculated the phylogenetic trait variance-covariance matrix between the pairs of focal variable and SVL using the function phyl.vcv() in phytools (Revell, 2012) with λ = 1 (i.e. Brownian motion), which we then scaled into a correlation matrix using cov2cor() in the stats package (R Core Team, 2022). Using the resulting correlation coefficients r_xy, r_xz, and r_yz, respectively, we then calculated the partial correlation coefficient r_xy.z between the x and y variables of interest while accounting for SVL (z) following Crawley’s (2007) equation: r_xy.z = with the associated t-statistics and 95% confidence intervals converted using standard conversion and the package effectsize (Ben-Shachar et al., 2020), respectively.

Multivariate allocation patterns

Since pairwise correlations do not necessarily capture more complex, multivariate allocation patterns, we used two additional approaches to explore how tissue sizes varied relative to others: A compositional analysis and a principal component analysis. In both analyses, we focused on those four tissues that covaried with brumation duration or deviated from proportionate scaling with body size: brain, body fat, testes, and hindlimb muscles.

Compositional analyses

In the first approach, we used the function acomp() in the R package compositions (van den Boogaart and Tolosana-Delgado, 2008) to partition total body mass of each species into a five-variable Aitchison composition in a logistic geometry (van den Boogaart and Tolosana-Delgado, 2013), consisting of the proportional representation of the four focal tissues and the remaining body mass combined. Since the focal tissues constituted a size-independent fraction of the total body, the closed composition of this combined mass should be unbiased relative to body size but can instead reveal differential contributions of the four tissues to their total in a multivariate context (Aitchison, 1982; Muldowney et al., 2001; van den Boogaart and Tolosana-Delgado, 2013). For phylogenetic correlations between these variables following the description above, we used centered log ratios obtained by the function clr() in the same package, which maintains the original variable structure. However, owing to the reliance on a full rank of the covariance in multivariate analyses, we used the ilr() function to project the D-part composition isometrically to a D−1 dimensional simplex (Aitchison, 1982), essentially representing the log ratios between the D parts. This multivariate object we subjected to a phylogenetic multivariate regression against brumation duration using the functions mvgls() and manova.gls() in the package mvMORPH (Clavel et al., 2015). For interpretation in the context of the original variable space, we back-transformed the coefficients using the ilrInv() function in compositions.

Phylogenetic principal component analysis

In addition to this compositional data analysis, we also peerformed a phylogenetically informed principal component analysis on the same focal tissues as log-transformed species means, using the phyl.pca() function of the package phytools (Revell, 2012). Here, we primarily focused on the directions of the loading vectors relative to the principal components and one another to glean information on the correlations between the original variables in the principal component space.

Phylogenetic confirmatory path analyses

We performed two phylogenetic confirmatory path analyses (Gonzalez-Voyer and von Hardenberg, 2014; von Hardenberg and Gonzalez-Voyer, 2013) based on pre-specified candidate structural equation models. In the first analysis, we explored the direct and indirect effects of climatic variables on the duration of the breeding season or formation of breeding aggregations. Using eight candidate structural equation models, we essentially tested if breeding aggregations were the result of variation in latitude as a proxy of seasonality, either directly or mediated by the duration of brumation and/or of the breeding season. To avoid overparameterization given that this analysis was based on only 43 species, we did not include testis mass (and so necessarily also SVL to control for body size) as additional variables in the same path models. Rather, to test for links between relative testis size and breeding parameters, we conducted separate directional tests of trait evolution (see below).

A second path analysis aimed to disentangle the different interrelationships between traits that could ultimately mediate the effect of brumation duration on brain and reproductive evolution. Here, we used 28 pre-specified structural equation models to test for direct and indirect links between the brumation period and the four main tissues. As brumation is determined primarily by the environment, we focused on models that explained variation in tissue investments rather than, say, brain size affecting brumation patterns. The effects of brumation on the four tissues were direct or indirect, for example mediated by the digestive tract (resource acquisition) and/or body fat (resource storage), thus providing context for the cognitive buffer, expensive brain, expensive tissue, energy trade-off and fat−brain trade-off hypotheses (Allman et al., 1993; Heldstab et al., 2016; Isler and van Schaik, 2006, 2009; Navarrete et al., 2011). Further, we allowed brain size to affect testis size and vice versa (i.e., expensive sexual tissue hypothesis; Pitnick et al. 2006), and included models in which these organs explained variation in the digestive tract or fat tissue instead of being affected by them (e.g., selection on brain or testis size might mediate selection for greater resource availability rather than resources influencing brain or testis evolution). Across the 28 candidate models, we tested different combinations of these predictions, with traits being explained by single or multiple predictors, or having individual or shared effects on other traits. Using the R package phylopath (van der Bijl, 2018), we examined the conditional independencies of each model, ranked all candidate models based on their C-statistic Information Criterion (CICc), and then averaged the coefficients of the models with ΔCICc ≤ 2 from the top model (von Hardenberg and Gonzalez-Voyer, 2013).

Directional tests of trait evolution

Since the first path analysis did not involve testis size to avoid overparametization, we separately tested for correlated evolution using directional tests of trait evolution (Pagel, 1994; Revell, 2012). One limitation of these models is that they rely on evolutionary transitions between binary states in each trait. Hence, we considered positive residuals of a log-log regression between testis mass and SVL as ‘relatively large testes’ and negative residuals as ‘relatively small testes.’ For the duration of the breeding season, we similarly split the distribution based on the mean duration, whereas aggregation formation was already coded as present or absent. Based on (the weight of) the Akaike Information Criterion (AIC), we then tested if changes in relative testis size and breeding parameters, respectively, were unilaterally dependent, mutually dependent, or independent (Pagel, 1994), using the fitPagel() function in the phytools package (Revell, 2012) with “fitDiscrete” as the optimization method and allowing all rates to differ (i.e., “ARD” model).

Acknowledgements

We thank K. Isler and C. van Schaik for helpful comments on an earlier draft of our manuscript, and J. Santin for insightful discussions. Financial support was provided by the National Natural Sciences Foundation of China (31772451 and 31970393 to WBL) and by the Swiss National Science Foundation (PP00P3_170669 and PP00P3_202662 to SL).

Data and code availability

All data and R scripts related to this article are deposited in the Figshare Data Repository (doi: 10.6084/m9.figshare.21078052).

Ethics statement

The specimens used in this study were collected with permission from the China West Normal University Ethical Committee for Animal Experiments (CWNU-20001), and the experimental protocols adhered to the current laws of China concerning animal experimentation.