Ecology

Establishing size-weight relationships for heavily targeted coral species is an important first step towards informing sustainable harvest limits19. Placing coral harvests into an ecological context is a core requirement for implementing a defensible stock assessment strategy, and this need is particularly critical given escalating disturbances and widespread reports of coral loss7,17,25. Using these relationships, managers can now easily sample and calculate biomass per unit area. It is important to point out that all sites sampled in our study represent fished locations, and there is no information available to test whether standing biomass has declined due to sustained coral harvesting at these locations. While these data may now provide a critical baseline for assessing the future effects of ongoing fishing, it is also important to sample at comparable locations where fishing is not permitted or has not occurred (where possible), to test for potential effects of recent and historical harvesting.Biomass per unit area data presented herein highlights the highly patchy abundance and biomass of targeted coral species14, which is evident based on the often vastly different mean and median values (Table 2). Examining biomass per unit area estimates for C. jardinei for example, which returned some of the highest biomass estimates, the 33.75 kg·m−2 maximum estimate from a transect stands as an extreme outlier, with 12 of the 16 other transects being below 0.2 kg·m−2. This indicates the challenges of managing species that occur in patchily distributed concentrations, particularly in a management area the size of the QCF. It is also important to note, these estimates are generated only on transects where the target species occurred, and therefore, should technically not be considered as an overall estimate of standing biomass. While the estimation of size-weight relationships is a step towards a standing biomass estimate, many challenges remain in terms of sampling or reliably predicting the occurrence of these patchily distributed species. Bruckner et al.14 attempted to overcome this management challenge in a major coral fishery region of Indonesia by categorising and sampling corals (in terms of coral numbers) in defined habitat types, and then extrapolating to estimated habitat area based on visual surveys and available data. This approach, utilising size-weight relationship derived biomass per unit area estimates (instead of coral numbers), may be a viable method for the QCF, however much more information is needed to understand the habitat associations (e.g., nearshore to offshore), and environmental gradients that influence the size and abundance of individual corals. Fundamentally, it is also clear that much more data is required to effectively assess the standing biomass of aquarium corals in the very large area of operation available to Australian coral fisheries.These corals are found in a range of environments, and it is important to consider available information on life history if attempting to use coral size-weight relationships to inform management strategies via standing biomass estimation. All corals in this study can be found as free living corals (at least post-settlement) in soft-sediment, inter-reefal habitats, from which they are typically harvested by commercial collectors19. However, only four of the 6 species are colonial (C. jardinei, D. axifuga, E. glabrescens, M. lordhowensis) while the remaining two species (H. cf. australis and T. geoffroyi) are more typically monostomatous or solitary. As indicated in previous work24, if larger colonial corals were to be fragmented during harvesting instead of removed entirely, fishery impacts would likely be lessened24. Given the power relationship between coral maximum diameter and weight, larger corals contribute disproportionately to the total available biomass of each species in a given area. The potential environmental benefit of leaving larger colonies (at least partially) intact is not limited to impacts on standing biomass, as this practice would likely be demographically beneficial given the greater reproductive potential (i.e., fecundity) of larger colonies, which also do not need to overcome barriers to replenishment of populations associated with new recruits (i.e., high mortality during and post-settlement26). This conclusion was drawn largely from data on branching taxa (e.g., Acropora), which are relatively resilient to fragmentation and commonly undergo fragmentation as a result of natural processes27,28,29. D. axifuga can be considered to exhibit a relatively similar branching growth form, however, the growth form of E. glabrescens and C. jardinei changes with size, moving from small discrete polyps to large phaceloid and flabello-meandroid colonies, respectively19. While larger colonies of E. glabrescens and C. jardinei may be relatively resilient to harvesting via fragmentation, the same may not be true for smaller colonies, or species with massive growth forms such as M. lordhowensis. Typically, for each species, the average reported weight was quite low, coinciding with the lower end of the sampled maximum diameter range. For colonial species, the harvested smaller maximum diameters (if fragments) are ideal from an ecological perspective as this will have the least impact possible on standing biomass, and may also leave a potentially mature breeding colony intact. Ultimately, in light of these considerations, the development of uniform and standardised industry-wide harvest guidelines to balance economic and ecological outcomes may be necessary. The development of these guidelines would require consultation with commercial harvesters, as well as considerable additional work in measuring ecological impacts and better understanding the cost of these impacts from an economic perspective. Conversely, if whole colonies are collected, which is necessarily the case for solitary species such as H. cf. australis and T. geoffroyi (and potentially smaller colonies of other species such as E. glabrescens and C. jardinei); smaller colonies may be collected before they reach sexual maturity, hindering their ability to contribute to population replenishment. Therefore, collection of small fragments should be encouraged for colonial species; while for monostomatous species where this is not possible, introduction of a minimum harvest size based on sexual maturity should be considered.Additionally, the need for further consideration of the selectivity of ornamental coral harvest fisheries3,4,30 when assessing standing biomass is evident. Due to various desirable traits, the majority of available biomass may not be targeted by collectors. As emphasised in this study, the focus on smaller corals is indicative of the trend towards collection of most of these species at the lower portion of their size range, at least compared to some of the maximum sizes recorded on transects (e.g., see Tables 1 and 2, section b). However, it is also important to consider that transects were conducted in areas subject to commercial collection and are likely to skew results and prevent clear conclusions relating to size selectivity. Sampling of unfished populations (i.e., any residing outside of permitted fishing zones) and/or spatial and temporal matching of catch data and transect data across a larger sample of operators will be required to properly address industry size selectivity trends. For instance, only 17.5% of C. jardinei corals measured on transects fell within the diameter range represented by data obtained from collectors, with 81.9% of corals measured on transects exceeding this range. If it is viable to collect fragments from larger colonies (which does appear to be the case for some corals such as C. jardinei), then a larger proportion of standing biomass outside of this size range could be targeted by fishers. As an additional consideration, only desirable colour morphs of these corals will be harvested, and due to lack of appropriate data, the prevalence of these morphs remains unclear. H. cf. australis and M. lordhowensis for example often occur in brown colour morphs, which are far less popular in markets where certain aesthetic qualities (e.g., specific, eye-catching colours or combinations of colours) are desired, such as the ornamental aquarium industry. Even without delving into further considerations such as heritability of phenotypic traits, management conclusions drawn from standing biomass estimates may be ineffective in the absence of efforts to account for selectivity in this fishery.The relationship between size and weight was found to differ between all corals, with the exception of C. jardinei and E. glabrescens. There can be some moderate similarity in skeletal structure between these two species, particularly between small colonies, reflecting the similar maximum diameter range of sampling in the current study. Subsequently, inherent physiological constraints may be imposed on corals that prevent the maintenance of growth rates between corals of smaller and larger sizes, for example, as the surface area to volume ratio declines with growth31. In the current study, all corals, with the exception of C. jardinei, showed evidence of allometric growth, as exhibited by an estimated exponent value different to 3. Sample size for C. jardinei was greatly limited, as this species typically forms extensive beds, and are rarely brought to facilities as whole colonies. Therefore, the lack of evidence for allometric growth may reflect higher error for the species coefficient parameter due to the comparatively small sample size for this species. This suggests that mass would not increase consistently with changes in colony size in 3 dimensions31, which seems likely considering the change in exhibited form described for E. glabrescens and C. jardinei previously. In the current context, this indicates that the estimated ‘a’ and ‘b’ constants are likely to vary as the sample range increases, reflecting the changes in the size-weight relationship between smaller and larger samples of these species. Therefore, ideally, these models should incorporate data that reflect the maximum diameter range of the species in the region of application to allow increased accuracy of biomass estimation. To achieve this will require additional fishery-independent sampling, as large colonies are rarely collected whole, though may be collected as fragments depending on the species. Sampling may be challenging for some species given the difficulty of physically collecting and replacing large whole colonies, particularly for inter-reefal species such as M. lordhowensis, which can occur in deep, soft sediment habitat, subject to strong currents. Importantly, obtaining ex situ or in situ growth rate data should be considered a priority for the management of heavily targeted species. This data is likely to be another necessary component (in conjunction with size-weight relationships) of any stock assessment model developed for LPS corals, and may also eliminate the need to collect large sample colonies to improve estimated size-weight relationships.The disproportionate focus on smaller corals (i.e., corals in the current study averaged between 4.28 and 11.48 cm in maximum diameter) is likely to lead to an underestimation of weight in corals at greater diameters when used as inputs for size-weight models. This may explain the apparent minor underestimation observed in some species (e.g., M. micromussa, T. geoffroyi). In the current context, this represents an added level of conservatism with estimates obtained from these equations. While the relationship between size and weight was particularly strong for some species, (mainly D. axifuga and T. geoffroyi), for other species, such as M. lordhowensis, growth curves tended towards underestimation at larger diameter values. As the mass of a coral is reflective of the amount of carbonate skeleton that has been deposited32, the coral skeleton may increase disproportionately to coral diameter if or when corals start growing vertically. For example, in massive corals such as M. lordhowensis, vertical growth (i.e., skeletal thickening) is often very negligible among smaller colonies, with thickening of the coral skeleton only becoming apparent once the coral has reached a threshold size in terms of horizontal planar area. Additional fisheries-independent sampling outside of the relatively narrow size range of harvested colonies will be required to address this source of error in future applications. Ecological context in the form of fishery independent data on stock size and structure is essential for effective management, especially in ensuring that exploitation levels are sustainable and appropriate limits are in place. Coral harvest fisheries offer managers an ecologically and biologically unique challenge, as the implementation of standard fisheries management techniques and frameworks is hampered by their coloniality and unique biology, as well as a general lack of relevant data for assessing standing biomass and population turnover, not to mention the evolving taxonomy of scleractinian corals33. Similarly, fishery-related management challenges such as the extreme selectivity in terms of targeted size-ranges and colour-morphs, plus the potentially vast difference in the impact of various collection strategies (i.e., whole colony collection vs fragmentation during collection) also complicates the application of typical fisheries stock assessment frameworks. The relationships and equations established in the current work offer an important first step for coral fisheries globally by laying the groundwork for a defensible, ecologically sound management strategy through estimation of standing biomass, thus bridging the gap between weight-based quotas and potential environmental impacts of ongoing harvesting. It is important to note that the species selected for the current work do not represent the extent of heavily targeted LPS corals. For example, Fimbriaphyllia ancora (Veron & Pichon, 1980), Fimbriaphyllia paraancora (Veron, 1990), Cycloseris cyclolites (Lamark, 1815), and Acanthophyllia deshayesiana (Michelin, 1850) are examples of other heavily targeted corals of potential environmental concern19, and management would also benefit from the estimation of size-weight relationships for these species. Moving forward, the next challenge for the coral harvest fisheries will be to comprehensively document and track the standing biomass of heavily targeted and highly vulnerable coral stocks, explicitly accounting for fisheries effects and also non-fisheries threats, especially global climate change. More

