Ecology

1.Luyster, F. S., Strollo, P. J., Zee, P. C. & Walsh, J. K. Sleep: A health imperative. Sleep 35, 727–734 (2012).PubMed 
PubMed Central 

Google Scholar 
2.Siegel, J. M. Clues to the functions of mammalian sleep. Nature 437, 1264–1271 (2005).CAS 
PubMed 
ADS 

Google Scholar 
3.Cirelli, C. & Tononi, G. Is sleep essential?. PLoS Biol. 6, 1605–1611 (2008).CAS 

Google Scholar 
4.Lesku, J. A., Roth, T. C. II., Amlaner, C. J. & Lima, S. L. A phylogenetic analysis of sleep architecture in mammals: The integration of anatomy, physiology, and ecology. Am. Nat. 168, 441–453 (2006).PubMed 

Google Scholar 
5.Banks, S. & Dinges, D. F. Behavioral and physiological consequences of sleep restriction. J. Clin. Sleep Med. 3, 519–528 (2007).PubMed 
PubMed Central 

Google Scholar 
6.Tougeron, K. & Abram, P. K. An ecological perspective on sleep disruption. Am. Nat. 190, E55–E66 (2017).PubMed 

Google Scholar 
7.Vyazovskiy, V. V. & Tobler, I. The temporal structure of behaviour and sleep homeostasis. PLoS ONE 7, 1–11 (2012).
Google Scholar 
8.Tobler, I. Is sleep fundamentally different between mammalian species?. Behav. Brain Res. 69, 35–41 (1995).CAS 
PubMed 

Google Scholar 
9.Vyazovskiy, V. Sleep, recovery, and metaregulation: Explaining the benefits of sleep. Nat. Sci. Sleep 7, 171 (2015).PubMed 
PubMed Central 

Google Scholar 
10.Vyazovskiy, V. V. & Delogu, A. NREM and REM sleep. Neuroscience 20, 203–219 (2014).
Google Scholar 
11.Orzeł-Gryglewska, J. Consequences of sleep deprivation. Int. J. Occup. Med. Environ. Health 23, 95–114 (2010).PubMed 

Google Scholar 
12.Meerlo, P., Sgoifo, A. & Suchecki, D. Restricted and disrupted sleep: Effects on autonomic function, neuroendocrine stress systems and stress responsivity. Sleep Med. Rev. 12, 197–210 (2008).PubMed 

Google Scholar 
13.Touitou, Y., Reinberg, A. & Touitou, D. Association between light at night, melatonin secretion, sleep deprivation, and the internal clock: Health impacts and mechanisms of circadian disruption. Life Sci. 173, 94–106 (2017).CAS 
PubMed 

Google Scholar 
14.Pires, G. N., Bezerra, A. G., Tufik, S. & Andersen, M. L. Effects of acute sleep deprivation on state anxiety levels: A systematic review and meta-analysis. Sleep Med. 68, 575–589 (2016).
Google Scholar 
15.Hudson, A. N., Van Dongen, H. P. A. A. & Honn, K. A. Sleep deprivation, vigilant attention, and brain function: A review. Neuropsychopharmacology https://doi.org/10.1038/s41386-019-0432-6 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
16.Tobler, I. & Sigg, H. Long-term motor activity recording of dogs and the effect of sleep deprivation. Experientia 42, 987–991 (1986).CAS 
PubMed 

Google Scholar 
17.Hänninen, L. Sleep and Rest in Calves Relationship To Welfare, Housing and Hormonal Activity (University of Helsinki, Helsinki, 2007).
Google Scholar 
18.Hsieh, W.-H.H. et al. Simulated shift work in rats perturbs multiscale regulation of locomotor activity. J. R. Soc. Interface 11, 20140318 (2014).PubMed 
PubMed Central 

Google Scholar 
19.Storch, C., Höhne, A., Holsboer, F. & Ohl, F. Activity patterns as a correlate for sleep–wake behaviour in mice. J. Neurosci. Methods 133, 173–179 (2004).PubMed 

Google Scholar 
20.Hicks, R. A., Moore, J. D., Hayes, C., Phillips, N. & Hawkins, J. REM sleep deprivation increases aggressiveness in male rats. Physiol. Behav. 22, 1097–1100 (1979).CAS 
PubMed 

Google Scholar 
21.Pires, G. N., Tufik, S. & Andersen, M. L. Grooming analysis algorithm: Use in the relationship between sleep deprivation and anxiety-like behavior. Prog. Neuro-Psychopharmacol. Biol. Psychiatry 41, 6–10 (2013).
Google Scholar 
22.Vorster, A. P. & Born, J. Sleep and memory in mammals, birds and invertebrates. Neurosci. Biobehav. Rev. 50, 103–119 (2015).PubMed 

Google Scholar 
23.Greives, T. J. et al. Costs of sleeping in: Circadian rhythms influence cuckoldry risk in a songbird. Funct. Ecol. 29, 1300–1307 (2015).
Google Scholar 
24.Fuchs, T., Haney, A., Jechura, T. J., Moore, F. R. & Bingman, V. P. Daytime naps in night-migrating birds: Behavioural adaptation to seasonal sleep deprivation in the Swainson’s thrush, Catharus ustulatus. Anim. Behav. 72, 951–958 (2006).
Google Scholar 
25.Klein, B. A., Klein, A., Wray, M. K., Mueller, U. G. & Seeley, T. D. Sleep deprivation impairs precision of waggle dance signaling in honey bees. Proc. Natl. Acad. Sci. U. S. A. 107, 22705–22709 (2010).CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
26.Owczarczak-Garstecka, S. C. & Burman, O. H. P. P. Can sleep and resting behaviours be used as indicators of welfare in shelter dogs (Canis lupus familiaris)?. PLoS ONE 11, e0163620 (2016).PubMed 
PubMed Central 

Google Scholar 
27.Langford, F. M. & Cockram, M. S. Is sleep in animals affected by prior waking experiences?. Anim. Welf. 19, 215–222 (2010).CAS 

Google Scholar 
28.Toth, L. A. & Bhargava, P. Animal models of sleep disorders. Comp. Med. 63, 91–104 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
29.Zanghi, B. M., Kerr, W., DeRivera, C., Araujo, J. & Milgram, B. Sleep and biorhythms as a function of age in the dog. J. Vet. Behav. 5, 159 (2010).
Google Scholar 
30.Zanghi, B. M. Circadian Biorhythms of Sleep/Wake Dogs, Activity/Rest Cycles In Adult and Aged dogs. in Nestle Purina Companion Aninmal Nutrition Summit vol. 31 114–117 (Purina Institute, 2010).31.Zanghi, B. M. et al. Characterizing behavioral sleep using actigraphy in adult dogs of various ages fed once or twice daily. J. Vet. Behav. 8, 195–203 (2013).
Google Scholar 
32.Bódizs, R., Kis, A., Gácsi, M. & Topál, J. Sleep in the dog: Comparative, behavioral and translational relevance. Curr. Opin. Behav. Sci. 33, 25–33 (2020).
Google Scholar 
33.Kis, A. et al. Development of a non-invasive polysomnography technique for dogs (Canis familiaris). Physiol. Behav. 130, 149–156 (2014).CAS 
PubMed 

Google Scholar 
34.Takeuchi, T. & Harada, E. Age-related changes in sleep-wake rhythm in dog. Behav. Brain Res. 136, 193–199 (2002).PubMed 

Google Scholar 
35.Adams, G. J. & Johnson, K. G. Sleep-wake cycles and other night-time behaviours of the domestic dog Canis familiaris. Appl. Anim. Behav. Sci. 36, 233–248 (1993).
Google Scholar 
36.Adams, G. J. & Johnson, K. G. Behavioural responses to barking and other auditory stimuli during night-time sleeping and waking in the domestic dog (Canis familiaris). Appl. Anim. Behav. Sci. 39, 151–162 (1994).
Google Scholar 
37.Bunford, N. et al. Differences in pre-sleep activity and sleep location are associated with variability in daytime/nighttime sleep electrophysiology in the domestic dog. Sci. Rep. 8, 7109 (2018).PubMed 
PubMed Central 
ADS 

Google Scholar 
38.Iotchev, I. B. et al. Age-related differences and sexual dimorphism in canine sleep spindles. Sci. Rep. 9, 10092 (2019).PubMed 
PubMed Central 
ADS 

Google Scholar 
39.Kis, A. et al. Sleep macrostructure is modulated by positive and negative social experience in adult pet dogs. Proc. R. Soc. B Biol. Sci. 284, 20171883 (2017).
Google Scholar 
40.Mong, J. A. & Cusmano, D. M. Sex differences in sleep: Impact of biological sex and sex steroids. Philos. Trans. R. Soc. B Biol. Sci. 371, 20150110 (2016).
Google Scholar 
41.Andersen, M. L. et al. Effects of sleep loss on sleep architecture in Wistar rats: Gender-specific rebound sleep. Prog. Neuro-Psychopharmacol. Biol. Psychiatry 32, 975–983 (2008).CAS 

Google Scholar 
42.McKillop, L. E. et al. Effects of aging on cortical neural dynamics and local sleep homeostasis in mice. J. Neurosci. 38, 3911–3928 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
43.Nakamura, T. J., Takasu, N. N. & Nakamura, W. The suprachiasmatic nucleus: Age-related decline in biological rhythms. J. Physiol. Sci. 66, 367–374 (2016).CAS 
PubMed 

Google Scholar 
44.Abou-Ismail, U. A., Burman, O. H. P., Nicol, C. J. & Mendl, M. Can sleep behaviour be used as an indicator of stress in group-housed rats (Rattus norvegicus)?. Anim. Welf. 16, 185–188 (2007).CAS 

Google Scholar 
45.Sadeh, A., Keinan, G. & Daon, K. Effects of stress on sleep: The moderating role of coping style. Heal. Psychol. 23, 542–545 (2004).
Google Scholar 
46.Okun, M. L. Biological consequences of disturbed sleep: Important mediators of health?. Jpn. Psychol. Res. 53, 163–176 (2011).PubMed 
PubMed Central 

Google Scholar 
47.Van Reeth, O. et al. Interactions between stress and sleep: From basic research to clinical situations. Sleep Med. Rev. 4, 201–219 (2000).
Google Scholar 
48.Novati, A. et al. Chronically restricted sleep leads to depression-like changes in neurotransmitter receptor sensitivity and neuroendocrine stress reactivity in rats. Sleep 31, 1579–1585 (2008).PubMed 
PubMed Central 

Google Scholar 
49.Taylor, K. D. & Mills, D. S. The effect of the kennel environment on canine welfare: A critical review of experimental studies. Anim. Welf. 16, 435–447 (2007).CAS 

Google Scholar 
50.Sales, G., Hubrecht, R., Peyvandi, A., Milligan, S. & Shield, B. Noise in dog kennelling: Is barking a welfare problem for dogs?. Appl. Anim. Behav. Sci. 52, 321–329 (1997).
Google Scholar 
51.Hewison, L. F., Wright, H. F., Zulch, H. E. & Ellis, S. L. H. Short term consequences of preventing visitor access to kennels on noise and the behaviour and physiology of dogs housed in a rescue shelter. Physiol. Behav. 133, 1–7 (2014).CAS 
PubMed 

Google Scholar 
52.Schwartz, J. R. L. & Roth, T. Neurophysiology of sleep and wakefulness: Basic science and clinical implications. Curr. Neuropharmacol. 6, 367–378 (2008).CAS 
PubMed 
PubMed Central 

Google Scholar 
53.Reicher, V. et al. Repeated afternoon sleep recordings indicate first-night-effect-like adaptation process in family dogs. J. Sleep Res. 29, 1–10 (2020).
Google Scholar 
54.Piccione, G. et al. Comparison of daily distribution of rest/activity in companion cats and dogs. Biol. Rhythm Res. 45, 615–623 (2014).
Google Scholar 
55.Kredlow, M. A., Capozzoli, M. C., Hearon, B. A., Calkins, A. W. & Otto, M. W. The effects of physical activity on sleep: A meta-analytic review. J. Behav. Med. 38, 427–449 (2015).PubMed 

Google Scholar 
56.Barber, N. Play and energy regulation in mammals. Q. Rev. Biol. 66, 129–147 (1991).CAS 
PubMed 

Google Scholar 
57.Ahloy-Dallaire, J., Espinosa, J. & Mason, G. Play and optimal welfare: Does play indicate the presence of positive affective states?. Behav. Processes 156, 3–15 (2018).PubMed 

Google Scholar 
58.Penev, P. D. Sleep deprivation and energy metabolism: To sleep, perchance to eat?. Curr. Opin. Endocrinol. Diabetes. Obes. 14, 374–381 (2007).PubMed 

Google Scholar 
59.Mavanji, V., Billington, C. J., Kotz, C. & Teske, J. A. Sleep and obesity: A focus on animal models. Neurosci. Biobehav. Rev. 36, 1015–1029 (2012).PubMed 
PubMed Central 

Google Scholar 
60.Zanghi, B. M., Kerr, W., de Rivera, C., Araujo, J. A. & Milgram, N. W. Effect of age and feeding schedule on diurnal rest/activity rhythms in dogs. J. Vet. Behav. 7, 339–347 (2012).
Google Scholar 
61.Knutson, K. L., Spiegel, K., Penev, P. & Van Cauter, E. The metabolic consequences of sleep deprivation. Sleep Med. Rev. 11, 163–178 (2007).PubMed 
PubMed Central 

Google Scholar 
62.Kiddie, J. & Collins, L. Identifying environmental and management factors that may be associated with the quality of life of kennelled dogs (Canis familiaris). Appl. Anim. Behav. Sci. 167, 43–55 (2015).
Google Scholar 
63.Taylor, K. D. & Mills, D. S. The effect of the kennel environment on canine welfare: A critical review of experimental studies (2007).
64.CONCEA. Resolução normativa no 12. Diretriz brasileira para o cuidado e a utilização de animais para fins científicos e didáticos. Diário Oficial da União No 186 (2013).65.Bateson, M. & Martin, P. Measuring Behaviour (Cambridge University Press, 2021) https://doi.org/10.1017/CBO9780511810893.Book 

Google Scholar 
66.Broom, D. M. & Fraser, A. F. Domestic Animal Behaviour and Welfare (CABI, 2015).
Google Scholar 
67.Luescher, U. A., McKeown, D. B. & Halip, J. Stereotypic or obsessive-compulsive disorders in dogs and cats. Vet. Clin. North Am. Small Anim. Pract. 21, 401–413 (1991).CAS 
PubMed 

