Ecology

Speakman, J. R. The cost of living: Field metabolic rates of small mammals. Adv. Ecol. Res. 30, 177–297 (1999).Article 

Google Scholar 
Brown, J. H., Gillooly, J. F., Allen, A. P., Savage, V. M. & West, G. B. Toward a metaboolic theory of ecology. Ecology 85(7), 1771–1789. https://doi.org/10.1890/03-9000 (2004).Article 

Google Scholar 
Larivée, M. L., Boutin, S., Speakman, J. R., McAdam, A. G. & Humphries, M. M. Associations between over-winter survival and resting metabolic rate in juvenile North American red squirrels. Funct. Ecol. 24(3), 597–607. https://doi.org/10.1111/j.1365-2435.2009.01680.x (2010).Article 

Google Scholar 
Corp, N., Gorman, M. L. & Speakman, J. R. Seasonal variation in the resting metabolic rate of male wood mice Apodemus sylvaticus from two contrasting habitats 15 km apart. J. Comp. Physiol. B 167(3), 229–239. https://doi.org/10.1007/s003600050069 (1997).Article 
CAS 

Google Scholar 
Lehto Hürlimann, M., Martin, J. G. A. & Bize, P. Evidence of phenotypic correlation between exploration activity and resting metabolic rate among populations across an elevation gradient in a small rodent species. Behav. Ecol. Sociobiol. 73(9), 131. https://doi.org/10.1007/s00265-019-2740-6 (2019).Article 

Google Scholar 
Reher, S., Rabarison, H., Montero, B. K., Turner, J. M. & Dausmann, K. H. Disparate roost sites drive intraspecific physiological variation in a Malagasy bat. Oecologia 198(1), 35–52. https://doi.org/10.1007/s00442-021-05088-2 (2022).Article 
ADS 

Google Scholar 
McDonald, R. I. et al. Research gaps in knowledge of the impact of urban growth on biodiversity. Nat. Sustain. https://doi.org/10.1038/s41893-019-0436-6 (2019).Article 

Google Scholar 
Shochat, E., Warren, P. S., Faeth, S. H., McIntyre, N. E. & Hope, D. From patterns to emerging processes in mechanistic urban ecology. Trends Ecol. Evol. 21(4), 186–191. https://doi.org/10.1016/j.tree.2005.11.019 (2006).Article 

Google Scholar 
United Nations, Department of Economic and Social Affairs, Population Division. World Urbanization Prospects 2018: Highlights. https://population.un.org/wup/Publications/ (2018).Alberti, M. et al. The complexity of urban eco-evolutionary dynamics. Bioscience 70(9), 772–793. https://doi.org/10.1093/biosci/biaa079 (2020).Article 

Google Scholar 
Birnie-Gauvin, K., Peiman, K. S., Gallagher, A. J., de Bruijn, R. & Cooke, S. J. Sublethal consequences of urban life for wild vertebrates. Environ. Rev. 24(4), 416–425. https://doi.org/10.1139/er-2016-0029 (2016).Article 

Google Scholar 
Diamond, S. E. & Martin, R. A. Physiological adaptation to cities as a proxy to forecast global-scale responses to climate change. J. Exp. Biol. 224((Suppl_1)), jeb22336. https://doi.org/10.1242/jeb.229336 (2021).Article 

Google Scholar 
Grimm, N. B. et al. Global change and the ecology of cities. Science 319(5864), 756–760. https://doi.org/10.1126/science.1150195 (2008).Article 
ADS 
CAS 

Google Scholar 
McDonnell, M. J. & Pickett, S. T. Ecosystem structure and function along urban-rural gradients: An unexploited opportunity for ecology. Ecology 71(4), 1232–1237. https://doi.org/10.2307/1938259 (1990).Article 

Google Scholar 
Francis, R. A. & Chadwick, M. A. What makes a species synurbic?. Appl. Geogr. 32(2), 514–521. https://doi.org/10.1016/j.apgeog.2011.06.013 (2012).Article 

Google Scholar 
Luniak, M. Synurbization–adaptation of animal wildlife to urban development. In Proc. 4th Int. Symposium Urban Wildl. Conserv (Tucson, University of Arizona, 2004).Coogan, S. C. P., Raubenheimer, D., Zantis, S. P. & Machovsky-Capuska, G. E. Multidimensional nutritional ecology and urban birds. Ecosphere 9(4), e02177. https://doi.org/10.1002/ecs2.2177 (2018).Article 

Google Scholar 
Lowry, H., Lill, A. & Wong, B. B. Behavioural responses of wildlife to urban environments. Biol. Rev. Camb. Philos. Soc. 88(3), 537–549. https://doi.org/10.1111/brv.12012 (2013).Article 

Google Scholar 
Łopucki, R., Klich, D., Ścibior, A. & Gołębiowska, D. Hormonal adjustments to urban conditions: Stress hormone levels in urban and rural populations of Apodemus agrarius. Urban Ecosyst. 22(3), 435–442. https://doi.org/10.1007/s11252-019-0832-8 (2019).Article 

Google Scholar 
McCleery, R. in Urban mammals in Urban Ecosystem Ecology (eds. Aitkenhead-Peterson, J., Volder, A.) 87–102 (American Society of Agronomy, 2010). https://doi.org/10.2134/agronmonogr55.c52010Uchida, K., Suzuki, K., Shimamoto, T., Yanagawa, H. & Koizumi, I. Seasonal variation of flight initiation distance in Eurasian red squirrels in urban versus rural habitat. J. Zool. 298(3), 225–231. https://doi.org/10.1111/jzo.12306 (2016).Article 

Google Scholar 
Kleerekoper, L., van Esch, M. & Salcedo, T. B. How to make a city climate-proof, addressing the urban heat island effect. Resour. Conserv. Recyl. 64, 30–38. https://doi.org/10.1016/j.resconrec.2011.06.004 (2012).Article 

Google Scholar 
Pickett, S. T. et al. Urban ecological systems: Scientific foundations and a decade of progress. J. Environ. Manag. 92(3), 331–362. https://doi.org/10.1016/j.jenvman.2010.08.022 (2011).Article 
CAS 

Google Scholar 
Rizwan, A. M., Dennis, L. Y. & Chunho, L. A review on the generation, determination and mitigation of Urban Heat Island. J. Environ. Sci. 20(1), 120–128 (2008).Article 
CAS 

Google Scholar 
Isaksson, C. Urban ecophysiology: Beyond costs, stress and biomarkers. J. Exp. Biol. 223(22), jeb203794. https://doi.org/10.1242/jeb.203794 (2020).Article 

Google Scholar 
Miles, L. S., Carlen, E. J., Winchell, K. M. & Johnson, M. T. J. Urban evolution comes into its own: Emerging themes and future directions of a burgeoning field. Evol. Appl. 14(1), 3–11. https://doi.org/10.1111/eva.13165 (2020).Article 

Google Scholar 
Gavett, A. P. & Wakeley, J. S. Blood constituents and their relation to diet in urban and rural house sparrows. Condor 88(3), 279–284. https://doi.org/10.2307/1368873 (1986).Article 

Google Scholar 
Murray, M. et al. Greater consumption of protein-poor anthropogenic food by urban relative to rural coyotes increases diet breadth and potential for human-wildlife conflict. Ecography 38(12), 1235–1242. https://doi.org/10.1111/ecog.01128 (2015).Article 

Google Scholar 
Pollock, C. J., Capilla-Lasheras, P., McGill, R. A. R., Helm, B. & Dominoni, D. M. Integrated behavioural and stable isotope data reveal altered diet linked to low breeding success in urban-dwelling blue tits (Cyanistes caeruleus). Sci. Rep. 7(1), 5014. https://doi.org/10.1038/s41598-017-04575-y (2017).Article 
ADS 
CAS 

Google Scholar 
Schulte-Hostedde, A. I., Mazal, Z., Jardine, C. M. & Gagnon, J. Enhanced access to anthropogenic food waste is related to hyperglycemia in raccoons (Procyon lotor). Conserv. Physiol. 6(1), coy026. https://doi.org/10.1093/conphys/coy026 (2018).Article 
CAS 

Google Scholar 
Fingland, K., Ward, S. J., Bates, A. J. & Bremner-Harrison, S. A systematic review into the suitability of urban refugia for the Eurasian red squirrel Sciurus vulgaris. Mamm. Rev. 52(1), 26–38. https://doi.org/10.1111/mam.12264 (2021).Article 

Google Scholar 
Jokimäki, J., Selonen, V., Lehikoinen, A. & Kaisanlahti-Jokimäki, M.-L. The role of urban habitats in the abundance of red squirrels (Sciurus vulgaris, L.) in Finland. Urban For. Urban Green. 27, 100–108. https://doi.org/10.1016/j.ufug.2017.06.021 (2017).Article 

Google Scholar 
Dausmann, K. H., Wein, J., Turner, J. M. & Glos, J. Absence of heterothermy in the European red squirrel (Sciurus vulgaris). Mammal. Biol. 78(5), 332–335. https://doi.org/10.1016/j.mambio.2013.01.004 (2013).Article 

Google Scholar 
Turner, J. M., Reher, S., Warnecke, L. & Dausmann, K. H. Eurasian red squirrels show little seasonal variation in metabolism in food-enriched habitat. Physiol. Biochem. Zool. 90(6), 655–662. https://doi.org/10.1086/694847 (2017).Article 

Google Scholar 
McNab, B. K. On the comparative ecological and evolutionary significance of total and mass-specific rates of metabolism. Physiol. Biochem. Zool. 72(5), 642–644 (1999).Article 
CAS 

Google Scholar 
Menzies, A. K. et al. Body temperature, heart rate, and activity patterns of two boreal homeotherms in winter: Homeostasis, allostasis, and ecological coexistence. Funct. Ecol. 34(11), 2292–2301. https://doi.org/10.1111/1365-2435.13640 (2020).Article 

Google Scholar 
Wauters, L. & Dhondt, A. Activity budget and foraging behaviour of the red squirrel (Sciurus vulgaris Linnaeus, 1758) in a coniferous habitat. Z. Säugetierkd. 52(6), 341–353 (1987).
Google Scholar 
Wauters, L., Swinnen, C. & Dhondt, A. A. Activity budget and foraging behaviour of red squirrels (Sciurus vulgaris) in coniferous and deciduous habitats. J. Zool. 227(1), 71–86. https://doi.org/10.1111/j.1469-7998.1992.tb04345.x (1992).Article 

Google Scholar 
Reher, S., Dausmann, K. H., Warnecke, L. & Turner, J. M. Food availability affects habitat use of Eurasian red squirrels (Sciurus vulgaris) in a semi-urban environment. J. Mammal. 97(6), 1543–1554. https://doi.org/10.1093/jmammal/gyw105 (2016).Article 

Google Scholar 
Moller, H. Foods and foraging behavior of red (Sciurus vulgaris) and grey (Sciurus carolinensis) squirrels. Mammal. Rev. 13(2–4), 81–98. https://doi.org/10.1111/j.1365-2907.1983.tb00270.x (1983).Article 

Google Scholar 
Krauze-Gryz, D. & Gryz, J. in A review of the diet of the red squirrel (Sciurus vulgaris) in different types of habitats in Red squirrels: Ecology, conservation & management in Europe (eds. Shuttleworth, C. M., Lurz, P. W. W., Hayward, M. W.) 39–50 (European Squirrel Initiative, London, 2015)Shuttleworth, C. M. in The effect of supplemental feeding on the red squirrel (Sciurus vulgaris), Doctoral dissertation (University of London, London, 1996).Birnie-Gauvin, K., Peiman, K. S., Raubenheimer, D. & Cooke, S. J. Nutritional physiology and ecology of wildlife in a changing world. Conserv. Physiol. https://doi.org/10.1093/conphys/cox030 (2017).Article 

Google Scholar 
Wist, B., Stolter, C. & Dausmann, K. H. Sugar addicted in the city: Impact of urbanisation on food choice and diet composition of the Eurasian red squirrel (Sciurus vulgaris). J. Urban Ecol. 8(1), juac012. https://doi.org/10.1093/jue/juac012 (2022).Article 

Google Scholar 
Burton, T., Killen, S. S., Armstrong, J. D. & Metcalfe, N. B. What causes intraspecific variation in resting metabolic rate and what are its ecological consequences?. Proc. Biol. Sci. 278(1724), 3465–3473. https://doi.org/10.1098/rspb.2011.1778 (2011).Article 
CAS 

Google Scholar 
Clarke, A. Costs and consequences of evolutionary temperature adaptation. Trends Ecol. Evol. 18(11), 573–581. https://doi.org/10.1016/j.tree.2003.08.007 (2003).Article 

Google Scholar 
Lovegrove, B. G. The influence of climate on the basal metabolic rate of small mammals: A slow-fast metabolic continuum. J. Comp. Physiol. B 173(2), 87–112. https://doi.org/10.1007/s00360-002-0309-5 (2003).Article 
CAS 

Google Scholar 
McNab, B. K. The energetics of endotherms. Ohio J. Sci. 74(6), 370–380 (1974).
Google Scholar 
Tattersall, G. J. et al. Coping with thermal challenges: Physiological adaptations to environmental temperatures. Compr. Physiol. 2(3), 2151–2202 (2012).Article 

Google Scholar 
Broggi, J. et al. Sources of variation in winter basal metabolic rate in the great tit. Funct. Ecol. 21(3), 528–533. https://doi.org/10.1111/j.1365-2435.2007.01255.x (2007).Article 

Google Scholar 
Schlünzen, K. H., Hoffmann, P., Rosenhagen, G. & Riecke, W. Long-term changes and regional differences in temperature and precipitation in the metropolitan area of Hamburg. Int. J. Climatol. 30(8), 1121–1136. https://doi.org/10.1002/joc.1968 (2010).Article 

Google Scholar 
Reher, S. & Dausmann, K. H. Tropical bats counter heat by combining torpor with adaptive hyperthermia. Proc. R. Soc. B Biol. Sci. 288(1942), 20202059. https://doi.org/10.1098/rspb.2020.2059 (2021).Article 

Google Scholar 
Rezende, E. L. & Bacigalupe, L. D. Thermoregulation in endotherms: Physiological principles and ecological consequences. J. Comp. Physiol. B 185(7), 709–727. https://doi.org/10.1007/s00360-015-0909-5 (2015).Article 
CAS 

Google Scholar 
Scholander, P. F., Hock, R., Walters, V., Johnson, F. & Irving, L. Heat regulation in some arctic and tropical mammals and birds. Biol. Bull. 99(2), 237–258. https://doi.org/10.2307/1538741 (1950).Article 
CAS 

Google Scholar 
Terblanche, J. S., Clusella-Trullas, S., Deere, J. A., Van Vuuren, B. J. & Chown, S. L. Directional evolution of the slope of the metabolic rate-temperature relationship is correlated with climate. Physiol. Biochem. Zool. 82(5), 495–503. https://doi.org/10.1086/605361 (2009).Article 

Google Scholar 
Gallo, K. P., Easterling, D. R. & Peterson, T. C. The influence of land use/land cover on climatological values of the diurnal temperature range. J. Clim. 9(11), 2941–2944. https://doi.org/10.1175/1520-0442(1996)009%3c2941:TIOLUC%3e2.0.CO;2 (1996).Article 
ADS 

Google Scholar 
Wang, K. et al. Urbanization effect on the diurnal temperature range: Different roles under solar dimming and brightening. J. Clim. 25(3), 1022–1027. https://doi.org/10.1175/jcli-d-10-05030.1 (2012).Article 
ADS 