Koonin, E. V. & Dolja, V. V. Virus world as an evolutionary network of viruses and capsidless selfish elements. Microbiol. Mol. Biol. Rev. 78, 278–303 (2014).Article 
CAS 

Google Scholar 
Pritham, E. J., Putliwala, T. & Feschotte, C. Mavericks, a novel class of giant transposable elements widespread in eukaryotes and related to DNA viruses. Gene 390, 3–17 (2007).Article 
CAS 

Google Scholar 
Kapitonov, V. V. & Jurka, J. Self-synthesizing DNA transposons in eukaryotes. Proc. Natl Acad. Sci. USA 103, 4540–4545 (2006).Article 
CAS 

Google Scholar 
Krupovic, M. & Koonin, E. V. Polintons: a hotbed of eukaryotic virus, transposon and plasmid evolution. Nat. Rev. Microbiol. 13, 105–115 (2015).Article 
CAS 

Google Scholar 
Koonin, E. V., Krupovic, M. & Yutin, N. Evolution of double-stranded DNA viruses of eukaryotes: from bacteriophages to transposons to giant viruses. Ann. N. Y. Acad. Sci. 1341, 10–24 (2015).Article 
CAS 

Google Scholar 
Yutin, N., Raoult, D. & Koonin, E. V. Virophages, polintons, and transpovirons: a complex evolutionary network of diverse selfish genetic elements with different reproduction strategies. Virol. J. 10, 158 (2013).Article 
CAS 

Google Scholar 
Krupovic, M., Bamford, D. H. & Koonin, E. V. Conservation of major and minor jelly-roll capsid proteins in Polinton (Maverick) transposons suggests that they are bona fide viruses. Biol. Direct 9, 6 (2014).Article 

Google Scholar 
Yutin, N., Shevchenko, S., Kapitonov, V., Krupovic, M. & Koonin, E. V. A novel group of diverse Polinton-like viruses discovered by metagenome analysis. BMC Biol. 13, 95 (2015).Article 

Google Scholar 
Bellas, C. M. & Sommaruga, R. Polinton-like viruses are abundant in aquatic ecosystems. Microbiome 9, 13 (2021).Article 
CAS 

Google Scholar 
Pagarete, A., Grébert, T., Stepanova, O., Sandaa, R.-A. & Bratbak, G. Tsv-N1: a novel DNA algal virus that infects Tetraselmis striata. Viruses 7, 3937–3953 (2015).Article 
CAS 

Google Scholar 
Bekliz, M., Colson, P. & La Scola, B. The expanding family of virophages. Viruses 8, 317 (2016).Article 

Google Scholar 
Fischer, M. G. The virophage family Lavidaviridae. Curr. Issues Mol. Biol. https://doi.org/10.21775/cimb.040.001 (2021).Desnues, C. et al. Provirophages and transpovirons as the diverse mobilome of giant viruses. Proc. Natl Acad. Sci. USA 109, 18078–18083 (2012).Article 
CAS 

Google Scholar 
Campos, R. K. et al. Samba virus: a novel mimivirus from a giant rain forest, the Brazilian Amazon. Virol. J. 11, 95 (2014).Article 

Google Scholar 
Gaia, M. et al. Broad spectrum of mimiviridae virophage allows its isolation using a mimivirus reporter. PLoS ONE 8, e61912 (2013).Article 
CAS 

Google Scholar 
Hackl, T., Duponchel, S., Barenhoff, K., Weinmann, A. & Fischer, M. G. Virophages and retrotransposons colonize the genomes of a heterotrophic flagellate. eLife 10, e72674 (2021).Article 
CAS 

Google Scholar 
Yau, S. et al. Virophage control of Antarctic algal host-virus dynamics. Proc. Natl Acad. Sci. USA 108, 6163–6168 (2011).Article 
CAS 

Google Scholar 
Gong, C. et al. Novel virophages discovered in a freshwater lake in China. Front. Microbiol. 7, 5 (2016).Article 

Google Scholar 
Zhou, J. et al. Three novel virophage genomes discovered from Yellowstone Lake metagenomes. J. Virol. 89, 1278–1285 (2014).Article 

Google Scholar 
Yutin, N., Kapitonov, V. V. & Koonin, E. V. A new family of hybrid virophages from an animal gut metagenome. Biol. Direct 10, 19 (2015).Article 

Google Scholar 
Stough, J. M. A. et al. Genome and environmental activity of a Chrysochromulina parva virus and its virophages. Front. Microbiol. 10, 703 (2019).Article 

Google Scholar 
La Scola, B. et al. The virophage as a unique parasite of the giant mimivirus. Nature 455, 100–104 (2008).Article 

Google Scholar 
Fischer, M. G. & Suttle, C. A. A virophage at the origin of large DNA transposons. Science 332, 231–234 (2011).Article 
CAS 

Google Scholar 
Gaia, M. et al. Zamilon, a novel virophage with Mimiviridae host specificity. PLoS ONE 9, e94923 (2014).Article 

Google Scholar 
Mougari, S. et al. Guarani virophage, a new Sputnik-like isolate from a Brazilian lake. Front. Microbiol. 10, 1003 (2019).Article 

Google Scholar 
Sheng, Y., Wu, Z., Xu, S. & Wang, Y. Isolation and identification of a large green alga virus (Chlorella Virus XW01) of Mimiviridae and its virophage (Chlorella Virus Virophage SW01) by using unicellular green algal cultures. J. Virol. 96, e02114–e02121 (2022).Article 

Google Scholar 
Baudoux, A. C. & Brussaard, C. P. D. Characterization of different viruses infecting the marine harmful algal bloom species Phaeocystis globosa. Virology 341, 80–90 (2005).Article 
CAS 

Google Scholar 
Santini, S. et al. Genome of Phaeocystis globosa virus PgV-16T highlights the common ancestry of the largest known DNA viruses infecting eukaryotes. Proc. Natl Acad. Sci. USA 110, 10800–10805 (2013).Article 
CAS 

Google Scholar 
Tarutani, K., Nagasaki, K. & Yamaguchi, M. Virus adsorption process determines virus susceptibility in Heterosigma akashiwo (Raphidophyceae). Aquat. Microb. Ecol. 42, 209–213 (2006).Article 

Google Scholar 
Gann, E. R., Gainer, P. J., Reynolds, T. B. & Wilhelm, S. W. Influence of light on the infection of Aureococcus anophagefferens CCMP 1984 by a ‘giant virus’. PLoS ONE 15, e0226758 (2020).Article 
CAS 

Google Scholar 
Van Etten, J. L., Burbank, D. E., Xia, Y. & Meints, R. H. Growth cycle of a virus, PBCV-1, that infects Chlorella-like algae. Virology 126, 117–125 (1983).Article 

Google Scholar 
Boyer, M. et al. Mimivirus shows dramatic genome reduction after intraamoebal culture. Proc. Natl Acad. Sci. USA 108, 10296–10301 (2011).Article 
CAS 

Google Scholar 
Desnues, C. & Raoult, D. Inside the lifestyle of the virophage. Intervirology 53, 293–303 (2010).Article 
CAS 

Google Scholar 
Sobhy, H., Scola, B. L., Pagnier, I., Raoult, D. & Colson, P. Identification of giant Mimivirus protein functions using RNA interference. Front. Microbiol. 6, 345 (2015).Article 

Google Scholar 
Fischer, M. G. & Hackl, T. Host genome integration and giant virus-induced reactivation of the virophage mavirus. Nature 540, 288–291 (2016).Article 
CAS 

Google Scholar 
Wodarz, D. Evolutionary dynamics of giant viruses and their virophages. Ecol. Evol. 3, 2103–2115 (2013).Article 

Google Scholar 
Farr, G. A., Zhang, L. & Tattersall, P. Parvoviral virions deploy a capsid-tethered lipolytic enzyme to breach the endosomal membrane during cell entry. Proc. Natl Acad. Sci. USA 102, 17148–17153 (2005).Article 
CAS 

Google Scholar 
Suhre, K., Audic, S. & Claverie, J.-M. Mimivirus gene promoters exhibit an unprecedented conservation among all eukaryotes. Proc. Natl Acad. Sci. USA 102, 14689–14693 (2005).Article 
CAS 

Google Scholar 
Legendre, M. et al. mRNA deep sequencing reveals 75 new genes and a complex transcriptional landscape in Mimivirus. Genome Res. 20, 664–674 (2010).Article 
CAS 

Google Scholar 
Smith, D. R., Arrigo, K. R., Alderkamp, A.-C. & Allen, A. E. Massive difference in synonymous substitution rates among mitochondrial, plastid, and nuclear genes of Phaeocystis algae. Mol. Phylogenet. Evol. 71, 36–40 (2014).Article 
CAS 

Google Scholar 
Krupovic, M., Kuhn, J. H. & Fischer, M. G. A classification system for virophages and satellite viruses. Arch. Virol. 161, 233–247 (2016).Article 
CAS 

Google Scholar 
Suplatov, D. A., Besenmatter, W., Svedas, V. K. & Svendsen, A. Bioinformatic analysis of alpha/beta-hydrolase fold enzymes reveals subfamily-specific positions responsible for discrimination of amidase and lipase activities. Protein Eng. Des. Sel. 25, 689–697 (2012).Article 
CAS 

Google Scholar 
Burt, A. & Koufopanou, V. Homing endonuclease genes: the rise and fall and rise again of a selfish element. Curr. Opin. Genet. Dev. 14, 609–615 (2004).Article 
CAS 

Google Scholar 
Sullivan, M. B. DNA extraction of cesium chloride-purified viruses using wizard prep columns. Protocols https://doi.org/10.17504/protocols.io.c26yhd (2016).González-Domínguez, J. & Schmidt, B. ParDRe: faster parallel duplicated reads removal tool for sequencing studies. Bioinformatics 32, 1562–1564 (2016).Article 

Google Scholar 
Guillard, R. R. L. Culture of phytoplankton for feeding marine invertebrates. In Culture of Marine Invertebrate Animals: Proceedings—1st Conference on Culture of Marine Invertebrate Animals Greenport (eds Smith, W. L., & Chanley, M. H.) 29– 60 (Springer, 1975).Cottrell, M. & Suttle, C. Wide-spread occurrence and clonal variation in viruses which cause lysis of a cosmopolitan, eukaryotic marine phytoplankter Micromonas pusilla. Mar. Ecol. Prog. Ser. 78, 1–9 (1991).Article 