Google Scholar 
68.Friard, O. & Gamba, M. BORIS: A free, versatile open-source event-logging software for video/audio coding and live observations. Methods Ecol. Evol. 7, 1325–1330 (2016).
Google Scholar 
69.Crawley, M. J. The R Book (Wiley, 2007). https://doi.org/10.1002/9780470515075.Book 
MATH 

Google Scholar 
70.Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting Linear Mixed-Effects Models Using lme4. J. Stat. Softw. 67, 1–48 (2015).
Google Scholar 
71.Becker, R. A., Chambers, J. M. & Wilks, A. R. The New S Language (Chapman and Hall/CRC, 1988).MATH 

Google Scholar  More

Study areaAll simulations were run at a 100 × 100 m resolution for the entire European Alps, which cover ~200,000 km². Elevations reach 4,810 m above sea level at the highest peak (Mont Blanc, elevational data were obtained from ref. 44). Mean annual temperature ranges from approximately −13 up to 16 °C and annual precipitation sums reach up to ~3,600 mm (climatic conditions were obtained from WorldClim45).Species dataTrue presences/absences were derived from complete species lists of 14,040 localized plots covering subalpine and alpine non-forest vegetation of the Alps, compiled from published46 and unpublished data sources. For more information see the supplementary information in ref. 21.Data on demographic rates as well as dispersal parameters were taken from ref. 21, Supplementary Table 1. Detailed values are provided in Supplementary Table 1.Environmental variablesCurrent climate dataMaps of current climatic conditions were generated on the basis of mean, minimum and maximum monthly temperature obtained from WorldClim45 and monthly precipitation sums derived from ref. 47 at a resolutions of 0.5 arcmin and 5 km, respectively. Temperature and precipitation data were downscaled to 100 m as described in ref. 21 and using ordinary kriging with elevation as covariable. As the reference periods of these two datasets do not match (temperature values represent averages for 1950–2000, while precipitation covers 1970–2005) temperature values were adapted using the E-OBS climate grids available online (www.ecad.eu/download/ensembles/download.php). We used these spatially refined temperature and precipitation grids to derive maps of mean annual temperature and mean annual precipitation sum. We chose only two climatic variables to keep models simple and, therefore, interpretation of results more straightforward. In fact, the climatic drivers of population performance and distribution can be more complex48 and vary among species, life history stages and vital rates49. However, as correlations between different descriptors of temperature (such as coldest month or warmest month temperature, Pearson correlation of 0.84) as well as between precipitation variables are high in the European Alps, we chose mean annual temperature and mean annual precipitation sum as they give the most basic description of how climatic conditions change over geographical and elevational gradients.Future climate dataMonthly time series of mean temperature as well as precipitation sums predicted for the twenty-first century were extracted from the Cordex data portal (http://cordex.org). We chose two IPCC5 scenarios from the RCP families representing mild (RCP 2.6) and severe (RCP 8.5) climate change to consider the uncertainty in the future climate predictions. Both scenarios were generated by Météo-France/Centre National de Recherches Météorologiques using the CNRM-ALADIN53 model, fed by output from the global circulation model CNRM-CM5 (ref. 50). The RCP 2.6 scenario assumes that radiative forcing reaches nearly 3 W m−2 (equal to 490 ppm CO2 equivalent) mid-century and will decrease to 2.6 W m−2 by 2100. In the RCP 8.5 scenario, radiative forcing continues to rise throughout the twenty-first century and reaches >8.5 W m−2 (equal to 1,370 ppm CO2 equivalent) in 210024.These time series were statistically downscaled (delta method) by (1) calculating differences (deltas) between monthly temperature and precipitation values hindcasted for the current climatic conditions (mean 1970–2005) and forecasted future values at the original spatial resolution of 11′; (2) spatially interpolating these differences to a resolution of 100 × 100 m using cubic splines and (3) adding them to the downscaled current climate data separately for each climatic variable21,36. Subsequently, we calculated running means (averaged over 35 years) for every tenth year (2030, 2040 through to 2080) for each climatic variable and finally derived, on the basis of the monthly data, mean annual temperature and mean annual precipitation sums for these decadal time steps. The application of 35-yr running means ensures that we use the same time interval for parameterization and prediction. Furthermore, for long-lived species such as alpine plants, running means over long time intervals appear most appropriate to characterize climatic niches33.Soil dataIn addition to the climatic data we used a map of the percentage of calcareous substrate within a cell (5′ longitude × 3′ latitude dissolved to 100 m resolution; further referred to as soil) as described in the supplementary information of ref. 21.Environmental suitability modellingWe parameterized logistic regression models (LRMs) with a logit link function using species presence/absence data as response and the three environmental (two bioclimatic and one soil) variables as predictors. All predictor variables entered the model as second-order polynomials in agreement with the standard unimodal niche concept.From the models, we also derived a threshold value to use for translating predicted probability of occurrence (as a measure of site suitability) into predicted presence or absence of each species at a site (called occurrence threshold, OT, henceforth). The threshold was defined such that it optimized the true skills statistic (TSS), a measure of predictive accuracy derived from comparing observed and predicted presence–absence maps51.Genetic model and niche partitioningSpecies-specific suitability curves for the annual mean temperature gradient derived from the LRMs were partitioned into suitability curves of ecologically different haplotypes of a species as defined by allelic variation in up to three diploid loci (Extended Data Fig. 6). The number of alleles was varied (n = 1, 2, 3, 5 and 10 alleles) as was the proportion of the entire species’ (temperature) niche covered by each haplotype. Models with more than one locus were run with diallelic loci, as otherwise computational efforts would have increased excessively (for each haplotype the number of seeds, juveniles and adults has to be stored and all seeds have to be distributed separately). Each combination of haplotype number and allelic niche size used in a particular simulation is further referred to as setting. Species-specific suitability curves along the other two dimensions (precipitation and soil) remained unpartitioned to ease interpretation of results. The implications of relaxing this assumption by modelling niche partitioning along different environmental gradients are hard to predict. Outcomes would probably depend on the direction and strength of individual specialization along these gradients, whether they are positively or negatively correlated, as well as on how both temperature and precipitation patterns will be affected by climate change. As a consequence, the patterns we found could be re-enforced, for example when the replacement of cold-adapted haplotypes within the species elevational range is further delayed, for example, because those haplotypes adapted to warmer conditions can cope less well with higher precipitation at higher elevations. Vice versa, maladaptation to the warming temperatures might be attenuated, for example, if temperature increase is paralleled by precipitation decrease and emerging drought stress. If, in this case, haplotypes from lower elevations can better cope with both higher temperatures and less water availability than those of median elevations, they may replace the latter faster at these median elevations and hence shorten the phase of maladaptation.Allelic effects were assumed to define the temperature optimum additively. Hence, the heterozygotes’ optimum is always exactly between the optima of the two corresponding homozygotes, corresponding to a codominant genetic model. Finally, all haplotypes corresponding to one setting were assumed to have constant (temperature) niche size, defined as a proportion (ω = 50%, 75% and 100%, for one haplotype only 100%) of the entire species’ (temperature) niche. The temperature niche was computed as the difference between the upper and lower temperature values at which the LRM-derived suitability curve predicted a suitability equal to OT (with precipitation and soil set to the respective optima of the species, also derived from the LRMs). To derive the same geographic distribution under current climatic conditions for each setting, the union of the niches of all haplotypes of one set has to approximate the niche of the single-species model (Extended Data Fig. 6). Note, however, that, the aspired equality of niches is impossible, as the niches resulting from a logistic regression with quadratic terms are always elliptic in shape. Therefore, the optima of the two extreme homozygotes (that is, those carrying haplotypes adapted to the coldest or warmest margins of the entire species’ niche) are fixed at: niche limits ± 1/2 × ω × niche size and all other optima are distributed between them at equal distances (Extended Data Fig. 6 gives a schematic illustration). As a consequence, the performance of the extreme haplotypes, both coldest and warmest, is modelled to be somewhat higher towards the cold and warm margins of the temperature niche, respectively, compared with the performance modelled for the species without intraspecific niche partitioning (compare the overlap of the black with the red and blue curves in Extended Data Fig. 6a). However, as haplotype number did not affect modelled range loss (compare Table 1), this marginal mismatch does not apparently impact our results. Homozygotes were ordered from the cold-adapted A1A1 up to the warm-adapted AnAn.Finally, the suitability curves of the genotypes were assumed to have the same value at their optimum as the species-specific suitability curve at this point (Extended Data Fig. 6).Artificial landscapesArtificial landscapes were defined as a raster of 50 × 112 cells (of 100 × 100 m). These rasters were homogeneous with respect to precipitation and soil carbon conditions which were set to the values optimal for each species according to the LRMs. With respect to temperature, by contrast, we implemented a gradient across the raster running from the minimum –9.1 °C to the maximum +3.8 °C temperature for which the LRM predicts values >OT across all six species. Buffering by 1 °C at both limits was done to avoid truncating simulation results. Further 4 °C at the lower limit were added to consider simulated temperature increase (below). The final temperature range represented by the raster ran from −14.1 to +4.8 °C. Temperature increased linearly along this axis by an increment of 0.171 °C per cell, derived from a 20° slope and a temperature decrease of 0.5 °C per 100 m of elevational change. Along the 50-cell breadth of the landscape, temperature values were kept constant. On the basis of these grids, we implemented a moderate and a severe climate change scenario, characterized by temperature increases of 2 and 4 °C, respectively, over 80 yr. Therefore, temperature of each raster cell increased annually by 0.025 and 0.05 °C, respectively.Simulations of spatiotemporal range dynamicsCATS21 is a spatially and temporarily explicit model operating on a two-dimensional grid (of 100 m mesh size in this case). CATS combines simulations of local species’ demography with species’ distribution models by scaling demographic rates relative to occurrence probabilities (suitabilities) predicted for the respective site or grid cell by the latter. Dispersal among grid cells is modelled as a combination of wind, exozoochoric and endozoochoric (that is, animal dispersal via attachment to the fur or ingestion and defecation, respectively) dispersal of plant seeds. Time proceeds in annual steps. The annually changing occurrence probabilities at each site affect the demography of local populations and hence, eventually, the number of seeds that are produced in each grid cell in the respective year. As a consequence, local populations grow or decrease, become extinct or establish anew and hence the species’ distribution across the whole study area changes as a function of the changing climate.Demographic modellingClimate-dependence of local demography was modelled by linking demographic rates (seed persistence, germination, survival, flowering frequency, seed yield and clonal reproduction) and carrying capacity to occurrence probabilities predicted by LRMs by means of sigmoidal functions. Furthermore, all rates were fixed at OT at a value ensuring stable population sizes; for more information see refs. 21,33. Demographic rates were confined by zero and a species-specific maximum value (Supplementary Table 1), which was assumed to be the same for all genotypes of a species. As a corollary, the demographic rates of all genotypes of a species are the same under optimal environmental conditions but their actual values for a particular site in a particular year differ due to different temperature optima of genotypes. In addition, germination, survival and clonal reproduction were modelled as density-dependent processes to account for intraspecific competition between genotypes. In our application, for all density-dependent functions, the species abundance is defined as the sum of all adult individuals in a given cell, regardless of their genotypes. Density dependence is commonly achieved by multiplying rates with (C – N)/C, where N is the species abundance and C is the (site- and genotype-specific) carrying capacity. This correction for density dependence causes the functions to drop to zero when N approximates C. To avoid the subsequent collapse of population sizes, we defined density-dependent rates as (C – N)/C × rate() + N/C × rate(OT), which ensures stable population sizes at densely populated sites occupied by only one genotype. To account for uncertainty in parameters of demographic rates, we assigned each species two value sets representing the upper and lower end of a plausible range of values on the basis of information derived from databases and own measurements (Supplementary Table 1).The simulations allowed for cross-pollination between genotypes. We used the relative amount of flowers (genotype-specific flowering frequency as defined by the sigmoid curve for the given suitability in the given raster cell for the given year × number of adults of that genotype in the population of that cell) to derive an estimate of the haplotype frequencies in the total pollen produced by the population within a grid cell. For the multiallelic case we allowed for recombination between loci with a recombination rate of 0.1%. These frequencies were set equal to the probability that particular haplotypes are transmitted to each year’s seed yield by pollination. Spatial pollen dispersal was accounted for in the following way: in each cell, 5% of the pollen involved in producing the annual seed yield, was assumed to stem from outside the respective raster cell. The proportions of different haplotypes in this 5% were derived from the overall pollen frequencies in all cells within a 700 m radius around the target cell (following estimates in ref. 52). Subsequently, produced seeds of each genotype were divided into resulting genotypes regarding the adult’s haplotype composition and the haplotype frequencies in the cells’ entire pollen load.Dispersal modellingFor wind dispersal of plant species we parameterized the analytical WALD kernel53 on the basis of measured seed traits and wind speed data from a meteorological station in the Central Alps of Austria. Exozoochorous and endozoochorous plant kernels were parameterized on the basis of correlated random walk simulations for the most frequent mammal dispersal vector in the study area, the chamois (Rupicapra rupicapra L.). For more details, see ref. 33. To account for uncertainties in species-specific dispersal rates, the proportion of seeds dispersed by the more far-reaching zoochorous kernels was assumed either as high (1–5%) or low (0.1–0.5%), setting upper and lower boundaries of a credible range of the dispersal ability of species.Simulation set up and simulation initializationTo test for the effects of climate change on genetic diversity in 2080, we ran CATS over the period 2000 to 2080 for each of the six study species across the entire Alps under a full factorial combination of (1) three niche sizes (50%, 75% and 100%); (2) six numbers of haplotypes (equal to two, three, five and ten alleles for one locus and four and eight for the diallelic two- and three-locus models, respectively); (3) three climate scenarios (current climate, RCP 2.6 and RCP 8.5); and (4) two sets of demographic and dispersal parameters. As a ‘control’ we also ran simulations for all climate scenarios and the two demographic and dispersal parameter sets for a setting with one genotype filling the whole niche of the species. To account for stochastic elements in CATS four replications were run for each combination of ‘treatments’.For simulations in artificial landscapes we used, instead of RCP 2.6 and RCP 8.5, ‘artificial’ climatic scenarios with an assumed warming of 2 and 4 °C, respectively, and no change in precipitation.All simulation runs were started with homozygotic individuals only. As initial distribution, for each simulation run each cell predicted to be environmentally suitable to the species (that is, occurrence probability of species >OT)—and within the real distribution range of the species28 (not relevant for simulations in artificial landscapes, of course)—was assumed to be occupied by an equal number of adults of each (homozygotic) genotype, with the total sum equal to the carrying capacity of the site. To accommodate this arbitrary within-cell genetic mixture of homozygotes (and missing heterozygotes) to actual local conditions we started simulations of range dynamics with a burn-in phase of 200 yr, run under constant current climatic conditions. To have a smooth transition from the burn-in phase under current climate (corresponding to the climate of the years 1970–2005; see current climate data) to future climate projections (starting with 2030) and to derive annual climate series, climate data were linearly interpolated between these two time intervals.Statistical analysisWe evaluated the contribution of climate scenario, haplotype number and haplotype niche size to overall species’ range change as well as the change in the frequency of the warm-adapted haplotype by means of linear models. In these models, log(range size in 2080/range size in 2000) and log(percentage of warm-adapted haplotype in 2080/percentage of warm-adapted haplotype in 2000), averaged over the four replicates and the two demographic and dispersal parameter sets, were the response variables. For the analysis of the change of the warm-adapted haplotype simulation settings with 100% niche size were ignored, as in this setting all haplotypes have the same temperature optimum (that is, neutral genetic variation). Approximate normality of residuals was confirmed by visual inspection.As indicator of the ‘topographic opportunity’ remaining to the species in the real world we calculated the area colonizable at elevations higher than those already occupied at the end of the simulation period. We therefore drew a buffer of 1 km around each cell predicted to be occupied in 2080 and then summed the area within these buffers at a higher elevation than the focal, occupied cell. Overlapping buffer areas were only counted once. This calculation was done for simulations conducted with one full-niche genotype per species only.Sensitivity analysisWe interpret the simulated relative decrease of warm-adapted haplotypes mainly as an effect of (1) the unrestricted expansion of cold-adapted haplotypes at the leading edge and (2) the resistance of the locally predominating haplotype that becomes increasingly maladapted with progressive climate warming, to recruitment by better-adapted haplotypes from below that are either rare or not present at all initially. However, the results achieved, and our conclusions, may be sensitive to assumptions about particular parameter values. Parameters that control the longevity of adult plants, and indirectly the rate of recruitment of new individuals, as well as those controlling gene flow via pollen (instead of seeds) may be particularly influential in this respect. We additionally ran simulations on artificial landscapes under alternative values of these parameters. In particular, we set the maximum age of plants to 10 yr instead of 100 yr and raised the proportion of locally produced pollen assumed to be transported up to 700 m to 10%. Both of these values are thus probably set to rather extreme levels: a maximum age of 10 yr is much shorter than the 30–50 yr assumed to be the standard age of (non-clonal) alpine plants31; and a cross-pollination rate between cells of 10% is high given that among the most important alpine pollinators only bumblebees are assumed to transport pollen >100 m regularly54,55. We add that we ran these additional simulations only in combination with the demographic species parameters set to high values because a short life span combined with low-level demographic parameters did not allow for stable populations of some species, even under current climatic conditions.For individual species, adapting plant age and cross-pollination rate between cells (Extended Data Fig. 7), did change the magnitude of loss of the warm-adapted haplotype. Nevertheless, for all of them the warm-adapted haplotype still became rarer when climate warmed and this effect increased with the level of warming. We are confident that our conclusions are qualitatively insensitive to variation of these parameters within a realistic range.Finally, in the simulations where we assumed a multilocus structure of the temperature niche, the recombination rate may also affect simulation results because it determines the rate by which new haplotypes can emerge. We also tested sensitivity of our simulations to doubling the recombination rate to 0.2%. Again, we found that a higher recombination rate had little qualitative effect on the results. Quantitatively, it resulted in a slightly pronounced relative decrease of the warmth-adapted haplotype in most species (Extended Data Fig. 8).Reporting SummaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Concentrations and sources of organic compoundsMid- and long-chain n-alkanesThroughout most sections, the biomarker composition of the Hala Hu record is dominated by long-chain n-alkanes (Fig. 2a). In most samples, concentrations were highest for nC31 (ca 5–35 µg/g d.w.) and decreasing with chain-length (Supplementary Fig. S3–1). Generally, the concentrations were low in the glacial period, and increased between ca. 10 and 8 ka cal BP, while gradually decreasing after ca. 5 ka cal BP. These alkanes can be attributed to vascular higher plants from the lake catchment27, which is characterized by alpine meadows. Mid-chain nC23 and nC25-alkanes are frequently attributed to aquatic macrophytes28. However, n-alkane patterns of aquatic and terrestrial plants overlap. Moreover, we consider the contribution from macrophyte to the sedimentary n-alkane pool at the coring location to be minor, because of the specific n-alkane pattern of the samples, the overall low concentrations of mid-chain n-alkanes, and the deep water depth of the core site. An exception is a section in the core dated to ca. 15–14 ka cal BP, that has concentrations from ca. 20 up to 50 µg/g d.w. for individual mid-chain n-alkanes (Fig. 2a). This can be explained by enhanced contribution from submerged aquatic plants29 when the lake level was much lower than present9,23. Alternatively, a contribution from other mid-chain producers, such as Sphagnum species30,31, is possible during phases of lake regression and the potential formation of peatlands around the lake shore.Fig. 2: Summary of measured proxy parameters in cores H7 and H11 (overlap 7.4−9.1 kyrs BP).a δD values and concentrations of mid- and long-chain n-alkanes. Arrows (ASM) indicate the maximum strength of the summer monsoon over Asia. b Concentrations of aquatic biomarkers (C20-HBI, alkenones, n-alkenes) and of microbial derived PMI; alkenone index (Uk´3738) and C37-alkenone δD values. c Concentrations of firemarkers (ΣM, L, G). d Titanium, sulfur, strontium, and calcium contents from XRF-scanning22,23. e Lake level reconstructed from ostracod assemblages9. Dark dashed interval 8.4−8.0 kyrs BP indicates mass flow layer. Light and medium grey shaded areas mark episodes of late glacial and mid-Holocene regime shifts within the aquatic ecosystem.Full size imageAlgal biomarkersAside from n-alkanes, a range of compounds of mostly aquatic origin can be identified in the aliphatic and ketone fraction of the sediment extracts. Unsaturated mid-chain alkenes nC21:1, nC23:1, nC25:1, and nC27:1 were abundant in traces in large parts of the core, but exhibit very high concentrations ( >100 µg/g d.w.) in the glacial period in a single sample at 17.5 ka cal BP and between 15 and 14 ka cal BP. Relatively high concentrations up to 25 µg/d.w. were also observed in the mid-Holocene sequence from ca. 9 to 5 ka cal BP (Fig. 2b). The origin of those compounds is not fully resolved and so far n-alkenes have not been detected in samples from aquatic and terrestrial vascular plants on the TP. However, algae such as eustigmatophytes and chlorophytes have been suggested as possible precursors in the open freshwater Lake Lugu, southeastern TP32, and in Lake Challa, eastern Africa33,34. Therefore, we consider phytoplankton as the most likely source of those compounds in Hala Hu.It is notable that different n-alkene distribution patterns are visible throughout the core, with a predominance of nC27:1 and nC25:1 in the mid-Holocene section and nC23:1 and nC25:1 in other core sections (Supplementary Fig. S3–1). This indicates either different source organisms for at least nC23:1 compared to nC27:1 (supported by more depleted δD-values, see below and Supplementary Fig. S3–3) or a change towards the synthesis of longer chain-lengths by the same source organisms, triggered by changing environmental conditions.The C20-highly branched isoprenoid (HBI) compound is another observed phytoplankton biomarker, which is widely absent in the glacial sequences but shows traces throughout the Holocene with peak abundances (up to 30 µg/g d.w) in the mid-Holocene (8−5 ka cal BP; Fig. 2b). As often been found in cyanobacterial and algal mats (e.g.,35,36,37), it has recently been assigned as a trophic indicator derived from diatoms in lake systems38.Other observed biomarkers of algal origin are alkenones, which are primarily produced by haptophytes, even though a variety of species are possible precursors39,40 The summed concentrations of C37-, C38-, and C39-compounds were low (300 ng/g d.w.; ca. 7.8–6.3 ka cal BP) (Fig. 2b). Here, the di- and tri-unsaturated homologues of all chain lengths appear to elute as peaks undisturbed by contamination, while the tetra-unsaturated alkenones show co-elution with a fatty acid ethyl ester (FAEE) at least for the C37-compound.Pentamethylicosane (PMI) is a compound that was detected in most samples in minor concentrations (400 ng/g d.w.) during an episode in the glacial period (ca 16.6−14 ka cal BP) and during the late Holocene (4.5 ka cal BP to present). This compound has been assigned to microorganisms, i.e bacteria and archaea, often related to the methane cycle41,42,43,44.Fire markersAnother analyzed biomarker group, anhydrous sugars levoglucosan (L), galactosan (G), and mannosan (M) are generated by combustion and pyrolysis of cellulose and hemicellulose, thus, are often referenced as “pyromarkers” or “firemarkers”45,46,47. They occur in low concentrations ( More

1.Science Plan and Implementation Strategy IGBP Report No. 53/IHDP Report No. 19 (Global Land Project, 2005).2.Millennium Ecosystem Assessment—Ecosystems and Human Well-Being: Desertification Synthesis Encyclopedia of the Anthropocene vols 1–5 (MEA, 2017).3.Huang, J., Yu, H., Guan, X., Wang, G. & Guo, R. Accelerated dryland expansion under climate change. Nat. Clim. Change 6, 166–171 (2015).Article 

Google Scholar 
4.Pimm, S. L. The complexity and stability of ecosystems. Nature 307, 321–326 (1984).Article 

Google Scholar 
5.Tilman, D. & Downing, J. A. Biodiversity and stability in grasslands. Nature 367, 363–365 (1994).Article 

Google Scholar 
6.Belnap, J. Surface disturbances: their role in acceleration desertification. Environ. Monit. Assess. 37, 38–57 (1995).Article 

Google Scholar 
7.Zhao, Y., Jia, R. L. & Wang, J. Towards stopping land degradation in drylands: water-saving techniques for cultivating biocrusts in situ. Land Degrad. Dev. 30, 2336–2346 (2019).Article 

Google Scholar 
8.Rodriguez-Caballero, E. et al. Dryland photoautotrophic soil surface communities endangered by global change. Nat. Geosci. 11, 185–189 (2018).CAS 
Article 

Google Scholar 
9.Coe, K. K. & Sparks, J. P. Physiology-based prognostic modeling of the influence of changes in precipitation on a keystone dryland plant species. Oecologia 176, 933–942 (2014).Article 

Google Scholar 
10.Ferrenberg, S., Tucker, C. L. & Reed, S. C. Biological soil crusts: diminutive communities of potential global importance. Front. Ecol. Environ. 15, 160–167 (2017).Article 

Google Scholar 
11.Belnap, J. & Gillette, D. A. Soil surface disturbance: impacts on potential wind erodibility of sand desert soils in SE Utah, USA. Land Degrad. Dev. 8, 355–362 (1997).Article 

Google Scholar 
12.Rutherford, W. A. et al. Albedo feedbacks to future climate via climate change impacts on dryland biocrusts. Sci. Rep. 7, 44188 (2017).13.Duniway, M. C. et al. Wind erosion and dust from US drylands: a review of causes, consequences, and solutions in a changing world. Ecosphere 10, e02650 (2019).14.Ferrenberg, S., Faist, A. M., Howell, A. & Reed, S. C. Biocrusts enhance soil fertility and Bromus tectorum growth, and interact with warming to influence germination. Plant Soil 429, 77–90 (2018).CAS 
Article 

Google Scholar 
15.Eldridge, D. J. et al. The pervasive and multifaceted influence of biocrusts on water in the world’s drylands. Glob. Change Biol. 26, 6003–6014 (2020).16.Ferrenberg, S., Reed, S. C. & Belnap, J. Climate change and physical disturbance cause similar community shifts in biological soil crusts. Proc. Natl Acad. Sci. USA 112, 12116–12121 (2015).CAS 
Article 

Google Scholar 
17.Reed, S. C. et al. Changes to dryland rainfall result in rapid moss mortality and altered soil fertility. Nat. Clim. Change 2, 752–755 (2012).CAS 
Article 

Google Scholar 
18.Concostrina-Zubiri, L. et al. Biological soil crusts across disturbance-recovery scenarios: effect of grazing regime on community dynamics. Ecol. Appl. 24, 1863–1877 (2014).CAS 
Article 

Google Scholar 
19.Weber, B., Bowker, M., Zhang, Y. & Belnap, J. in Biological Soil Crusts: An Organizing Principle in Drylands (eds Weber, B., Büdel, B. & Belnap, J.) 479–498 (Springer, 2016).20.Reynolds, J. F. et al. Global desertification: building a science for dryland development. Science 316, 847–851 (2007).CAS 
Article 

Google Scholar 
21.Berdugo, M. et al. Global ecosystem thresholds driven by aridity. Science 367, 787–790 (2020).CAS 
Article 

Google Scholar 
22.Bestelmeyer, B. T. et al. Analysis of abrupt transitions in ecological systems. Ecosphere 2, e03360 (2011).Article 

Google Scholar 
23.Bestelmeyer, B. T. et al. Desertification, land use, and the transformation of global drylands. Front. Ecol. Environ. 13, 28–36 (2015).Article 

Google Scholar 
24.Beisner, B., Haydon, D. & Cuddington, K. Alternative stable states in ecology. Front. Ecol. Environ. 1, 376–382 (2003).25.Fukami, T. & Nakajima, M. Community assembly: alternative stable states or alternative transient states? Ecol. Lett. 14, 973–984 (2011).Article 

Google Scholar 
26.Belnap, J. & Büdel, B. in Biological Soil Crusts: An Organizing Principle in Drylands (eds Weber, B., Büdel, B. & Belnap, J.) 305–320 (Springer, 2016).27.Belnap, J. & Warren, S. D. Measuring restoration success: a lesson from Patton’s tank tracks. Ecol. Bull. 79, 33 (1998).28.Belnap, J. & Elderidge, D. in Biological Soil Crusts: Structure, Function and Management (eds Belnap, J. & Lange, O. L.) 363–383 (Springer, 2001).29.Turner, M. G. Disturbance and landscape dynamics in a changing world. Ecology 91, 2833–2849 (2010).Article 

Google Scholar 
30.Scheffer, M. & Carpenter, S. R. Catastrophic regime shifts in ecosystems: linking theory to observation. Trends Ecol. Evol. 18, 648–656 (2003).Article 

Google Scholar 
31.Sala, O. E. & Lauenroth, W. K. Small rainfall events: an ecological role in semiarid regions. Oecologia 53, 301–304 (1982).32.Cayan, D. R. et al. Future dryness in the Southwest US and the hydrology of the early 21st century drought. Proc. Natl Acad. Sci. USA 107, 21271–21276 (2010).CAS 
Article 

Google Scholar 
33.Christensen, N. S., Wood, A. W., Nathalie, V., Lettenmaier, D. P. & Palmer, R. N. The effects of climate change on the hydrology and water resources of the Colorado river basin. Clim. Change 62, 337 (2004).Article 

Google Scholar 
34.Herrick, J. et al. Field soil aggregate stability kit for soil quality and rangeland health evaluations. Catena 44, 27–35 (2001).Article 

Google Scholar 
35.Escolar, C., Martínez, I., Bowker, M. A. & Maestre, F. T. Warming reduces the growth and diversity of biological soil crusts in a semi-arid environment: implications for ecosystem structure and functioning. Phil. Trans. R. Soc. B 367, 3087–3099 (2012).Article 

Google Scholar 
36.Scheffer, M. et al. Creating a safe operating space for iconic ecosystems: manage local stressors to promote resilience to global change. Science 347, 1317–1319 (2015).CAS 
Article 

Google Scholar 
37.Collins, S. L., Micheli, F. & Hartt, L. A method to determine rates and patterns of variability in ecological communities. Oikos 91, 285–293 (2000).Article 

Google Scholar 
38.Rillig, M. C. et al. The role of multiple global change factors in driving soil functions and microbial biodiversity. Science 366, 886–890 (2019).CAS 
Article 

Google Scholar 
39.IPCC. Climate Change 2014: Impacts, Adaptations, and Vulnerability (eds Field, C. B. et al.) (Cambridge Univ. Press, 2014).40.Mirzabaev, A. et al. in IPCC Special Report on Global Warming of 1.5 °C (eds Masson-Delmotte, V. et al.) Ch. 3 (WMO, 2018).41.Torres-Cruz, T. J. et al. Species-specific nitrogenase activity in lichen-dominated biological soil crusts from the Colorado Plateau, USA. Plant Soil 429, 113–125 (2018).CAS 
Article 

Google Scholar 
42.Omernik, J. M. & Griffith, G. E. Ecoregions of the conterminous United States: evolution of a hierarchical spatial framework. Environ. Manag. 54, 1249–1266 (2014).Article 

Google Scholar 
43.Kuske, C. R., Yeager, C. M., Johnson, S., Ticknor, L. O. & Belnap, J. Response and resilience of soil biocrust bacterial communities to chronic physical disturbance in arid shrublands. ISME J. 6, 886–897 (2011).Article 

Google Scholar 
44.Tucker, C. L., Ferrenberg, S. & Reed, S. C. Climatic sensitivity of dryland soil CO2 fluxes differs dramatically with biological soil crust successional state. Ecosystems 22, 15–32 (2018). https://doi.org/10.1007/s10021-018-0250-445.Cable, J. M. & Huxman, T. E. Precipitation pulse size effects on Sonoran Desert soil microbial crusts. Oecologia 141, 317–324 (2004).Article 

Google Scholar 
46.Karl, T. R., Knight, R. W. & Plummer, N. Trends in high-frequency climate variability in the twentieth century. Nature 377, 217–220 (1995).CAS 
Article 

Google Scholar 
47.Kunkel, K. E., Easterling, D. R., Redmond, K. & Hubbard, K. Temporal variations of extreme precipitation events in the United States: 1895–2000. Geophys. Res. Lett. 30, 1895–2000 (2003).Article 

Google Scholar 
48.Kim, J. A projection of the effects of the climate change induced by increased CO2 on extreme hydrologic events in the Western US. Clim. Change 68, 153–168 (2005).CAS 
Article 

Google Scholar 
49.Smith, S. J. et al. Climate change impacts for the conterminous USA: an integrated assessment part 1. Scenarios and context. Clim. Change 69, 7–25 (2005). https://doi.org/10.1007/1-4020-3876-3_250.Schwinning, S., Belnap, J., Bowling, D. R. & Ehleringer, J. R. Sensitivity of the Colorado Plateau to change: climate, ecosystems, and society. Ecol. Soc. 13, 28 (2008).51.Daly, C. et al. Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States. Int. J. Climatol. 28, 2031–2064 (2008).Article 

Google Scholar 
52.Jonasson, S. The point intercept method for non-destructive estimation of biomass. Phytocoenologia 11, 385–388 (1983).Article 

Google Scholar 
53.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2019).54.Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag, 2016).55.Oksanen, A. J. et al. Vegan: Community Ecology Package. Rpackage version 2.5-7 https://CRAN.R-project.org/package=vegan (2020).56.Anderson, M. J. A new method for non-parametric multivariate analysis of variance. Austral Ecol. 26, 32–46 (2008).Article 

