How hibernation in frogs drives brain and reproductive evolution in opposite directions

Environmental seasonality can promote the evolution of larger brains through cognitive and behavioral flexibility but can also hamper it when temporary food shortage is buffered by stored energy. Multiple hypotheses linking brain evolution with resource acquisition and allocation have been proposed for warm-blooded organisms, but it remains unclear how these extend to cold-blooded taxa whose metabolism is tightly linked to ambient temperature. Here, we integrated these hypotheses across frogs and toads in the context of varying brumation (hibernation) durations and their environmental correlates. We showed that protracted brumation covaried negatively with brain size but positively with reproductive investment, likely in response to brumation-dependent changes in the socio-ecological context and associated selection on different tissues. Our results provide novel insights into resource allocation strategies and possible constraints in trait diversification, which may have important implications for the adaptability of species under sustained environmental change.


Introduction 35
Seasonal food scarcity challenges animal energy budgets. A positive or less negative energy 36 balance across seasons can be achieved by a more constant net energy intake than predicted 37 solely by food abundance (Sol, 2009), or by investing less in costly organs (Heldstab et al.,38 8 brain size appears to be more constrained than that of other organs when selection favors 162 larger body size. 163 Given these distinct allocation patterns between brain and other tissues, we next tested whether 164 extended brumation constrains brain size evolution, akin to suggestions for mammals based on 165 presence/absence of hibernation (Heldstab et al., 2018). In PGLS models, absolute brain size  Table S7). This trend was 171 distinct from body size changes in response to brumation, as SVL remained independent of 172 brumation duration (Table S5). On average, brumating species tended to possess relatively 173 smaller brains compared to those less likely to experience prolonged brumation (Table S8), 174 aligning with mammalian studies utilizing hibernation presence/absence (Heldstab et al., 2018). 175 This trend extended to those 91 species categorized as brumating for some period (Table S9), 176 reinforcing the link between brumation and brain evolution beyond coarse binary classification. 177 This pattern further persisted when recalculating brumation periods using conservative 178 thresholds of 2°C or 4°C below their experimentally derived thresholds, simulating shelter 179 buffering (details and validation in Material and Methods; Table S10). These patterns generally held true when using the presence/absence of brumation (Table S8), 9 buffered temperature fluctuations (Tables S10), or excluding the 25 species unlikely to brumate 186 (except for the non-significant effect on body fat; Table S9). Other tissue sizes remained 187 independent of brumation (Tables S7 to S10) and showed not significant changes between 188 sampling periods (Fig. 2). 189 Comparing pre-and post-brumation males, averaged across all species, we found a roughly 190 50% reduction in fat tissue and a 100% increase in testis size, indicating resource depletion and 191 testicular regrowth during brumation, respectively (Fig. 1B). Only brain size deviated from zero 192 among the remaining tissues, but this change was minimal compared to the body fat and testes, 193 and within the range of many unchanged tissues. Thus, the biological significance of this 194 putative increase in brain size during brumation remains uncertain, possibly reflecting general 195 differences between sampled individuals. We will thus refrain from further interpretation. Among To examine whether the increase in relative testis size with prolonged brumation might be 200 mediated by a shorter, more synchronized mating season (Wells, 2007), we tested for links 201 between brumation duration and different breeding parameters. Prolonged brumation notably 202 shortened the breeding season (r = −0.57, t41 = −4.47, P < 0.001, λ = 0.00 [0.00, 0.38]; Fig.  203 S2A), exerting a stronger effect than climatic variables (Table S11), particularly when 204 considered together (Table S12). Hence, the effect of these climatic variables may be mediated 205 by brumation. Furthermore, a shorter breeding season increased the probability of dense 206 breeding aggregations (phylogenetic logistic regression: N = 42, z = −3.03, P = 0.002, α = 0.02; 207  Table S13). Finally, when combining these data with published data on the 211 density of breeding populations (Lüpold et al., 2017; N = 8 species overlapping), a trend 212 emerged toward higher mean population densities in species with shorter breeding seasons (r = 213 −0.69, t6 = −2.37, P = 0.06, λ = 0.00 [0.00, 1.00]), albeit with a small sample size (Fig. S2C). 214 To explore possible causal links between breeding parameters and relative testis size, we 215 employed directional tests of trait evolution (Pagel, 1994;Revell, 2012), assessing whether 216 changes in two binary traits are unilaterally or mutually dependent, or independent (Pagel, 217 1994). To this end, we converted our continuous to binary variables (see Material and Methods). 218 Relating small/large testes to short/long breeding seasons, the best-supported scenario based 219 on the Akaike Information Criterion (AIC) was independent evolution. Yet, the model with 220 changes in relative testis size dependent on those in the breeding season found similar support 221 (ΔAIC = 0.83, wAIC = 0.33 compared to independent model with wAIC = 0.50; Fig. S4A), differing 222 from the remaining two models with clearly higher AIC scores (ΔAIC ≥ 3.40, wAIC ≤ 0.09). In 223 another directional test, between relative testis size and aggregation formation, the independent 224 model was the best-supported (wAIC = 0.66), followed by the scenario of increased relative testis 225 size in response to aggregation formation (Fig. S4B), albeit above a ΔAIC cut-off of 2 (ΔAIC = 226 2.38, wAIC = 0.20), with the remaining models being clearly less supported (ΔAIC ≥ 3.64, wAIC ≤ 227 0.11). Hence, it is at least possible that breeding conditions mediate the positive relationship 228 between brumation duration and relative testis size in our relatively small sample of species, 229 whilst a response of breeding conditions to variation in relative testis size is clearly rejected by 230 our analyses. 231 232 Covariation between tissues. As all tissues rely on the same finite resources, their responses 233 to brumation are likely interconnected. Pairwise partial correlations, controlling for SVL and 234 phylogeny, revealed positive covariation or non-significant associations among all tissue 235 masses (Fig. S5). As such, our data do not support the expensive tissue (Aiello and Wheeler, 236 1995) or the more general energy trade-off hypotheses (Isler and van Schaik, 2006), which 237 predict brain size trade-offs with the digestive tract or other costly organs, respectively. 238 However, since brain size differed from fat, hindlimb muscles and testes in allometric However, the resource distribution among these tissues varied considerably between species. 247 Pairwise correlations between the four focal tissues, transformed to centered log ratios (van den 248 Boogaart and Tolosana-Delgado, 2008) and controlling for phylogeny, showed brain mass to 249 covary negatively with fat and testis mass, whilst testis mass covaried negatively with hindlimb 250 muscle mass but not with body fat (Fig. 3A-D; Table S14). 251 To further examine the effect of brumation duration on all five variables simultaneously, we 252 the expected coefficient of 0.20 if brumation had no effect, whilst brain mass, hindlimb muscles 258 and the rest of the body had lower coefficients (0.16, 0.17, and 0.18, respectively; also see Fig.  259 3). 260 Finally, a phylogenetically informed principal component analysis (Revell, 2012) confirmed the 261 negative associations of brain size with body fat and testis mass, and that of testis mass with 262 hindlimb muscles (Fig. 3E,F, Table S15). Here, the first three principal components (PC1 to 263 PC3) explained 84.7%, 8.0% and 5.1% of total variance, respectively. Although PC2 and PC3 264 explained a relatively small proportion of the total variance, they separated the different tissues.  Table S16). The averaged model (Fig. 4)  Our study on anurans, with validated brumation periods and direct measures of expensive 285 tissues, provides novel insights into brain and reproductive evolution in 'cold-blooded' 286 organisms exposed to environmental seasonality. We found that species with longer brumation 287 exhibited relatively smaller brains and allocated greater fat reserves primarily in reproduction, 288 possibly due to the shorter breeding season with its socio-ecological implications. 289 290