Google Scholar 
Krueger, F., James, F., Ewels, P., Afyounian, E. & Schuster-Boeckler, B. FelixKrueger/TrimGalore: v0.6.7 – DOI via Zenodo. https://doi.org/10.5281/zenodo.5127899 (2021).Bankevich, A. et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J. Comput. Biol. 19, 455–477 (2012).Article 
CAS 

Google Scholar 
Schneider, C. A., Rasband, W. S. & Eliceiri, K. W. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671–675 (2012).Article 
CAS 

Google Scholar 
Patel, A. et al. Virus and prokaryote enumeration from planktonic aquatic environments by epifluorescence microscopy with SYBR Green I. Nat. Protoc. 2, 269–276 (2007).Article 
CAS 

Google Scholar 
Bolte, S. & Cordelières, F. P. A guided tour into subcellular colocalization analysis in light microscopy. J. Microsc. 224, 213–232 (2006).Article 
CAS 

Google Scholar 
Brussaard, C. P. D. Optimization of procedures for counting viruses by flow cytometry. Appl. Environ. Microbiol. 70, 1506–1513 (2004).Article 
CAS 

Google Scholar 
Kirzner, S., Barak, E. & Lindell, D. Variability in progeny production and virulence of cyanophages determined at the single-cell level. Environ. Microbiol. Rep. 8, 605–613 (2016).Article 

Google Scholar 
Ziv, I. et al. A perturbed ubiquitin landscape distinguishes between ubiquitin in trafficking and in proteolysis. Mol. Cell. Proteomics 10, M111.009753 (2011).HaileMariam, M. et al. S-Trap, an ultrafast sample-preparation approach for shotgun proteomics. J. Proteome Res. 17, 2917–2924 (2018).Article 
CAS 

Google Scholar 
Rappsilber, J., Mann, M. & Ishihama, Y. Protocol for micro-purification, enrichment, pre-fractionation and storage of peptides for proteomics using StageTips. Nat. Protoc. 2, 1896–1906 (2007).Article 
CAS 

Google Scholar 
Kong, A. T., Leprevost, F. V., Avtonomov, D. M., Mellacheruvu, D. & Nesvizhskii, A. I. MSFragger: ultrafast and comprehensive peptide identification in mass spectrometry-based proteomics. Nat. Methods 14, 513–520 (2017).Article 
CAS 

Google Scholar 
Li, W. & Godzik, A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22, 1658–1659 (2006).Article 
CAS 

Google Scholar 
Lechner, M. et al. Proteinortho: detection of (Co-)orthologs in large-scale analysis. BMC Bioinformatics 12, 124 (2011).Article 

Google Scholar 
Buchfink, B., Reuter, K. & Drost, H.-G. Sensitive protein alignments at tree-of-life scale using DIAMOND. Nat. Methods 18, 366–368 (2021).Article 
CAS 

Google Scholar 
Cox, J. et al. Accurate proteome-wide label-free quantification by delayed normalization and maximal peptide ratio extraction termed MaxLFQ. Mol. Cell. Proteomics 13, 2513–2526 (2014).Article 
CAS 

Google Scholar 
Mirdita, M. et al. ColabFold: making protein folding accessible to all. Nat. Methods 19, 679–682 (2022).Jumper, J. et al. Highly accurate protein structure prediction with AlphaFold. Nature 596, 583–589 (2021).Article 
CAS 

Google Scholar 
O’Connell, J. et al. NxTrim: optimized trimming of Illumina mate pair reads. Bioinformatics 31, 2035–2037 (2015).Article 

Google Scholar 
Li, D., Liu, C.-M., Luo, R., Sadakane, K. & Lam, T.-W. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics 31, 1674–1676 (2015).Article 
CAS 

Google Scholar 
Luo, R. et al. SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. GigaScience 1, 2047-217X-1–18 (2012).Chevreux, B., Wetter, T. & Suhai, S. Genome sequence assembly using trace signals and additional sequence information. In Proc. German Conference on Bioinformatics 45–56 (Fachgruppe Bioinformatik, 1999).Deng, Z. & Delwart, E. ContigExtender: a new approach to improving de novo sequence assembly for viral metagenomics data. BMC Bioinformatics 22, 119 (2021).Article 
CAS 

Google Scholar 
Walker, B. J. et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS ONE 9, e112963 (2014).Article 

Google Scholar 
Eddy, S. R. Accelerated profile HMM searches. PLoS Comput. Biol. 7, e1002195 (2011).Article 
CAS 

Google Scholar 
Katoh, K., Misawa, K., Kuma, K. & Miyata, T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 30, 3059–3066 (2002).Article 
CAS 

Google Scholar 
Minh, B. Q. et al. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol. Biol. Evol. 37, 1530–1534 (2020).Article 
CAS 

Google Scholar 
Hoang, D. T., Chernomor, O., von Haeseler, A., Minh, B. Q. & Vinh, L. S. UFBoot2: improving the ultrafast bootstrap approximation. Mol. Biol. Evol. 35, 518–522 (2018).Article 
CAS 

Google Scholar 
Kozlov, A. M., Darriba, D., Flouri, T., Morel, B. & Stamatakis, A. RAxML-NG: a fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference. Bioinformatics 35, 4453–4455 (2019).Article 
CAS 

Google Scholar 
Barbera, P. et al. EPA-ng: massively parallel evolutionary placement of genetic sequences. Syst. Biol. 68, 365–369 (2019).Article 

Google Scholar 
Capella-Gutiérrez, S., Silla-Martínez, J. M. & Gabaldón, T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25, 1972–1973 (2009).Article 

Google Scholar 
Steinegger, M. & Söding, J. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nat. Biotechnol. 35, 1026–1028 (2017).Article 
CAS 

Google Scholar 
Steinegger, M. et al. HH-suite3 for fast remote homology detection and deep protein annotation. BMC Bioinformatics 20, 473 (2019).Article 

Google Scholar 
Enright, A. J., Van Dongen, S. & Ouzounis, C. A. An efficient algorithm for large-scale detection of protein families. Nucleic Acids Res. 30, 1575–1584 (2002).Article 
CAS 

Google Scholar 
Bolduc, B. et al. vConTACT: an iVirus tool to classify double-stranded DNA viruses that infect Archaea and Bacteria. PeerJ 5, e3243 (2017).Article 

Google Scholar 
Zimmermann, L. et al. A completely reimplemented MPI bioinformatics toolkit with a new HHpred server at its core. J. Mol. Biol. 430, 2237–2243 (2018).Article 
CAS 

Google Scholar 
Heger, A. & Holm, L. Rapid automatic detection and alignment of repeats in protein sequences. Proteins 41, 224–237 (2000).Article 
CAS 

Google Scholar 
Chase, E., Desnues, C. & Blanc, G. Integrated viral elements unveil the dual lifestyle of Tetraselmis spp. polinton-like viruses. Virus Evol. 8, veac068 (2022).Egge, E. S., Eikrem, W. & Edvardsen, B. Deep-branching novel lineages and high diversity of haptophytes in the Skagerrak (Norway) uncovered by 454 pyrosequencing. J. Eukaryot. Microbiol. 62, 121–140 (2015).Article 
CAS 

Google Scholar 
Hovde, B. T. et al. Chrysochromulina: genomic assessment and taxonomic diagnosis of the type species for an oleaginous algal clade. Algal Res. 37, 307–319 (2019).Article 

Google Scholar 
Andersen, R. A., Bailey, J. C., Decelle, J. & Probert, I. Phaeocystis rex sp. nov. (Phaeocystales, Prymnesiophyceae): a new solitary species that produces a multilayered scale cell covering. Eur. J. Phycol. 50, 207–222 (2015).Article 

Google Scholar 
Stepanova, O. A. Black Sea algal viruses. Russ. J. Mar. Biol. 42, 123–127 (2016).Article 

Google Scholar 
Alarcón-Schumacher, T., Guajardo-Leiva, S., Antón, J. & Díez, B. Elucidating viral communities during a phytoplankton bloom on the West Antarctic Peninsula. Front. Microbiol. 10, 1014 (2019).Article 

Google Scholar  More