Google Scholar 
57.Venables, W. & Ripley, B. Modern Applied Statistics with S. (Springer, 2002).58.Lenth, R., Singmann, H., Love, J., Buerkner, P. & Herve, M. Package ‘emmeans’ https://github.com/rvlenth/emmeans (2018).59.Signorell, A. DescTools: Tools for Descriptive Statistics (2021).60.Hallett, L. M. et al. codyn: an R package of community dynamics metrics. Methods Ecol. Evol. 7, 1146–1151 (2016).61.Wood, S. N. Generalized Additive Models: An Introduction with R 2nd edn 1–476 (CRC/Taylor & Francis, 2017).62.Bürkner, P. C. brms: an R package for Bayesian multilevel models using Stan. J. Stat. Softw. 80, 1–28 (2017).63.Gelman, A. & Rubin, D. B. Inference from iterative simulation using multiple sequences. Stat. Sci. 7, 547–511 (1992).
Google Scholar 
64.Vehtari, A., Gelman, A. & Gabry, J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat. Comput. 27, 1413–1432 (2017).Article 

Google Scholar 
65.Modrák, M., Barrett, M., Weber, F. & Coronado, E. bayesplot: Plotting for Bayesian Models. R package version 1.8.0 https://mc-stan.org/bayesplot/ (2021).66.Makowski, D., Ben-Shachar, M. & Lüdecke, D. bayestestR: describing effects and their uncertainty, existence and significance within the Bayesian framework. J. Open Source Softw. 4, 1541 (2019).Article 