Environmental and resource considerations in brain evolution 291
We demonstrated that species in cooler and more seasonal climates were more cold tolerant, 292 thereby likely optimizing their active period. Yet, low temperatures impair foraging and digestion 293 in ectotherms (e.g., Fontaine et al., 2018;Riddle, 1909), such that high seasonality may lead to 294 longer brumation and smaller brains. These results confirm that, unlike birds (Sol, 2009)  Supporting a large brain may not be sustainable without continued resource intake, or larger 298 brains could be less tolerant to hypoxic conditions during brumation (Sukhum et al., 2016). 299 However, selection for relatively larger brains may also simply be stronger in species with longer 300 active (and short brumation) periods owing to extended cognitive benefits such as predator In pairwise comparisons, the relative sizes of the tissues examined here, including the brain, 304 were generally positively correlated. These results reject both the expensive tissue and energy 305 trade-off hypotheses (Aiello and Wheeler, 1995;Isler and van Schaik, 2006), which predict 306 trade-offs of brain size with the size of the digestive tract or other costly organs, respectively. 307 This lack of support in anurans aligns with a previous report in mammals (Navarrete et al.,308 2011) despite their smaller brains and vastly different ecology and physiology, including a lower 309 metabolic rate and largely lacking physiological thermoregulation. When focusing jointly on the 310 four tissues (brain, body fat, testes, hindlimb muscles) that covaried with brumation duration, 311 however, relative brain size covaried negatively with the relative mass of both fat tissue and this notion, arboreal species tended to be leaner compared to (semi)aquatic or terrestrial 323 species (Table S17), controlling for brumation duration and relative brain size, both of which we 324 had shown to covary with body fat (Figs. 2 and 3). 325 326