DatasetsImage acquisitionFor this study, we assembled three datasets: one for training of the RootDetector Convolutional Neural Network (Training-Set), one for a performance comparison between humans and RootDetector in segmenting roots in minirhizotron images (Comparison-Set), and one for the validation of the algorithm (Validation-Set). The Training-Set contained 129 images comprised of 17 randomly selected minirhizotron images sampled in a mesocosm experiment (see “Mesocosm sampling” Section), 47 randomly selected minirhizotron images sampled in a field study (see “Field sampling” Section) as well as the 65 minirhizotron images of soy roots published by Wang et al.15. The Comparison-Set contained 25 randomly selected minirhizotron images from the field-study which all were not part of the images included in the Training- and Validation-Sets. The Validation-Set contained 10 randomly selected minirhizotron images from the same field study, which had not been used in the Training-Set. All images were recorded with 2550 ✕ 2273 pixels at 300 dpi with a CI-600 In-Situ Root Imager (CID Bio-Science Inc., Camas, WA, USA) and stored as .tiff files to reduce compression loss. For all training and evaluation purposes we used raw, unprocessed output images from the CI-600.Mesocosm samplingThe mesocosm experiment was established in 2018 on the premises of the Institute for Botany and Landscape Ecology of the University of Greifswald (Fig. S1). It features 108 heavy duty plastic buckets of 100 l each, filled to two thirds of their height with moderately decomposed sedge fen peat. Each mesocosm contained one minirhizotron (inner diameter: 64 mm, outer diameter: 70 mm, length: 650 mm) installed at a 45°angle and capped in order to avoid penetration by light. The mesocosms were planted with varying compositions of plant species that typically occur in north-east German sedge fens (Carex rostrata, Carex acutiformis, Glyceria maxima, Equisetum fluviatile, Juncus inflexus, Mentha aquatica, Acorus calamus and Lycopus europaeus). The mesocosms were subjected to three different water table regimes: stable at soil surface level, stable at 20 cm below soil surface and fluctuating between the two levels every two weeks. The minirhizotrons were scanned weekly at two levels of soil depth (0–20 cm and 15–35 cm) between April 2019 and December 2021, resulting in roughly 9500 minirhizotron images of 216 × 196 mm. Manual quantification of root length would, based on own experience, take approximately three hours per image, resulting in approximately 28,500 h of manual processing for the complete dataset. Specimens planted were identified by author Dr. Blume-Werry, however no voucher specimen were deposited. All methods were carried out in accordance with relevant institutional, national, and international guidelines and legislation.Field samplingThe field study was established as part of the Wetscapes project in 201716. The study sites were located in Mecklenburg-Vorpommern, Germany, in three of the most common wetland types of the region: alder forest, percolation fen and coastal fen (Fig. S2). For each wetland type, a pair of drained versus rewetted study sites was established. A detailed description of the study sites and the experimental setup can be found in Jurasinski et al.16. At each site, 15 minirhizotrons (same diameter as above, length: 1500 mm) were installed at 45° angle along a central boardwalk. The minirhizotrons have been scanned biweekly since April 2018, then monthly since January 2019 at two to four levels of soil depth (0–20 cm, 20–40 cm, 40–60 cm and 60–80 cm), resulting in roughly 12,000 minirhizotron images of 216 × 196 cm, i.e. an estimated 36,000 h of manual processing for the complete dataset. Permission for the study was obtained from the all field owners. Figure 1Overview of the RootDetector system. The main component is a semantic segmentation network based on the U-Net architecture. The root length is estimated by skeletonizing the segmentation output and applying the formula introduced by Kimura et al.17. During training only, a weight map puts more emphasis on fine roots.Full size imageThe CNN RootDetectorImage annotationFor the generation of training data for the CNN, human analysts manually masked all root pixels in the 74 images of the Training-Set using GIMP 2.10.12. The resulting ground truth data are binary, black-and-white images in Portable Network Graphics (.png) format, where white pixels represent root structures and black pixels represent non-root objects and soil (Fig. 2b). All training data were checked and, if required, corrected by an expert (see “Selection of participants” for definition). The Validation-Set was created in the same way but exclusively by experts.Figure 2Example of segmentation and result of skeletonization. A 1000 by 1000 pixel input image (a), the manually annotated ground truth image (b), the RootDetector estimation image (c), the combined representation image (error map, d with green indicating true positives, red indicating false positive, blue indicating false negatives), the skeletonized RootDetector estimation image (e), and the skeletonized ground truth image (f).Full size imageArchitectureRootDetector’s core consists of a Deep Neural Network (DNN) based on the U-Net image segmentation architecture[27]nd is implemented in TensorFlow and Keras frameworks18. Although U-Net was originally developed for biomedical applications, it has since been successfully applied to other domains due to its generic design.RootDetector is built up of four down-sampling blocks, four up-sampling blocks and a final output block (Fig. 1). Every block contains two 3 × 3 convolutional layers, each followed by rectified linear units (ReLU). The last output layer instead utilizes Sigmoid activation. Starting from initial 64 feature channels, this number is doubled in every down-block and the resolution is halved via 2 × 2 max-pooling. Every up-block again doubles the resolution via bilinear interpolation and a 1 × 1 convolution which halves the number of channels. Importantly, after each up-sampling step, the feature map is concatenated with the corresponding feature map from the down-sampling path. This is crucial to preserve fine spatial details.Our modifications from the original architecture include BatchNormalization19 after each convolutional layer which significantly helps to speed up the training process and zero-padding instead of cropping as suggested by Ronneberger, Fischer, & Brox20 to preserve the original image size.In addition to the root segmentation network, we trained a second network to detect foreign objects, specifically the adhesive tape that is used as a light barrier on the aboveground part of the minirhizotrons. We used the same network architecture as above and trained in a supervised fashion with the binary cross-entropy loss. During inference, the result is thresholded (predefined threshold value: 0.5) and used without post-processing.TrainingWe pre-trained RootDetector on the COCO dataset21 to generate a starting point. Although the COCO dataset contains a wide variety of image types and classes not specifically related to minirhizotron images, Majurski et al.22 showed, that for small annotation counts, transfer-learning even from unrelated datasets may improve a CNNs performance by up to 20%. We fine-tuned for our dataset with the Adam optimizer23 for 15 epochs and trained on a total of 129 images from the Training-Set (17 mesocosm images, 47 field-experiment images, 65 soy root images). To enhance the dataset size and reduce over-fitting effects, we performed a series of augmentation operations as described by Shorten & Khoshgoftaar24. In many images, relatively coarse roots ( > 3 mm) occupied a major part of the positive (white) pixel space, which might have caused RootDetector to underestimate fine root details overall. Similarly, negative space (black pixels) between tightly packed, parallel roots was often very small and might have impacted the training process to a lesser extent when compared to large areas with few or no roots (Fig. 2). To mitigate both effects, we multiplied the result of the cross-entropy loss map with a weight map which emphasizes positive–negative transitions. This weight map is generated by applying the following formula to the annotated ground truth images:$$omega left( x right) = 1 – left( {tanh left( {2tilde{x} – 1} right)} right)^{2}$$
(1)
where ω(x) is the average pixel value of the annotated weight map in a 5 × 5 neighborhood around pixel x. Ronneberger, Fischer, & Brox20 implemented a similar weight map, however with stronger emphasis on space between objects. As this requires computation of distances between two comparatively large sets of points, we adapted and simplified their formula to be computable in a single 5 × 5 convolution.For the loss function we applied a combination of cross-entropy and Dice loss 25:$${mathcal{L}} = {mathcal{L}}_{CE} + lambda {mathcal{L}}_{Dice} = – frac{1}{N}sumnolimits_{i} {wleft( {x_{i} } right)y_{i} log left( {x_{i} } right) + lambda frac{{2sumnolimits_{i} {x_{i} y_{i} } }}{{sumnolimits_{i} {x_{i}^{2} sumnolimits_{i} {y_{i}^{2} } } }}}$$
(2)

where x are the predicted pixels, y the corresponding ground truth labels, N the number of pixels in an image and λ a balancing factor which we set to 0.01. This value was derived empirically. The Dice loss is applied per-image to counteract the usually high positive-to-negative pixel imbalance. Since this may produce overly confident outputs and restrict the application of weight maps, we used a relatively low value for λ.Output and post-processingRootDetector generates two types of output. The first type of output are greyscale .png files in which white pixels represent pixels associated with root structures and black pixels represent non-root structures and soil (Fig. 2c). The advantage of .png images is their lossless ad artifact-free compression at relatively small file sizes. RootDetector further skeletonizes the output images and reduces root-structures to single-pixel representations using the skeletonize function of scikit-image v. 0.17.1 (26; Fig. 2e,f). This helps to reduce the impact of large diameter roots or root-like structures such as rhizomes in subsequent analyses and is directly comparable to estimations of root length. The second type of output is a Comma-separated values (.csv) file, with numerical values indicating the number of identified root pixels, the number of root pixels after skeletonization, the number of orthogonal and diagonal connections between pixels after skeletonization and an estimation of the physical combined length of all roots for each processed image. The latter is a metric commonly used in root research as in many species, fine roots provide most vital functions such as nutrient and water transport3. Therefore, the combined length of all roots in a given space puts an emphasis on fine roots as they typically occupy a relatively smaller fraction of the area in a 2D image compared to often much thicker coarse roots. To derive physical length estimates from skeletonized images, RootDetector counts orthogonal- and diagonal connections between pixels of skeletonized images and employs the formula proposed by Kimura et al.17 (Eq. 3).$$L = left[ {N_{d}^{2} + left( {N_{d} + N_{o} /2} right)^{2} } right]^{{1/2}} + N_{o} /2$$
(3)
where Nd is the number of diagonally connected and No the number of orthogonally connected skeleton pixels. To compute Nd we convolve the skeletonized image with two 2 × 2 binary kernels, one for top-left-to-bottom-right connections and another for bottom-left-to-top-right connections and count the number of pixels with maximum response in the convolution result. Similarly, No is computed with a 1 × 2 and a 2 × 1 convolutional kernels.Performance comparisonSelection of participantsFor the performance comparison, we selected 10 human analysts and divided them into three groups of different expertise levels in plant physiology and with the usage of digital root measuring tools. The novice group consisted of 3 ecology students (2 bachelor’s, 1 master’s) who had taken or were taking courses in plant physiology but had no prior experience with minirhizotron images or digital root measuring tools. This group represents undergraduate students producing data for a Bachelor thesis or student assistants employed to process data. The advanced group consisted of 3 ecology students (1 bachelor’s, 2 master’s) who had already taken courses in plant physiology and had at least 100 h of experience with minirhizotron images and digital root measuring tools. The expert group consisted of 4 scientists (2 PhD, 2 PhD candidates) who had extensive experience in root science and at least 250 h of experience with digital root measuring tools. All methods were carried out in accordance with relevant institutional, national, and international guidelines and legislation and informed consent was obtained from all participants.Instruction and root tracingAll three groups were instructed by showing them a 60 min live demo of an expert tracing roots in minirhizotron images, during which commonly encountered challenges and pitfalls were thoroughly discussed. Additionally, all participants were provided with a previously generated, in-depth manual containing guidelines on the identification of root structures, the correct operation of the root tracing program and examples of often encountered challenges and suggested solutions. Before working on the Comparison-Set, all participants traced roots in one smaller-size sample image and received feedback from one expert.Image preparation and root tracingBecause the minirhizotron images acquired in the field covered a variety of different substrates, roots of different plant species, variance in image quality, and because tracing roots is very time consuming, we decided to maximize the number of images by tracing roots only in small sections, in order to cover the largest number of cases possible. To do this, we placed a box of 1000 × 1000 pixels (8.47 × 8.47 cm) at a random location in each of the images in the Comparison-Set and instructed participants to trace only roots within that box. Similarly, we provided RootDetector images where the parts of the image outside the rectangle were occluded. All groups used RootSnap! 1.3.2.25 (CID Bio-Science Inc., Camas, WA, USA;27), a vector based tool to manually trace roots in each of the 25 images in the comparison set. We decided on RootSnap! due to our previous good experience with the software and its’ relative ease of use. The combined length of all roots was then exported as a csv file for each person and image and compared to RootDetector’s output of the Kimura root length.ValidationWe tested the accuracy of RootDetector on a set of 10 image segments of 1000 by 1000 pixels cropped from random locations of the 10 images of the Validation-Set. These images were annotated by a human expert without knowledge of the estimations by the algorithm and were exempted from the training process. As commonly applied in binary classification, we use the F1 score as a metric to evaluate the performance RootDetector. F1 is calculated from precision (Eq. 4) and recall (Eq. 5) and represents their harmonic mean (Eq. 6). Ranging from 0 to 1, higher values indicate high classification (segmentation) performance. As one of the 10 image sections contained no roots and thus no F1 Score was calculable, it was excluded from the validation. We calculated the F1 score for each of the nine remaining image sections and averaged the values as a metric for overall segmentation performance.$$Precision;(P) = frac{{tp}}{{tp + fp}}$$
(4)
$$Recall;(R) = frac{{tp}}{{tp + fn}}$$
(5)
$$F1 = 2*frac{{P*R}}{{P + R}}$$
(6)
where P = precision, R = recall, tp = true positives; fp = false positives, fn = false negatives.Statistical analysisWe used R Version 4.1.2 (R Core Team, 2021) for all statistical analyses and R package ggplot2 Version 3.2.128 for visualizations. Pixel identification-performance comparisons were based on least-squares fit and the Pearson method. Root length estimation-performance comparisons between groups of human analysts (novice, advanced, expert) and RootDetector were based on the respective estimates of total root length plotted over the minirhizotron images in increasing order of total root length. Linear models were calculated using the lm function for each group of analysts. To determine significant differences between the groups and the algorithm, 95% CIs as well as 83% CIs were displayed and RootDetector root length outside the 95% CI were considered significantly different from the group estimate at α = 0.0529. The groups of human analysts were considered significantly different if their 83% CIs did not overlap, as the comparison of two 83% CIs approximates an alpha level of 5%30,31.This study is approved by Ethikkommission der Universitätsmedizin Greifswald, University of Greifswald, Germany. More