Google Scholar 
Fristoe, T. S. et al. Metabolic heat production and thermal conductance are mass-independent adaptations to thermal environment in birds and mammals. Proc. Natl. Acad. Sci. USA 112(52), 15934–15939. https://doi.org/10.1073/pnas.1521662112 (2015).Article 
ADS 
CAS 

Google Scholar 
Sándor, K. et al. Urban nestlings have reduced number of feathers in Great Tits (Parus major). Ibis 163(4), 1369–1378. https://doi.org/10.1111/ibi.12948 (2021).Article 

Google Scholar 
Beliniak, A., Krauze-Gryz, D., Jasińska, K., Jankowska, K. & Gryz, J. Contrast in daily activity patterns of red squirrels inhabiting urban park and urban forest. Hystrix https://doi.org/10.4404/hystrix-00476-2021 (2021).Article 

Google Scholar 
Thomas, L. S., Teich, E., Dausmann, K., Reher, S. & Turner, J. M. Degree of urbanisation affects Eurasian red squirrel activity patterns. Hystrix 29(2), 175–180. https://doi.org/10.4404/hystrix-00065-2018 (2018).Article 

Google Scholar 
Krauze-Gryz, D., Gryz, J. & Brach, M. Spatial organization, behaviour and feeding habits of red squirrels: Differences between an urban park and an urban forest. J. Zool. 315(1), 69–78. https://doi.org/10.1111/jzo.12905 (2021).Article 

Google Scholar 
Jarman, T. E., Gartrell, B. D. & Battley, P. F. Differences in body composition between urban and rural mallards Anas platyrhynchos. J. Urban Ecol. 6(1), juaa011. https://doi.org/10.1093/jue/juaa011 (2020).Article 

Google Scholar 
Cruz-Neto, A. P. & Bozinovic, F. The relationship between diet quality and basal metabolic rate in endotherms: Insights from intraspecific analysis. Physiol. Biochem. Zool. 77(6), 877–889 (2004).Article 

Google Scholar 
Geluso, K. & Hayes, J. P. Effects of dietary quality on basal metabolic rate and internal morphology of European starlings (Sturnus vulgaris). Physiol. Biochem. Zool. 72(2), 189–197 (1999).Article 
CAS 

Google Scholar 
Seebacher, F. Is endothermy an evolutionary by-product?. Trends Ecol. Evol. 35(6), 503–511. https://doi.org/10.1016/j.tree.2020.02.006 (2020).Article 

Google Scholar 
Perissinotti, P. P., Antenucci, C. D., Zenuto, R. & Luna, F. Effect of diet quality and soil hardness on metabolic rate in the subterranean rodent Ctenomys talarum. Comp. Biochem. Physiol. Mol. Integr. Physiol. 154(3), 298–307. https://doi.org/10.1016/j.cbpa.2009.05.013 (2009).Article 
CAS 

Google Scholar 
Thorp, C. R., Ram, P. K. & Florant, G. L. Diet alters metabolic rate in the yellow-bellied marmot (Marmota flaviventris) during hibernation. Physiol. Zool. 67(5), 1213–1229. https://doi.org/10.1086/physzool.67.5.30163890 (1994).Article 

Google Scholar 
Silva, S. I., Jaksic, F. M. & Bozinovic, F. Interplay between metabolic rate and diet quality in the South American fox Pseudalopex culpaeus. Comp. Biochem. Physiol. Mol Integr. Physiol. 137(1), 33–38. https://doi.org/10.1016/j.cbpb.2003.09.022 (2004).Article 
CAS 

Google Scholar 
Rewkiewicz-Dziarska, A., Wielopolska, A. & Gill, J. Hematological indices of Apodemus agrarius (Pallas, 1771) from different urban environments. Bull. Acad. Polon. Sci. Ser. Sci. Biol. 25(4), 261–268 (1977).CAS 

Google Scholar 
Ohrnberger, S. A., Hambly, C., Speakman, J. R. & Valencak, T. G. Limits to sustained energy intake XXXII: Hot again: Dorsal shaving increases energy intake and milk output in golden hamsters (Mesocricetus auratus). J Exp. Biol. https://doi.org/10.1242/jeb.230383 (2020).Article 

Google Scholar 
Speakman, J. R. & Król, E. The heat dissipation limit theory and evolution of life histories in endotherms—Time to dispose of the disposable soma theory?. Integr. Comp. Biol. 50(5), 793–807. https://doi.org/10.1093/icb/icq049 (2010).Article 

Google Scholar 
Diamond, S. E., Chick, L. D., Perez, A., Strickler, S. A. & Martin, R. A. Evolution of thermal tolerance and its fitness consequences: Parallel and non-parallel responses to urban heat islands across three cities. Proc. R. Soc. B Biol. Sci. 285(1882), 20180036. https://doi.org/10.1098/rspb.2018.0036 (2018).Article 

Google Scholar 
Isaksson, C. & Hahs, A. Urbanization, oxidative stress and inflammation: A question of evolving, acclimatizing or coping with urban environmental stress. Funct. Ecol. 29(7), 913–923. https://doi.org/10.1111/1365-2435.12477 (2015).Article 

Google Scholar 
Sokolova, I. M. & Lannig, G. Interactive effects of metal pollution and temperature on metabolism in aquatic ectotherms: Implications of global climate change. Clim. Res. 37(2–3), 181–201 (2008).Article 

Google Scholar 
Carey, H. V., Andrews, M. T. & Martin, S. L. Mammalian hibernation: Cellular and molecular responses to depressed metabolism and low temperature. Physiol. Rev. 83(4), 1153–1181 (2003).Article 
CAS 

Google Scholar 
Pereira, M. E., Aines, J. & Scheckter, J. L. Tactics of heterothermy in eastern gray squirrels (Sciurus carolinensis). J. Mammal. 83(2), 467–477 (2002).Article 

Google Scholar 
Breuner, C. W., Wingfield, J. C. & Romero, L. M. Diel rhythms of basal and stress-induced corticosterone in a wild, seasonal vertebrate. Gambel’s white-crowned sparrow. J Exp. Zool. 284(3), 334–342. https://doi.org/10.1002/(SICI)1097-010X(19990801)284:3%3c334::AID-JEZ11%3e3.0.CO;2-# (1999).Article 
CAS 

Google Scholar 
Careau, V., Thomas, D., Humphries, M. M. & Réale, D. Energy metabolism and animal personality. Oikos 117(5), 641–653. https://doi.org/10.1111/j.0030-1299.2008.16513.x (2008).Article 

Google Scholar 
Fletcher, Q. E. et al. Seasonal stage differences overwhelm environmental and individual factors as determinants of energy expenditure in free-ranging red squirrels. Funct. Ecol. 26(3), 677–687. https://doi.org/10.1111/j.1365-2435.2012.01975.x (2012).Article 

Google Scholar 
Barthel, L. & Berger, A. Unexpected gene-flow in urban environments: The example of the European Hedgehog. Animals 10(12), 2315. https://doi.org/10.3390/ani10122315 (2020).Article 

Google Scholar 
Fusco, N. A., Carlen, E. J. & Munshi-South, J. Urban landscape genetics: are biologists keeping up with the pace of urbanization?. Current Landsc. Ecol. Rep. 6(2), 35–45. https://doi.org/10.1007/s40823-021-00062-3 (2021).Article 

Google Scholar 
Ziege, M. et al. Population genetics of the European rabbit along a rural-to-urban gradient. Sci. Rep. 10(1), 2448. https://doi.org/10.1038/s41598-020-57962-3 (2020).Article 
ADS 
CAS 

Google Scholar 
Morash, A. J., Neufeld, C., MacCormack, T. J. & Currie, S. The importance of incorporating natural thermal variation when evaluating physiological performance in wild species. J. Exp. Biol. 221(14), jeb164673. https://doi.org/10.1242/jeb.164673 (2018).Article 

Google Scholar 
Pörtner, H.-O., et al. Climate change 2022: Impacts, adaptation and vulnerability. IPCC Sixth Assessment Report (2022).Anderies, J. M., Katti, M. & Shochat, E. Living in the city: Resource availability, predation, and bird population dynamics in urban areas. J. Theor. Biol. 247(1), 36–49. https://doi.org/10.1016/j.jtbi.2007.01.030 (2007).Article 
ADS 
MATH 

Google Scholar 
Shochat, E. Credit or debit? Resource input changes population dynamics of city-slicker birds. Oikos 106(3), 622–626. https://doi.org/10.1111/j.0030-1299.2004.13159.x (2004).Article 

Google Scholar 
Koprowski, J. L. Handling tree squirrels with a safe and efficient restraint. Wildl. Soc. B 30(1), 101–103. https://doi.org/10.2307/3784642 (2002).Article 

Google Scholar 
Magris, L. & Gurnell, J. Population ecology of the red squirrel (Sciurus vulgaris) in a fragmented woodland ecosystem on the Island of Jersey Channel Islands. J. Zool. 256(1), 99–112. https://doi.org/10.1017/s0952836902000134 (2002).Article 

Google Scholar 
Bethge, J., Wist, B., Stalenberg, E. & Dausmann, K. Seasonal adaptations in energy budgeting in the primate Lepilemur leucopus. J Comp. Physiol. B 187(5–6), 827–834. https://doi.org/10.1007/s00360-017-1082-9 (2017).Article 

Google Scholar 
Dausmann, K. H., Glos, J. & Heldmaier, G. Energetics of tropical hibernation. J Comp. Physiol. B 179(3), 345–357. https://doi.org/10.1007/s00360-008-0318-0 (2009).Article 
CAS 

Google Scholar 
Kobbe, S., Nowack, J. & Dausmann, K. H. Torpor is not the only option: Seasonal variations of the thermoneutral zone in a small primate. J. Comp. Physiol. B 184(6), 789–797. https://doi.org/10.1007/s00360-014-0834-z (2014).Article 

Google Scholar 
Lighton, J. R. Measuring Metabolic Rates: A Manual for Scientists (Oxford University Press, 2018).Book 

Google Scholar 
Bethge, J., Razafimampiandra, J. C., Wulff, A. & Dausmann, K. H. Sportive lemurs elevate their metabolic rate during challenging seasons and do not enter regular heterothermy. Conserv. Physiol. 9(1), coab075. https://doi.org/10.1093/conphys/coab075 (2021).Article 

Google Scholar 
Reher, S., Ehlers, J., Rabarison, H. & Dausmann, K. H. Short and hyperthermic torpor responses in the Malagasy bat Macronycteris commersoni reveal a broader hypometabolic scope in heterotherms. J. Comp. Physiol. B 188(6), 1015–1027. https://doi.org/10.1007/s00360-018-1171-4 (2018).Article 
CAS 

Google Scholar 
Grolemund, G. & Wickham, H. Dates and times made easy with lubridate. J Stat. Softw. 40(3), 1–25 (2011).Article 

Google Scholar 
Wickham, H., François, R., Henry, L. & Müller, K. RStudio. dplyr: A Grammar of Data Manipulation (1.0. 7) (2021).Zeileis, A. & Grothendieck, G. zoo: S3 infrastructure for regular and irregular time series. J. Stat. Softw. 14(6), 1–27. https://doi.org/10.18637/jss.v014.i06 (2005).Article 

Google Scholar 
Sarkar, D. Lattice: Multivariate Data Visualization with R (Springer Science & Business Media, New York, 2008).Book 
MATH 

Google Scholar 
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. https://doi.org/10.18637/jss.v067.i01 (2015).Article 

Google Scholar 
Kuznetsova, A., Brockhoff, P. B. & Christensen, R. H. lmerTest package: Tests in linear mixed effects models. J. Stat. Softw. 82(13), 1–26 (2017).Article 

Google Scholar 
Wickham, H. ggplot2: Elegant graphics for data analysis (Springer, 2016).Book 
MATH 

Google Scholar 
Fox, J. Effect displays in R for generalised linear models. J. Stat. Softw. 8(15), 1–27 (2003).Article 

Google Scholar 
Garamszegi, L. Z. et al. Changing philosophies and tools for statistical inferences in behavioral ecology. Behav. Ecol. 20(6), 1363–1375. https://doi.org/10.1093/beheco/arp137 (2009).Article 

Google Scholar 
Symonds, M. R. E. & Moussalli, A. A brief guide to model selection, multimodel inference and model averaging in behavioural ecology using Akaike’s information criterion. Behav. Ecol. Sociobiol. 65(1), 13–21. https://doi.org/10.1007/s00265-010-1037-6 (2010).Article 

Google Scholar 
Whittingham, M. J., Stephens, P. A., Bradbury, R. B. & Freckleton, R. P. Why do we still use stepwise modelling in ecology and behaviour?. J. Anim. Ecol. 75(5), 1182–1189. https://doi.org/10.1111/j.1365-2656.2006.01141.x (2006).Article 

Google Scholar 
Barton, K. & Barton, M. K. MuMIn: Multi-Model Inference. R package version 1.43.17; https://CRAN.R-project.org/package=MuMIn (2020).Zuur, A., Ieno, E. N., Walker, N., Saveliev, A. A. & Smith, G. M. Mixed effects models and extensions in ecology with R ( Springer Science & Business Media 2009).Burnham, K. P. & Anderson, D. R. Multimodel inference: Understanding AIC and BIC in model selection. Soc. Method. Res. 33(2), 261–304 (2004).Article 

Google Scholar 
Johnson, J. B. & Omland, K. S. Model selection in ecology and evolution. Trends Ecol. Evol. 19(2), 101–108. https://doi.org/10.1016/j.tree.2003.10.013 (2004).Article 

Google Scholar 
Lorah, J. Effect size measures for multilevel models: Definition, interpretation, and TIMSS example. Large-scale Assess. Educ. 6(1), 8. https://doi.org/10.1186/s40536-018-0061-2 (2018).Article 

Google Scholar 
Selya, A. S., Rose, J. S., Dierker, L. C., Hedeker, D. & Mermelstein, R. J. A practical guide to calculating cohen’s f2, a measure of local effect size, from PROC MIXED. Front. Psychol. 3, 111–111. https://doi.org/10.3389/fpsyg.2012.00111 (2012).Article 

Google Scholar 
Lüdecke, D. sjPlot: Data visualization for statistics in social science. R package version 2.8.5 2020; https://CRAN.R-project.org/package=sjPlot (2020).Nakagawa, S., Johnson, P. C. & Schielzeth, H. The coefficient of determination R2 and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded. J. R. Soc. Interface 14(134), 20170213 (2017).Article 

Google Scholar 
Nakagawa, S. & Schielzeth, H. Repeatability for Gaussian and non-Gaussian data: A practical guide for biologists. Biol. Rev. Camb. Philos. Soc. 85(4), 935–956. https://doi.org/10.1111/j.1469-185X.2010.00141.x (2010).Article 

Google Scholar 
Stoffel, M. A., Nakagawa, S. & Schielzeth, H. rptR: Repeatability estimation and variance decomposition by generalized linear mixed-effects models. Methods Ecol. Evol. 8(11), 1639–1644. https://doi.org/10.1111/2041-210X.12797 (2017).Article 

Google Scholar  More