Brumation and testis evolution 327
Species with prolonged brumation also had relatively smaller hindleg muscles and larger testes. 328 The negative relationship between hindleg muscle mass and brumation duration may be linked 329 to more movement during a longer active period, including predator evasion 330 Marchisin and Anderson, 1978). Larger testes may result from a shorter breeding season, 331 leading to denser and more synchronous mating activity (Wells, 2007), as suggested by our 332 In addition to the average size of the testes, their seasonal change also varied with the 339 brumation period. Seasonally breeding anurans regress and regrow their testes between mating 340 seasons (Ogielska and Bartmańska, 2009). Non-brumating species can use energy uptake to 341 compensate for testicular recrudescence, whilst those with a short breeding season after a 342 prolonged inactive period depend on the stored fat to regrow their testes before or immediately 343 after emergence from their hibernaculum. Hence, resources are diverted away from the brain 344 and other organs, especially in species such as Brachytarsophrys spp., in which the fully 345 developed testes combined weigh 12−14 times more than the brain (Data S1). 346 347 Brain-testis trade-off 348 brain size and revealed its direct and indirect positive effects on relative testis size, mediated by 350 the amount of adipose tissue, which responded to variation in the inactive period (energetic 351 demand) and the size of the digestive tract (energy uptake). That body fat did not contribute to 352 brain size evolution in this more comprehensive analysis compared to pairwise correlations 353 suggests that the fat−brain trade-off may not be direct. Rather, longer brumation, and thus a 354 short active period, may enhance selection on fat storage for testicular investments in addition 355 to starvation avoidance, while reducing selection for larger brains due to a shifted balance 356 between cognitive benefits and energetic costs (Fig. 3). acquisition and allocation patterns, and it remains to be seen to what extent variation in the 388 adaptability, and thus resilience, between species exposed to environmental change is 389 attributable to such competing needs between investments and species-specific constraints. shortly before entering their hibernacula (Data S1 and S2). For each species, we sampled all 397 males at a single location in southern and western China with known longitude, latitude, and 398 elevation (Data S3). Upon transfer to the laboratory, we sacrificed the individuals by single-399 pithing, measured their snout-vent length (SVL) to the nearest 0.01 mm with calipers and then 400 preserved them in 4% phosphate-buffered formalin for tissue fixation. 401 After two months of preservation, we weighed each complete specimen to the nearest 0.1 mg 402 using an electronic balance to obtain body mass before dissecting them following a strict 403 protocol. We separately extracted the brain, heart, liver, lungs, kidneys, spleen, digestive tract, 404 testes, limb muscles, and fat stores, cleaned these tissues and immediately weighed them to 405 the nearest 0.1 mg with an electronic balance. We additionally measured the length of the 406 digestive tract to the nearest 0.01 mm using calipers. We excluded emaciated individuals or 407 those exhibiting visible organ pathologies from our analyses. 408 409

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