Lauchlan, S. S. & Nagelkerken, I. Species range shifts along multistressor mosaics in estuarine environments under future climate. Fish Fish. 21, 32–46 (2020).Article 

Google Scholar 
Gao, G., Zhao, X., Jiang, M. & Gao, L. Impacts of marine heatwaves on algal structure and carbon sequestration in conjunction with ocean warming and acidification. Front. Mar. Sci. 8, 758651 (2021).Article 

Google Scholar 
Asch, R. G. Climate change and decadal shifts in the phenology of larval fishes in the California Current ecosystem. Proc. Natl. Acad. Sci. 112, E4065–E4074 (2015).Article 
ADS 
CAS 

Google Scholar 
Lonhart, S. I., Jeppesen, R., Beas-Luna, R., Crooks, J. A. & Lorda, J. Shifts in the distribution and abundance of coastal marine species along the eastern Pacific Ocean during marine heatwaves from 2013 to 2018. Mar. Biodivers. Rec. 12, 13 (2019).Article 

Google Scholar 
Morley, J. W. et al. Projecting shifts in thermal habitat for 686 species on the North American continental shelf. PLoS ONE 13, e0196127 (2018).Article 

Google Scholar 
Vergés, A. et al. The tropicalization of temperate marine ecosystems: Climate-mediated changes in herbivory and community phase shifts. Proc. R. Soc. B 281, 20140846 (2014).Article 

Google Scholar 
Wernberg, T. et al. Climate-driven regime shift of a temperate marine ecosystem. Science 353, 169–172 (2016).Article 
ADS 
CAS 

Google Scholar 
Cheung, W. W. L. et al. Marine high temperature extremes amplify the impacts of climate change on fish and fisheries. Sci. Adv. https://doi.org/10.1126/sciadv.abh0895 (2021).Article 

Google Scholar 
Ling, S. D., Johnson, C. R., Frusher, S. D. & Ridgway, K. R. Overfishing reduces resilience of kelp beds to climate-driven catastrophic phase shift. Proc. Natl. Acad. Sci. 106, 22341–22345 (2009).Article 
ADS 
CAS 

Google Scholar 
Pessarrodona, A. et al. Tropicalization unlocks novel trophic pathways and enhances secondary productivity in temperate reefs. Funct. Ecol. 36, 659–673 (2022).Article 

Google Scholar 
Hobday, A. J. et al. A hierarchical approach to defining marine heatwaves. Prog. Oceanogr. 141, 227–238 (2016).Article 
ADS 

Google Scholar 
Holbrook, N. J. et al. A global assessment of marine heatwaves and their drivers. Nat. Commun. 10, 2624 (2019).Article 
ADS 

Google Scholar 
Smale, D. A. et al. Marine heatwaves threaten global biodiversity and the provision of ecosystem services. Nat. Clim. Chang. 9, 306–312 (2019).Article 
ADS 

Google Scholar 
Cheung, W. W. L. & Frölicher, T. L. Marine heatwaves exacerbate climate change impacts for fisheries in the northeast Pacific. Sci. Rep. 10, 6678 (2020).Article 
ADS 
CAS 

Google Scholar 
Garrabou, J. et al. Marine heatwaves drive recurrent mass mortalities in the Mediterranean Sea. Glob. Change Biol. 28, 5708–5725 (2022).Article 
CAS 

Google Scholar 
Wernberg, T. et al. An extreme climatic event alters marine ecosystem structure in a global biodiversity hotspot. Nat. Clim. Change 3, 78–82 (2013).Article 
ADS 

Google Scholar 
Cure, K. et al. Distributional responses to marine heat waves: insights from length frequencies across the geographic range of the endemic reef fish Choerodon rubescens. Mar. Biol. 165, 1 (2018).Article 

Google Scholar 
Jacox, M. G., Tommasi, D., Alexander, M. A., Hervieux, G. & Stock, C. A. Predicting the evolution of the 2014–2016 California current system marine heatwave from an ensemble of coupled global climate forecasts. Front. Mar. Sci. https://doi.org/10.3389/fmars.2019.00497 (2019).Article 

Google Scholar 
Gentemann, C. L., Fewings, M. R. & García-Reyes, M. Satellite sea surface temperatures along the West Coast of the United States during the 2014–2016 northeast Pacific marine heat wave. Geophys. Res. Lett. 44, 312–319 (2017).Article 
ADS 

Google Scholar 
Cavanaugh, K. C., Reed, D. C., Bell, T. W., Castorani, M. C. N. & Beas-Luna, R. Spatial variability in the resistance and resilience of giant kelp in southern and baja California to a multiyear heatwave. Front. Mar. Sci. https://doi.org/10.3389/fmars.2019.00413 (2019).Article 

Google Scholar 
Cavole, L. M. et al. Biological impacts of the 2013–2015 warm-water anomaly in the Northeast Pacific: Winners, losers, and the future. Oceanography 29, 273–285 (2016).Article 

Google Scholar 
Sen Gupta, A. et al. Drivers and impacts of the most extreme marine heatwave events. Sci. Rep. 10, 19359 (2020).Article 
ADS 
CAS 

Google Scholar 
Rykaczewski, R. R. & Checkley, D. M. Influence of ocean winds on the pelagic ecosystem in upwelling regions. Proc. Natl. Acad. Sci. 105, 1965–1970 (2008).Article 
ADS 
CAS 

Google Scholar 
Thompson, A. R. et al. Putting the Pacific marine heatwave into perspective: The response of larval fish off southern California to unprecedented warming in 2014–2016 relative to the previous 65 years. Glob. Change Biol. 28, 1766–1785 (2022).Article 
CAS 

Google Scholar 
Suryan, R. M. et al. Ecosystem response persists after a prolonged marine heatwave. Sci. Rep. 11, 6235 (2021).Article 
ADS 
CAS 

Google Scholar 
Bates, A. E. et al. Resilience and signatures of tropicalization in protected reef fish communities. Nat. Clim. Change 4, 62–67 (2014).Article 
ADS 

Google Scholar 
Behrens, M. & Lafferty, K. Effects of marine reserves and urchin disease on southern Californian rocky reef communities. Mar. Ecol. Prog. Ser. 279, 129–139 (2004).Article 
ADS 

Google Scholar 
Bernhardt, J. R. & Leslie, H. M. Resilience to climate change in coastal marine ecosystems. Ann. Rev. Mar. Sci. 5, 371–392 (2013).Article 

Google Scholar 
Caselle, J. E., Davis, K. & Marks, L. M. Marine management affects the invasion success of a non-native species in a temperate reef system in California, USA. Ecol. Lett. 21, 43–53 (2018).Article 

Google Scholar 
Micheli, F. et al. Evidence that marine reserves enhance resilience to climatic impacts. PLoS ONE 7, e40832 (2012).Article 
ADS 
CAS 

Google Scholar 
Olds, A. D. et al. Marine reserves help coastal ecosystems cope with extreme weather. Glob. Change Biol. 20, 3050–3058 (2014).Article 
ADS 

Google Scholar 
Freedman, R. M., Brown, J. A., Caldow, C. & Caselle, J. E. Marine protected areas do not prevent marine heatwave-induced fish community structure changes in a temperate transition zone. Sci. Rep. 10, 21081 (2020).Article 
ADS 
CAS 

Google Scholar 
Bates, A. E. et al. Climate resilience in marine protected areas and the ‘Protection Paradox’. Biol. Cons. 236, 305–314 (2019).Article 

Google Scholar 
Kirlin, J. et al. California’s Marine Life Protection Act Initiative: Supporting implementation of legislation establishing a statewide network of marine protected areas. Ocean Coast. Manag. 74, 3–13 (2013).Article 

Google Scholar 
Saarman, E. T. et al. An ecological framework for informing permitting decisions on scientific activities in protected areas. PLoS ONE 13, e0199126 (2018).Article 

Google Scholar 
Caselle, J. E., Rassweiler, A., Hamilton, S. L. & Warner, R. R. Recovery trajectories of kelp forest animals are rapid yet spatially variable across a network of temperate marine protected areas. Sci. Rep. 5, 14102 (2015).Article 
ADS 
CAS 

Google Scholar 
Hamilton, S. L., Caselle, J. E., Malone, D. P. & Carr, M. H. Incorporating biogeography into evaluations of the Channel Islands marine reserve network. PNAS 107, 18272–18277 (2010).Article 
ADS 
CAS 

Google Scholar 
Wendt, D. E. & Starr, R. M. Collaborative research: An effective way to collect data for stock assessments and evaluate marine protected areas in California. Mar. Coast. Fish. 1, 315–324 (2009).Article 

Google Scholar 
Côté, I. M. & Darling, E. S. Rethinking ecosystem resilience in the face of climate change. PLoS Biol. 8, e1000438 (2010).Article 

Google Scholar 
Holling, C. S. Resilience and stability of ecological systems. Annu. Rev. Ecol. Syst. 4, 1–23 (1973).Article 

Google Scholar 
Li, L. et al. Subregional differences in groundfish distributional responses to anomalous ocean bottom temperatures in the northeast Pacific. Glob. Change Biol. 25, 2560–2575 (2019).Article 
ADS 

Google Scholar 
Dawson, M. N. Phylogeography in coastal marine animals: A solution from California?. J. Biogeogr. 28, 723–736 (2001).Article 

Google Scholar 
Horn, M. H., Allen, L. G. & Lea, R. N. Biogeography. In The Ecology of Marine Fishes: California and Adjacent Waters (ed. Allen, L.) 3–25 (University of California Press, 2006). https://doi.org/10.1525/california/9780520246539.003.0001.Chapter 

Google Scholar 
Horn, M. H. & Allen, L. G. A distributional analysis of California coastal marine fishes. J. Biogeogr. 5, 23–42 (1978).Article 

Google Scholar 
Garrabou, J. et al. Mass mortality in Northwestern Mediterranean rocky benthic communities: Effects of the 2003 heat wave. Glob. Change Biol. 15, 1090–1103 (2009).Article 
ADS 