Google Scholar 
67.Phillips, M. L., Howell, A., Lauria, C. M., Belnap, J. & Reed, S. C. Data and software code from two long-term experiments (1996–2011 and 2005–2018) at three sites on the Colorado Plateau of North America (US Geological Survey, 2021); https://doi.org/10.5066/P9RUN1TP More

Study areaThe North Sea is a continental sea connected to the Atlantic Ocean through the English Channel in the southwest and between northern Shetland along the 61° latitude parallel to Norway in the north (Fig. 1). It is bordered by Norway, Denmark, Germany, the Netherlands, Belgium, France and the UK, and has a surface of 570,000 km2. The North Sea has an average depth of 95 m, yet maximum depths of ca. 700 m are found in the Norwegian Trench. The maximum tidal amplitude of the North Sea can reach up to 8 m, average winter sea surface temperatures are ca. 6 °C and average summer temperatures reach ca. 17°C33. The English Channel encompasses the marine strait between the UK and France. It covers 75,000 km2, has an average depth of 63 m, a maximum depth of 174 m and can reach a maximum tidal amplitude up to 12 m. The average winter and summer sea surface temperatures in the English Channel are ca. 5 and 20 °C, respectively54.TaggingIn total, 320 silver eels were tagged with pop-off archival tags (Table 1; Supplementary Table S2). In Belgium, 238 eels were caught and tagged at a drainage system upstream of the Yser Estuary (hereafter referred to as the Belgian eels) in 2018–2020 via nets that were attached to gravitational discharge sluice gates (coordinates: 51.127 N, 2.761 E) in October, November and December (n2018 = 102, n2019 = 60 and n2020 = 76). In Germany, 82 eels were tagged in 2011 and 2012. In early December 2011, seven eels were caught at Lake Plön (coordinates: 54.137 N, 10.334 E) with fyke nets. During September, October and November 2012, eels were caught in the Rivers Eider (n = 30; coordinates: 54.190 N, 9.093 E) and Havel (n = 45; coordinates: 52.419 N, 12.571 E) with fyke and stow nets, respectively.Upon capture, the eels were anaesthetized with 0.3 ml/L clove oil (Belgium), 0.4 ml/L ethylene glycol monophenyl ether (Germany 2011) or 120 mg/L MS-222 (Germany 2012), and various morphometric characteristics were measured to identify the life stage55: total length (to the nearest mm), weight (to the nearest g), horizontal and vertical eye diameter (to the nearest 0.01 mm in Belgium and to the nearest 0.1 mm in Germany) and pectoral fin length (to the nearest 0.01 mm and 0.1 mm in Belgium and Germany, respectively). Given that their total body length was  > 450 mm, all eels were considered female55. According to the morphometrics, five Belgian eels could be considered in the premigratory stage (FIII); however, based on visual inspection, they were considered silver eels (i.e. silver-coloured abdomen, dark grey on the dorsal side, jaw hinge not proceeding beyond the eye, enlarged eyes and dark coloured pectoral fins). The other 315 eels identified as silver eels based on both morphometry and visual inspection (201 FIV stage and 114 FV stage).Eels weighing ≥ 550 g were externally fitted with a G5 PDST (CEFAS Technology Ltd, UK), which log temperature and pressure (providing information on depth). They were attached applying the three-point Westerberg attachment method56. Two tag types were used: one with a separate tag and pop-off mechanism (Germany) and one where both mechanisms were integrated (Belgium). The flotation collar of the PDSTs was painted bright red, contained contact information and a cash reward to stimulate retrieval by the general public (e.g. beach combers and fishermen). The seven eels caught in 2011 in Germany (minimum 1220 g) were fitted with PSATs (X-Tag, Microwave Telemetry Inc., USA), also using the Westerberg-method56. Like the PDSTs, the PSATs record temperature and pressure. After release, they drift to the surface and transmit the data to the user via the ARGOS satellite system (www.argos-system.org). For the specifications of the different tags, we refer to Supplementary Table S3.Upon recovery from the anaesthetic, eels tagged with PDSTs were released close to their capture locations in the rivers Eider (coordinates 2011: 54.381 N, 9.009 E; coordinates 2012: 54.379 N, 9.013 E), Elbe (coordinates 1: 53.793 N, 9.402 E; coordinates 2: 53.569 N, 9.700 E; coordinates 3: 53.396 N, 10.171 E) and Yser (coordinates: 51.135 N, 2.757 E) (Table 1). The seven eels captured for PSAT tagging in 2011 were held for several weeks in the Thünen Institute of Fisheries Ecology, then tagged and released the same day; others were tagged in the field.PreprocessingOnce downloaded, the temperature and pressure data obtained from the PDSTs was subsampled to 1-min (Belgian eels) or 2-min (German eels) intervals to reduce the datasets and improve geolocation calculation time; this discrepancy is due to the minimum logging rate of the tags (Supplementary Table S3). Linear regression was applied to correct for pressure sensor drift over time. Indeed, pressure values increased over time even if the tag was kept at atmospheric pressure level. The regression was applied between 15 min before release and the moment the tag popped off and reached the surface, since the tag was then considered at sea level and hence to be under zero pressure.The PSAT data were retrieved through the ARGOS satellite system as a subset with 15-min intervals and converted to values of pressure and temperature. Contemporaneous values of temperature and depth were not always transmitted due to the transmission method. As a consequence of the tag release programming, the transmission of the first position for one of the tags was only received five days after the tag reached the sea surface.GeolocationThe daily movements of each electronically tagged European eel were reconstructed using an adapted version of the tidal geolocation model of Pedersen et al.57. The geolocation model uses a novel Fokker–Planck based method to combine the tidal location method of Metcalfe and Arnold58 with a hidden Markov model (HMM), such that an individual’s daily location d is modelled conditionally on its previous location (d − 1), its inferred behavioural state ds, where behaviour is defined by a single diffusivity parameter (i.e. the maximum amount of movement permitted in a given day), and the observations made between d and d − 1. In this case, observations consisted of the recorded depth (m; D1, …, Dn) and temperature (°C; T1, …, Tn), where n is the number of measurements made per day (the HMM down-samples to 10-min intervals, hence 144 measurements per day), and any hydrostatic (tidal) data which are derived from the sinusoidal pressure cycle recorded in the depth data when a fish is at rest on the seafloor. In addition to bathymetry and tidal amplitude with phase, the model was developed to include sea surface temperature (SST), which can provide additional validation when fish are swimming at or near the surface (i.e. depth ≤ 20 m)59,60, and temperature at depth, which can provide additional validation when fish remain at depths well below the sea surface61,62.The model was run in three different configurations for each recovered dataset: (i) using the tidal location model only (as for Pedersen et al.57), hereafter termed TLM geolocation; (ii) using the TLM plus sea surface temperature (as for Wright et al.60), hereafter termed SST geolocation, and (iii) using temperature at the surface and sub-surface, hereafter termed 3D geolocation (Supplemental Fig. S3). The final trajectory output for the PDST Belgian eels and PSAT German eels was obtained via 3D geolocation, while SST geolocation was used for the PDST German eels. The reason for this discrepancy is that the German PDST eels stayed closer to the coast and in shallower water. Consequently, the 3D geolocation results were more prone to error due to coastal influences on water temperature. As a result, we used the SST geolocation method for these datasets to obtain more reliable results.Data for the model were derived from publicly available resources. Gridded global bathymetry data were obtained from the general bathymetric chart of the oceans (Gebco; British Oceanographic Data Centre, Liverpool, United Kingdom, 2009). Tidal constituents were obtained from the Oregon State University Tidal Prediction model, as described in Egbert and Erofeeva63. Sea surface temperature data were sourced from OSTIA64, while temperature at depth data were sourced from the operational Mercator global ocean analysis and forecast system65. These datasets were downloaded from the Copernicus Marine Environmental Model Service (CMEMS: documented here http://resources.marine.copernicus.eu/documents/PUM/CMEMS-GLO-PUM-001-024.pdf). Data were sourced so as to fit the spatial scale of the model (30°N to 80°N and from 110°W to 60°E) and coarsened to reduce model run-time by modifying the spatial grid to a 1/10th of a degree resolution. The output of the model is a nonparametric probability distribution of the geographical position from which a most probable location, for each day at liberty, and a most probable movement path can be estimated.Prior to running the model, a number of constraints and input parameters were defined to ensure that the model ran effectively. The recapture information was either set as (a) the latitude and longitude where the tag was recaptured, with a high confidence ( 200 m) and hence did not exhibit diel vertical migrations, the input estimates of longitude were based on a simple linear interpolation from release to estimated pop-up. However, for eels that did reach oceanic depths, the time of local noon was estimated (based on the timing of significant diel vertical migrations, as for Righton et al.16), and used to estimate longitude. Geolocation was conducted with MATLAB software66.Migration routesOnly datasets containing ≥ 100 km of net tracking distance were included for further analysis, leading to 54 datasets from the 96 retrieved tags and 320 tagged eels. The net tracking distance was identified as the distance along the reconstructed trajectory between the release of the tagged eel and the pop-off event. When an eel was ingested by a predator, leading to the tag tracking the predator rather than the eel, the data were excluded from the day the eel was predated. The 100 km cut-off point was arbitrarily chosen to select migration paths of sufficient length for further analysis (e.g. migration direction); tracks had a minimum deployment duration of 4 days.Migration speedTo exclude a size-effect, we first applied an independent two-sample t-test to confirm eel sizes (i.e. weight) did not differ between Belgian and German eels. The assumptions of normality (Shapiro–Wilk test), homogeneity of variances (F-test) and independence were met (weight measurements are individual-specific and therefore independent).Next, an independent two-sample t-test was conducted to test if the total migration speeds (i.e. the ground speed along the reconstructed trajectory between the release of the tagged eel and the pop-off or predation event) differed between Belgian and German eels. The assumptions were tested and met as described above.Finally, we tested if the daily migration speed (i.e. the ground speed along the reconstructed trajectory per day) differed according to the eel’s position (i.e. modelled latitude and longitude) via a linear mixed effects model. The tag IDs were implemented as a random effect to account for autocorrelation. Since the two-sample t-test showed a significant difference between Belgian and German total migration speeds, we performed a separate analysis on eels from both countries. Assumptions of normality, homogeneity of variances and independence were tested and met.The migration speed analyses were conducted in R (version 3.6.3)67. The packages ‘lme4’ and ‘nlme’ were used to conduct the linear mixed effects model.Ethical statementEels were tagged using approved protocols by trained and individually licensed scientists working under national project authority in accordance with institutional and national guides for the care and use of laboratory animals. These guidelines are consistent with Institutional Review Board/Institutional Animal Care and Use Committee guidelines. Tagging in Belgium was carried out in accordance with the Belgian national and regional regulations for animal welfare and treatment (Permit ID: EC INBO-011). Tagging in Germany followed German legislation concerning care and use of laboratory animals, and ethical permission for the experiments was given by the Ministry of Energy, Agriculture, the Environment, and Rural Areas of the federal state Schleswig–Holstein (reference numbers V312-72241.123-34 (90-8/11) and V311-7224.123.3 (93-6/12) for tagging in 2011 and 2012 respectively). More

Plant defence traitsWe compiled species level data for five plant traits: wood density (WD), leaf and stem spinescence, latex production, and leaf size, for tropical and extra-tropical South and Central American woody species (i.e., the Neotropical biogeographic realm). WD was obtained for 2577 species from ref. 44. We only used wood density data from Zanne et al.44, because this study used WD measured in stems, whereas most other studies with available data used WD measured in branches. Leaf size data were obtained for 2660 woody species from Wright et al.37. We did not include leaf size from herbaceous species because herbaceous and woody species are influence by different megafauna guilds, suggesting distinct mechanisms, and because this dataset37 only included data for 253 Neotropical herbaceous species. The presence or absence of stem (and/or branch) spines (mostly thorns, but also prickles) were obtained from Dantas and Pausas45 for Neotropical savanna and forest species (1004 species) and complemented with other literature sources for other ecoregions (listed in the supplementary materials) using the names of the species for which we had WD and Leaf Size data. Our final stem spines dataset included 2843 woody species. We also compiled data on the presence of latex in plant stems and leaves for all the species for which we had data on other traits (3160 species; references in the supplementary materials). Finally, we also compiled data on leaf spines. While we managed to find leaf spine data for a total of 2173 woody species, we found spinescence in leaves to be especially concentrated in the palm Family (Arecaceae; 198 out of 221 species with leaf spines). Moreover, out of the non-palm species, all but three species also presented stem spines, indicating that, for other taxa, leaf spines might be dependent on the presence of stem spines at the region (in palms, 51% have stem spines). Thus, we only used leaf spinescence data of palm species (694 species) from the global Palm Traits Database 1.046.For wood density and leaf size, we often had more than one trait value per species (1005 and 831 species with more than one trait value, respectively). Thus, we computed the species mean trait value. This rarely occurred for binary traits (spinescence and latex) and, when occurred, the maximum value was used (0 for absence and 1 for presence). This later decision was based in the assumption that omitting the presence of spines or latex is more likely than incorrectly reporting the presence when it is absent. Moreover, some of these traits can be plastic18.From species to ecoregionsWe searched for geographical distribution data (coordinates) from the Global Biodiversity Information Facility (GBIF) for all of the species in each species-trait dataset (Data available from GBIF using the following doi: WD: https://doi.org/10.15468/dl.3vua3x; Stem spines: https://doi.org/10.15468/dl.ar5ddj; Latex: https://doi.org/10.15468/dl.m8dzjd; Leaf spines: https://doi.org/10.15468/dl.vv8gw4; Leaf size: https://doi.org/10.15468/dl.k98nxc). For this search, we used tools provided by the “rgbif” package for R in which species names are updated to the most recent classification and the returned occurrences also include those associated with synonyms (i.e., the “backbone” method). We labelled the obtained geographical coordinates according to their ecoregion and biogeographical realm (following Dinerstein et al.47) and cropped out occurrences falling outside of the Neotropical realm. Since occurrence data was not available to all the species in our initial trait dataset, the number of species used to calculate ecoregion level means was reduced to 2110 species, for wood density, 2133, for leaf size, 2629, for stem spines, 2714, for latex, and 657, for leaf spines. A detailed evaluation of the representativity of this data in relation to ecoregion- and Neotropical- level patterns can be found in the Supplementary Methods. Based on the occurrence data and their ecoregion label, we built a species abundance (columns) by ecoregion (rows) matrix for each trait.We obtained ecoregion scale abundance-weighted means for continuous traits (WD and Leaf Size) by: (1) Multiplying species abundance in each grid cell of the ecoregion by the mean species value; (2) Summing up the row values; (3) dividing the resulting row sum by the total species abundance (row sum prior to trait multiplication), and (4) calculating the ecoregions’ means (across all of the grid cells). For Stem Spines and Latex (binary traits), we used a similar procedure, but the maximum (0 for absence and 1 for presence) value was used instead of the mean in step (1), and step (2) was directly used to calculate the number of presences (i.e., 1 s). Moreover, instead of the steps (3) and (4), we calculated the number of absences as the difference between the total abundance (row sums before trait multiplication) and the values obtained in step (2). This process resulted in weighted means for WD and stem spinescence for 173 ecoregions, and Leaf Size and Latex for 174, out of the 179 Neotropical ecoregions. For leaf spinescence, we used a similar approach, although, because of the fewer species, the abundance estimate from GBIF was less reliable. Thus, we transformed the ecoregion species abundance to presence/absence before multiplying the trait values (0/1 for absence/presence). We obtained leaf spinescence data for 159 out of the 179 Neotropical ecoregions. The species- and ecoregion- level data is provided in the Supplementary Data and in ref. 47.Historical megafauna distributionWe obtained data on historical distribution of megafauna species from the MegaPast2Future/PHYLACINE_1.2 dataset24, a dataset containing distribution maps (96.5 km of spatial resolution) and functional traits for mammal species of the last 130,000 years. From this dataset, we obtained the probable past distribution of extinct large mammal herbivore (hereafter, “megafauna”) species, if these species were still alive today (“Present Natural” scenario; see details below). The “Present Natural” distribution of extinct species in this database is based on the estimated historical distribution (i.e., preceding anthropogenic range modifications) of extant species that are known (from the fossil record) to have coexisted with the extinct species. In this approach, an extinct species is considered to have been present in a given grid cell if at least 50% of the extant species that were found coexisting with the extinct species in the fossil (and subfossil) record was predicted to have occurred in the same cell prior to anthropogenic range modifications24,48. This approach assumes that, since extant and extinct species coexisted in the same locations, they must have had similar ecological requirements. It also assumes that megafauna extinction had anthropogenic causes, instead of causes related to climate change49, which is largely accepted in the literature50.We extracted the “Present Natural” distribution of extinct mammal (coded “EP” for IUCN status; i.e., “extinct in prehistory”, meaning before 1500 CE) whose body mass was higher than 50 kg (megafauna), and for which at least 90% of their diet consisted of plants (i.e., strict herbivores). For each Ecoregion, we began by calculating two megafauna-related metrics: extinct megafauna species richness (Mrich) and their mean body mass (Mbm). For this, we cropped the distribution maps of the megafauna species (containing 1 for presence and 0 for absence of each species) to the Neotropical realm. To calculate Mrich, we (1) counted species presences within each of the grid cells in the global grid (i.e., calculated the cell’s megafauna richness); (2) assigned the corresponding ecoregion label to the resulting richness grid cells, subset the richness cell values corresponding to the Neotropical region; and (3) calculated the mean for each Neotropical ecoregion. For Mbm, we replaced the presences of the megafauna species in the initial raster object (grid cell map of each megafauna species) by their body masses and calculated the grid cell-level mean body mass, before calculating the ecoregion-level means. We also calculated megafauna density and secondary productivity based on allometric equations that relate these metrics to megafauna body mass. However, we did not used megafauna density and secondary productivity because they were strongly correlated to megafauna richness (Supplementary Fig. 3). More details on how these metrics were calculated can be found in the Supplementary Methods.We also obtained diet preference information from the literature for most megafauna species that occurred in the Neotropical region (details and references in the Supplementary Material). Based on these information, we calculated the richness of large browser (MBrich for megabrowser richness), grazer (MGrich for megagrazer richness), and mixed-feeder (MMfrich for mega mixed-feeder richness) species by sub setting the megafauna species by grid cell array before the richness calculation in order to select only species that were classified within the correspondent subgroup.Extant herbivore mammal distributionWe also compiled data on the distribution, body mass and diet of extant and recently extinct (i.e., extinct after 1500 CE) herbivore mammal species (for simplicity, called ‘extant’ species in this study). As with megafauna maps, the distributions used represented reconstructions for periods preceding anthropogenic reduction of extant herbivores ranges (“Present Natural” scenario), based on abiotic, biotic and geographic variables48, rather than the currently observed distribution. This scenario was used because modern anthropogenic range reductions are too recent to produce substantial geographic effects at this spatial scale. These data were obtained by sub setting the MegaPast2Future/PHYLACINE_1.2 dataset to exclude species that were coded “EP” for IUCN status and that were not strict herbivores (at least 90% of the diet constituting of plants). We subsequently associated diet information to these species using data from ref. 51 and excluded all species that did not feed mainly on aboveground vegetative plant tissues (i.e., species that fed mostly on fruits, seed, roots were excluded). This later filtering was because the number of herbivores that feed mostly on seed and fruit increase with decreasing size (and this dataset included small mammals). We subsequently calculated the same metrics as for the extinct megafauna species (except for the richness of mixed-feeders as our source for diets50 labelled species according to dominant feeding pattern). For this, we used the same approach described for extinct megafauna species. We did not use a size threshold for extant species because there were only 13 extant mammal herbivore species with over 50 kg in the Neotropical region, most of which were grazers (9 species; 4 species were mixed-feeders and none were browsers). Therefore, we relied on the mean body mass metric calculated for extant mammals to detect potential size-related effects.Climate, soil, fire, insularity, and hurricanesFor each Ecoregion, we obtained data on climate (mean annual precipitation and temperature, and rainfall seasonality) and soil (sand content, pH, and cation exchange capacity) variables. Climate data was obtained from WorldClim 2.1 (10 min spatial resolution) and was based on climate data from 1970 to 200052. Soil data were obtained from SoilGrids (5 km of spatial resolution)53, and consisted of mean values for two depths, 0.05 and 2 m. We calculated Ecoregion level means for all of the soil and climate variables after intersecting the climate and soil grid maps with the ecoregion map.We obtained the number (a proxy for frequency) and intensity of wildfires per ecoregion area using the MODIS active fire location product (MCD14ML)54. We only considered fires (i.e., hotspots) with detection confidence of 95% or higher occurring from November 2000 to December 2019 (both included). To ensure that only wildfires were considered, we associated each fire pixel with a land cover type (300 m of spatial resolution) from ref. 55 for a buffer area of 1000 m surrounding the fire pixel centroid. We excluded all of the fires occurring in areas in which more than 10% of the surrounding land cover pixels corresponded to agricultural, urban and water classes. We calculated the number of wildfires per ecoregion area by dividing the fire count of each Ecoregion by the ecoregion area, and multiplying the resulting value by the proportion of vegetated land cover pixels (same classes used to exclude fires in anthropogenic areas and water bodies above). Fire intensity was estimated as the average fire radiative power across all detected MODIS hotspots in the ecoregion. Ecoregions lacking large preserved vegetated areas (criteria above) were excluded from subsequent analyses.Using the ecoregion map, we also classified ecoregions into insular (1), when most of the ecoregion area was located in islands, vs. continental (0), otherwise. This was performed because island biogeography theory predicts that, in island, species richness should be low due to low colonization and high extinction rates. Insularity has also been shown to reduce megafauna body size (i.e., the island rule), even though the mechanisms are not fully understood56. We also compiled data on hurricane activity, as woody density was suggested to confer resistance against this disturbance57. We used data from 1990 to 2019 from the HURDAT2 dataset58, containing six-hourly information about the location of all of the known tropical and subtropical cyclones (0.1° latitude/longitude). We used the sum of hurricane occurrences per ecoregions divided by ecoregion area as an indicator of hurricane activity.Statistical analysesTo understand megafauna patterns, we began by fitting (multiple) regression models with habitat-related (fire, climate, soil) and geographical (insularity) variables as predictors. We expected that megafauna richness in general was higher under savanna conditions (arid nutrient-rich or mesic nutrient-poor environments with frequent fires)1,22. We also expected that megafauna richness and body mass were affected negatively by insularity (i.e., following the island biogeography theory and island rule). Before the analyses, we tested the correlations among all of the variables that would eventually be entered as predictors in the same model for both the megafauna and trait models (Supplementary Table 1), in order to avoid multicollinearity associated with highly correlated variables (here, r ≥ 0.60). Since mean annual precipitation and soil pH were strongly positively correlated (r = −0.78), for all of the analyses (including the analyses with functional traits, described below), model selection was performed separately for these two variables (i.e., two different model selection procedures, one containing each of the two variables among the initial set of predictors). We selected the best among the two resulting models as that with the lowest AIC (differences higher than two points in all of the cases). To make sure that no multicollinearity remained we also calculated the Variation Inflation Factor (VIF) for all of the predictor variables as 1/tolerance, where tolerance is calculated as 1 minus the R2 of all of the model regressing a predictive variable against all of the other predictors. In all of the models, VIF was 3.33 or smaller (i.e., a tolerance of 0.30 or higher), indicating absence of multicollinearity.Model simplification was carried interactively using stepwise (both forward and backward) searching for the model with the lowest AIC (using R’s “step” function) and subsequently retaining only the significant variables (p ≤ 0.05). We calculated the Pearson r statistics as a measure of effect size for the selected variables as well as the associated confidence intervals, using the packages “parameters” and “effectsize” for R. The average contribution of each predictor variable was also calculated, using the package “dominanceanalysis”, as the mean difference in R2 before and after removing the target variable from models containing all of the possible subset combinations of the selected predictor variables, including the full selected model.For testing whether the studied plant functional traits were related to our megafauna indicators, we fit linear models to WD and leaf size, and generalized linear models (GLM; binomial family) for spinescence and latescence, using ecoregion as the unit. For spinescence and latescence, we used the matrix containing the count of spiny/latex and non-spiny/non-latex plants (species abundance; for stem spines and latex) or number of species with or without spines (for leaf spines; see above) as response variables. The predictor variables included the animal indicators for extinct megafauna and extant herbivores, as well as climate, soil, and fire predictors (and, for WD, hurricane counts). Because total, as well as megagrazer, megabrowser, and mega mixed-feeder species richness were strongly positively correlated (Supplementary Table 1), we used the richness difference between grazers and browsers to evaluate the effect of diet (Supplementary Fig. 1). For consistency, we used the same diet variable for extant and extinct species. Since we did not identify strong correlations among extinct megafauna and extant herbivore indicators (Supplementary Table 1), these variables were all entered simultaneously in the same initial models. As with the analyses of the megafauna indices, we also used r as effect size and calculated the average predictor contribution in terms of R2 for these models. For the later, we used the MacFadden Pseudo-R2 in the GLM models as implemented in the “pscl” and “dominanceanalysis” packages for R, as this statistic is the most comparable with R2 from linear multiple regression (Maximum Likelihood and Cragg and Uhler’s Pseudo-R2 were also calculated for the logistic models), and adjusted R2 for continuous traits. Islands were not included in these models, as island plants were expected to respond differently due to the effects of insularity on animal species richness, precluding megafauna and extant mammal richness from being accurate proxies for consumer abundance. For stem spines, we always included a quadratic term to both megafauna and extant mammal herbivore body mass, as evidence suggest that medium-size herbivores (i.e., approximately 250 kg) are important selective drivers of this trait12. If a significant relationship with our herbivory indicators (both extant and extinct) were significant but not indicative of a selective effect by herbivores (for more defended plants), this relationship was discarded (along with related variables, such as diet); this happened only once, for leaf size, which increased with extant herbivore richness (Supplementary Table 8).For all of the general linear regression models, assumptions of normality, homoscesticity and lack of spatial autocorrelation in the residuals were checked using the Kolmogorov–Smirnov, Breusch–Pagan and Moran’s I tests, respectively. For the later, ecoregions were considered neighbours when they were adjacent and non-neighbour otherwise. In some cases, heteroscesticity was detected and, thus, the significance of the coefficients was tested using heteroskedasticity-consistent covariance matrix estimation. If one or more variable lost their significance they were stepwise removed from the final model, beginning by the least significant, until all remaining variables had a significant effect. Overdispersion in the generalized linear model was also detected and dealt with using overdispersed binomial logit models, as implemented in the “dispmod” package for R, in which weights are interactively calculated and used to maintain the residual deviance lower than the degrees of freedom. To confirm that the detected associations between megafauna indices and plant traits were robust, we also tested the coefficient significance using randomization of the plant species by ecoregion matrices (see Supplementary Methods for details).To test the prediction that Neotropical ecoregions could be broadly classified into the three hypothesised antiherbiomes, we used hierarchical clustering on principal component axes of the ecoregion by trait matrix (five plant traits, standardized to zero mean and unit variance). We selected the number of clusters associated with the highest loss of inertia (within group variability) when progressively increasing the number of clusters, using the R package “FactoMineR”. This procedure allowed the recognition of large regions characterised by specific patterns of defence strategies (‘antiherbiomes’). We subsequently tested for axes score, megafauna and environmental differences among the resulting antiherbiomes to verify whether and how trait, climate and soil patterns matched those described for African ecosystems, and to understand the megafaunal differences among the antiherbiomes. For these comparisons, we used Kruskal-Wallis and post-hoc pairwise Dunn tests, using the Benjamini & Hochberg59 (1995) correction of P-values for multiple comparisons in both cases, and exclusively included continental ecoregions. For spines, we used the proportion of spinescent plants/species (rather than the number of “yes” and “no” used on previous analyses) in the principal component analysis. Because palms were missing from 20 ecoregions, we completed the values for these ecoregions using predicted model probabilities. To better understand these associations between traits and the environmental and megafauna variables, we also regressed the PCA axes against the same predictors used for traits.We also developed a framework to identify forest ecoregions most likely to have experienced a biome shift after megafauna extinction using antiherbiome, biome and megafauna distribution data. Ecoregions likely to have experienced a savanna-to-forest shift since the Pleistocene are those that: (1) are currently forest-dominated; (2) are classified in antiherbiomes analogous to African arid nutrient-rich or mesic nutrient-poor savannas; and (3) were megafauna- and, especially, megagrazer- rich during the Pleistocene (richness equal or greater than the 0.75 quantile: 14 species for Mrich, and 3 for exclusively grazing species; MGrich). We validated the distribution of these areas with fossil evidence (22 sites) from the Last Glacial Maximum and mid-Holocene (see Supplementary Methods and Supplementary Table 9). For this, we also used information about the present dominant vegetation type in the fossil sites, extracted from the reference sources (see Supplementary Table 9), to segregate savanna-forest shifts from data coming from stable savanna patches within forest or long-term savanna regions. We also contrasted the predicted patterns with the present location of savanna patches within the Amazon Forest region from ref. 60.All statistical analyses and data handling were carried out in the R (v.4.0.2) environment, using the previously mentioned packages, in addition to FSA, gridExtra, grid, lattice, lmtest, latticeExtra, olsrr, raster, rgdal, rgeos, sandwich, spatialreg, spdep and vegan, using codes provided in ref. 47.Reporting SummaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.Gaston, K. J., Jackson, S. F., Cantú-Salazar, L. & Cruz-Piñón, G. The ecological performance of protected areas. Annu. Rev. Ecol. Evol. Syst. 39, 93–113 (2008).
Google Scholar 
2.Staver, A. C., Abraham, J. O., Hempson, G. P., Karp, A. T. & Faith, J. T. The past, present, and future of herbivore impacts on savanna vegetation. J. Ecol. 109, 2804–2822 (2021).
Google Scholar 
3.Bond, W. J. Keystone species. In Biodiversity and Ecosystem Function (eds Schulze, E. D. & Mooney, H. A.) 237–253 (Springer, 1994).
Google Scholar 
4.Cole, M. M. The influence of soils, geomorphology and geology on the distribution of plant communities in savanna ecosystems. In Ecology of Tropical Savannas (eds Huntley, B. J. & Walker, B. H.) 145–174 (Springer, 1982).
Google Scholar 
5.Huntley, B. J. & Walker, B. H. (eds) Ecology of Tropical Savannas Vol. 42 (Springer, 1982).
Google Scholar 
6.Frost, P. et al. Responses of Savannas to Stress and Disturbance: A Proposal for a Collaborative Programme of Research (Biology International 10, 1985).
Google Scholar 
7.Medina, E. & Silva, J. F. Savannas of northern South America: A steady state regulated by water–ire interactions on a background of low nutrient availability. J. Biogeogr. 17, 403–413 (1990).
Google Scholar 
8.Fensham, R. J., Fairfax, R. J. & Archer, S. R. Rainfall, land use and woody vegetation cover change in semi-arid Australian savanna. J. Ecol. 93, 596–606 (2005).
Google Scholar 
9.Mucina, L. & Rutherford, M. C. The Vegetation of South Africa, Lesotho and Swaziland (South African National Biodiversity Institute, 2006).
Google Scholar 
10.Staver, A. C., Botha, J. & Hedin, L. Soils and fire jointly determine vegetation structure in an African savanna. New Phytol. 216, 1151–1160 (2017).CAS 
PubMed 

Google Scholar 
11.Jakubka, D. et al. Effects of climate, habitat and land use on the cover and diversity of the savanna herbaceous layer in Burkina-Faso, West Africa. Folia Geobot. 52, 129–142 (2017).
Google Scholar 
12.Bond, W. J. Open Ecosystems: Ecology and Evolution Beyond the Forest Edge (Oxford University Press, 2019).
Google Scholar 
13.O’Connor, T. G. Composition and population responses of an African savanna grassland to rainfall and grazing. J. Appl. Ecol. 31, 155–171 (1994).
Google Scholar 
14.Walker, B. & Langridge, J. Predicting savanna vegetation structure on the basis of plant available moisture (PAM) and plant available nutrients (PAN): A case study from Australia. J. Biogeogr. 24, 813–825 (1997).
Google Scholar 
15.Sankaran, M. et al. Determinants of woody cover in African savannas. Nature 438, 846–849 (2005).CAS 
ADS 
PubMed 

Google Scholar 
16.Bucini, G. & Hanan, N. P. A continental-scale analysis of tree cover in African savannas. Glob. Ecol. Biogeogr. 16, 593–605 (2007).
Google Scholar 
17.Venter, F. J., Scholes, R. J. & Eckhardt, H. C. The abiotic template and its associated vegetation pattern. In The Kruger Experience: Ecology and Management of Savanna Heterogeneity (eds du Toit, J. T. et al.) 83–129 (Island Press, Washington, D.C., 2003).18.Valeix, M. et al. Vegetation structure and ungulate abundance over a period of increasing elephant abundance in Hwange National Park, Zimbabwe. J. Trop. Ecol. 23, 87–93 (2007).
Google Scholar 
19.Asner, G. P. et al. Large-scale impacts of herbivores on the structural diversity of African savannas. Proc. Natl Acad. Sci. USA 106, 4947–4952 (2009).CAS 
ADS 
PubMed 
PubMed Central 

Google Scholar 
20.Archibald, S. & Hempson, G. P. Competing consumers: Contrasting the patterns and impacts of fire and mammalian herbivory in Africa. Philos. Trans. R. Soc. B 371, 20150309 (2016).
Google Scholar 
21.Archibald, S., Bond, W. J., Stock, W. D. & Fairbanks, D. H. K. Shaping the landscape: Fire–grazer interactions in an African savanna. Ecol. Appl. 15, 96–109 (2005).
Google Scholar 
22.Staver, A. C. & Bond, W. J. Is there a ‘browse trap’? Dynamics of herbivore impacts on trees and grasses in an African savanna. J. Ecol. 102, 595–602 (2014).
Google Scholar 
23.Jeltsch, F., Weber, G. E. & Grimm, V. Ecological buffering mechanisms in savannas: A unifying theory of long-term tree-grass coexistence. Plant Ecol. 105, 161–171 (2000).
Google Scholar 
24.February, E. C., Higgins, S. I., Bond, W. J. & Swemmer, L. Influence of competition and rainfall manipulation on the growth responses of savanna trees and grasses. Ecology 94, 1155–1164 (2013).PubMed 

Google Scholar 
25.Savadogo, P., Tiveau, D., Sawadogo, L. & Tigabu, M. Herbaceous species responses to long-term effects of prescribed fire, grazing and selective tree cutting in the savanna-woodlands of West Africa. Perspect. Plant Ecol. Evol. Syst. 10, 179–195 (2008).
Google Scholar 
26.Smith, M. D. et al. Long-term effects of fire frequency and season on herbaceous vegetation in savannas of the Kruger National Park, South Africa. J. Plant Ecol. 6, 71–83 (2013).
Google Scholar 
27.Thrash, I. Impact of water provision on herbaceous vegetation in Kruger National Park, South Africa. J. Arid Environ. 38, 437–450 (1998).ADS 

Google Scholar 
28.Staver, A. C., Bond, W. J., Stock, W. D., van Rensburg, S. J. & Waldram, M. S. Browsing and fire interact to suppress tree density in an African savanna. Ecol. Appl. 19, 1909–1919 (2009).PubMed 

Google Scholar 
29.Smit, I. P. J. & Ferreira, S. M. Management intervention affects river-bound spatial dynamics of elephants. Biol. Conserv. 143, 2172–2181 (2010).
Google Scholar 
30.Loarie, S. R., van Aarde, R. J. & Pimm, S. L. Elephant seasonal vegetation preferences across dry and wet savannas. Biol. Conserv. 142, 3099–3107 (2009).
Google Scholar 
31.Young, K. D., Ferreira, S. M. & Van Aarde, R. J. Elephant spatial use in wet and dry savannas of southern Africa. J. Zool. 278, 189–205 (2009).
Google Scholar 
32.Codron, J. et al. Elephant (Loxodonta africana) diets in Kruger National Park, South Africa: spatial and landscape differences. J. Mammal. 87, 27–34 (2006).
Google Scholar 
33.Timberlake, J. Colophospermum mopane: Annotated Bibliography and Review (Zimbabwe Bulletin of Forestry Research No. 11, 1995).
Google Scholar 
34.Pyšek, P. et al. Into the great wide open: Do alien plants spread from rivers to dry savanna in the Kruger National Park?. NeoBiota 60, 61–77 (2020).
Google Scholar 
35.Eckhardt, H. C., van Wilgen, B. W. & Biggs, H. C. Trends in woody vegetation cover in the Kruger National Park, South Africa, between 1940 and 1998. Afr. J. Ecol. 38, 108–115 (2000).
Google Scholar 
36.Groen, T. A. Spatial Matters: How Spatial Patterns and Processes Affect Savanna Dynamics. PhD Thesis, Wageningen University, The Netherlands (2007).37.MacFadyen, S., Hui, C., Verburg, P. H. & Van Teeffelen, A. J. A. Quantifying spatiotemporal drivers of environmental heterogeneity in Kruger National Park, South Africa. Landsc. Ecol. 31, 2013–2029 (2016).
Google Scholar 
38.Munyati, C. & Ratshibvumo, T. Differentiating geological fertility derived vegetation zones in Kruger National Park, SouthAfrica, using Landsat and MODIS imagery. J. Nat. Conserv. 18, 169–179 (2010).
Google Scholar 
39.Colgan, M. S., Asner, G. P., Levick, S. R., Martin, R. E. & Chadwick, O. A. Topo-edaphic controls over woody plant biomass in South African savannas. Biogeosciences 9, 1809–1821 (2012).ADS 

Google Scholar 
40.Walter, H. & Burnett, J. H. Ecology of Tropical and Subtropical Vegetation Vol. 539 (Oliver and Boyd, 1971).
Google Scholar 
41.Scholes, R. J., Bond, W. J. & Eckhardt, H. C. Vegetation dynamics in the Kruger ecosystem. In The Kruger Experience: Ecology and Management of Savanna Heterogeneity (eds du Toit, J. T. et al.) 243–262 (Island Press, 2003).
Google Scholar 
42.Knoop, W. T. & Walker, B. H. Interactions of woody and herbaceous vegetation in southern African savanna. J. Ecol. 73, 235–253 (1985).
Google Scholar 
43.Tilman, D. The resource-ratio hypothesis of plant succession. Am. Nat. 125, 827–852 (1985).
Google Scholar 
44.Scholes, R. J. & Archer, S. R. Tree–grass interactions in savannas. Annu. Rev. Ecol. Syst. 28, 517–544 (1997).
Google Scholar 
45.Muvengwi, J., Davies, A. B., Parrini, F. & Witkowski, E. T. F. Contrasting termite diversity and assemblages on granitic and basaltic African savanna landscapes. Insect. Soc. 65, 25–35 (2018).
Google Scholar 
46.Trollope, W. S. W., Potgieter, A. L. F. & Zambatis, N. Assessing veld condition in the Kruger National Park using key grass species. Koedoe 32, 67–93 (1989).
Google Scholar 
47.Ehleringer, J. R. & Monson, R. K. Evolutionary and ecological aspects of photosynthetic pathway variation. Annu. Rev. Ecol. Syst. 24, 411–439 (1993).
Google Scholar 
48.Chytrý, M., Tichý, L. & Roleček, J. Local and regional patterns of species richness in Central European vegetation types along the pH/calcium gradient. Folia Geobot. 38, 429–442 (2003).
Google Scholar 
49.Owen-Smith, N., Page, B., Teren, G. & Druce, D. J. Megabrowser impacts on woody vegetation in savannas. In Savanna Woody Plants and Large Herbivores (eds Scogings, P. F. & Sankaran, M.) 585–611 (Wiley, 2019).
Google Scholar 
50.Nasseri, N. A., McBrayer, L. D. & Schulte, B. A. The impact of tree modification by African elephant (Loxodonta africana) on herpetofaunal species richness in northern Tanzania. Afr. J. Ecol. 49, 133–140 (2010).
Google Scholar 
51.Asner, G. P. & Levick, S. R. Landscape-scale effects of herbivores on treefall in African savannas. Ecol. Lett. 15, 1211–1217 (2012).PubMed 

Google Scholar 
52.Haynes, G. Elephants (and extinct relatives) as earth-movers and ecosystem engineers. Geomorphology 157, 99–107 (2012).ADS 

Google Scholar 
53.Kruger, L. M., Coetzee, J. A. & Vickers, K. The Impacts of Elephants on Woodlands and Associated Biodiversity (Summary report to South African National Parks, Organization for Tropical Studies, 2007).54.Guy, P. R. The influence of elephants and fire on a Brachystegia julbernardia woodland in Zimbabwe. J. Trop. Ecol. 5, 215–226 (1989).
Google Scholar 
55.Laws, R. M., Parker, I. S. C. & Johnstone, R. C. B. Elephants and Their Habitats: The Ecology of Elephants in North Bunyoro, Uganda (Clarendon Press, 1975).
Google Scholar 
56.Thompson, P. J. The role of elephants, fire and other agents in the decline of Brachystegia woodlands. J. S. Afr. Wildl. Manag. Assoc. 5, 11–18 (1975).
Google Scholar 
57.Barnes, M. E. Effects of large herbivores and fire on the regeneration of Acacia erioloba woodlands in Chobe National Park, Botswana. Afr. J. Ecol. 39, 340–350 (2001).
Google Scholar 
58.Scogings, P. F., Johansson, T., Hjältén, J. & Kruger, J. Responses of woody vegetation to exclusion of large herbivores in semi-arid savannas. Austral Ecol. 37, 56–66 (2012).
Google Scholar 
59.Verweij, R. J. T., Higgins, S. I., Bond, W. J. & February, E. C. Water sourcing by trees in a mesic savanna: Responses to severing deep and shallow roots. Environ. Exp. Bot. 74, 229–236 (2011).
Google Scholar 
60.Smit, I. et al. Effects of fire on woody vegetation structure in African savanna. Ecol. Appl. 20, 1865–1875 (2010).PubMed 

Google Scholar 
61.MacFadyen, S., Hui, C., Verburg, P. H. & Van Teeffelen, A. J. A. Spatiotemporal distribution dynamics of elephants in response to density, rainfall, rivers and fire in Kruger National Park, South Africa. Divers. Distrib. 25, 880–894 (2019).
Google Scholar 
62.O’Connor, T. G., Goodman, P. S. & Clegg, B. A functional hypothesis of the threat of local extirpation of woody plant species by elephant in Africa. Biol. Conserv. 136, 329–345 (2007).
Google Scholar 
63.Didan, K. MOD13Q1 MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid V006. NASA EOSDIS Land Processes DAAC. Accessed 2021-06-08 from https://doi.org/10.5067/MODIS/MOD13Q1.006 (2015).64.Kreft, H. & Jetz, W. Global patterns and determinants of vascular plant diversity. Proc. Natl Acad. Sci. 104, 5925–5930 (2007).CAS 
ADS 
PubMed 
PubMed Central 

Google Scholar 
65.Bohdalková, E., Toszogyova, A., Šímová, I. & Storch, D. Universality in biodiversity patterns: variation in species-temperature and species-productivity relationships reveals a prominent role of productivity in diversity gradients. Ecography 44, 1366–1378 (2021).
Google Scholar 
66.Knight, R. S., Crowe, T. M. & Siegfried, W. R. Distribution and species richness of trees in southern Africa. J. S. Afr. Bot. 48, 455–480 (1982).
Google Scholar 
67.Gaylard, A., Owen-Smith, N. & Redfern, J. Surface water availability: implications for heterogeneity and eco-system processes. In The Kruger Experience: Ecology and Management of Savanna Heterogeneity (eds du Toit, J. T. et al.) 171–188 (Island Press, 2003).
Google Scholar 
68.Venter, F. J. A Classification of Land for Management Planning in the Kruger National Park. PhD Thesis, University of South Africa (1990).69.du Toit, J. T. et al. (eds) The Kruger Experience: Ecology and Management of Savanna Heterogeneity (Island Press, 2003).
Google Scholar 
70.Smit, I. P. J., Smit, C. F., Govender, N., van der Linde, M. & MacFadyen, S. Rainfall, geology and landscape position generate large-scale spatiotemporal fire pattern heterogeneity in an African savanna. Ecography 36, 447–459 (2013).
Google Scholar 
71.Obermeijer, A. A. A preliminary list of plants found in the Kruger National Park. Ann. Transvaal Mus. 17, 185–227 (1937).
Google Scholar 
72.van der Schijff, H. P. ‘n Ekologiese Studie van die Flora van die Nasionale Krugerwildtuin. D.Sc. Thesis, Potchefstroom University (1957).73.van der Schijff, H. P. The affinities of the flora of the Kruger National Park. Kirkia 7, 109–120 (1968).
Google Scholar 
74.Coetzee, B. J. Phytosociology, Vegetation Structure and Landscapes of the Central District. Kruger National Park, South Africa. Dissertationes Botanicae, J. Cramer, Vaduz (1983).75.Gertenbach, W. P. D. Landscapes of the Kruger National Park. Koedoe 26, 9–121 (1983).
Google Scholar 
76.Siebert, F. & Eckhardt, H. C. The vegetation and floristics of the Nkhuhlu exclosures, Kruger National Park. Koedoe 50, 126–144 (2008).
Google Scholar 
77.Siebert, F., Eckhardt, H. C. & Siebert, S. J. The vegetation and floristics of the Letaba exclosures, Kruger National Park, South Africa. Koedoe 52, 1–12 (2010).
Google Scholar 
78.Wigley, B. J., Fritz, H., Coetsee, C. & Bond, W. J. Herbivores shape woody plant communities in the Kruger National Park: Lessons from three long-term exclosures. Koedoe 56, 1–12 (2014).
Google Scholar 
79.Enslin, B. W., Potgieter, A. L. F., Biggs, H. C. & Biggs, R. Long term effects of fire frequency and season on the woody vegetation dynamics of the Sclerocarya birrea/Acacia nigrescens savanna of the Kruger National Park. Koedoe 43, 27–37 (2000).
Google Scholar 
80.Brits, J., Van Rooyen, M. W. & Van Rooyen, N. Ecological impact of large herbivores on the woody vegetation at selected watering points on the eastern basaltic soils in the Kruger National Park. Afr. J. Ecol. 40, 53–60 (2002).
Google Scholar 
81.Todd, S. W. Gradients in vegetation cover, structure and species richness of Nama-Karoo shrublands in relation to distance from livestock watering points. J. Appl. Ecol. 43, 293–304 (2006).
Google Scholar 
82.Foxcroft, L. C., Henderson, L., Nichols, G. R. & Martin, B. W. A revised list of alien plants for the Kruger National Park. Koedoe 46, 21–44 (2003).
Google Scholar 
83.Foxcroft, L. C., Richardson, D. M. & Wilson, J. R. Ornamental plants as invasive aliens: Problems and solutions in Kruger National Park, South Africa. Environ. Manag. 41, 32–51 (2008).ADS 

Google Scholar 
84.Mueller-Dombois, D. & Ellenberg, H. Aims and Methods of Vegetation Ecology (Wiley, 1974).
Google Scholar 
85.van der Maarel, E. Transformation of cover-abundance values in phytosociology and its effects on community similarity. Vegetatio 38, 97–114 (1979).
Google Scholar 
86.Magurran, A. E. Ecological Diversity and Its Measurement (Croom Helm, 1988).
Google Scholar 
87.Pooley, E. Wildflowers of Kwazulu-Natal and the Eastern Region (Natal Flora Publications Trust, 1998).
Google Scholar 
88.Schmidt, E., Lötter, M. & McCleland, W. Trees and Shrubs of Mpumalanga and Kruger National Park (Jacana Media, 2002).
Google Scholar 
89.van der Walt, R. Wild Flowers of the Limpopo Valley (Retha van der Walt, 2009).
Google Scholar 
90.Oudtshoorn, F. V. Guide to Grasses of Southern Africa 3rd edn. (Briza Publications, 2018).
Google Scholar 
91.R Development Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, 2013).
Google Scholar 
92.Lenth, R. Emmeans: Estimated Marginal Means, aka Least-Square Means. R package version 1.2.1. https://CRAN.R-project.org/package=emmeans (2018).93.Lepš, J. & Šmilauer, P. Multivariate Analysis of Ecological Data Using CANOCO 5 (Cambridge University Press, 2014).MATH 

Google Scholar 
94.Dray, S., Legendre, P. & Peres-Neto, P. R. Spatial modelling: A comprehensive framework 562 for principal coordinate analysis of neighbour matrices (PCNM). Ecol. Modell. 196, 483–563 (2006).
Google Scholar 
95.Šmilauer, P. & Lepš, J. Multivariate Analysis of Ecological Data Using Canoco 5 2nd edn. (Cambridge University Press, 2014).MATH 

Google Scholar  More

1.Bastviken, D. Methane. in Encyclopedia of Inland Waters (ed. Likens, G. E.) 783–805 (Elsevier, 2009). https://doi.org/10.1016/B978-012370626-3.00117-42.Hanson, R. S. & Hanson, T. E. Methanotrophic bacteria. Microbiol. Rev. 60, 439–471 (1996).CAS 
PubMed 
PubMed Central 

Google Scholar 
3.Thottathil, S. D., Reis, P. C. J., del Giorgio, P. A. & Prairie, Y. T. The extent and regulation of summer methane oxidation in Northern Lakes. J. Geophys. Res. Biogeosciences 123, 3216–3230 (2018).CAS 
ADS 

Google Scholar 
4.Kankaala, P., Taipale, S., Nykänen, H. & Jones, R. I. Oxidation, efflux, and isotopic fractionation of methane during autumnal turnover in a polyhumic, boreal lake. J. Geophys. Res. Biogeosciences 112, 1–7 (2007).
Google Scholar 
5.Kankaala, P., Huotari, J., Peltomaa, E., Saloranta, T. & Ojala, A. Methanotrophic activity in relation to methane efflux and total heterotrophic bacterial production in a stratified, humic, boreal lake. Limnol. Oceanogr. 51, 1195–1204 (2006).CAS 
ADS 

Google Scholar 
6.Bastviken, D., Ejlertsson, J., Sundh, I. & Tranvik, L. Methane as a source of carbon and energy for lake pelagic food webs. Ecology 84, 969–981 (2003).
Google Scholar 
7.Kankaala, P., Lopez Bellido, J., Ojala, A., Tulonen, T. & Jones, R. I. Variable production by different pelagic energy mobilizers in Boreal Lakes. Ecosystems 16, 1152–1164 (2013).CAS 

Google Scholar 
8.Morana, C. et al. Methanotrophy within the water column of a large meromictic tropical lake (Lake Kivu, East Africa). Biogeosciences 12, 2077–2088 (2015).ADS 

Google Scholar 
9.Grey, J. The incredible lightness of being methane-fuelled: stable isotopes reveal alternative energy pathways in aquatic ecosystems and beyond. Front. Ecol. Evol. 4, 1–14 (2016).
Google Scholar 
10.Jones, R. I. & Grey, J. Biogenic methane in freshwater food webs. Freshw. Biol. 56, 213–229 (2011).CAS 

Google Scholar 
11.Kankaala, P., Taipale, S. & Grey, J. Experimental d13C evidence for a contribution of methane to pelagic food webs in lakes. Limnol. Oceanogr. 51, 2821–2827 (2006).CAS 
ADS 

Google Scholar 
12.Guérin, F. & Abril, G. Significance of pelagic aerobic methane oxidation in the methane and carbon budget of a tropical reservoir. J. Geophys. Res. 112, 1–14 (2007).
Google Scholar 
13.Rasilo, T., Hutchins, R. H. S., Ruiz-González, C. & del Giorgio, P. A. Transport and transformation of soil-derived CO2, CH4 and DOC sustain CO2 supersaturation in small boreal streams. Sci. Total Environ. 579, 902–912 (2017).CAS 
PubMed 
ADS 

Google Scholar 
14.Soued, C. & Prairie, Y. T. The carbon footprint of a Malaysian tropical reservoir: measured versus modeled estimates highlight the underestimated key role of downstream processes. Biogeosciences 17, 515–227 (2020).CAS 
ADS 

Google Scholar 
15.Del Giorgio, P. A. & Gasol, J. M. Physiological structure and single-cell activity in marine bacterioplankton. in Microbial Ecology of the Oceans: Second Edition (ed. Kirchman, D. L.) 243–298 (John Wiley & Sons, Inc, 2008). https://doi.org/10.1002/9780470281840.ch816.Reis, P. C. J., Ruiz-González, C., Soued, C., Crevecoeur, S. & Prairie, Y. T. Rapid shifts in methanotrophic bacterial communities mitigate methane emissions from a tropical hydropower reservoir and its downstream river. Sci. Total Environ. 748, 141374 (2020).CAS 
PubMed 
ADS 

Google Scholar 
17.Thottathil, S. D., Reis, P. C. J. & Prairie, Y. T. Methane oxidation kinetics in northern freshwater lakes. Biogeochemistry 143, 105–116 (2019).CAS 

Google Scholar 
18.Milucka, J. et al. Methane oxidation coupled to oxygenic photosynthesis in anoxic waters. ISME J. 9, 1991–2002 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
19.Zigah, P. K. et al. Methane oxidation pathways and associated methanotrophic communities in the water column of a tropical lake. Limnol. Oceanogr. 60, 553–572 (2015).ADS 

Google Scholar 
20.Mayr, M. J. et al. Growth and rapid succession of methanotrophs effectively limit methane release during lake overturn. Commun. Biol. 3, 1–9 (2020).
Google Scholar 
21.Bussmann, I., Rahalkar, M. & Schink, B. Cultivation of methanotrophic bacteria in opposing gradients of methane and oxygen. FEMS Microbiol. Ecol. 56, 331–344 (2006).CAS 
PubMed 

Google Scholar 
22.Kankaala, P., Eller, G. & Jones, R. I. Could bacterivorous zooplankton affect lake pelagic methanotrophic activity? Fundam. Appl. Limnol. / Arch. f.ür. Hydrobiol. 169, 203–209 (2007).
Google Scholar 
23.Khmelenina, V. N. et al. Structural and functional features of methanotrophs from hypersaline and alkaline lakes. Microbiology 79, 472–482 (2010).CAS 

Google Scholar 
24.Westfall, C. S. & Levin, P. A. Bacterial cell size: multifactorial and multifaceted. Annu. Rev. Microbiol. 71, 499–517 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
25.Chien, A. C., Hill, N. S. & Levin, P. A. Cell size control in bacteria. Curr. Biol. 22, R340–R349 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
26.Velimirov, B. Nanobacteria, ultramicrobacteria and starvation forms: a search for the smallest metabolizing bacterium. Microbes Environ. 16, 67–77 (2001).
Google Scholar 
27.Reis, P. C. J., Thottathil, S. D., Ruiz-González, C. & Prairie, Y. T. Niche separation within aerobic methanotrophic bacteria across lakes and its link to methane oxidation rates. Environ. Microbiol. 22, 738–751 (2020).CAS 
PubMed 

Google Scholar 
28.Garcia-Chaves, M. C., Cottrell, M. T., Kirchman, D. L., Ruiz-González, C. & del Giorgio, P. A. Single-cell activity of freshwater aerobic anoxygenic phototrophic bacteria and their contribution to biomass production. Isme J. 10, 1579–1588 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
29.Jürgens, K. & Matz, C. Predation as a shaping force for the phenotypic and genotypic composition of planktonic bacteria. Antonie van Leeuwenhoek. Int. J. Gen. Mol. Microbiol. 81, 413–434 (2002).
Google Scholar 
30.Rautio, M. & Vincent, W. F. Benthic and pelagic food resources for–zooplankton in shallow high-latitude lakes and ponds. Freshw. Biol. 51, 1038–1052 (2006).CAS 

Google Scholar 
31.Rissanen, A. J. et al. Gammaproteobacterial methanotrophs dominate methanotrophy in aerobic and anaerobic layers of boreal lake waters. Aquat. Microb. Ecol. 81, 257–276 (2018).
Google Scholar 
32.Zimmermann, M. et al. Microbial methane oxidation efficiency and robustness during lake overturn. Limnol. Oceanogr. Lett. 6, 320–328 (2021).CAS 

Google Scholar 
33.Puri, A. W. et al. Genetic tools for the industrially promising methanotroph Methylomicrobium buryatense. Appl. Environ. Microbiol. 81, 1775–1781 (2015).PubMed 
PubMed Central 
ADS 

Google Scholar 
34.Strong, P. J., Kalyuzhnaya, M., Silverman, J. & Clarke, W. P. A methanotroph-based biorefinery: potential scenarios for generating multiple products from a single fermentation. Bioresour. Technol. 215, 314–323 (2016).CAS 
PubMed 

Google Scholar 
35.Oswald, K. et al. Aerobic gammaproteobacterial methanotrophs mitigate methane emissions from oxic and anoxic lake waters. Limnol. Oceanogr. 61, S101–S118 (2016).
Google Scholar 
36.Smith, E. M. & Prairie, Y. T. Bacterial metabolism and growth efficiency in lakes: the importance of phosphorus availability. Limnol. Oceanogr. 49, 137–147 (2004).CAS 
ADS 

Google Scholar 
37.Del Giorgio, P. A., Cole, J. J., Caraco, N. F. & Peters, R. H. Linking planktonic biomass and metabolism to net gas fluxes in northern temperate lakes. Ecology 80, 1422–1431 (1999).
Google Scholar 
38.Kellerman, A. M., Dittmar, T., Kothawala, D. N. & Tranvik, L. J. Chemodiversity of dissolved organic matter in lakes driven by climate and hydrology. Nat. Commun. 5, 3804 (2014).CAS 
PubMed 
ADS 

Google Scholar 
39.Sun, L., Perdue, E. M., Meyer, J. L. & Weis, J. Use of elemental composition to predict bioavailability of dissolved organic matter in a Georgia river. Limnol. Oceanogr. 42, 714–721 (1997).CAS 
ADS 

Google Scholar 
40.Kellerman, A. M., Kothawala, D. N., Dittmar, T. & Tranvik, L. J. Persistence of dissolved organic matter in lakes related to its molecular characteristics. Nat. Geosci. 8, 454–457 (2015).CAS 
ADS 

Google Scholar 
41.Guillemette, F. & del Giorgio, P. A. Reconstructing the various facets of dissolved organic carbon bioavailability in freshwater ecosystems. Limnol. Oceanogr. 56, 734–748 (2011).CAS 
ADS 

Google Scholar 
42.Logue, J. B. et al. Experimental insights into the importance of aquatic bacterial community composition to the degradation of dissolved organic matter. ISME J. 10, 533–545 (2016).CAS 
PubMed 

Google Scholar 
43.Salcher, M. M., Posch, T. & Pernthaler, J. In situ substrate preferences of abundant bacterioplankton populations in a prealpine freshwater lake. ISME J. 7, 896–907 (2013).CAS 
PubMed 

Google Scholar 
44.Sobek, S., Tranvik, L. J., Prairie, Y., Kortelainen, P. & Cole, J. J. Patterns and regulation of dissolved organic carbon: an analysis of 7,500 widely distributed lakes. Limnol. Oceanogr. 52, 1208–1219 (2007).CAS 
ADS 

Google Scholar 
45.Kalyuzhnaya, M. G. et al. Highly efficient methane biocatalysis revealed in a methanotrophic bacterium. Nat. Commun. 4, 2785 (2013).CAS 
PubMed 
ADS 

Google Scholar 
46.Oshkin, I. Y. et al. Methane-fed microbial microcosms show differential community dynamics and pinpoint taxa involved in communal response. ISME J. 9, 1119–1129 (2014).PubMed 
PubMed Central 

Google Scholar 
47.Chistoserdova, L. & Kalyuzhnaya, M. G. Current trends in methylotrophy. Trends Microbiol 26, 703–714 (2018).CAS 
PubMed 

Google Scholar 
48.Martinez-Cruz, K. et al. Anaerobic oxidation of methane by aerobic methanotrophs in sub-Arctic lake sediments. Sci. Total Environ. 607–608, 23–31 (2017).PubMed 
ADS 

Google Scholar 
49.Samad, M. S. & Bertilsson, S. Seasonal variation in abundance and diversity of bacterial methanotrophs in five temperate lakes. Front. Microbiol. 8, 1–12 (2017).
Google Scholar 
50.Ricão Canelhas, M., Denfeld, B. A., Weyhenmeyer, G. A., Bastviken, D. & Bertilsson, S. Methane oxidation at the water-ice interface of an ice-covered lake. Limnol. Oceanogr. 61, S78–S90 (2016).ADS 

Google Scholar 
51.Houser, J. N. Water color affects the stratification, surface temperature, heat content, and mean epilimnetic irradiance of small lakes. Can. J. Fish. Aquat. Sci. 63, 2447–2455 (2006).
Google Scholar 
52.Caplanne, S. & Laurion, I. Effect of chromophoric dissolved organic matter on epilimnetic stratification in lakes. Aquat. Sci. 70, 123–133 (2008).CAS 

Google Scholar 
53.Oswald, K. et al. Light-dependent aerobic methane oxidation reduces methane emissions from seasonally stratified lakes. PLoS ONE 10, e0132574 (2015).PubMed 
PubMed Central 

Google Scholar 
54.Savvichev, A. S. et al. Light-dependent methane oxidation is the major process of the methane cycle in the water column of the Bol’shie Khruslomeny Polar Lake. Microbiology 88, 370–374 (2019).CAS 

Google Scholar 
55.Baines, S. B. & Pace, M. L. The production of dissolved organic matter by phytoplankton and its importance to bacteria: patterns across marine and freshwater systems. Limnol. Oceanogr. 36, 1078–1090 (1991).ADS 

Google Scholar 
56.Cole, J. J., Likens, G. E. & Strayer, D. L. Photosynthetically produced dissolved organic carbon: an important carbon source for planktonic bacteria. Limnol. Oceanogr. 27, 1080–1090 (1982).CAS 
ADS 

Google Scholar 
57.Dumestre, J. et al. Influence of light intensity on methanotrophic bacterial activity in Petit Saut reservoir, French Guiana. Appl. Environ. Microbiol. 65, 534–539 (1999).CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
58.Murase, J. & Sugimoto, A. Inhibitory effect of light on methane oxidation in the pelagic water column of a mesotrophic lake (Lake Biwa, Japan). Limnol. Oceanogr. 50, 1339–1343 (2005).CAS 
ADS 

Google Scholar 
59.Moran, M. A. & Hodson, R. E. Bacterial production on humic and nonhumic components of dissolved organic carbon. Limnol. Oceanogr. 35, 1744–1756 (1990).CAS 
ADS 

Google Scholar 
60.Azam, F. et al. The ecological role of water-column microbes in the sea. Mar. Ecol. Prog. Ser. 10, 257–263 (1983).ADS 

Google Scholar 
61.Roulet, N. & Moore, T. R. Browning the waters. Nature 444, 283–284 (2006).CAS 
PubMed 
ADS 

Google Scholar 
62.Weyhenmeyer, G. A., Prairie, Y. T. & Tranvik, L. J. Browning of boreal freshwaters coupled to carbon-iron interactions along the aquatic continuum. PLoS ONE 9, e88104 (2014).PubMed 
PubMed Central 
ADS 

Google Scholar 
63.O’Reilly, C. M. et al. Rapid and highly variable warming of lake surface waters around the globe. Geophys. Res. Lett. 42, 10773–10781 (2015).ADS 

Google Scholar 
64.Yamamoto, S., Alcauskas, J. B. & Crozier, T. E. Solubility of methane in distilled water and seawater. J. Chem. Eng. Data 21, 78–80 (1976).CAS 

Google Scholar 
65.Cantin, A., Beisner, B. E., Gunn, J. M., Prairie, Y. T. & Winter, J. G. Effects of thermocline deepening on lake plankton communities. Can. J. Fish. Aquat. Sci. 68, 260–276 (2011).
Google Scholar 
66.Smith, D. C. & Azam, F. A simple, economical method for measuring bacterial protein synthesis rates in seawater using 3H-leucine. Mar. Microb. Food Webs 6, 107–114 (1992).
Google Scholar 
67.Del Giorgio, P. A., Pace, M. L. & Fischer, D. Relationship of bacterial growth efficiency to spatial variation in bacterial activity in the Hudson River. Aquat. Microb. Ecol. 45, 55–67 (2006).
Google Scholar 
68.Eller, G., Stubner, S. & Frenzel, P. Group-specific 16S rRNA targeted probes for the detection of type I and type II methanotrophs by fluorescence in situ hybridisation. FEMS Microbiol. Lett. 198, 91–97 (2001).CAS 
PubMed 

Google Scholar 
69.Zeder, M. ACME tool3. (2014).70.Callieri, C. et al. Bacteria, Archaea, and Crenarchaeota in the epilimnion and hypolimnion of a deep holo-oligomictic lake. Appl. Environ. Microbiol. 75, 7298–7300 (2009).CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
71.Lew, S. & Glińska-Lewczuk, K. Environmental controls on the abundance of methanotrophs and methanogens in peat bog lakes. Sci. Total Environ. 645, 1201–1211 (2018).CAS 
PubMed 
ADS 

Google Scholar 
72.Fagerbakke, K. M., Heldal, M. & Norland, S. Content of carbon, nitrogen, oxygen, sulfur and phosphorus in native aquatic and cultured bacteria. Aquat. Microb. Ecol. 10, 15–27 (1996).
Google Scholar 
73.Read J. S. et al. Derivation of lake mixing and stratification indices from high-resolution lake buoy data. Environmental Modelling and Software. 26, 1325–1336 (2011).
Google Scholar 
74.R Core Team. R: a language and environment for statistical computing. (2019). https://www.R-project.org/.75.RStudio Team. RStudio: Integrated Development for R. RStudio, Inc., Boston, MAURL. (2018). http://www.rstudio.com/.76.Wickham, H., François, R., Henry, L. & Müller, K. dplyr: A Grammar of Data Manipulation. R package version 0.8.3. (2019). https://cran.r-project.org/package=dplyr.77.Sievert, C. Interactive Web-Based Data Visualization with R, plotly, and shiny. Chapman and Hall/CRC Florida. (2020).78.Wickham, H. ggplot2: Elegant Graphics for Data Analysis. (Springer-Verlag, New York, 2016). https://ggplot2.tidyverse.org.79.Wilke, C.O. cowplot: Streamlined Plot Theme and Plot Annotations for ‘ggplot2’. R package version 1.0.0. (2019). https://CRAN.R-project.org/package=cowplot.80.Auguie, B. gridExtra: Miscellaneous Functions for ‘Grid’ Graphics. R package version 2.3. (2017). https://CRAN.R-project.org/package=gridExtra.81.Reis, P. C. J., Thottathil, S. D. & Prairie, Y. T. Dataset: the role of methanotrophy in the microbial carbon metabolism of temperate lakes. (1.0.0) [Data set]. Zenodo. (2021). https://doi.org/10.5281/zenodo.5737277. More