Brumation period 418
One way that anurans can physiologically respond to seasonality is by adjusting their thermal 419 sensitivity and thus brumation period (Wells, 2007), which in turn could directly or indirectly 420 affect the evolution of brain size (Heldstab et al., 2018). Hence, we estimated the brumation We defined the brumation period as the number of consecutive days in each year that remained 444 below this threshold. For simplicity we determined the active rather than brumation period, 445 starting with the first day that the mean daily temperature rose above the activity threshold and  Table S1). 457 Based on this data validation, we used for each species the mean brumation period predicted 458 from our experimentally simulated temperature thresholds. However, to test for potential 459 buffering effects of burrowing in the soil relative to the air temperatures reported by the 460 meteorological stations, we also repeated these estimates by using more conservative thermal 461 thresholds. Here, we restricted the putative brumation days to those with a reported air 462 temperature of either 2°C or 4°C below the experimentally derived inactivity thresholds, 463 simulating prolonged activity by seeking shelter in burrows. The 2°C threshold was based on a 464 pilot study comparing direct measurements of air and burrow temperatures for four different 465 burrows in each of five of our study species (burrow depths: 32.0 ± 3.2 to 121.0 ± 17.8 cm; Fig.  466 S9). Across these species, the burrow-to-air temperature difference reached 1.03 ± 0.35°C to 467 2.45 ± 0.60°C in measurements around the peak of the brumation period (i.e., early January; 468 Fig. S9). However, since these temporal snapshots were based on sites at relatively low 469 elevation (≤320 m a.s.l.) due to accessibility of burrows during winter, we also used a second, 470 more conservative buffer (4°C below activity range) for comparison. These temperature buffers 471 shortened the predicted brumation periods to a varying degree between species (Fig. S8B); yet 472 the predicted periods covaried strongly between the different temperature thresholds (all r > 473 0.90, t114 > 21.96, P < 0.0001, all λ < 0.01). 474 475

Phylogeny reconstruction 476
To reconstruct the phylogeny, we obtained the sequences of three nuclear and six mitochondrial studies on species-specific life histories. These data were available for 43 of our species (Data 504 S3). We used dates when the first and last clutches were observed in focal ponds as a proxy of 505 mating activity, given that males release their sperm during oviposition in these external 506 fertilizers. For each species, dates from at least two years were combined and averaged to 507 obtain the mean duration of the breeding season. 508 We further recorded whether dense mating aggregations are typically observed in these 509 species. We have previously shown that larger mating clusters, with multiple males clasping the 510 same females, have a significant effect on the evolution of testis size due to the resulting 511 competition among sperm for fertilization (Lüpold et al., 2017). Here, we had no detailed data on 512 the sizes of aggregations and so were only able to code the typical presence or absence of 513 aggregations as a binary variable (Data S3). 514 Finally, we used our direct estimates of species-specific population densities from our previous 515 study (Lüpold et al., 2017) to test whether a shorter breeding season results in denser breeding 516 populations. Although population density is a more direct measure than the occurrence of 517 aggregations, such data were available for only eight of our species, each based on multiple 518 populations per species (Lüpold et al., 2017). All these data were not necessarily derived from 519 the same years or populations of our main dataset, but given the within-species repeatability in 520  Table S1), these differences should be relatively small compared to the 522 interspecific variation and mostly introduce random noise. 523 524

Data analyses 525
General methods. We conducted all statistical analyses in R v.4.2.0 (R Core Team, 2022), 526 using log-transformed data for all phenotypic traits, and for the CV in temperature among the 527 ecological variables. To account for non-independence of data due to common ancestry 528 (Freckleton et al., 2002;Pagel, 1999), we conducted phylogenetic generalized least-squares 529 (PGLS) or phylogenetic logistic regressions (e.g., for occurrence of breeding aggregations), 530 using the R package phylolm (Ho and Ané, 2014) and our reconstructed phylogeny. To account 531 for variation around the species means, we bootstrapped for each model (at 100 fitted 532 replicates) the standardized regression coefficients along with the phylogenetic scaling 533 parameter λ and calculated their corresponding 95% confidence intervals. The λ values indicate 534 phylogenetic independence near zero and strong phylogenetic dependence near one 535 (Freckleton et al., 2002). 536 Unless stated otherwise, all PGLS models focusing on the relative mass of tissues as the 537 response included snout-vent length (SVL) as a covariate in addition to the focal predictor 538 variable(s). We chose SVL instead of body mass because it is the commonly used measure of 539 body size in anurans and independent of seasonal fluctuations in tissues such as body fat, 540 testes, or limb muscles. One exception, however, was the analysis of phylogenetically informed 541 allometric relationships, for which we cubed SVL such that a slope of 1 equaled unity (isometry). 542 For these allometric relationships we calculated ordinary (generalized) least-squares rather than 543 reduced major-axis regressions, because their greater sensitivity to changes in the steepness, 544 but lower sensitivity to changes in scatter, capture allometric slopes more adequately (Kilmer 545 and Rodrí guez, 2017). 546 Pairwise correlations between tissues. To examine the covariation between different tissues 547 across species, we first calculated pairwise partial correlations controlling for SVL and 548 phylogeny. To this end, we calculated the phylogenetic trait variance-covariance matrix between 549 the pairs of focal variable and SVL using the function phyl.vcv() in phytools (Revell, 2012) with λ 550 = 1 (i.e. Brownian motion), which we then scaled into a correlation matrix using cov2cor() in the 551 stats package (R Core Team, 2022). Using the resulting correlation coefficients rxy, rxz, and ryz, 552 respectively, we then calculated the partial correlation coefficient rxy.z between the x and y 553 , with the associated t-statistics and 95% confidence intervals converted using 555 standard conversion ( = √ 1− 2 ) and the package effectsize (Ben-Shachar et al., 2020), 556