Google Scholar 
Smale, D. A. & Wernberg, T. Extreme climatic event drives range contraction of a habitat-forming species. Proc. R. Soc. B 280, 20122829 (2013).Article 

Google Scholar 
O’Leary, B. C. et al. Addressing criticisms of large-scale marine protected areas. Bioscience 68, 359–370 (2018).Article 

Google Scholar 
California Department of Fish and Wildlife. California Sheephead, Bodianus (formerly Semicossyphus) pulcher, Enhanced Status Report. (2021).Pinsky, M. L., Selden, R. L. & Kitchel, Z. J. Climate-driven shifts in marine species ranges: Scaling from organisms to communities. Ann. Rev. Mar. Sci. 12, 153–179 (2020).Article 

Google Scholar 
Francour, P., Mangialajo, L. & Pastor, J. Mediterranean marine protected areas and non-indigenous fish spreading. In Fish Invasions of the Mediterranean Sea: Change and Renewal (eds Golani, D. & Appelbaum-Golani, B.) 127–144 (Pensoft Publisher, 2010).
Google Scholar 
Couce, E., Ridgwell, A. & Hendy, E. J. Future habitat suitability for coral reef ecosystems under global warming and ocean acidification. Glob. Change Biol. 19, 3592–3606 (2013).Article 
ADS 

Google Scholar 
Bennett, S., Wernberg, T., Harvey, E. S., Santana-Garcon, J. & Saunders, B. J. Tropical herbivores provide resilience to a climate-mediated phase shift on temperate reefs. Ecol. Lett. 18, 714–723 (2015).Article 

Google Scholar 
Trainer, V. L. et al. Pelagic harmful algal blooms and climate change: Lessons from nature’s experiments with extremes. Harmful Algae 91, 101591 (2020).Article 

Google Scholar 
Gliwicz, Z. M., Babkiewicz, E., Kumar, R., Kunjiappan, S. & Leniowski, K. Warming increases the number of apparent prey in reaction field volume of zooplanktivorous fish. Limnol. Oceanogr. 63, S30–S43 (2018).Article 
ADS 

Google Scholar 
Nielsen, J. M. et al. Responses of ichthyoplankton assemblages to the recent marine heatwave and previous climate fluctuations in several Northeast Pacific marine ecosystems. Glob. Change Biol. 27, 506–520 (2021).Article 
ADS 

Google Scholar 
du Pontavice, H., Gascuel, D., Reygondeau, G., Stock, C. & Cheung, W. W. L. Climate-induced decrease in biomass flow in marine food webs may severely affect predators and ecosystem production. Glob. Change Biol. 27, 2608–2622 (2021).Article 
ADS 

Google Scholar 
Arimitsu, M. L. et al. Heatwave-induced synchrony within forage fish portfolio disrupts energy flow to top pelagic predators. Glob. Change Biol. 27, 1859–1878 (2021).Article 
ADS 
CAS 

Google Scholar 
Oken, K. L., Essington, T. E. & Fu, C. Variability and stability in predation landscapes: A cross-ecosystem comparison on the potential for predator control in temperate marine ecosystems. Fish Fish. 19, 489–501 (2018).Article 

Google Scholar 
Baum, J. K. & Worm, B. Cascading top-down effects of changing oceanic predator abundances. J. Anim. Ecol. 78, 699–714 (2009).Article 

Google Scholar 
Jacox, M. G. et al. Impacts of the 2015–2016 El Niño on the California current system: Early assessment and comparison to past events. Geophys. Res. Lett. 43, 7072–7080 (2016).Article 
ADS 

Google Scholar 
Brodeur, R. D., Auth, T. D. & Phillips, A. J. Major shifts in pelagic micronekton and macrozooplankton community structure in an upwelling ecosystem related to an unprecedented marine heatwave. Front. Mar. Sci. https://doi.org/10.3389/fmars.2019.00212 (2019).Article 

Google Scholar 
Field, J. C. et al. Spatiotemporal patterns of variability in the abundance and distribution of winter-spawned pelagic juvenile rockfish in the California Current. PLoS ONE 16, e0251638 (2021).Article 
CAS 

Google Scholar 
Schroeder, I. D. et al. Source water variability as a driver of rockfish recruitment in the California current ecosystem: Implications for climate change and fisheries management. Can. J. Fish. Aquat. Sci. 76, 950–960 (2019).Article 
CAS 

Google Scholar 
Echeverria, T. W. Thirty-four species of California rockfishes: Maturity and seasonality of reproduction. Fish. Bull. 85, 229–250 (1987).
Google Scholar 
Miller, A. & Sydeman, W. Rockfish response to low-frequency ocean climate change as revealed by the diet of a marine bird over multiple time scales. Mar. Ecol. Prog. Ser. 281, 207–216 (2004).Article 
ADS 

Google Scholar 
Johnson, K. F. et al. Status of lingcod (Ophiodon elongatus) along the southern U.S. west coast in 2021. 195 p. (2021).Winemiller, K. O. & Rose, K. A. Patterns of life-history diversification in North American fishes: Implications for population regulation. Can. J. Fish. Aquat. Sci. 49, 2196–2218 (1992).Article 

Google Scholar 
Stuart-Smith, R. D., Brown, C. J., Ceccarelli, D. M. & Edgar, G. J. Ecosystem restructuring along the great barrier reef following mass coral bleaching. Nature 560, 92–96 (2018).Article 
ADS 
CAS 

Google Scholar 
Starr, R. M. et al. Variation in responses of fishes across multiple reserves within a network of marine protected areas in temperate waters. PLoS ONE 10, e0118502 (2015).Article 

Google Scholar 
Ziegler, S. L. et al. External fishing effort regulates positive effects of no-take marine protected areas. Biol. Cons. 269, 109546 (2022).Article 

Google Scholar 
Jarvis, E. T. & Lowe, C. G. The effects of barotrauma on the catch-and-release survival of southern California nearshore and shelf rockfish (Scorpaenidae, Sebastes spp.). Can. J. Fish. Aquat. Sci. 65, 1286–1296 (2008).Article 

Google Scholar 
Brooks, R. et al. Nearshore Fishes Abundance and Distribution Data, California Collaborative Fisheries Research Program (CCFRP). (2022).García-Reyes, M. & Sydeman, W. J. California multivariate ocean climate indicator (MOCI) and marine ecosystem dynamics. Ecol. Ind. 72, 521–529 (2017).Article 

Google Scholar 
R Development Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing. (2021).Oksanen, J. et al. vegan: Community Ecology Package. (2020).Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D. & R Core Team. nlme: Linear and Nonlinear Mixed Effects Models. (2021). More