Real-world dataThe dataset consisted of 2986 observations from 902 freshwater shallow lakes in Denmark and North America (data extracted from the LAGOSNE database on 22 February 2022 via R LAGOSNE package version 2.0.2)56 (Supplementary Fig. 9). The Danish lakes were sampled for one or several years from 1984 to 2020 (data extracted in October 2021 from https://odaforalle.au.dk/main.aspx) (Supplementary Fig. 10). Prerequisites for inclusion in the analysis were that lakes had been sampled for physical and chemical variables at least four times or at least three times over the growing season (May to September) for the Danish or North American lakes, respectively, had a mean depth of less than 3 m and were freshwater. Water chemistry samples were analysed using standard methods and data for total phosphorus (TP), total nitrogen (TN) and chlorophyll-a are included here57. The mean and range of TP, TN and chlorophyll-a for the combined sites is given in Table 1, along with the values for each region separately.To gain a longer-term perspective on the relationship between nutrients and chlorophyll-a, we calculated the across-year averages of the summer means of TP, TN and chlorophyll-a, sequentially increasing numbers of years included in the mean up to a total of a five-year mean, at which point there were only 99 lakes left in the dataset. In calculating the multi-year means we allowed a maximum gap of 2 years between observations (i.e. two observations could cover 3 years) to avoid including time series with too many missing years in between. Hence, only lakes with sufficient numbers of sequential data were included, resulting in a large drop in lake numbers as the length of the multi-year mean increased (Table 2).Numerical methodsDiagnostic tests or proxies of alternative equilibriaWe modelled the response of chlorophyll-a to TP and TN using generalised-linear models58 with Gamma distribution and an identity link on untransformed data for single-year and multiple-year means up to 5-year means. We used the Gamma distribution, as chlorophyll fit this significantly better than a normal or log-normal distribution. We used psuedo R2 of the model along with the patterns of residuals, and finally, we plotted the kernel density of the chlorophyll-a values as diagnostics of the presence, absence or prevalence of alternative equilibria in the simulated and real work data.To test how appropriate these diagnostics or proxies of alternative stable states in terms of how well they identify the existence of alternative stable states in randomly sampled multi-year data, we

1.

Simulated two scenarios for the main manuscript, with and without alternative stable states in the data, which were as close to the real-life data as possible. The results of these scenarios appear in the main text (please see details below in the “Data Simulation” section).

2.

We provide multiple scenarios with different degrees, or prevalence, of alternative stable states in the data, see simulations of alternative stable state scenarios. The results of these scenarios appear in Supplementary note 2.

Hierarchical bootstrap approachThere are a large number of permutations of data, both real-world and simulated, that can provide a mean of the two to five sequential years from each lake in the time series data. It was vital to have a method that selects the data for analysis that provides a valid and comparable representation of both real work and simulated data and the models’ errors. In order to provide this we used a non-parametric hierarchical bootstrap procedure38. The flowchart shows the data preparation and data analysis steps of the hierarchical bootstrap procedure (Fig. 4). In the first step (step 1 in Fig. 4), all possible longer-term means are calculated for each lake. To keep as much data as possible, we decided to allow for up to 2 years of gap in the data between years. Taking the 5-year mean data as an example, if data from a lake existed for the years 1991 and 1994−1997, a 5-year mean would be calculated for the years 1991, 1994, 1995, 1996 and 1997. However, if the time series would contain a larger gap, e.g. data would only exist for the years 1991 and 1995–1998, no 5-year mean could be calculated. After the selection procedure, all the 2-year, 3-year and 5-year means are transferred into a new table (step 2 in Fig. 4).Fig. 4Data preparation and analysis steps of the hierarchical bootstrap procedure.Full size imageThe procedure is the same for each temporal scale from 2-year means to 5-year means. For the example of 5 mean years, lakes are randomly sampled from the full 5-year mean dataset in step 2 (Fig. 4) with replacement up to the number of lakes as in the original dataset, for the 5-year means 99 (step 3a). Here, the same lake can appear multiple times or not at all. This step is common for every bootstrap procedure59. However, since we have nested data (5-year means within lakes), we need a second step, in which for every resampled lake in step 3a, one 5-year mean is chosen (step 3b in Fig. 4). Then the three GLM models are produced from the randomly selected data in step 3c (Fig. 4). These steps are then repeated 1000 times to get a good representation of the uncertainties of the model. To ensure a fair comparison between single-year data and their equivalent multi-year mean data, we repeated the bootstrap procedure with single years only using only the lakes for which we also calculated multi-year means. To take the five-year mean as an example, there were 99 lakes where we could calculate at least one 5-year mean observation. First, we ran the bootstrap procedure to calculate 5-year mean values of TP, TN and chlorophyll-a (1000 times) and then took single years’ values of TP, TN and chlorophyll-a (1000 times) from exactly the same 99 lakes. With this approach, exactly the same datasets with the same lakes and observations within lakes are used for the calculation of the multi-year means and their single-year counterparts, making for a robust analysis. The GLM models did not always converge. If either the TP, TN or TP*TN model with interaction did not converge, the iteration was not used in further analysis. The number of converging models equal for each iteration of random samples is given in the results.The described hierarchical approach is the best way to reflect the structure of the original data. A simple, non-hierarchical bootstrap would favour lakes with more five-year means over lakes with fewer five-year means, simply because these make up a larger part of the data. Furthermore, sampling without replacement at the lake level would result in five-year means from lakes with few data dominating the produced random dataset, as every lake would be sampled every time, which then would result in high model leverage of 5-year means from lakes with less data. In contrast, the hierarchical procedure ensures that every lake has the same chance to end up in the randomly sampled bootstrap, in the second step, it ensures that of each sampled lake, every 5-year mean has the same chance to end up in the random dataset. These notions are in agreement with the findings of an assessment on how to properly resample hierarchical data by non-parametric bootstrap38.Data simulationGeneral approach of simulation assumptions and proceduresWe generated random scatter for the generalised-linear model based on Gamma distributions for two different “populations” of lakes with two different intercepts and slopes. At first, we calculated the linear equations for the two populations:For each population i and j, 99 samples (equalling the number of lakes in real-life data with 5-year means, nlake) were generated with a specific number of data points depending on the scenario (nyear) each, hence nlake = 1−99 for each population of lakes, e.g. with 20 years (nyear = 20) each.We found the real nutrient data to be normally distributed, with total nitrogen (TN) having a range between 0.33 and 4.93 mg/L and a constant coefficient of variation (CV, with a mean CV of 0.35) across this range (the same is true for total phosphorus (TP) at a shorter range). Hence, for each nlake, the x for the nyear = 20 were generated based on the mean range (mean per lake of the real-life data) and CV (0.35) from the real-life TN concentration data, hence with a range of 0.33 to 4.33 mg/L. Therefore the values and random variability of x in the simulations are close to the true values of the TN concentrations. The x is then fed into the linear equations above.To the resulting yi and yj we added random noise based on the Gamma distribution (using the rgamma function in R). We used a Gamma distribution because the Chlorophyll-a concentration also follows a Gamma distribution. The variability of a Gamma distribution is expressed by the shape variable. The variability of chlorophyll-a, its shape value, equals 2.63. This shape value was used in the Gamma distribution of yi and yj. The final calculated yi and yj had therefore a random rate calculated as shape/yi or shape/yj. Hence, their variability in the y dimension was close to the true chlorophyll-a variability.The data from both lake populations were then pooled and randomly sampled using the same hierarchical bootstrap procedure with 500 iterations for the scenarios in the supplementary materials and with 1000 iterations for main text simulation scenarios, which is identical to what was done for the real-world data.Simulation scenarios based on characteristics of real-world dataThe real-world 5-year mean data consisted of 99 lakes with 5–20 years of data for each lake. For the simulation scenario in the main text, we therefore randomly sampled between 5 and 20 data points for each of the 99 simulated lakes based on the x distribution described above. Intercepts and slopes of the simulation, resembled the range of the true data (see scatter plots in Fig. 2 of the main manuscript).In the alternative stable state scenario, we chose two slopes and intercepts for different populations of lakes:

Population i: ai = 0, bi = 40

Population j: aj = 50, bj = 120

We based the slopes and intercepts of the ASS scenarios on the diagnostic combination defined by Scheffer and Carpenter7 which propose an abrupt shift in (a) the time series, (b) the multimodal distribution of states and (c) the dual relationship to a controlling factor. Here, the idea is that an ecosystem will jump from one state to the next at the same (nutrient) conditions (different intercept and/or slope, condition a within ref. 7), where any change in the nutrient will have different effects on algae or macrophytes (best represented by different slope, condition c), resulting in a multimodal distribution of the response (condition b). Hence, simulations are in line with what is predicted for ASS, but we took great lengths to also show other possibilities with the simulations in the Supplementary information to ensure we did not overlook any occasional occurrence of alternative equilibria.Here, the appearance of alternative stable states in the data could happen at any point in the time series of a single lake, or the entire time series could include only one of the two alternative stable states. To accommodate these alternative stable state constellations (for each of which we made a separate simulation scenario, (see Supplementary Note 2, “Simulations of alternative stable state constellations”), we forced the alternative stable state scenario to be constructed of 1/3 of data with one state, 1/3 of data with the second state and 1/3 of data where both alternative states could occur. In the latter case, the alternative stable state appeared after the first 20% but before the last 20% of the time series. Since the variability and range of x (nutrient) and y (chlorophyll -a response) is simulated as close as possible to the real-world data in all scenarios, the measures taken here (variable time series and combination of different alternative stable state scenario constellations) produce a simulation as close to the real-world data as possible. Specifically, we found the real-world nutrient data to be normally distributed, with total nitrogen (TN) having a range between 0.33 and 4.93 mg/L and a constant coefficient of variation (CV, with a mean CV of 0.35) across this range (the same is true for total phosphorus (TP) at a shorter range). Hence, for each simulated lake, the x were generated based on this mean range and CV. Furthermore, the resulting yi and yj were randomised by using a Gamma distribution (using the rgamma function in R). We used a Gamma distribution because the chlorophyll-a concentration also follows a Gamma distribution. The variability of a Gamma distribution is expressed by the shape variable. The variability of chlorophyll-a, its shape value, equals 2.63. This shape value was used in the Gamma distribution of yi and y. The final calculated yi and yj had, therefore a random rate calculated as shape/yi or shape/y. Hence, their variability in the y dimension was close to the true chlorophyll-a variability.For the scenario without alternative stable states, both populations of data had the same intercept and slope:

Population i: ai = 0, bi = 40

Population j: aj = 0, bj = 40.

Please see Supplementary Note 2 for further simulations of different potential constellations of alternative states. There we show that our approach finds alternative stable states in response to nutrient concentration, even if they appear in time series from different lakes.Assessment of diagnostic tests or proxies of alternative equilibriaWe modelled the response of chlorophyll-a to TP and TN using generalised-linear models3 with Gamma distribution and an identity link on untransformed data for single-year and multiple-year means up to 5-year means. We used the Gamma distribution, as chlorophyll fit this significantly better than a normal or log-normal distribution. We used R2 of the model along with the patterns of residuals, and finally, we plotted the kernel density of the chlorophyll-a values as diagnostics of the presence, absence or prevalence of alternative equilibria in the simulated and real work data.The comparison of how the diagnostics/proxies of alternative stable states respond to the variation in the prevalence of alternative equilibria in the simulated datasets provides a robust assessment of their ability to identify both the presence and absence of alternative equilibria. It is the response of these diagnostic tests over time, with the increase in the temporal perspective as more years are added to the mean values of TP, TN and chlorophyll-a, that are key to the identification of the presence and or absence of alternative equilibria in a given dataset. The simulations show that a dataset which contains alternative equilibria will show (1) no improvement in R2 as the temporal perspective of the data increases (more years in the multi-year mean); (2) an increased bimodality in the residuals of the models of nutrients predicting chlorophyll-a will increase as more years are added to the multi-year mean and (3) the kernel density function of chlorophyll-a will display increasingly bimodality as more years are added to the mean. In the absence of alternative equilibria, the patterns differ with an R2, and increase in unimodality of residuals and a consistent unimodal pattern in the kernel density function. Thus, the diagnostic tests provide a robust test of both the presence and absence of alternative equilibria in a given dataset.Alternative stable state assessment for real data with limited data rangeIt could be the case that alternative stable states do not appear in the full dataset but only in a limited TN and TP concentration range. We filtered and re-analyzed the data, only keeping data points within the following two ranges: – TN concentration = 0.5−2 mg/L–TP concentration = 0.05−0.4 mg/L. In the filtered data, 1329 out of the original 2876 single-year data points, 289 out of 1028 3-year mean data points and 212 out of the 864 five-mean year data points remained. The filtered data consisted of data points from 550, 48 and 27 lakes for the single-year data, 3-year means and 5-year means, respectively. The smaller range resulted in lower R² of the models, yet the pattern that multi-year means result in higher R² compared to single-year data was largely consistent, apart from the 5-year mean TN models for which both, the single-year and mean data resulted in very low R² (Supplementary Fig. 6). Furthermore, due to the lower number of samples, the errors of all proxies are higher, making conclusions more difficult than for the full data. Still, we do not see any clear indication of alternative stable states in the scatter plots (two groups of dots are not appearing (Supplementary Fig. 5), the Kernel density plots (or model residuals (Supplementary Fig. 6)). i.e. no signs of bimodality in residuals or Kernel density plots. Please see details on this analysis in the supplementary material.Details and the R code for the steps for the random multi-year sampling can be found in the supplementary materials.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More

Game theory is a framework for analysing the outcome of the strategic interaction between decision makers1. The fundamental concept is that of a Nash equilibrium where no player can improve her payoff by a unilateral strategy change. Typically, the Nash equilibrium is considered to be the optimal outcome of a game, however in social dilemmas the individual optimal outcome is at odds with the collective optimal outcome2. This means that one player can improve her payoff at the expense of the other by unilaterally deviating, but if both deviate, they end up with lower payoffs. In this type of games, the mutually beneficial, but non-Nash equilibrium strategy is called cooperation. However, in this context cooperation should not be interpreted as an interest in the welfare of others, as players only aim to secure a high payoff for themselves.In this framework, payoff maximisation is considered to be rational, but when such rational players then seize every opportunity to gain at the opponent’s expense, they may counterintuitively both end up with low payoffs. A game that clearly exhibits this contradiction is the Traveler’s Dilemma. Since its formulation in 1994 by the economist Kaushik Basu3, the game has become one of the most studied in the economics literature. Additionally, it has been discussed in theoretical biology in the context of evolutionary game theory.In general, the dilemma relies on the individuals’ incentive to undercut the opponent. To be more specific, players are motivated to claim a lower value than their opponent to reach a higher payoff at the opponent’s expense. Such incentive leads players to a systematic mutual undercutting until the lowest possible payoff is reached, which is the unique Nash equilibrium. It seems paradoxical that players defined as rational in a game theoretical sense end up with such a poor outcome. Therefore, the question that naturally arises is how can this poor outcome be prevented and how cooperation can be achieved.To address these questions, it can be helpful to better understand price wars, which consist in the mutual undercutting of prices to gain market share. In addition, it can provide information about human behaviour, because economic experiments have shown that individuals prefer to choose the cooperative high payoff action, instead of the Nash equilibrium4.Our analysis focuses on showing that the Traveler’s Dilemma can be decomposed into a local and a global game. If the payoff optimisation is constrained to the local game, then players will inevitably end up in the Nash equilibrium. However, if players escape the local maximisation and optimise their payoff for the global game, they can reach the cooperative high payoff equilibrium.Here, we show that the cooperative equilibrium can be reached in a game like the Traveler’s Dilemma due to diversity, which we define as the presence of suboptimal strategies. The appearance of strategies far from those of the residents allows for the local maximisation process to be escaped, such that an optimisation at a global level takes place. Overall this can lead to cooperation because by considering “suboptimal strategies” that play against each other it is possible to reach higher payoffs, both collectively and individually.GameThe Traveler’s Dilemma is a two-player game. Player i has to choose a claim, (n_i), from the action space, consisting of all integers on the interval [L, U], where (0 le L < U). The payoffs are determined as follows: If both players, i and j, choose the same value ((n_i = n_j)), both get paid that value. There is a reward parameter (R >1), such that if (n_i < n_j), then i receives (n_i + R) and j gets (n_i- R) Thus, the payoff of player i playing against player j is$$begin{aligned} pi _{ij} = {left{ begin{array}{ll} n_i& text { if } n_i = n_j\ n_i + R& text { if } n_i < n_j\ n_j - R& text { if } n_i > n_j end{array}right. } end{aligned}$$
(1)
Thus, a player is better off by choosing a slightly lower value than the opponent: when j plays (n_j), then it is best for i to play (n_j-1). The iteration of this reasoning, which we will call the stairway to hell, leads to the only Nash equilibrium of the game, ({L,L}), where both players choose the lowest possible claim. The classical game theory method to arrive to this equilibrium is called iterative elimination of dominated strategies5.The game can be visualised through its payoff matrix (Fig. 1). For simplicity, we use the values from the original formulation: (L=2), (U=100) and (R=2). The payoff matrix shows that the Traveler’s Dilemma can be decomposed into a local and a global game. Let us begin with the local game. When the action space of the game is reduced to two adjacent actions n and (n+1) (black boxes in Fig. 1), the Traveler’s Dilemma with (R=2) is equivalent to the Prisoner’s Dilemma6. In general, for any value of R, the Traveler’s Dilemma becomes a Prisoner’s Dilemma for any pair of actions n and (n+s), where ( 1 le s le R-1 ). For example, for (R=4) the pair of actions n and (n+1), n and (n+2), n and (n+3) follow the same game structure as the Prisoner’s Dilemma. Therefore, the Traveler’s Dilemma consists of many embedded Prisoners’ Dilemmas. This means that at a local level the game is a Prisoner’s Dilemma.If we now consider actions that are distant from each other in the action space, e.g. 2 and 100, we can observe a coordination game structure (gray boxes in Fig. 1), where ({100,100}) is payoff and risk dominant7,8. In general, any pair of actions n and (n+s), where ( R le s le U-n), construct a coordination game. As a result, the Traveler’s Dilemma becomes a coordination game at a global level, which has different equilibria than the local game.Figure 1Payoff matrix of the Traveler’s Dilemma. Visualisation of the payoff scheme described by Eq. (1). For simplicity, the action space is ( {n_i in {mathbb {N}} mid 2 le n_i le 100}) and the reward parameter is (R=2). The Traveler’s Dilemma can be decomposed into a local and a global game. At a local level the game is a Prisoner’s Dilemma (black boxes). This happens when the action space is reduced to any pair of actions n and (n+s), where ( 1 le s le R-1 ). While at the global level, we can observe a coordination game (gray boxes). This level is defined as any pair of actions n and (n+s), where ( R le s le U-n).Full size imageParadoxSocial dilemmas appear paradoxical in the sense that self-interested competing players, when rationally playing the Nash equilibrium, end up with a payoff that clearly goes against their self-interest. But with the Traveler’s Dilemma, the paradox goes further, as suggested in its original formulation3. Classical game theory proposes ({L,L}) as the Nash equilibrium of the game. However, it seems unlikely and implausible that, with R being moderately low, say (R=2), for individuals to play the Nash equilibrium. This has been confirmed in economic experiments where individuals rather choose values close to the upper bound of the interval. Such experiments have also shown that the chosen value depends on the reward parameter (R), where an increasing value of R shifts players’ decision towards ({L,L})4. Nonetheless, classical game theory states that the Nash equilibrium of the game is independent of R.Consequently, the aim of this paper is to seek and explore simple mechanisms through which the apparent non-rational cooperative behaviour can come about. We also examine the effect of the reward parameter on the game’s outcome. Given that the Traveler’s Dilemma paradox emerges in the classical game theory framework, we analyse the game using evolutionary game theory tools5,9,10. This dynamical approach allows us to explore adaptive behaviour outside of the stationary classical game theory framework. To be more precise, for this approach individuals dynamically adjust their actions according to their payoffs.The key point of course is to understand how the system can converge to high claims. We show that this behaviour is possible because the Traveler’s Dilemma can be decomposed into a local and a global game. If the payoff maximisation is constrained to the local level, then the stairway to hell leads the system to the Nash equilibrium; given that locally the game is a Prisoner’s Dilemma. On the other hand, at a global level the game follows a coordination game structure, where the high claim actions are payoff dominant. Thus, for the system to reach a high claim equilibrium the maximisation process needs to jump from the local to the global level.Our analysis led us to identify the mechanism of diversity as responsible for enabling this jump and preventing players from going down the stairway to hell. This mechanism works on the idea that to reach a high claim equilibrium, players have to benefit from playing a high claim. For a population setting, it means that players need to have the chance to encounter opponents also playing high. From a learning model point of view, it refers to the belief that the opponent will also play high, at least with a certain probability. If the belief is shared by both players, they should both play high and reach the cooperative equilibrium. Here, we explore these two types of models to unveil the mechanism leading to cooperation.Population based models unveil diversity as the cooperative mechanism via the effect of mutations on the game’s outcome. This is shown for the replicator-mutator equation and the Wright–Fisher model. Similarly, a two-player learning model approach, more in line with human reasoning, shows that if players are free to adopt a higher payoff action from a diverse action set during their introspection process, they can reach the cooperative equilibrium. This result is obtained using introspection dynamics.Finally, we explain how diversity is the underlying mechanism that enables the convergence to high claims in previously proposed models. To be more precise, we show that diversity is required because it allows for the maximisation process to jump from the local to the global level. More

Pardo, J. D., Lennie, K. & Anderson, J. S. Can we reliably calibrate deep nodes in the tetrapod tree? Case studies in deep tetrapod divergences. Front. Genet. 11, 1159 (2020).Article 

Google Scholar 
Rage, J.-C. & Roček, Z. Redescription of Triadobatrachus massinoti (Piveteau, 1936) an anuran amphibian from the early Triassic. Palaeontographica A 206, 1–16 (1989).
Google Scholar 
Evans, S. E. & Borsuk-Białynicka, M. A stem-group frog from the Early Triassic of Poland. Acta Palaeontol. Pol. 43, 573–580 (1998).Article 

Google Scholar 
Heckert, A. B., Mitchell, J. S., Schneider, V. P. & Olsen, P. E. Diverse new microvertebrate assemblage from the Upper Triassic Cumnock Formation, Sanford Subbasin, North Carolina, USA. J. Paleontol. 86, 368–390 (2012).Article 

Google Scholar 
Stocker, M. R. et al. The earliest equatorial record of frogs from the Late Triassic of Arizona. Biol. Lett. 15, 20180922 (2019).Article 

Google Scholar 
Schoch, R. R., Werneburg, R. & Voigt, S. A Triassic stem-salamander from Kyrgyzstan and the origin of salamanders. Proc. Natl Acad. Sci. USA 117, 11584–11588 (2020).Article 
ADS 
CAS 

Google Scholar 
Anderson, J. S., Reisz, R. R., Scott, D., Fröbisch, N. B. & Sumida, S. S. A stem batrachian from the Early Permian of Texas and the origin of frogs and salamanders. Nature 453, 515–518 (2008).Article 
ADS 
CAS 

Google Scholar 
Anderson, J. S. Focal review: the origin(s) of modern amphibians. Evol. Biol. 35, 231–247 (2008).Article 

Google Scholar 
Sigurdsen, T. & Bolt, J. R. The Lower Permian amphibamid Doleserpeton (Temnospondyli: Dissorophoidea), the interrelationships of amphibamids, and the origin of modern amphibians. J. Vertebr. Paleontol. 30, 1360–1377 (2010).Article 

Google Scholar 
Schoch, R. R. The putative lissamphibian stem-group: phylogeny and evolution of the dissorophoid temnospondyls. J. Paleontol. 93, 137–156 (2019).Article 

Google Scholar 
Jenkins, P. A. & Walsh, D. M. An Early Jurassic caecilian with limbs. Nature 365, 246–250 (1993).Article 
ADS 

Google Scholar 
Jenkins, F. A., Walsh, D. M. & Carroll, R. L. Anatomy of Eocaecilia micropodia, a limbed caecilian of the Early Jurassic. Bull. Mus. Comp. Zool. 158, 285–365 (2007).Article 

Google Scholar 
Maddin, H. C., Jenkins, F. A. Jr & Anderson, J. S. The braincase of Eocaecilia micropodia (Lissamphibia, Gymnophiona) and the origin of caecilians. PLoS ONE 7, e50743 (2012).Article 
ADS 
CAS 

Google Scholar 
Pardo, J. D., Small, B. J. & Huttenlocker, A. K. Stem caecilian from the Triassic of Colorado sheds light on the origins of Lissamphibia. Proc. Natl Acad. Sci. USA 114, E5389–E5395 (2017).Article 
ADS 
CAS 

Google Scholar 
Nussbaum, R. A. The evolution of a unique dual jaw‐closing mechanism in caecilians: (Amphibia: Gymnophiona) and its bearing on caecilian ancestry. J. Zool. 199, 545–554 (1983).Article 

Google Scholar 
Kleinteich, T., Haas, A. & Summers, A. P. Caecilian jaw-closing mechanics: integrating two muscle systems. J. R. Soc. Interface 5, 1491–1504 (2008).Article 

Google Scholar 
Sherratt, E., Gower, D. J., Klingenberg, C. P. & Wilkinson, M. Evolution of cranial shape in caecilians (Amphibia: Gymnophiona). Evol. Biol. 41, 528–545 (2014).Article 

Google Scholar 
Schmidt, A. & Wake, M. H. Olfactory and vomeronasal systems of caecilians (Amphibia: Gymnophiona). J. Morphol. 205, 255–268 (1990).Article 

Google Scholar 
Pincheira‐Donoso, D., Meiri, S., Jara, M., Olalla‐Tárraga, M. Á. & Hodgson, D. J. Global patterns of body size evolution are driven by precipitation in legless amphibians. Ecography 42, 1682–1690 (2019).Article 

Google Scholar 
San Mauro, D., Vences, M., Alcobendas, M., Zardoya, R. & Meyer, A. Initial diversification of living amphibians predated the breakup of Pangaea. Am. Nat. 165, 590–599 (2005).Article 

Google Scholar 
Padian, K. & Sues, H.-D. in Great Transformations in Vertebrate Evolution (eds Dial, K. P., Shubin, N. & Brainerd, E. L.) 351–374 (Univ. Chicago Press, 2021).Santos, R. O., Laurin, M. & Zaher, H. A review of the fossil record of caecilians (Lissamphibia: Gymnophionomorpha) with comments on its use to calibrate molecular timetrees. Biol. J. Linn. Soc. 131, 737–755 (2020).Article 

Google Scholar 
Evans, S. E. & Sigogneau‐Russell, D. A stem‐group caecilian (Lissamphibia: Gymnophiona) from the Lower Cretaceous of North Africa. Palaeontology 44, 259–273 (2001).Article 

Google Scholar 
Ramezani, J. et al. High-precision U-Pb zircon geochronology of the Late Triassic Chinle Formation, Petrified Forest National Park (Arizona, USA): temporal constraints on the early evolution of dinosaurs. GSA Bull. 123, 2142–2159 (2011).Article 
CAS 

Google Scholar 
Rasmussen, C. et al. U-Pb zircon geochronology and depositional age models for the Upper Triassic Chinle Formation (Petrified Forest National Park, Arizona, USA): implications for Late Triassic paleoecological and paleoenvironmental change. GSA Bull. 133, 539–558 (2021).Article 
CAS 

Google Scholar 
Nordt, L., Atchley, S. & Dworkin, S. Collapse of the Late Triassic megamonsoon in western equatorial Pangea, present-day American Southwest. GSA Bull. 127, 1798–1815 (2015).Article 
CAS 

Google Scholar 
Martz, J. W. & Parker, W. G. in Terrestrial Depositional Systems (eds Zeigler, K. E. & Parker, W. G.) 39–125 (Elsevier, 2017).Daza, J. D. et al. Enigmatic amphibians in mid-Cretaceous amber were chameleon-like ballistic feeders. Science 370, 687–691 (2020).Article 
ADS 
CAS 

Google Scholar 
Gardner, J. D. Monophyly and affinities of albanerpetontid amphibians (Temnospondyli; Lissamphibia). Zool. J. Linn. Soc. 131, 309–352 (2001).Article 

Google Scholar 
Bolt, J. R. Lissamphibian origins: possible protolissamphibian from the Lower Permian of Oklahoma. Science 166, 888–891 (1969).Article 
ADS 
CAS 

Google Scholar 
Gardner, J. D. & Averianov, A. O. Albanerpetontid amphibians from the Upper Cretaceous of Middle Asia. Acta Palaeontol. Pol. 43, 453–476 (1998).
Google Scholar 
Carroll, R. L. The Palaeozoic ancestry of salamanders, frogs and caecilians. Zool. J. Linn. Soc. 150, 1–140 (2007).Article 

Google Scholar 
Müller, H., Oommen, O. V. & Bartsch, P. Skeletal development of the direct-developing caecilian Gegeneophis ramaswamii (Amphibia: Gymnophiona: Caeciliidae). Zoomorphology 124, 171–188 (2005).Article 

Google Scholar 
Ahlberg, P. E. & Clack, J. A. Lower jaws, lower tetrapods—a review based on the Devonian genus Acanthostega. Earth Environ. Sci. Trans. R. Soc. Edinb. 89, 11–46 (1998).Article 

Google Scholar 
Bolt, J. R. & Lombard, R. E. The mandible of the primitive tetrapod Greererpeton, and the early evolution of the tetrapod lower jaw. J. Paleontol. 75, 1016–1042 (2001).Article 

Google Scholar 
Shishkin, M. A. & Sulej, T. The Early Triassic temnospondyls of the Czatkowice 1 tetrapod assemblage. Acta Palaeontol. Pol. 65, 31–77 (2009).
Google Scholar 
Anderson, J. S., Scott, D. & Reisz, R. R. The anatomy of the dermatocranium and mandible of Cacops aspidephorus Williston, 1910 (Temnospondyli: Dissorophidae), from the Lower Permian of Texas. J. Vertebr. Paleontol. 40, e1776720 (2020).Article 

Google Scholar 
Wilkinson, M., San Mauro, D., Sherratt, E. & Gower, D. J. A nine-family classification of caecilians (Amphibia: Gymnophiona). Zootaxa 2874, 41–64 (2011).Article 

Google Scholar 
Jared, C. et al. Skin gland concentrations adapted to different evolutionary pressures in the head and posterior regions of the caecilian Siphonops annulatus. Sci. Rep. 8, 3576 (2018).Article 
ADS 

Google Scholar 
O’Reilly, J. C., Ritter, D. A. & Carrier, D. R. Hydrostatic locomotion in a limbless tetrapod. Nature 386, 269–272 (1997).Article 
ADS 

Google Scholar 
Muttoni, G. & Kent, D. V. Jurassic monster polar shift confirmed by sequential paleopoles from Adria, promontory of Africa. J. Geophys. Res. 124, 3288–3306 (2019).Article 
ADS 

Google Scholar 
Parsons, T. S. & Williams, E. E. The relationships of the modern Amphibia: a re-examination. Q. Rev. Biol. 38, 26–53 (1963).Article 

Google Scholar 
Marjanović, D. & Laurin, M. A reevaluation of the evidence supporting an unorthodox hypothesis on the origin of extant amphibians. Contrib. Zool. 77, 149–199 (2008).Article 

Google Scholar 
Jenkins, X. A. et al. Using manual ungual morphology to predict substrate use in the Drepanosauromorpha and the description of a new species. J. Vertebr. Paleontol. 40, e1810058 (2020).Article 

Google Scholar 
Kligman, B. T., Marsh, A. D., Nesbitt, S. J., Parker, W. G. & Stocker, M. R. New trilophosaurid species demonstrates a decline in allokotosaur diversity across the Adamanian–Revueltian boundary in the Late Triassic of western North America. Palaeodiversity 13, 25–37 (2020).Article 

Google Scholar 
Marsh, A. D., Smith, M. E., Parker, W. G., Irmis, R. B. & Kligman, B. T. Skeletal anatomy of Acaenasuchus geoffreyi Long and Murry, 1995 (Archosauria: Pseudosuchia) and its implications for the origin of the aetosaurian carapace. J. Vertebr. Paleontol. 40, e1794885 (2020).Article 

Google Scholar 
Marsh, A. D. & Parker, W. G. New dinosauromorph specimens from Petrified Forest National Park and a global biostratigraphic review of Triassic dinosauromorph body fossils. PaleoBios https://doi.org/10.5070/P9371050859 (2020).Kligman, B. T., Marsh, A. D., Sues, H.-D. & Sidor, C. A. A new non-mammalian eucynodont from the Chinle Formation (Triassic: Norian), and implications for the early Mesozoic equatorial cynodont record. Biol. Lett. 16, 20200631 (2020).Article 

Google Scholar 
Huttenlocker, A. K., Pardo, J. D., Small, B. J. & Anderson, J. S. Cranial morphology of recumbirostrans (Lepospondyli) from the Permian of Kansas and Nebraska, and early morphological evolution inferred by micro-computed tomography. J. Vertebr. Paleontol. 33, 540–552 (2013).Article 

Google Scholar 
Pardo, J. D., Szostakiwskyj, M., Ahlberg, P. E. & Anderson, J. S. Hidden morphological diversity among early tetrapods. Nature 546, 642–645 (2017).Article 
ADS 
CAS 

Google Scholar 
Marjanović, D. & Laurin, M. Phylogeny of Paleozoic limbed vertebrates reassessed through revision and expansion of the largest published relevant data matrix. PeerJ 6, e5565 (2019).Article 

Google Scholar 
Goloboff, P. A. & Catalano, S. A. TNT version 1.5, including a full implementation of phylogenetic morphometrics. Cladistics 32, 221–238 (2016).Article 

Google Scholar 
Huelsenbeck, J. P. & Ronquist, F. MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics 17, 754–755 (2001).Article 
CAS 

Google Scholar 
Lewis, P. O. A likelihood approach to estimating phylogeny from discrete morphological character data. Syst. Biol. 50, 913–925 (2001).Article 
CAS 

Google Scholar 
Eltink, E., Schoch, R. R. & Langer, M. C. Interrelationships, palaeobiogeography and early evolution of Stereospondylomorpha (Tetrapoda: Temnospondyli). J. Iber. Geol. 45, 251–267 (2019).Article 

Google Scholar 
Bystrow, A. Dvinosaurus als neotenische Form der Stegocephalen. Acta Zool. 19, 209–295 (1938).Article 

Google Scholar 
Dutuit, J.-M. Introduction à l’étude paléontologique du Trias continental Marocain. Description des premiers stegocephales recueillis dans le couloir d’Argana (Atlas Occidental). Mémoires du Muséum National d’Histoire 36, 1–253 (1976).
Google Scholar 
Dias, E. V., Dias-da-Silva, S. & Schultz, C. L. A new short-snouted rhinesuchid from the Permian of southern Brazil. Revista Brasileira de Paleontologia 23, 98–122 (2020).Article 

Google Scholar 
Damiani, R. J. & Kitching, J. W. A new brachyopid temnospondyl from the Cynognathus Assemblage Zone, Upper Beaufort Group, South Africa. J. Vertebr. Paleontol. 23, 67–78 (2003).Article 

Google Scholar 
Schoch, R. R. & Witzmann, F. Cranial morphology of the plagiosaurid Gerrothorax pulcherrimus as an extreme example of evolutionary stasis. Lethaia 45, 371–385 (2012).Article 

Google Scholar 
Schoch, R. R. Studies on braincases of early tetrapods: Structure, morphological diversity, and phylogeny-1 Trimerorhacis and other prmitive temnospondyls. Neues Jahrbuch für Geologie und Paläontologie-Abhandlungen 213, 233–259 (1999).Article 

Google Scholar 
Ruta, M. & Bolt, J. R. The brachyopoid Hadrokkosaurus bradyi from the early Middle Triassic of Arizona, and a phylogenetic analysis of lower jaw characters in temnospondyl amphibians. Acta Palaeontol. Pol. 53, 579–592 (2008).Article 

Google Scholar 
Bystrow, A. & Efremov, J. Benthosuchus sushkini Efr.—a labyrinthodont from the Eotriassic of Sharzhenga River. Trudy Paleontol. Inst. 10, 1–152 (1940).
Google Scholar 
Warren, A. Karoo tupilakosaurid: a relict from Gondwana. Earth Environ. Sci. Trans. R. Soc. Edinb. 89, 145–160 (1998).Article 

Google Scholar 
Holmes, R. B., Carroll, R. L. & Reisz, R. R. The first articulated skeleton of Dendrerpeton acadianum (Temnospondyli, Dendrerpetontidae) from the Lower Pennsylvanian locality of Joggins, Nova Scotia, and a review of its relationships. J. Vertebr. Paleontol. 18, 64–79 (1998).Article 

Google Scholar 
Steyer, J. S. The first articulated trematosaur ‘amphibian’ from the Lower Triassic of Madagascar: implications for the phylogeny of the group. Palaeontol. 45, 771–793 (2002).Article 

Google Scholar 
Englehorn, J., Small, B. J. & Huttenlocker, A. A redescription of Acroplous vorax (Temnospondyli: Dvinosauria) based on new specimens from the Early Permian of Nebraska and Kansas, USA. J. Vertebr. Paleontol. 28, 291–305 (2008).Article 

Google Scholar 
Warren, A. Laidleria uncovered: a redescription of Laidleria gracilis Kitching (1957), a temnospondyl from the Cynognathus Zone of South Africa. Zool. J. Linn. Soc. 122, 167–185 (1998).Article 

Google Scholar 
Bolt, J. R. & Chatterjee, S. A new temnospondyl amphibian from the Late Triassic of Texas. J. Paleontol. 74, 670–683 (2000).Article 

Google Scholar 
Milner, A. & Sequeira, S. The temnospondyl amphibians from the Viséan of east Kirkton, West Lothian, Scotland. Earth Environ. Sci. Trans. R. Soc. Edinb. 84, 331–361 (1993).
Google Scholar 
Schoch, R. R. & Milner, A. R. Encyclopedia of Paleoherpetology, Part 3A. Temnospondyli (Verlag Dr. Friedrich Pfeil, 2014).Damiani, R., Schoch, R. R., Hellrung, H., Werneburg, R. & Gastou, S. The plagiosaurid temnospondyl Plagiosuchus pustuliferus (Amphibia: Temnospondyli) from the Middle Triassic of Germany: anatomy and functional morphology of the skull. Zool. J. Linn. Soc. 155, 348–373 (2009).Article 

Google Scholar 
Chernin, S. A new brachyopid, Batrachosuchus concordi sp. nov. from the Upper Luangwa Valley, Zambia with a redescription of Batrachosuchus browni Broom, 1903. Palaeontol. Afr. 20, 87–109 (1977).
Google Scholar 
Sulej, T. Osteology, variability, and evolution of Metoposaurus, a temnospondyl from the Late Triassic of Poland. Acta Palaeontol. Pol. 64, 29–139 (2007).
Google Scholar  More

Bradford, M. A. et al. Soil carbon science for policy and practice. Nat. Sustain. 2, 1070–1072 (2019).Article 

Google Scholar 
Lehmann, J. & Kleber, M. The contentious nature of soil organic matter. Nature 528, 60–68 (2015).Article 

Google Scholar 
Liang, C., Schimel, J. P. & Jastrow, J. D. The importance of anabolism in microbial control over soil carbon storage. Nat. Microbiol. 2, 17105 (2017).Article 

Google Scholar 
Liang, C., Amelung, W., Lehmann, J. & Kästner, M. Quantitative assessment of microbial necromass contribution to soil organic matter. Glob. Change Biol. 25, 3578–3590 (2019).Article 

Google Scholar 
Wang, B. R., An, S. S., Liang, C., Liu, Y. & Kuzyakov, Y. Microbial necromass as the source of soil organic carbon in global ecosystems. Soil Biol. Biochem. 162, 108422 (2021).Article 

Google Scholar 
Kästner, M. & Miltner, A. in The Future of Soil Carbon (eds Garcia, C. et al.) Ch. 5 (Academic Press, 2018).Buckeridge, K. M. et al. Sticky dead microbes: rapid abiotic retention of microbial necromass in soil. Soil Biol. Biochem. 149, 107929 (2020).Article 

Google Scholar 
Kallenbach, C. M., Grandy, A. S., Frey, S. D. & Diefendorf, A. F. Microbial physiology and necromass regulate agricultural soil carbon accumulation. Soil Biol. Biochem. 91, 279–290 (2015).Article 

Google Scholar 
Kallenbach, C. M., Frey, S. D. & Grandy, A. S. Direct evidence for microbial-derived soil organic matter formation and its ecophysiological controls. Nat. Commun. 7, 13630 (2016).Article 

Google Scholar 
Emerson, J. B. et al. Schrödinger’s microbes: tools for distinguishing the living from the dead in microbial ecosystems. Microbiome 5, 86 (2017).Article 

Google Scholar 
Zhang, Y. et al. Simulating measurable ecosystem carbon and nitrogen dynamics with the mechanistically defined MEMS 2.0 model. Biogeosciences 18, 3147–3171 (2021).Article 

Google Scholar 
Ackermann, M., Stearns Stephen, C. & Jenal, U. Senescence in a bacterium with asymmetric division. Science 300, 1920–1920 (2003).Article 

Google Scholar 
Aguilaniu, H., Gustafsson, L., Rigoulet, M. & Nyström, T. Asymmetric inheritance of oxidatively damaged proteins during cytokinesis. Science 299, 1751–1753 (2003).Article 

Google Scholar 
Maheshwari, R. & Navaraj, A. Senescence in fungi: the view from Neurospora. FEMS Microbiol. Lett. 280, 135–143 (2008).Article 

Google Scholar 
See, C. R. et al. Hyphae move matter and microbes to mineral microsites: integrating the hyphosphere into conceptual models of soil organic matter stabilization. Glob. Change Biol. 28, 2527–2540 (2022).Article 

Google Scholar 
Pusztahelyi, T. et al. Comparative studies of differential expression of chitinolytic enzymes encoded by chiA, chiB, chiC and nagA genes in Aspergillus nidulans. Folia Microbiologica 51, 547–554 (2006).Article 

Google Scholar 
Bartoszewska, M. & Kiel, J. A. The role of macroautophagy in development of filamentous fungi. Antioxid. Redox Signal. 14, 2271–2287 (2011).Article 

Google Scholar 
Josefsen, L. et al. Autophagy provides nutrients for nonassimilating fungal structures and is necessary for plant colonization but not for infection in the necrotrophic plant pathogen Fusarium graminearum. Autophagy 8, 326–337 (2012).Article 

Google Scholar 
Heaton, L. L., Jones, N. S. & Fricker, M. D. Energetic constraints on fungal growth. Am. Nat. 187, E27–E40 (2016).Article 

Google Scholar 
Taiz, L. & Zeiger, E. Plant Physiology 4th edn (Spektrum Akademischer Verlag, 2008).Bowman, E. J. & Bowman, B. J. in Cellular and Molecular Biology of Filamentous Fungi (eds Borkovich, K. & Ebbole, D.) 179–190 (ASM Press, 2010).Voigt, O. & Pöggeler, S. Self-eating to grow and kill: autophagy in filamentous ascomycetes. Appl. Microbiol. Biotechnol. 97, 9277–9290 (2013).Article 

Google Scholar 
Grimmett, I. J., Shipp, K. N., Macneil, A. & Barlocher, F. Does the growth rate hypothesis apply to aquatic hyphomycetes? Fungal Ecol. 6, 493–500 (2013).Article 

Google Scholar 
Camenzind, T., Philipp Grenz, K., Lehmann, J. & Rillig, M. C. Soil fungal mycelia have unexpectedly flexible stoichiometric C:N and C:P ratios. Ecol. Lett. 24, 208–218 (2021).Article 

Google Scholar 
Mason-Jones, K., Robinson, S. L., Veen, G. F., Manzoni, S. & van der Putten, W. H. Microbial storage and its implications for soil ecology. ISME J. 16, 617–629 (2022).Article 

Google Scholar 
Gow, N. A. R., Latge, J. P. & Munro, C. A. The fungal cell wall: structure, biosynthesis, and function. Microbiol. Spectr. 5, FUNK-0035–2016 (2017).Article 

Google Scholar 
Steiner, U. K. Senescence in bacteria and its underlying mechanisms. Front. Cell Dev. Biol. 9, 668915 (2021).Article 

Google Scholar 
Allocati, N., Masulli, M., Di Ilio, C. & De Laurenzi, V. Die for the community: an overview of programmed cell death in bacteria. Cell Death Dis. 6, e1609 (2015).Article 

Google Scholar 
Peeters, S. H. & de Jonge, M. I. For the greater good: programmed cell death in bacterial communities. Microbiol. Res. 207, 161–169 (2018).Article 

Google Scholar 
Wang, J. & Bayles, K. W. Programmed cell death in plants: lessons from bacteria? Trends Plant Sci. 18, 133–139 (2013).Article 

Google Scholar 
Nagamalleswari, E., Rao, S., Vasu, K. & Nagaraja, V. Restriction endonuclease triggered bacterial apoptosis as a mechanism for long time survival. Nucleic Acids Res. 45, 8423–8434 (2017).Article 

Google Scholar 
Kysela, D. T., Brown, P. J. B., Huang, K. C. & Brun, Y. V. Biological consequences and advantages of asymmetric bacterial growth. Annu. Rev. Microbiol. 67, 417–435 (2013).Article 

Google Scholar 
Bayles, K. W. Bacterial programmed cell death: making sense of a paradox. Nat. Rev. Microbiol. 12, 63–69 (2014).Article 

Google Scholar 
Flemming, H.-C. & Wuertz, S. Bacteria and archaea on Earth and their abundance in biofilms. Nat. Rev. Microbiol. 17, 247–260 (2019).Article 

Google Scholar 
Coleman, D. C. & Wall, D. H. in Soil Microbiology, Ecology and Biochemistry 4th edn (ed. Paul, E. A.) Ch. 5 (Academic Press, 2015).Hungate, B. A. et al. The functional significance of bacterial predators. mBio 12, e00466-21 (2021).Article 

Google Scholar 
Kuzyakov, Y. & Mason-Jones, K. Viruses in soil: nano-scale undead drivers of microbial life, biogeochemical turnover and ecosystem functions. Soil Biol. Biochem. 127, 305–317 (2018).Article 

Google Scholar 
Williamson, K. E., Fuhrmann, J. J., Wommack, K. E. & Radosevich, M. Viruses in soil ecosystems: an unknown quantity within an unexplored territory. Annu. Rev. Virol. 4, 201–219 (2017).Article 

Google Scholar 
Sokol, N. W. et al. Life and death in the soil microbiome: how ecological processes influence biogeochemistry. Nat. Rev. Microbiol. 20, 415–430 (2022).Article 

Google Scholar 
Bonkowski, M. & Clarholm, M. J. A. P. Stimulation of plant growth through interactions of bacteria and protozoa: testing the auxiliary microbial loop hypothesis. Acta Protozool. 51, 237–247 (2012).
Google Scholar 
Potapov, A. M., Pollierer, M. M., Salmon, S., Šustr, V. & Chen, T.-W. Multidimensional trophic niche revealed by complementary approaches: gut content, digestive enzymes, fatty acids and stable isotopes in Collembola. J. Anim. Ecol. 90, 1919–1933 (2021).Article 

Google Scholar 
Esteban, G. F. & Fenchel, T. M. in Ecology of Protozoa: The Biology of Free-living Phagotrophic Protists (eds Esteban, G. F. & Fenchel, T. M.) 33–54 (Springer, 2020).Koksharova, O. A. Bacteria and phenoptosis. Biochemistry 78, 963–970 (2013).
Google Scholar 
Tilman, D. Resource Competition and Community Structure (Princeton Univ. Press, 1982).Boddy, L. Interspecific combative interactions between wood-decaying basidiomycetes. FEMS Microbiol. Ecol. 31, 185–194 (2000).Article 

Google Scholar 
Hibbing, M. E., Fuqua, C., Parsek, M. R. & Peterson, S. B. Bacterial competition: surviving and thriving in the microbial jungle. Nat. Rev. Microbiol. 8, 15–25 (2010).Article 

Google Scholar 
Müller, S. et al. Predation by Myxococcus xanthus induces Bacillus subtilis to form spore-filled megastructures. Appl. Environ. Microbiol. 81, 203–210 (2015).Article 

Google Scholar 
Laskowska, E. & Kuczynska-Wisnik, D. New insight into the mechanisms protecting bacteria during desiccation. Curr. Genet. 66, 313–318 (2020).Article 

Google Scholar 
Rillig, M. C., Ryo, M. & Lehmann, A. Classifying human influences on terrestrial ecosystems. Glob. Change Biol. 27, 2273–2278 (2021).Article 

Google Scholar 
Dörr, T., Moynihan, P. J. & Mayer, C. Bacterial cell wall structure and dynamics. Front. Microbiol. 10, 02051 (2019).Article 

Google Scholar 
Corredor, B., Lang, B. & Russell, D. Effects of nitrogen fertilization on soil fauna in a global meta-analysis. Preprint at Res. Sq. https://doi.org/10.21203/rs.3.rs-1438491/v1 (2022).Blankinship, J. C., Niklaus, P. A. & Hungate, B. A. A meta-analysis of responses of soil biota to global change. Oecologia 165, 553–565 (2011).Article 

Google Scholar 
Manzoni, S., Chakrawal, A., Spohn, M. & Lindahl, B. D. Modeling microbial adaptations to nutrient limitation during litter decomposition. Front. For. Glob. Change 4, 686945 (2021).Article 

Google Scholar 
Frank, D. et al. Effects of climate extremes on the terrestrial carbon cycle: concepts, processes and potential future impacts. Glob. Change Biol. 21, 2861–2880 (2015).Article 

Google Scholar 
Gunina, A. & Kuzyakov, Y. From energy to (soil organic) matter. Glob. Change Biol. 28, 2169–2182 (2022).Article 

Google Scholar 
Fernandez, C. W. & Koide, R. T. Initial melanin and nitrogen concentrations control the decomposition of ectomycorrhizal fungal litter. Soil Biol. Biochem. 77, 150–157 (2014).Article 

Google Scholar 
Kästner, M., Miltner, A., Thiele-Bruhn, S. & Liang, C. Microbial necromass in soils—linking microbes to soil processes and carbon turnover. Front. Environ. Sci. 9, 756378 (2021).Article 

Google Scholar 
Buckeridge, K. M., Creamer, C. & Whitaker, J. Deconstructing the microbial necromass continuum to inform soil carbon sequestration. Funct. Ecol. 36, 1396–1410 (2022).Article 

Google Scholar 
Lehmann, J. et al. Persistence of soil organic carbon caused by functional complexity. Nat. Geosci. 13, 529–534 (2020).Article 

Google Scholar 
Blazewicz, S. J. et al. Taxon-specific microbial growth and mortality patterns reveal distinct temporal population responses to rewetting in a California grassland soil. ISME J. 14, 1520–1532 (2020).Article 

Google Scholar 
Kallenbach, C. M., Wallenstein, M. D., Schipanksi, M. E. & Grandy, A. S. Managing agroecosystems for soil microbial carbon use efficiency: ecological unknowns, potential outcomes, and a path forward. Front. Microbiol. 10, 1146 (2019).Article 

Google Scholar 
Liang, C. Soil microbial carbon pump: mechanism and appraisal. Soil Ecol. Lett. 2, 241–254 (2020).Article 

Google Scholar 
Sinsabaugh, R. L., Manzoni, S., Moorhead, D. L. & Richter, A. Carbon use efficiency of microbial communities: stoichiometry, methodology and modelling. Ecol. Lett. 16, 930–939 (2013).Article 

Google Scholar 
van Groenigen, J. W. et al. Sequestering soil organic carbon: a nitrogen dilemma. Environ. Sci. Technol. 51, 4738–4739 (2017).Article 

Google Scholar 
Greenlon, A. et al. Quantitative stable-isotope probing (qSIP) with metagenomics links microbial physiology and activity to soil moisture in Mediterranean-climate grassland ecosystems (in the press).Mafla-Endara, P. M. et al. Microfluidic chips provide visual access to in situ soil ecology. Commun. Biol. 4, 889 (2021).Article 

Google Scholar 
Schaible, G. A., Kohtz, A. J., Cliff, J. & Hatzenpichler, R. Correlative SIP-FISH-Raman-SEM-NanoSIMS links identity, morphology, biochemistry, and physiology of environmental microbes. ISME Commun. 2, 52 (2022).Article 

Google Scholar 
See, C. R. et al. Distinct carbon fractions drive a generalisable two-pool model of fungal necromass decomposition. Funct. Ecol. 35, 796–806 (2021).Article 

Google Scholar 
Wang, C. et al. Stabilization of microbial residues in soil organic matter after two years of decomposition. Soil Biol. Biochem. 141, 107687 (2020).Article 

Google Scholar 
Veresoglou, S. D., Halley, J. M. & Rillig, M. C. Extinction risk of soil biota. Nat. Commun. 6, 8862 (2015).Article 

Google Scholar 
Potapov, A. M. et al. Feeding habits and multifunctional classification of soil-associated consumers from protists to vertebrates. Biol. Rev. 97, 1057–1117 (2022).Article 

Google Scholar 
Trap, J., Bonkowski, M., Plassard, C., Villenave, C. & Blanchart, E. Ecological importance of soil bacterivores for ecosystem functions. Plant Soil 398, 1–24 (2016).Article 

Google Scholar 
Dooley, S. R. & Treseder, K. K. The effect of fire on microbial biomass: a meta-analysis of field studies. Biogeochemistry 109, 49–61 (2012).Article 

Google Scholar 
Muñoz-Leoz, B., Ruiz-Romera, E., Antigüedad, I. & Garbisu, C. Tebuconazole application decreases soil microbial biomass and activity. Soil Biol. Biochem. 43, 2176–2183 (2011).Article 

Google Scholar 
Meyer, M., Diehl, D., Schaumann, G. E. & Muñoz, K. Agricultural mulching and fungicides—impacts on fungal biomass, mycotoxin occurrence, and soil organic matter decomposition. Environ. Sci. Pollut. Res. 28, 36535–36550 (2021).Article 

Google Scholar 
Thiery, S. & Kaimer, C. The predation strategy of Myxococcus xanthus. Front. Microbiol. 11, 2 (2020).Article 

Google Scholar 
Laloux, G. Shedding light on the cell biology of the predatory bacterium Bdellovibrio bacteriovorus. Front. Microbiol. 10, 3136 (2020).Article 

Google Scholar  More

Our study is focussed on population trends of large carnivores; a culturally important group32, essential for regulating ecosystem function33. Large carnivores represent an important study group as their population status is unclear, with reports of devastating declines33 contrasted with remarkable recoveries23. Further, as a well-studied taxa with abundant trend and trait datasets, large carnivores present a good system to evaluate important drivers of trends without being impacted by poor inference from missing data34. Finally, as large carnivores are considered indicator species of the overall status of biodiversity within an area35, our inference may provide insight beyond our focal taxa.Population trendsWe sourced population (defined by the authors of the original studies, who reported on population trends for one or more studied groups of individuals) trend information for species in the families Canidae, Felidae, Hyaenidae, and Ursidae of the order Carnivora from two large trend datasets: CaPTrends12 and the Living Planet Database13. The CaPTrends database is the product of a semi-systematic literature search for population trends of large carnivore species (from the families listed above); the dataset possesses trend information for 50 species from locations around the world, and trends are reported in a variety of ways. The Living Planet Database contains population abundance time-series for vertebrates from thousands of sites around the world and is one of the larger population trend datasets. Combined, these datasets produce a cumulative 1123 trends (after removing duplicates and records we deemed unreliable or unsuitable), derived from >10,000 individual population estimates. In the Living Planet Database, and for most records in CaPTrends, trends are reported as a time-series of abundance (or density) estimates. We modelled these time-series with log-linear regressions, where abundance (the response) was loge transformed, and year of abundance estimates was selected as the predictor. We included a continuous Ornstein-Uhlenbeck (OU) autoregressive process to control for temporal autocorrelation in these models. The OU process estimates covariance between abundance values, under the assumption that abundances in time point 1 will be more similar to abundances in time point 2, than time-point 3, 4, 5, etc. Accounting for covariance resolves non-independence within time-series. We extracted the slope coefficient which represents the annual instantaneous rate of change, sometimes called the population growth rate (rt). Alongside the abundance time-series, CaPTrends also has three other quantitative datatypes, all of which we converted into an annual instantaneous rate of change (rt): (1) a mean finite rate of change; (2) estimates of percentage abundance change between two points in time; and (3) time-series’ of population change estimates (e.g. in year 1 the population doubled and in year 2 it halved). We converted all annual instantaneous rates of change into an annual rate of change percentage to improve interpretability. These annual rates of change ranged from −75 to 68%, but the majority of values fell within −10 to 10% (Supplementary Fig. 1a).Alongside the quantitative records, 138 populations in the CaPTrends dataset were only described qualitatively with categories: Increase, Stable, and Decrease. These qualitative records were more common for populations located in traditionally poorer-sampled countries (e.g. with lower human development), so whilst they are less informative (only describing the direction and not the magnitude), we deem them important to reduce bias (Fig. 1). As a result, we used a combination of percentage annual rates of change (N = 985) and qualitative categories (N = 138) as our response in our model (see below), representing 50 large carnivore species.CovariatesFor each population, we extracted sixteen covariates (each z-transformed) that fell into four categories: land-use, climate, governance, and traits. Our covariates were designed to cover a diverse array of factors that could influence population trends in large carnivores (Supplementary Table 1). Each covariate is described briefly in Fig. 2 with full descriptions of how variables were derived in the Supplementary material: Covariates.One of the challenges in identifying how covariates—which can vary in space and time—impact population trends is matching the spatial and temporal scale of the covariate with the population i.e. how much of the population is affected by the covariate at a given point in time. To tackle the spatial element of this problem, we used data on the area of extent of each population (e.g. how large is the spatial extent of the population or monitoring zone) to generate a circular distribution zone around the population’s coordinate centroid. We refer to this as the ‘population area’ hereafter. We then sampled covariate values within each population area, with more sampling points in larger areas (range: 13–295 sampling points, Supplementary Fig. 2b). For covariates which varied over time, we extracted the covariates across the ‘population monitoring period’, which refers to the period (from start to end year) the population was monitored for. However, as evidence suggests there can be a lag period between impact or change and any detectable changes in population abundance3, we tested how 0-, 5-, and 10-year lags in covariates changed model fits and effect sizes. We implemented these lags by extending the start of the population monitoring period backwards for each given lag e.g. for a 10-year lag, a normal population monitoring period of 1990–2000, would then capture covariates between 1980–2000. Sensitivity analysis showed a 10-year lag had the greatest balance of improved model fit, with high taxonomic and spatial coverage (see Supplementary: Sensitivity analysis).ModellingAt its core, our model is a linear mixed effects model, regressing annual rates of change against a combined 23 covariates and interactions, using random intercepts to account for phylogenetic and spatial nesting. The model was written in BUGS language and implemented in JAGS 4.3.036 via R 4.0.337. The model structure is summarised below, with a full description in Supplementary: Modelling.ResponseWe modelled our annual rate of change response with a normal error prior. However, to allow the two different types of population trend data (quantitative rates of change and qualitative descriptions of change) to be included in the same model, we treated the qualitative records as partially known. Specifically, we censored the qualitative records to indicate that the true value is unknown, but it occurs within a specified range, with annual rate changes ranging from −50 to 0%, −5 to 5%, and 0 to 50% within the decrease, stable and increase categories, respectively. The overlapping nature of these thresholds is by design, as we wanted to acknowledge that there is likely a grey area between the different categories. For instance, in one study, a 2% trend could be called stable, whilst a different study would consider this as increasing, our overlapping thresholds address this grey area. Admittedly, our category thresholds were arbitrarily selected—this is as a consequence of there being no strict rules on what population change is needed to be assigned a given category. However, despite being arbitrary, they were still carefully selected. For instance, our censoring range thresholds are similar to the range of the observed change (−75 to 68%). Further, whilst we don’t have a clear definition for what an increasing or decreasing population looks like (is it 1% or 10%), we can be confident that increasing and decreasing populations will fall above and below 0%, respectively. The stable category is most vulnerable to subjectivity, and so without clear definitions, we set a large range e.g. the maximum and minimum value we considered could be plausibly called stable was 5% and −5%, respectively.Many of the qualitative and short-term (brief monitoring period) quantitative records address known data biases as they occur in less-well represented regions, species, and time-periods (Fig. 1). However, these lower quality records can be more prone to error. As a result, we developed a weighting term within the model to inflate uncertainty around trends derived over a short timeframe, with few abundance observations, and less robust methods—see Supplementary: Modelling—Weighted error.CovariatesPrior to modelling, we identified missing values in some covariates (e.g. some species were missing Maximum longevity values), which can be problematic for inference if ignored34. We used imputation approaches38,39 to predict these missing values and recorded the associated imputation uncertainty alongside these predictions. Within our model, we accounted for uncertainty in the imputed estimates by treating imputed values of the covariates as distributions instead of point estimates. Specifically, for each imputed value we assigned a normal distribution defined by the mean and standard deviation of the imputed estimates. This approach allowed us to capture imputation uncertainty and improve inference robustness.With 16 covariates and a further seven interactive effects (23 effects in total), we were conscious of overparameterizing the model. As a result, we split these parameters into three groups: (1) core parameters—which included main effects that were considered likely drivers of population change; (2) optional parameters—which included main effects we considered interesting but with little evidence to-date of any influence on trends; and (3) interaction parameters—which includes the seven proposed interaction terms. We included our core parameters (Change in human density, Primary land loss, Population area, Body mass, Change in extreme heat, Governance, and Protected area coverage) in every model, but used Kuo and Mallick variable selection40 to identify parameters from the optional and interaction groups that improve model fit whilst balancing the risk of overfitting.Random interceptsWe used a hierarchical model structure to account for phylogenetic and spatial non-independence in the data, including species as a random intercept nested with genus, and country as a random intercept nested within sub-regions, as defined by the United Nations (https://www.un.org/about-us/member-states).Model runningWe ran the full model through three chains, each with 150,000 iterations. The first 50,000 iterations in each chain were discarded, and we only stored every 10th iteration along the chain (thinning factor of 10). We opted for a large chain and burn-in due to the model complexity, and to allow a broad selection of parameter combinations to be tested under variable selection. We assessed convergence of the full model on all parameters monitored in the sensitivity analysis, as well as the model intercept, and all 23 main and interactive effect slope coefficients. We checked the standard assumptions of a mixed effect linear model (normal residuals and heterogeneity of variance), and tested the residuals to ensure there was no spatial (Moran’s test) or phylogenetic (Pagel’s lambda) autocorrelation. We also conducted posterior predictive checks to ensure independently simulated values were broadly reminiscent of model predicted values.We report the median slope coefficient and associated credible intervals for each of the main and interactive effects, and produce marginal effect plots for a selection of important parameters. These marginal effects hold all other covariates at zero (which is the equivalent of the mean, as covariates were z-transformed).LimitationsDeveloping macro-scale models of population change is challenging as response data are biased41 and hard to summarise42, and response-covariate relationships are likely complex and numerous2. Within our workflow, we attempted to address these challenges, and overall, this allowed us to achieve a moderate model fit (conditional R2 ~ 0.4). We minimised biases in the trend data by integrating qualitative trends with quantitative estimates, which allowed us to increase the taxonomic and spatial scale of the work. However, biases are likely still present to some extent. For instance, whilst we have population trend data covering the full parameter space of our most influential variable (change in human development), we have more population trends in high human development countries (Supplementary Fig. 20)—given these biases, caution should be used when interpreting results. While we could not avoid some biases, we found inference was similar across different fragments of the data and model structures (Supplementary results: Sensitivity analysis). We also attempted to capture complexity by covering a more comprehensive array of covariates than many previous analyses, but we still lack data on likely important aspects that are cryptic and difficult to measure (e.g. poaching, persecution, culling, and the conservation benefits of being flagship species). Further, there are temporal lags between disturbance-events and observable changes in the population10 and we tested several to incorporate the lag that maximised model fit. However, it is possible that responses to different types of disturbance (e.g. habitat loss and climate change) have different lags, although this has not been quantified. Long lags (the maximum lag we explored was 10-years) may also occur and be associated with slow recoveries, but an absence of longer temporal extents in the response and covariate data largely prohibits this analysis at global scales (long temporal extent data is less available outside of the global north).Counterfactual scenariosTo explore how observed changes in land-use, climate and human development have influenced population trends, we developed three counterfactual scenarios, where we compared observed population change to predicted population change if habitat, climate, and human development remained static. For instance, in the climate change counterfactual scenario, we predicted each population trend using the global model (all covariate parameters) with available covariate data (e.g. land-use, governance and trait covariates), as well as taxa and location data (to provide sensitivity to the models varying random intercepts), but set the climate change covariate data to zero (in this case, change in extreme heat and change in drought). We then subtracted these counterfactual predictions from the observed trends to define ‘Difference in annual rate of change (%)’, whereby a positive value indicates carnivore populations would be in better shape (fewer declines) under the counterfactual scenario, and vice-versa. We summarise counterfactual scenarios by reporting the median Difference in annual rate of change and 95% quantiles across the observed 1123 populations.Socioeconomic development and non-linearity in carnivore trendsGiven the large effect of human development change on carnivore population trends within our counterfactual scenarios, we further explore the potential impacts of human development change (i.e. changes in the socioeconomic standards of society) on the dynamics of potential carnivore abundance change. Specifically, we test how changing the rate of human development growth of a hypothetical low human development country could impact carnivore abundances. We test this by simulating time series of human development change between the years 1960 and 2020 along three common development pathways for low human development countries, each given: (1) a mean rate of change in human development (%) defined as Slow (1.25%), Moderate (1.5%) and Fast (1.75%); (2) a shared deceleration rate set to −0.02% per year—a key feature of the human development data is that as human development grows, its growth rate decreases; and (3) a shared initial human development value which we set as 0.2 (a hypothetical low human development country) at year 1960 (Fig. 4a). All our selected parameter values are representative of the human development data (Supplementary Fig. 2), with the Moderate pathway being largely typical for a country with an initial human development value of 0.2, while Slow and Fast represent plausible extremes.We then used our fitted model (Fig. 2) to evaluate how the three pathways of Change in Human development would affect annual abundance of a hypothetical carnivore. This involved predicting the annual rate of change in abundance using the Change in human development pathways and the marginal effect of the Change in human development parameter from the fitted model—setting all other covariates in the model to zero, which in our z-transformed variables represents the mean. We then used the predicted annual rates of change in abundance to project carnivore abundance up to the year 2020, from an arbitrary baseline abundance of 100 in the year 1960 (Fig. 4c). These projections capture the 95% credible intervals around the human development change model coefficient, and assume constant and average values for all other effects (e.g. primary habitat loss or climate change).Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More

Study site, field experiment, and soil samplingThis study relied on an experimental field station that simulates atmospheric CO2 enrichment and warming (Fig. 1a). The experimental site is located at Kangbo village (31°30′N, 120°33′E), Guli Township, Changshu municipality in Jiangsu Province, China. The climate type is a subtropical monsoon climate with a mean annual precipitation between 1100–1200 mm and an annual average temperature of approximately 16 °C. The high rainfall and temperature mainly occur from May through September25. The soils are Gleyic Stagnic Anthrosols derived from clayey lacustrine deposits. The main properties of the topsoil (0–15 cm) before the field experiment were as follows: pH (H2O) 7.0, bulk density 1.2 g cm−3, and contents of organic C and total N of 16.0 g kg−1 and 1.9 g kg−1, respectively.The FACE system has been well described in previous publications26,53. The climate change treatments were set according to the representative concentration pathway (RCP) scenario that modeled CO2 atmospheric concentration and temperature elevation in approximately 30–40 years. Elevated atmospheric CO2 concentration and warming of crop canopy air were maintained 24 h a day during the crop-growing period. Each treatment was replicated in three rings with the same infrastructure, and the rings were arranged in a split row design (Fig. 1b). All rings were buffered by adjacent open fields to avoid any treatment cross-over. Rice and wheat cultivations of all plots were managed with local conventional practices. Soil samples for qSIP incubation were collected in June 2020.Evaluation of nutrient pulses in soil after water additionTo determine if there is a nutrient pulse after soil rewetting, we estimated the dynamics of several biochemical characteristics (i.e., CO2 fluxes, hydrolase activity, DOC, DN, TC, and TN) in soils before and during the incubation. Soil samples for nutrient flux measurements were collected from the free-air CO2 enrichment and warming experimental station in July 2022.To measure the CO2 production rate, soil (2.00 g) was set in a 29 mm diameter glass vial (Volume 23 mL). Water (400 μL) was added evenly to soil in each vial. Gas samples were then collected from each incubation vial at multiple time points (i.e., 3 h, 6 h, 9 h, 12 h, 24 h, 72 h, and 144 h after water addition), and CO2 concentrations were determined in all gas samples with a Trace Gas Chromatograph (Agilent 7890, Santa Clara, CA, USA). Other biochemical characteristics were measured in parallel incubations, with destructive sampling. Briefly, soil (40.00 g) was set in a plastic jar, and 8 mL water was added evenly to the soil. All the incubation vials and jars were placed in the dark at 25 °C and soil moisture content was adjusted twice a day. Destructive sampling was conducted after 1, 3, and 6 days of incubation. Including prewet soil (the soil before water addition), a total of 48 soil samples were obtained. Soil samples were weighed before and after being oven dried at 105 °C, and soil moisture content was measured using the gravimetric method. Fluorescein Diacetate (FDA) hydrolase activity of soil was measured by using soil FDA hydrolase activity assay kit (Solarbio®, Beijing, China) according to the manufacturer’s protocol. Soil total carbon and total nitrogen were measured with a C/N elemental analyzer (multi EA® 5000, analytikjena, Germany). Dissolved organic carbon and dissolved nitrogen were analyzed using a TOC/TN analyzer (multi N/C® 3100, analytikjena, Germany) after extraction with distilled water.Simulating nutrient pulse after rewetting dry soil by adding 18O-waterTo determine whether and how microbial growth varied in response to pulse events after long-term acclimation to warming and CO2 enrichment, we estimated the population-specific growth rates of active microbes by conducting a 18O-water incubation experiment combined with DNA quantitative stable isotope probing (DNA-qSIP) (Fig. 1c). The incubation conditions were similar to those reported in a previous study6. In brief, approximately 60 g fresh soil of each treatment were sieved (2 mm) and air-dried (24 h at room temperature) immediately after transport to the laboratory. Then, triplicate samples of dry soils (2.00 g) were incubated in the dark at room temperature in sterile plastic aerobic culture tubes (17 × 100 mm) with 400 μL of 98 atom% H218O or natural abundance water (H216O) for 6 days, with harvests at four time points (T = 0, 1, 3, 6 d after water addition). DNA was extracted from the soils of 0 d incubation treatment immediately after water addition (~30 s interval), representing the prewet treatment. At each harvest, soils were destructively sampled and immediately stored at −80 °C. A total of 84 soil samples (4 climate change treatments × 3 replicates at 0 h with H216O addition +  4 treatments × 3 subsequent time points × 3 replicates × 2 types of H2O addition) were collected.DNA extraction and isopycnic centrifugationTotal DNA from all the collected soil samples was extracted using the FastDNA™ SPIN Kit for Soil (MP Biomedicals, Cleveland, OH, USA) according to the manufacturer’s instructions. The concentration of extracted DNA was determined fluorometrically using Qubit® DNA HS (High Sensitivity) Assay Kits (Yeasen Biotechnology, Shanghai, China) on a Qubit® 4 fluorometer (Thermo Scientific™, Waltham, MA, USA). The DNA samples of day 1, day 3, and day 6 were used for isopycnic centrifugation, and the detailed pipeline was described previously with minor modifications6. Briefly, 3 μg DNA were added into 1.85 g mL−1 CsCl gradient buffer (0.1 M Tris-HCl, 0.1 M KCl, 1 mM EDTA, pH = 8.0) with a final buoyant density of 1.718 g mL−1. Approximately 5.1 mL of the solution was transferred to an ultracentrifuge tube (Beckman Coulter QuickSeal, 13 mm × 51 mm) and heat-sealed. All tubes were spun in an Optima XPN-100 ultracentrifuge (Beckman Coulter) using a VTi 65.2 rotor at 177,000 × g at 18 °C for 72 h with minimum acceleration and braking.Immediately after centrifugation, the contents of each ultracentrifuge tube were separated into 20 fractions (~250 μL each fraction) by displacing the gradient medium with sterile water at the top of the tube using a syringe pump (Longer Pump, LSP01‐2A, China). The buoyant density of each fraction was measured using a digital hand-held refractometer (Reichert, Inc., Buffalo, NY, USA) from 10 μL volumes. Fractionated DNA was precipitated from CsCl by adding 500 μL 30% polyethylene glycol (PEG) 6000 and 1.6 M NaCl solution, incubated at 37 °C for 1 h, and then washed twice with 70% ethanol. The DNA of each fraction was then dissolved in 30 μL of Tris‐EDTA buffer. Detailed information (company names and catalog numbers) of the reagents and consumables for qSIP experiment is provided in Supplementary Data 2.Quantitative PCR and sequencingTotal 16S rRNA gene copies for DNA samples of all the fractions and the day-0 soils were quantified using the primers for V4-V5 regions: 515F (5′‐GTG CCA GCM GCC GCG G‐3′) and 907R (5′‐CCG TCA ATT CMT TTR AGT TT‐3′)54. Plasmid standards were prepared by inserting a copy of purified PCR product from soil DNA into Escherichia coli. The E. coli was then cultured, followed by plasmid extraction and purification. The concentration of plasmid was measured using Qubit DNA HS Assay Kits. Standard curves were generated using 10‐fold serial dilutions of the plasmid. Each reaction was performed in a 25-μL volume containing 12.5 μL SYBR Premix Ex Taq (TaKaRa Biotechnology, Otsu, Shiga, Japan), 0.5 μL of forward and reverse primers (10 μM), 0.5 μL of ROX Reference Dye II (50×), 1 μL of template DNA, and 10 μL of sterile water. A two-step thermocycling procedure was performed, which consisted of 30 s at 95 °C, followed by 40 cycles of 5 s at 95 °C, 34 s at 60 °C (at which time the fluorescence signal was collected). Following qPCR cycling, melting curves were obtained to ensure that the results were representative of the target gene. Average PCR efficiency was 96% and the average slope was −3.38, with all standard curves having R2 ≥ 0.99.The DNA of day-0 samples (unfractionated) and the fractionated DNA of fractions with buoyant density between 1.695 and 1.735 g/mL were selected for 16S rRNA amplicon sequencing by using the same primers of qPCR (i.e., 515F/907R). Eleven out of 20 fractions from each ultracentrifuge tube (density between 1.695 and 1.735 g/mL) were selected because they contained more than 99% gene copy numbers of the 20 fractions. A total of 804 DNA samples (12 unfractionated DNA samples + 72 × 11 fractionated DNA samples) were sequenced using the NovaSeq6000 platform (Genesky Biotechnologies, Shanghai, China).The sequences were quality-filtered using the USEARCH v.11.055. In brief, sequences  0.5 were removed. Chimeras were identified and removed. Subsequently, high-quality sequences were clustered into operational taxonomic units (OTUs) using the UPARSE algorithm at a 97% identity threshold, and the most abundant sequence from each OTU was selected as a representative sequence. The taxonomic affiliation of the representative sequence was determined using the RDP classifier (version 16)56. In total, 51,127,459 reads of the bacterial 16S rRNA gene and 11,898 OTUs were obtained. The 16S rRNA amplicon sequences were uploaded to the National Genomics Data Center (NGDC) Genome Sequence Archive (GSA) with accession number CRA006507.Quantitative stable isotope probing calculationsWe used the amount of 18O incorporated into DNA to estimate the growth rates of active taxa24,57. The density shifts of OTUs between 16O and 18O treatments were calculated following the qSIP procedures24. Briefly, the number of 16S rRNA gene copies per OTU in each density fraction was calculated by multiplying the OTU’s relative abundance (acquisition by sequencing) by the total number of 16S rRNA gene copies (acquisition by qPCR). Then, the GC content and molecular weight of a particular taxon were calculated. Further, the change in 18O isotopic composition of 16S rRNA genes for each taxon was estimated. We assumed an exponential growth model over the course of the incubations, and absolute population growth rates were estimated over each of the three time intervals of the incubation: 0–1, 0–3, and 0–6 d corresponding to the samplings at days 1, 3, and 6. The absolute growth rate is a function of the rate of appearance of 18O-labeled 16S rRNA genes. Therefore, the growth rate (g) of taxon i was calculated as:$${g}_{i}={{{{{rm{ln}}}}}}left(frac{{N}_{{{{{{rm{TOTAL}}}}}}{it}}}{{N}_{{{{{{rm{LIGHT}}}}}}{it}}}right)times frac{1}{t}$$
(1)
Where NTOTALit is the number of total gene copies for taxon i and NLIGHTit represents the unlabeled 16S rRNA gene abundances of taxon i at the end of the incubation period (time t). NLIGHTit is calculated by a function with four variables: NTOTALit, molecular weights of DNA (taxon i) in the 16O treatment (MLIGHTi) and in the 18O treatment (MLABi), and the maximum molecular weight of DNA that could result from assimilation of H218O (MHEAVYi)24. We further calculated the average growth rates (represented by the production of new 16S rRNA gene copies of each taxon per g dry soil per day) over each of the three time intervals of the incubation: 0–1, 0–3, and 0–6 d, using the following equation39:$$frac{d{N}_{i}}{{dt}}={N}_{{{{{{rm{TOTAL}}}}}}{it}}left(1-{e}^{-{g}_{i}t}right)times frac{1}{t}$$
(2)
Where t is the incubation time (d). All data calculations were performed using the qSIP package (https://github.com/bramstone/qsip) in R (v. 3.6.2).Grouping of taxa into growth strategiesWe compared the average growth rates of taxa at three time intervals (n = 3 in each time interval) and classified the species into rapid, intermediate, and slow growth strategies based on the timing of the maximum growth rate (Fig. 1d): (1) Rapid responders: Species had the highest growth rates by 1 day of the incubation; (2) Intermediate responders: Species had the highest growth rates at the 3-day incubation; (3) Slow responders: Species had the highest growth rates at the 6-day incubation. The taxa with growth rates significantly greater than zero can be divided into one of three strategies in each treatment.Analyses of phylogenetic conservationPhylogenetic tree analyses were performed in Galaxy/DengLab (http://mem.rcees.ac.cn:8080/) with PyNAST Alignment and FastTree functions58,59. The trees were visualized and edited using iTOL60. To estimate the phylogenetic patterns of growth strategies, phylogenetic dispersion (D) was calculated by the function ‘phylo.d’ of package “caper” in R (v. 3.6.2). The D statistic equal to 1 means the observed trait has a phylogenetically random distribution across the tips of the phylogeny, and 0 means the observed trait is as clumped as if it had evolved by Brownian motion29. Increasing phylogenetic clumping in the binary trait is indicated by values of D decreasing from 1. The p values were obtained by performing 1000 permutations to test D for significant departure from 1 (random distribution). Furthermore, the phylogenetic signal metrics Blomberg’s K and Pagel’s λ, were used to test for significantly nonrandom phylogenetic distributions using the package “phylosignal” in R. The p values of both indices were used to test the significance of phylogenetic signals.The nearest taxon index (NTI) was calculated to determine the degree of phylogenetic clustering as described previously5. The mean nearest taxon distance (MNTD) were calculated in the “picante” package of R61. The values of NTI are equivalent to the negative output of the standardized effect size (SES) of the observed MNTD distances, which test whether the distribution of growth strategies across different phylogenetic groups is random or nonrandom. NTI values > 0 and their p values < 0.05 represent phylogenetic clustering, while NTI values < 0 and p values > 0.95 indicate phylogenetic over-dispersion. The p values of NTI between 0.05 and 0.95 represent random phylogenetic distributions. The data were converted into binary matrices (1 s and 0 s) before all phylogenetic analyses. For the OTUs that were present in more than one treatment and had differing responses (in the analyses when all climate treatments are combined, i.e., in Table 1 and Supplementary Fig. 6), the OTU was classified to the respective growth responder when that taxa exhibited that specific growth strategy in any treatment.Quantification of explained variance in the distribution patterns of strategiesThe phylogenetic distance matrix obtained from the phylogenetic tree and three binary matrices (including four climate treatments) for three growth strategies were decomposed by the package ‘FactoMineR’ in R. The distribution matrices (binary) of three growth responders in the four climate treatments were used to represent the impact of CO2 and temperature on bacterial growth. To predict the microbial growth strategy (i.e., “y”, a 1 or 0 indicating the membership in a response category), variance partitioning was performed by using the first four phylogenetic principal components (PCs) (accounting for over 80% of phylogenetic variance) and two environment PCs (accounting for over 65% of habitat variance). The fraction of the variance explained by phylogenetic constraints and environmental acclimation was calculated by the package ‘car’ in R.Statistical analysesUncertainty of growth rates (95% confidence interval) was estimated using a bootstrapping procedure with 1000 iterations6. The cumulative growth rates at phylum-level were estimated as the sum of taxon-specific growth rates of those OTUs affiliated to the same phylum. The per capita growth rates of each OTU were calculated by dividing absolute growth rates by the total 16S rRNA gene abundance of taxon i. Ecological processes of active populations (e.g., density-dependence) after rewetting dry soil were estimated by correlating species-specific fitness (refers to per capita growth rates of each OTU) with initial population size using linear regression analyses (R v.3.6.2, ‘lm’ function). The beta diversity of the bacterial community at 0 d incubation was visualized by unconstrained principal coordinate analysis (PCoA) and tested by permutational multivariate analysis of variance (PERMANOVA) based on Bray-Curtis distance, using the vegan package in R (v. 3.6.2)62. Comparisons between climate scenarios were tested by two-way ANOVA and LSD post hoc tests (SPSS 19 for Windows, IBM Corp., Armonk, NY, USA).Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More

StatementsWe report our study in accordance with ARRIVE guidelines.Structure of artificial turf of ICUA schematic illustration of a ground plan of the artificial turf sports field of the ICU is shown in Fig. 1. This artificial turf was installed in 2013 by Japanese company B. The field is surrounded by ditches, and there are three drains that connect to sewer pipes. The artificial turf field of TGU was installed in 2011 by Japanese company C.Characterization of rubber tips of artificial turf field of ICU and TGURT were collected from the artificial fields of ICU and TGU. RT for the artificial turf field of ICU were made of residual part of rubber for making tires, window frames and windshields of automobiles. RT of ICU consists of a mixture of EPDM (ethylene-propylene-diene) and SBR (styrene-butadiene rubber) (personal communication from a Japanese company B). The RT of TGU was made of rubber of the residual part of rubber for making tires, window frames, etc. (personal communication from a Japanese company C). Information on raw material of the RT was not manifested.RT collected from the fields (ICU and TGU) was sieved to estimate the particle sizes. The RT of the ICU varied from 0.053 to 3.35 mm, and that of TGU varied from 0.212 to 3.35 mm. The specific gravity of the RT was obtained as follows: A certain amount of RT was weighed and poured into a 10 ml graduated cylinder containing some water. The total volume of the RT was obtained by measuring the rise in the meniscus of the water. The specific gravities of the tested RT were 1.28 (ICU) and 1.28 (TGU). Elemental analyses of RT (ICU and TGU) were conducted using micro-PIXE line analysis47, and calcium, sulfur, zinc, and iron were detected, but lead was under the detection limit from the RT of both ICU and TGU.Sampling of sediments in the ditches of the fieldTo examine the migration of RT from the field to the ditches, approximately 200 g of sediments in the ditches was sampled at four different sites, D1–D4 (Fig. 1), in the ICU. The ditch surrounding the field is made by connecting U-shaped concrete blocks and concrete lids. The inner width, length, and depth of the block are 24, 60, and 24 cm, respectively. The size of the lid is 33 × 60 × 4.5 cm with 1.5 × 9.0 cm snicks at short sides, which make an opening of the ditch of 3.0 × 9.0 cm size between two lids.Each sample was weighed (wet weight) and washed with water using a fine sieve to remove the soil. After the removal of the soil, the sediment was dried, and RT was collected manually. The collected RT was weighed, and the percentages of RT in the sediments were calculated (weight/weight).Goldfish and crucian carpA common variety of goldfish C. auratus of different sizes were obtained from a fish merchant in Saitama Prefecture and from a pet shop in Tokyo and then kept in the ICU. Approximately 200 fish of four different sizes (large, body weight (BW) ~ 100 g; medium, BW, ~ 30 g; small-medium, BW, ~ 15 g; small, BW ~ 4.0 g) were kept in three 800-L stock tanks maintained at 20 °C under a 16-h light/8-h dark (16 L/8 D) photoperiod (lights on at 06:00). Small body size fish were kept in a floating cage in one of the stock tanks. The fish were fed commercial floating goldfish feed (Itosui) once a day ad libitum. The fish stock tanks had circulation filtration systems equipped with sand filters. The filter was cleaned every week to maintain the water quality. The health condition of the fish was judged by their appetite. All the experimental fish (mixed sex) in the present study were kept in stock tanks for over two weeks before they were used for experiments. A total of 127 goldfish were used for the present study. The sample size of each experiment was determined by the results of preliminary experiments. Our preliminary survey confirmed that the fish feed we used did not contain RT-like substances. Therefore, the sample sizes of the control groups (goldfish) were smaller than those of the experimental groups to sacrifice fewer fish. All goldfish and crucian carp experiments were conducted in the ICU.Approximately 30 wild juvenile crucian carp C. auratus subsp. 2 weighing 1.4–4.6 g were obtained from a fish merchant in Saitama Prefecture and kept in an 800-L stock tank in the same conditions as that for goldfish. A total of 16 crucian carp were used for the present study.For the experiments, fish were transferred from the stock tanks to experimental 60-L glass aquaria, which were maintained at 20 °C under a 16-h light/8-h dark (16 L/8 D) photoperiod (lights on at 06:00). The experimental aquaria had a running water system, and dechlorinated tap water was added at 20 ml/min. Plastic box filters were also set to each experimental aquarium to maintain water quality. When stock fish were transferred to experimental aquaria, fish were randomly allocated to the aquaria. All the methods for using goldfish and crucian carp were performed in accordance with the guidelines of the Animal Experimentation Committee of International Christian University. The conduct of the present study was approved by the Animal Experimentation Committee of International Christian University.Co-ingestion of feed and RT by goldfish of three different body sizesWe examined whether RT are ingested by goldfish with feed and whether the body size of fish affects the ingestion of RT using three different body sizes of fish, large, medium, and small. First, we conducted an experiment using large body size fish (N = 24; BW, 91.9 ± 21.6 g, mean and SD). Three goldfish of large body size were transferred from the stock tank to the experimental 60-L glass aquarium and kept for three days for acclimation of the environment and sinking fish feed. Fish were fed 3.0 g of large-size feed (Japan Pet Design Co. Ltd.) once a day. On the fourth day, fish were fed a mixture of RT collected from the field (ICU, 300 mg) and large feed (3.0 g). Control fish were fed only fish feed. At 90 min after feeding, the fish were transferred to a pail containing 0.05% 2-phenoxyethanol solution and deeply anesthetized. After body weight measurement, fish were dissected. We observed the intestine to determine whether RT was ingested. When RT was observed in the intestine, we collected the tips and counted the number of tips in each fish. The experimental tests were repeated eight times, and the data were combined.Second, we conducted an experiment using medium body size fish (N = 24; BW, 30.4 ± 12.4 g). Three goldfish of medium body size were transferred from the stock tank to the experimental 60-L glass aquarium and kept for three days for acclimation. Fish were fed 1.0 g of medium-size feed (Kyorin) once a day. On the fourth day, fish were fed a mixture of RT (ICU, 300 mg) and medium feed (1.0 g). Control fish were fed only fish feed. At 90 min after feeding, fish were anesthetized and dissected, and then the intestine was observed as described above. The experimental tests were repeated eight times, and the data were combined.Third, we conducted an experiment using small body size fish (N = 40; BW, 4.4 ± 1.5 g). Four goldfish of small body size were transferred from the stock tank to the experimental 60-L glass aquarium and kept for three days for acclimation. Fish were fed 0.5 g of small-size feed (Kyorin) once a day. On the fourth day, fish were fed a mixture of RT (ICU, 300 mg) and small feed (0.5 g). Because of the small size of fish, RT of small size particles (212–500 µm) were collected with sieves and used for the tests. Control fish were fed only fish feed. At 90 min after feeding, fish were anesthetized and dissected, and then the intestine was observed as described above. The experimental tests were repeated ten times, and the data were combined.In the first three experiments, all three control groups showed no ingestion of RT. From the results of the three experiments, it was clear that our experimental system was not contaminated with RT. Therefore, we omitted making control groups for further experiments to decrease the number of fish sacrificed from the standpoint of fish welfare.Fourth, we examined whether RT collected from TGU was ingested by goldfish. We conducted an experiment using large body size fish (N = 12; BW, 140.3 ± 27.0 g). Three goldfish of large body size were transferred from the stock tank to the experimental 60-L glass aquarium and kept for three days for acclimation. Fish were fed 3.0 g of large-size feed once a day. On the fourth day, fish were fed a mixture of RT (TGU, 300 mg) and large feed (3.0 g). At 90 min after feeding, fish were anesthetized and dissected, and then the intestine was observed as described above. The experimental tests were repeated four times, and the data were combined.We conducted an additional experiment with a similar design to those of the four experiments to take photographs of the fish and RT using fish of small-medium body size (N = 9; BW, 12.8 ± 2.7 g). Three fish were transferred from the stock tank to the experimental 60-L glass aquarium and kept for two days for acclimation. Fish were given 0.5 g of medium-size feed once a day. On the third day, fish were given a mixture of RT of ICU (30 pieces; size 0.5–1.0 mm) and medium feed (0.5 g). At 60 min after feeding, fish were anesthetized and dissected, and photographs of RT in the intestine were taken. The experimental tests were repeated three times, and the data were combined.Active ingestion of RT by goldfishWe examined whether goldfish actively ingest RT when RT are given without fish feed using large body size fish (N = 9; BW, 122.4 ± 20.8 g). Three fish were transferred from the stock tank to the experimental 60-L glass aquarium and kept for three days for acclimation. Fish were fed 3.0 g of large-size feed once a day. On the fourth day, fish were given 300 mg of RT (ICU) on the bottom of the aquarium. At 90 min after the placement of RT, fish were anesthetized and dissected, and then the intestine was observed as described above. The experimental tests were repeated three times, and the data were combined.Retention and elimination of ingested RT in the intestine of goldfishWe examined how long RT was retained in the intestine using large body size goldfish (N = 9; BW, 101.6 ± 11.4 g). Three goldfish were transferred from the stock tank to the experimental 60-L glass aquarium and kept for three days for acclimation. Fish were fed 3.0 g of large-size feed once a day. On Day 4, fish were given 1.0 g of RT (ICU). At 90 min after the placement of RT, each fish was individually transferred to three experimental 60-L glass aquaria. Then, each fish was fed 1.0 g of the feed. At 24 and 48 h (Day 5 and Day 6) after the transfer, we collected feces from fish and some water from the bottom of the aquaria. We observed whether RT was eliminated from the fish into the aquaria. When RT was observed in the feces and the bottom of the aquarium, we collected the RT and counted the number of RT. On Day 5, after the RT observation, each fish was fed 1.0 g of the feed. On Day 6, after RT observation in feces and water, the fish were anesthetized and dissected. We observed whether the intestine retained RT. The experimental tests were repeated three times, and the data were combined.Ingestion of RT by wild crucian carpWe examined whether wild Japanese crucian carp ingest RT. The experiment was conducted using juvenile crucian carp (N = 16, BW, 2.8 ± 0.9 g). Sixteen fish were transferred from the stock tank to three experimental 60-L glass aquaria (5 or 6 fish per aquarium) and kept for six days for acclimation. Fish were fed with 0.2 g of small-size feed once a day. On the seventh day, fish were fed a mixture of RT (ICU, 30 mg) and the small feed (0.2 g) or RT alone (30 mg). Because of the small size of fish, RT of small size particles (212–500 µm) were collected with sieves and used for the test. Control fish were fed only fish feed (0.2 g). At 6 h after feeding, fish were anesthetized and dissected, and then the intestine was observed as described above. More