respectively. 557
Multivariate allocation patterns. Since pairwise correlations do not necessarily capture more 558 complex, multivariate allocation patterns, we used two additional approaches to explore how 559 tissue sizes varied relative to others: A compositional analysis and a principal component 560 analysis. In both analyses, we focused on those four tissues that covaried with brumation 561 duration or deviated from proportionate scaling with body size: brain, body fat, testes, and 562 hindlimb muscles. tissues and the remaining body mass combined. Since the focal tissues constituted a size-568 independent fraction of the total body, the closed composition of this combined mass should be 569 unbiased relative to body size but can instead reveal differential contributions of the four tissues 570 to their total in a multivariate context (Aitchison, 1982;Muldowney et al., 2001; van den 571 Boogaart and Tolosana-Delgado, 2013). For phylogenetic correlations between these variables 572 following the description above, we used centered log ratios obtained by the function clr() in the 573 same package, which maintains the original variable structure. However, owing to the reliance 574 on a full rank of the covariance in multivariate analyses, we used the ilr() function to project the 575 D-part composition isometrically to a D−1 dimensional simplex (Aitchison, 1982), essentially 576 representing the log ratios between the D parts. This multivariate object we subjected to a 577 phylogenetic multivariate regression against brumation duration using the functions mvgls() and species, we did not include testis mass (and so necessarily also SVL to control for body size) as 596 additional variables in the same path models. Rather, to test for links between relative testis size 597 and breeding parameters, we conducted separate directional tests of trait evolution (see below). 598 A second path analysis aimed to disentangle the different interrelationships between traits that 599 could ultimately mediate the effect of brumation duration on brain and reproductive evolution. 600 Here, we used 28 pre-specified structural equation models to test for direct and indirect links 601 between the brumation period and the four main tissues. As brumation is determined primarily 602 by the environment, we focused on models that explained variation in tissue investments rather 603 than, say, brain size affecting brumation patterns. The effects of brumation on the four tissues 604 were direct or indirect, for example mediated by the digestive tract (resource acquisition) and/or 605 body fat (resource storage), thus providing context for the cognitive buffer, expensive brain, 606 expensive tissue, energy trade-off and fat−brain trade-off hypotheses (Allman et al., 1993;607 Heldstab et al., 2016;Isler andvan Schaik, 2006, 2009;Navarrete et al., 2011). Further, we 608 allowed brain size to affect testis size and vice versa (i.e., expensive sexual tissue hypothesis; 609 Pitnick et al. 2006), and included models in which these organs explained variation in the 610 digestive tract or fat tissue instead of being affected by them (e.g., selection on brain or testis 611 size might mediate selection for greater resource availability rather than resources influencing 612 brain or testis evolution). Across the 28 candidate models, we tested different combinations of 613 these predictions, with traits being explained by single or multiple predictors, or having individual 614 or shared effects on other traits. Using the R package phylopath (van der Bijl, 2018), we 615 examined the conditional independencies of each model, ranked all candidate models based on 616 their C-statistic Information Criterion (CICc), and then averaged the coefficients of the models 617 with ΔCICc ≤ 2 from the top model (von Hardenberg and Gonzalez-Voyer, 2013). 618 Directional tests of trait evolution. Since the first path analysis did not involve testis size to 619 avoid overparametization, we separately tested for correlated evolution using directional tests of 620 trait evolution (Pagel, 1994;Revell, 2012). One limitation of these models is that they rely on 621 evolutionary transitions between binary states in each trait. Hence, we considered positive 622 residuals of a log-log regression between testis mass and SVL as 'relatively large testes' and 623 negative residuals as 'relatively small testes.' For the duration of the breeding season, we 624 similarly split the distribution based on the mean duration, whereas aggregation formation was 625 already coded as present or absent. Based on (the weight of) the Akaike Information Criterion 626 (AIC), we then tested if changes in relative testis size and breeding parameters, respectively, 627 were unilaterally dependent, mutually dependent, or independent (Pagel, 1994), using the 628