Our findings highlight the critical role of polyP in Sodalinema stali-formed cyanobacterial mats, as it was dynamically accumulated and recycled during acclimation to P fluctuations.Cellular response to progressive P starvationAnalogous to planktonic cyanobacteria, growth under low P availability could be sustained by recycling polyP, which acted as a primary P source (Fig. 2a) [16, 23, 24]. We further attribute the rapid reduction of easily dispensable cellular P-containing compounds to the substitution of cellular phospholipids with S- or N-containing membrane lipids to maintain growth at the onset of P stress (Fig. 2a) [15, 23]. However, the exhaustion of this easily dispensable P pool limited proliferation in Phase 2, and the metabolic strategy switched from a focus on growth towards maintenance (Fig. 5). The interpretation of prevailing cellular processes based on our results is graphically summarized and explained in detail below (Fig. 5).Fig. 5: Schematic interpretation of cellular phosphorus (P) cycling in a cyanobacterial mat, based on significant changes of the monitored parameters (arbitrary units).a At low P availability, initially contained polyphosphate (polyP) was recycled simultaneously with phosphate uptake to sustain growth at constant C:N:P ratios. Further proliferation at the onset of P stress in Phase 1 was sustained by mobilization of cellular P, e.g. phospholipids, which led to rapidly increasing C:N:P ratios. Severe P stress in Phase 2, indicated by increasing APase activity, prevented proliferation and photosynthesis, indicated by a loss of green chlorophyll pigments. PolyP accumulation by deficiency response occurs with severely increasing P stress, whereby globular DNA accumulation indicates the allocation of P contained in DNA into polyP. P re-addition to the P-stressed cells in Phase 3 triggered overplus uptake and narrow C:N:P ratios, transitioning to luxury uptake at higher C:N:P ratios following polyP recycling. b At high P availability, polyP in Phase 1 was accumulated by overplus uptake at narrow C:N:P ratios, transitioning to luxury uptake at higher C:N:P ratios during polyP recycling in Phase 2. P-deprivation in Phase 3 did not affect the cells, which we attributed to a sufficient amount of phosphate in the residual medium or within the biofilm matrix. Arrows indicate phosphorus transformation processes, whereby arrows pointing towards DNA represent cell growth. Yellow granules = polyP, blue granules = globular DNA spheres, P = phospholipids, S = substitute lipids.Full size imageSevere P stress in Phase 2 was indicated by the colour change from green towards yellow-green (Fig. S1) and increasing APase activity (Fig. 2a). The colour change suggested the loss of photosynthetic pigments [40], but we could not clarify whether this occurred through active cellular pigment reduction or degradation of available chlorophyll e.g., by oxidation. The increasing APase activity (Fig. 2a) suggested that Sodalinema stali is capable of hydrolysing organic P [14]. Even though APase expression did not trigger proliferation, it likely hydrolysed a potentially available organic P pool, as increasing DIC, NH4 and decreasing pH indicated progressive decay and remineralisation of organic matter (Fig. 1a). This suggests that in analogous oligotrophic environments with often fluctuating conditions, the strategy has to be maximizing the utilization of external P sources contained in organic and inorganic sediment particles that get trapped in the EPS [41]. The sediment can contain large amounts of organic P [42] and the fluctuating physico-chemical gradients in the EPS matrix due to high daytime pH and low oxygen conditions at night, facilitate P desorption from metal oxides, leading to higher dissolved phosphate concentrations within the mat, compared to the overlying water body [3]. However, alternating redox conditions at the SWI could also trigger polyP release from benthic microorganisms to the sediments, where it could act as a P source for the benthic food-chain, or ultimately trigger the formation of mineral P phases [32], to sustainably remove P from the aquatic cycle. Either way, we suggest that polyP-containing cyanobacterial mats critically impact P fluxes at the SWI.With persisting severe P stress and increasing APase activity in Phase 2, polyP accumulation as a deficiency response was observed (Fig. 2a), which has been reported from planktonic cyanobacteria of different habitats [24, 29, 23], as well as stream periphyton [28]. However, the reasons causing this deficiency response remain unresolved. In marine phytoplankton of the oligotrophic Sargasso Sea, Martin et al. [23] excluded that polyP-rich cells were in a perpetual overplus state with ‘undetectable’ pulses driving this state and suggested that polyP accumulation occurred as a cellular stress response. In other studies, reduced biosynthesis of P-rich rRNA coincided with deficiency responses [26, 28] and led to the suggestion that polyP accumulation at P concentrations below a certain threshold required for growth occurs because of P allocation changes away from growth and towards storage. Further, APase can hydrolytically cleave phosphate groups from nucleic acids and convert DNA-lipid-P to DNA-lipids, which were shown to self-assemble into globular lipid-based DNA micelles [43]. These preferentially anchor on cell membranes [44], and indeed, such DNA spheres were found to accumulate at the cell’s polar membranes in our experiments adjacent to polyP during deficiency response (Fig. 4a: Phase 2,c). Therefore, we suggest that intracellular P recovery by cleavage from P-rich DNA and reallocation to polyP, and potentially reduced rRNA synthesis [31], is also a strategy in benthic mats of Sodalinema stali as a response to severe P stress when P availability is too low to sustain growth. This supports the theory of a reallocation of resources away from growth towards flexibly available P and energy storage. Such direct intracellular P cycling could be beneficial to help retain P within the cyanobacterial population; while external P moieties such as dissolved organic P within the matrix can act as an additional P source, they are also likely to be subject to nutrient competition between cyanobacteria and other organisms inhabiting the matrix.Such effects of potential interactions in terms of nutrient competition or provision between cyanobacteria and mutualistic microorganisms contained within the same EPS matrix are difficult to assess and we cannot exclude some potential effects on our results. However, mutualistic microorganisms that are naturally contained in many cyanobacterial or algal cultures are often critical for metacommunity functioning and hence, working with axenic mat-forming strains may even further falsify any obtained results. Furthermore, microscopic analyses revealed that Sodalinema always dominated the biomass and hence, it is here considered reasonable to work with a non-axenic culture.Cellular response to a simulated P pulseIn P-deficient cells, the affinity of the P uptake system is typically increased to maximize P uptake for future pulses [13, 45]. The simulated P pulse to the P-stressed cells in Phase 3 led to a rapid increase of the cellular P content by 1260% relative to C within 3 days (Fig. 2a), whereby P was accumulated to a significant part as polyP, which is characteristic for overplus uptake [25]. Many different types of oligotrophic aquatic habitats experience only temporal P pulses, e.g., from redox changes at the benthic interface leading to P release from the sediment [32], storm run-off [28], upwelling [46], or excretions of aquatic animals [47]. The capability of microorganisms to immediately take up, store, and efficiently re-use this P by overplus uptake is hence of critical importance for a population to sustain a potential subsequent period of low P availability. Overplus uptake is typically accompanied by the overall slow growth of the population and cellular recovery from P starvation, including ultrastructural organization and recovery of the photosynthetic apparatus [48]. This took one week after re-feeding of P-starved Nostoc sp. PCC 7118 cells [48]—a timeframe very similar to the delayed onset of photosynthesis observed in our study, indicated by the elevated pH at day 9 (Fig. 1a). Regarding overplus-triggering mechanisms following P pulses, Solovochenko et al. [48] suggested that overplus uptake occurs due to a delayed down-regulation of high-affinity Pi transporters, which are active during P starvation, and emphasized the simultaneous advantage of osmotically inert polyP accumulation as a response to dramatically high phosphate concentrations in the cells. Even though APase levels declined following our experimental P re-addition, they were significantly elevated for at least 9 days (Fig. 2a). As our experimental design involved replacing the medium with APase-free, BG11 + medium after Phase 2, we assume that the APase detected in Phase 3 was actively produced, and we conclude that previously relevant, low-P response mechanisms are slowly disengaged with some sort of lag, even when ambient P is repleted. Following cellular recovery, Sodalinema now recycled stored polyP instead of further accumulating it during the transition from overplus-to luxury uptake, which was reflected in the increasing C:N:P molar ratios and decreasing polyP levels without significant additional phosphate uptake (Figs. 1a, 2a, 5).Qualitative observations on polyP distributionMost methods applied to analyse polyP in microorganisms are quantitative and do not contain information on its spatial distribution within a population. The here observed variable distribution of polyP between the cells during luxury uptake and deficiency response, as well as the retention of polyP in few individual filaments during polyP recycling in Phase 1 of the low P experiment (Fig. 4) suggests strategies of either slow growth with a retention of polyP, or of high growth with polyP recycling. This was also suggested for cells of a unicellular Synechocystis sp. PCC 6803 population during overplus uptake [47]. In contrast, polyP in our experiment was distributed homogeneously between all cyanobacterial cells during overplus uptake (Fig. 4a: Phase 3, Fig. 4b: Phase 1). Yet, we are unaware of any polyP distribution study in multicellular or mat-forming cyanobacteria and hence, further mechanisms of interactions, e.g., cell-to-cell communication [49, 50], might also contribute to purposeful differentiation of cells or filaments within a common matrix.In summary, our study shows that the mat-forming Sodalinema stali (1) is capable of luxury uptake, overplus uptake and deficiency response with a heterogenous polyP distribution during polyP recycling, luxury uptake and deficiency response, while (2) dynamically adjusting cellular P content to changing phosphate concentrations. (3) Proliferation is sustained under the expense of polyP, followed by P acquisition from other easily dispensable cellular P-containing compounds under the onset of P stress. (4) Further, biosynthetic allocation changes away from growth towards maintenance with relative polyP accumulation at the expense of P-rich DNA are conducted under severe P stress. Our findings demonstrate the extraordinary capabilities of mat-forming cyanobacteria to adapt their P acquisition strategies to strong P fluctuations. While lasting proliferation under P limitation requires the mobilization of additional P sources through regeneration of P from particulate matter, the transition to net P accumulation under excess ambient P is rapid and effective. Since current projections of climate and land use change include intensified pulses of P load to aquatic ecosystems [50], e.g., through external input from surplus of agriculture fertilizer, inefficient wastewater treatment plants, and internal loads via the mobilization of legacy P, these P ‘bioaccumulators’ could form an important component in P remediation by temporarily accumulating P within the mat, and synthesizing polyP that could ultimately stimulate the formation of mineral P phases to sustainably remove P from the aquatic cycle. More

To calculate total root biomass C colonized by AM and EcM fungi, we developed a workflow that combines multiple publicly available datasets to ultimately link fine root stocks to mycorrhizal colonization estimates (Fig. 1). These estimates were individually derived for 881 different spatial units that were constructed by combining 28 different ecoregions, 15 land cover types and six continents. In a given spatial unit, the relationship between the proportion of AM- and EcM-plants aboveground biomass and the proportion of AM- and EcM-associated root biomass depends on the prevalence of distinct growth forms. Therefore, to increase the accuracy of our estimates, calculations were made separately for woody and herbaceous vegetation and combined in the final step and subsequently mapped. Below we detail the specific methodologies we followed within the workflow and the main assumptions and uncertainties associated.Fig. 1Workflow used to create maps of mycorrhizal fine root biomass carbon. The workflow consists of two main steps: (1) Estimation of total fine root stock capable to form mycorrhizal associations with AM and EcM fungi and (2) estimation of the proportion of fine roots colonized by AM and EcM fungi.Full size imageDefinition of spatial unitsAs a basis for mapping mycorrhizal root abundances at a global scale, we defined spatial units based on a coarse division of Bailey’s ecoregions23 After removing regions of permanent ice and water bodies, we included 28 ecoregions defined according to differences in climatic regimes and elevation (deposited at Dryad-Table S1). A map of Bailey’s ecoregions was provided by the Oak Ridge National Laboratory Distributed Active Archive Center24 at 10 arcmin spatial resolution. Due to potential considerable differences in plant species identities, ecoregions that extended across multiple continents were split for each continent. The continent division was based upon the FAO Global Administrative Unit Layers (http://www.fao.org/geonetwork/srv/en/). Finally, each ecoregion-continent combination was further divided according to differences in land cover types using the 2015 Land Cover Initiative map developed by the European Space Agency at 300 m spatial resolution (https://www.esa-landcover-cci.org/). To ensure reliability, non-natural areas (croplands and urban areas), bare areas and water bodies were discarded (Table 1). In summary, a combination of 28 ecoregions, 15 land cover types and six continents were combined to define a total of 881 different spatial units (deposited at Dryad-Table S2). The use of ecoregion/land cover/continent combination provided a much greater resolution than using a traditional biome classification and allowed to account for human-driven transformations of vegetation, the latter based on the land cover data.Table 1 List of land cover categories within the ESA CCI Land Cover dataset, used to assemble maps of mycorrhizal root biomass.Full size tableMycorrhizal fine root stocksTotal root C stocksEstimation of the total root C stock in each of the spatial units was obtained from the harmonized belowground biomass C density maps of Spawn et al.20. These maps are based on continental-to-global scale remote sensing data of aboveground biomass C density and land cover-specific root-to-shoot relationships to generate matching belowground biomass C maps. This product is the best up-to-date estimation of live root stock available. For subsequent steps in our workflow, we distinguished woody and herbaceous belowground biomass C as provided by Spawn et al.20. As the tundra belowground biomass C map was provided without growth form distinction, it was assessed following a slightly different workflow (see Section 2.2.3 for more details). To match the resolution of other input maps in the workflow, all three belowground biomass C maps were scaled up from the original spatial resolution of 10-arc seconds (approximately 300 m at the equator) to 10 arc‐minutes resolution (approximately 18.5 km at the equator) using the mean location of the raster cells as aggregation criterion.As the root biomass C maps do not distinguish between fine and coarse roots and mycorrhizal fungi colonize only the fine fractions of the roots, we considered the fine root fraction to be 88,5% and 14,1% of the total root biomass for herbaceous and woody plants, respectively. These constants represent the mean value of coarse/fine root mass ratios of herbaceous and woody plants provided by the Fine-Root Ecology Database (FRED) (https://roots.ornl.gov/)25 (deposited at Dryad-Table S3). Due to the non-normality of coarse/fine root mass ratios, mean values were obtained from log-transformed data and then back-transformed for inclusion into the workflow.Finally, the belowground biomass C maps consider the whole root system, but mycorrhizal colonization occurs mainly in the upper 30 cm of the soil18. Therefore, we estimated the total fine root stocks in the upper 30 cm by applying the asymptotic equation of vertical root distribution developed by Gale & Grigal26:$$y=1-{beta }^{d}$$where y is the cumulative root fraction from the soil surface to depth d (cm), and β is the fitted coefficient of extension. β values of trees (β = 0.970), shrubs (β = 0.978) and herbs (β = 0.952) were obtained from Jackson et al.27. A mean value was then calculated for trees and shrubs to obtain a woody vegetation β value of 0.974. As a result, we estimated that 54.6% of the total live root of woody vegetation and 77.1% of herbaceous vegetation is stored in the upper 30 cm of the soil. In combination, this allowed deriving fine root C stocks in the upper 30 cm of woody and herbaceous vegetation.The proportion of root stocks colonized by AM and EcMThe proportion of root stock that forms associations with AM or EcM fungi was obtained from the global maps of aboveground biomass distribution of dominant mycorrhizal types published by Soudzilovskaia et al.14. These maps provide the relative abundance of EcM and AM plants based on information about the biomass of grass, shrub and tree vegetation at 10arcmin resolution. To match with belowground root woody plants biomass data, proportions of AM trees and shrubs underlying the maps of Soudzilovskaia et al.14 were summed up to obtain the proportion of AM woody vegetation. The same was done for EcM trees and shrubs.Our calculations are subjected to the main assumption that, within each growth form, the proportion of aboveground biomass associated with AM and EcM fungi reflects the proportional association of AM and EM fungi to belowground biomass. We tested whether root:shoot ratios were significantly different between AM and EcM woody plants (the number of EcM herbaceous plants is extremely small17). Genera were linked to growth form based on the TRY database (https://www.try-db.org/)19 and the mycorrhizal type association based on the FungalRoots database17. Subsequently, it was tested whether root:shoot ratios of genera from the TRY database (https://www.try-db.org/)19 were significantly different for AM vs EcM woody plants. No statistically significant differences (ANOVA-tests p-value = 0.595) were found (Fig. 2).Fig. 2Mean and standar error of root to shoot ratios of AM and EcM woody plant species.Full size imageEstimation of mycorrhizal fine root stocksWe calculated the total biomass C of fine roots that can potentially be colonized by AM or EcM fungi by multiplying the total woody and herbaceous fine root C biomass in the upper 30 cm of the soil by the proportion of AM and EcM of woody and herbaceous vegetation. In the case of tundra vegetation, fine root C stocks were multiplied by the relative abundance of AM and EcM vegetation without distinction of growth forms (for simplicity, this path was not included in Fig. 1, but can be seen in Fig. 3. As tundra vegetation consists mainly of herbs and small shrubs, the distinction between woody and herbaceous vegetation is not essential in this case.Fig. 3Workflow used to create mycorrhizal fine root biomass C maps specific for tundra areas.Full size imageFinally, we obtained the mean value of mycorrhiza growth form fine root C stocks in each of the defined spatial units. These resulted in six independent estimations: AM woody, AM herbaceous, EcM woody, EcM herbaceous, AM tundra and EcM tundra total fine root biomass C (Fig. 4).Fig. 4Fine root biomass stocks capable to form association with AM (a) and EcM (b) fungi for woody, herbaceous and tundra vegetation. Final AM and EcM stock result from the sum of the growth form individual maps. There were no records of fine root biomass of EcM herbaceous vegetation.Full size imageThe intensity of root colonization by mycorrhizal fungiColonization databaseThe FungalRoot database is the largest up-to-date compilation of intensity of root colonization data, providing 36303 species observations for 14870 plant species. Colonization data was filtered to remove occurrences from non-natural conditions (i.e., from plantations, nurseries, greenhouses, pots, etc.) and data collected outside growing seasons. Records without explicit information about habitat naturalness and growing season were maintained as colonization intensity is generally recorded in the growing season of natural habitats. When the intensity of colonization occurrences was expressed in categorical levels, they were converted to percentages following the transformation methods stated in the original publications. Finally, plant species were distinguished between woody and herbaceous species using the publicly available data from TRY (https://www.trydb.org/)19. As a result, 9905 AM colonization observations of 4494 species and 521 EcM colonization observations of 201 species were used for the final calculations (Fig. 5).Fig. 5Number of AM (a) and EcM (b) herbaceous and woody plant species and total observations obtained from FungalRoot database.Full size imageThe use of the mean of mycorrhizal colonization intensity per plant species is based on two main assumptions:

1)

The intensity of root colonization is a plant trait: It is known that the intensity of mycorrhizal infections of a given plant species varies under different climatic and soil conditions28,29, plant age30 and the identity of colonizing fungal species31. However, Soudzilovskaia et al.9 showed that under natural growth conditions the intraspecific variation of root mycorrhizal colonization is lower than interspecific variation, and is within the range of variations in other plant eco-physiological traits. Moreover, recent literature reported a positive correlation between root morphological traits and mycorrhizal colonization, with a strong phylogenetic signature of these correlations32,33. These findings provide support for the use of mycorrhizal root colonization of plants grown in natural conditions as a species-specific trait.

2)

The percentage of root length or root tips colonized can be translated to the percentage of biomass colonized: intensity of root colonization is generally expressed as the proportion of root length colonized by AM fungi or proportion of root tips colonized by EcM fungi (as EcM infection is restricted to fine root tips). Coupling this data with total root biomass C stocks requires assuming that the proportion of root length or proportion of root tips colonized is equivalent to the proportion of root biomass colonized. While for AM colonization this equivalence can be straightforward, EcM colonization can be more problematic as the number of root tips varies between tree species. However, given that root tips represent the terminal ends of a root network34, the proportion of root tips colonized by EcM fungi can be seen as a measurement of mycorrhizal infection of the root system and translated to biomass independently of the number of root tips of each individual. Yet, it is important to stress that estimations of fine root biomass colonized by AM and EcM as provided in this paper might not be directly comparable.

sPlot databaseThe sPlotOpen database21 holds information about the relative abundance of vascular plant species in 95104 different vegetation plots spanning 114 countries. In addition, sPlotOpen provides three partially overlapping resampled subset of 50000 plots each that has been geographically and environmentally balanced to cover the highest plant species variability while avoiding rare communities. From these three available subsets, we selected the one that maximizes the number of spatial units that have at least one vegetation plot. We further checked if any empty spatial unit could be filled by including sPlot data from other resampling subsets.Plant species in the selected subset were classified as AM and EcM according to genus-based mycorrhizal types assignments, provided in the FungalRoot database17. Plant species that could not be assigned to any mycorrhizal type were excluded. Facultative AM species were not distinguished from obligated AM species, and all were considered AM species. The relative abundance of species with dual colonization was treated as 50% AM and 50% ECM. Plant species were further classified into woody and herbaceous species using the TRY database.Estimation of the intensity of mycorrhizal colonizationThe percentage of AM and EcM root biomass colonized per plant species was spatially upscaled by inferring the relative abundance of AM and EcM plant species in each plot. For each mycorrhizal-growth form and each vegetation plot, the relative abundance of plant species was determined to include only the plant species for which information on the intensity of root colonization was available. Then, a weighted mean intensity of colonization per mycorrhizal-growth form was calculated according to the relative abundance of the species featuring that mycorrhizal-growth form in the vegetation plot. Lastly, the final intensity of colonization per spatial unit was calculated by taking the mean value of colonization across all plots within that spatial unit. These calculations are based on 38127 vegetation plots that hold colonization information, spanning 384 spatial units.The use of vegetation plots as the main entity to estimate the relative abundance of AM and EcM plant species in each spatial unit assumes that the plant species occurrences and their relative abundances in the selected plots are representative of the total spatial unit. This is likely to be true for spatial units that are represented by a high number of plots. However, in those spatial units where the number of plots is low, certain vegetation types or plant species may be misrepresented. We addressed this issue in our uncertainty analysis. Details are provided in the Quality index maps section.Final calculation and maps assemblyThe fraction of total fine root C stocks that is colonized by AM and EcM fungi was estimated by multiplying fine root C stocks by the mean root colonization intensity in each spatial unit. This calculation was made separately for tundra, woody and herbaceous vegetation.To generate raster maps based on the resulting AM and EcM fine root biomass C data, we first created a 10 arcmin raster map of the spatial units. To do this, we overlaid the raster map of Bailey ecoregions (10 arcmin resolution)24, the raster of ESA CCI land cover data at 300 m resolution aggregated to 10 arcmin using a nearest neighbour approach (https://www.esa-landcover-cci.org/) and the FAO polygon map of continents (http://www.fao.org/geonetwork/srv/en/), rasterized at 10 arcmin. Finally, we assigned to each pixel the corresponding biomass of fine root colonized by mycorrhiza, considering the prevailing spatial unit. Those spatial units that remained empty due to lack of vegetation plots or colonization data were filled with the mean value of the ecoregion x continent combination.Quality index mapsAs our workflow comprises many different data sources and the extracted data acts in distinct hierarchical levels (i.e plant species, plots or spatial unit level), providing a unified uncertainty estimation for our maps is particularly challenging. Estimates of mycorrhizal fine root C stocks are related mainly to belowground biomass C density maps and mycorrhizal aboveground biomass maps, which have associated uncertainties maps provided by the original publications. In contrast, estimates of the intensity of root colonization in each spatial unit have been associated with three main sources of uncertainties:

1)

The number of observations in the FungalRoot database. The mean species-level intensity of mycorrhizal colonization in the vegetation plots has been associated with a number of independent observations of root colonization for each plant species. We calculated the mean number of observations of each plant species for each of the vegetation plots and, subsequently the mean number of observations (per plant species) from all vegetation plots in each spatial unit. These spatial unit averaged number of observations ranged from 1 to 14 in AM and from 1 to 26 in EcM. A higher number of observations would indicate that the intraspecific variation in the intensity of colonization is better captured and, therefore, the species-specific colonization estimates are more robust.

2)

The relative plant coverage that was associated with colonization data. From the selected vegetation plots, only a certain proportion of plant species could be associated with the intensity of colonization data in FungalRoot database. The relative abundance of the plant species with colonization data was summed up in each vegetation plot. Then, we calculated the average values for each spatial unit. Mean abundance values ranged from 0.3 to 100% in both AM and EcM spatial units. A high number indicates that the dominant plant species of the vegetation plots have colonization data associated and, consequently, the community-averaged intensity of colonization estimates are more robust.

3)

The number of vegetation plots in each spatial unit. Each of the spatial units differs in the number of plots used to calculate the mean intensity of colonization, ranging from 1 to 1583 and from 1 to 768 plots in AM and EcM estimations, respectively. A higher number of plots is associated with a better representation of the vegetation variability in the spatial units, although this will ultimately depend on plot size and intrinsic heterogeneity (i.e., a big but homogeneous spatial unit may need fewer vegetation plots for a good representation than a small but very heterogeneous spatial unit).

We provide independent quality index maps of the spatial unit average of these three sources of uncertainty. These quality index maps can be used to locate areas where our estimates have higher or lower robustness. More

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Google Scholar  More

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

1.

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

2.

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

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

Population i: ai = 0, bi = 40

Population j: aj = 50, bj = 120

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

Population i: ai = 0, bi = 40

Population j: aj = 0, bj = 40.

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