Philippa Kaur

1.Adam PS, Borrel G, Gribaldo S. Evolutionary history of carbon monoxide dehydrogenase/acetyl-CoA synthase, one of the oldest enzymatic complexes. Proc Natl Acad Sci USA. 2018;115:E1166–73.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Borrel G, Adam PS, McKay LJ, Chen LX, Sierra-García IN, Sieber CM, et al. Wide diversity of methane and short-chain alkane metabolisms in uncultured archaea. Nat Microbiol. 2019;4:603–13.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
3.Ragsdale SW, Pierce E. Acetogenesis and the Wood–Ljungdahl pathway of CO2 fixation. BBA-Proteins Proteom. 2008;1784:1873–98.CAS 
Article 

Google Scholar 
4.Müller B, Sun L, Schnürer A. First insights into the syntrophic acetate‐oxidizing bacteria–a genetic study. MicrobiologyOpen. 2013;2:35–53.PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
5.Schuchmann K, Müller V. Energetics and application of heterotrophy in acetogenic bacteria. Appl Environ Microbiol. 2016;82:4056–69.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
6.Schuchmann K, Müller V. Autotrophy at the thermodynamic limit of life: a model for energy conservation in acetogenic bacteria. Nat Rev Microbiol. 2014;12:809–21.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
7.Evans PN, Boyd JA, Leu AO, Woodcroft BJ, Parks DH, Hugenholtz P, et al. An evolving view of methane metabolism in the Archaea. Nat Rev Microbiol. 2019;17:219–32.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Vavourakis CD, Andrei AS, Mehrshad M, Ghai R, Sorokin DY, Muyzer G. A metagenomics roadmap to the uncultured genome diversity in hypersaline soda lake sediments. Microbiome. 2018;6:1–18.Article 

Google Scholar 
9.Liu YF, Chen J, Liu ZL, Shou LB, Lin DD, Zhou L, et al. Anaerobic degradation of paraffins by thermophilic Actinobacteria under methanogenic conditions. Environ Sci Technol. 2020;54:10610–20.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
10.Zhang Y, Wu G, Jiang H, Yang J, She W, Khan I, et al. Abundant and rare microbial biospheres respond differently to environmental and spatial factors in tibetan hot springs. Front Microbiol. 2018;9:2096.PubMed 
PubMed Central 
Article 

Google Scholar 
11.Hua ZS, Han YJ, Chen LX, Liu J, Hu M, Li SJ, et al. Ecological roles of dominant and rare prokaryotes in acid mine drainage revealed by metagenomics and metatranscriptomics. ISME J. 2015;9:1280–94.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
12.Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol. 2012;19:455–77.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
13.Kang DD, Froula J, Egan R, Wang Z. MetaBAT, an efficient tool for accurately reconstructing single genomes from complex microbial communities. PeerJ. 2015;3:e1165PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
14.Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 2015;25:1043–55.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Chaumeil PA, Mussig AJ, Hugenholtz P, Parks DH. GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database. Bioinformatics. 2020;36:1925–7.CAS 

Google Scholar 
16.Hyatt D, Chen GL, LoCascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinform. 2010;11:119.Article 
CAS 

Google Scholar 
17.Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12:59.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
18.Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007;35:W182–5.PubMed 
PubMed Central 
Article 

Google Scholar 
19.Lowe TM, Eddy SR. tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 1997;25:955–64.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
20.Lagesen K, Hallin P, Rødland EA, Stærfeldt HH, Rognes T, Ussery DW. RNAmmer: consistent and rapid annotation of ribosomal RNA genes. Nucleic Acids Res. 2007;35:3100–8.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
21.Nguyen LT, Schmidt HA, von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32:268–74.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
22.Kalyaanamoorthy S, Minh BQ, Wong TK, von Haeseler A, Jermiin LS. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat Methods. 2017;14:587.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
23.Greening C, Biswas A, Carere CR, Jackson CJ, Taylor MC, Stott MB, et al. Genomic and metagenomic surveys of hydrogenase distribution indicate H2 is a widely utilised energy source for microbial growth and survival. ISME J. 2016;10:761–77.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
24.Carnevali PBM, Schulz F, Castelle CJ, Kantor RS, Shih PM, Sharon I, et al. Hydrogen-based metabolism as an ancestral trait in lineages sibling to the Cyanobacteria. Nat Commun. 2019;10:1–15.Article 
CAS 

Google Scholar 
25.Jaffe AL, Castelle CJ, Dupont CL, Banfield JF. Lateral gene transfer shapes the distribution of RuBisCO among candidate phyla radiation bacteria and DPANN archaea. Mol Biol Evol. 2019;36:435–46.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
26.Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792–7.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
27.Capella-Gutiérrez S, Silla-Martínez JM, Gabaldón T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics. 2009;25:1972–3.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
28.Søndergaard D, Pedersen CNS, Greening C. HydDB: a web tool for hydrogenase classification and analysis. Sci Rep. 2016;6:1–8.Article 
CAS 

Google Scholar 
29.El-Gebali S, Mistry J, Bateman A, Eddy SR, Luciani A, Potter SC, et al. The Pfam protein families database in 2019. Nucleic Acids Res. 2019;47:D427–32.CAS 
Article 

Google Scholar 
30.Letunic I, Bork P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic Acids Res. 2016;44:W242–5.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
31.Emms DM, Kelly S. OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy. Genome Biol. 2015;16:157.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
32.Ronquist F, Teslenko M, Van Der Mark P, Ayres DL, Darling A, Höhna S, et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012;61:539–42.PubMed 
PubMed Central 
Article 

Google Scholar 
33.Csűös M. Count: evolutionary analysis of phylogenetic profiles with parsimony and likelihood. Bioinformatics. 2010;26:1910–2.Article 
CAS 

Google Scholar 
34.Hua ZS, Wang YL, Evans PN, Qu YN, Goh KM, Rao YZ, et al. Insights into the ecological roles and evolution of methyl-coenzyme M reductase-containing hot spring Archaea. Nat Commun. 2019;10:1–11.Article 
CAS 

Google Scholar 
35.Jiao JY, Liu L, Hua ZS, Fang BZ, Zhou EM, Salam N, et al. Microbial dark matter coming to light: challenges and opportunities. Natl Sci Rev. 2021;8:nwaa280.36.Ludwig W, Strunk O, Westram R, Richter L, Meier H, Yadhukumar. ARB: a software environment for sequence data. Nucleic Acids Res. 2004;32:1363–71.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Cole JR, Wang Q, Fish JA, Chai B, McGarrell DM, Sun Y, et al. Ribosomal Database Project: data and tools for high throughput rRNA analysis. Nucleic Acids Res. 2014;42:D633–42.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
38.Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7:335–6.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
39.Rognes T, Flouri T, Nichols B, Quince C, Mahé F. VSEARCH: a versatile open source tool for metagenomics. PeerJ. 2016;4:e2584.PubMed 
PubMed Central 
Article 

Google Scholar 
40.Menzel P, Ng KL, Krogh A. Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nat Commun. 2016;7:1–9.Article 
CAS 

Google Scholar 
41.Brown CT, Olm MR, Thomas BC, Banfield JF. Measurement of bacterial replication rates in microbial communities. Nat Biotechnol. 2016;34:1256–63.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
42.Jain C, Rodriguez-R LM, Phillippy AM, Konstantinidis KT, Aluru S. High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat Commun. 2018;9:1–8.Article 
CAS 

Google Scholar 
43.Langmead B, Wilks C, Antonescu V, Charles R. Scaling read aligners to hundreds of threads on general-purpose processors. Bioinformatics. 2019;35:421–32.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
44.Salam N, Jiao JY, Zhang XT, Li WJ. Update on the classification of higher ranks in the phylum Actinobacteria. Int J Syst Evol Microbiol. 2020;70:1331–55.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
45.Richter M, Rosselló-Móra R. Shifting the genomic gold standard for the prokaryotic species definition. Proc Natl Acad Sci USA. 2009;106:19126–31.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.Chun J, Oren A, Ventosa A, Christensen H, Arahal DR, da Costa MS, et al. Proposed minimal standards for the use of genome data for the taxonomy of prokaryotes. Int J Syst Evol Micr. 2018;68:461–6.CAS 
Article 

Google Scholar 
47.Merino N, Kawai M, Boyd ES, Colman DR, McGlynn SE, Nealson KH, et al. Single-cell genomics of novel actinobacteria with the Wood–Ljungdahl pathway discovered in a serpentinizing system. Front Microbiol. 2020;11:1031.PubMed 
PubMed Central 
Article 

Google Scholar 
48.Coleman GA, Davín AA, Mahendrarajah T, Spang AA, Hugenholtz P, Szöllősi GJ, et al. A rooted phylogeny resolves early bacterial evolution. 2020. https://doi.org/10.1101/2020.07.15.205187.49.Spear JR, Walker JJ, McCollom TM, Pace NR. Hydrogen and bioenergetics in the Yellowstone geothermal ecosystem. Proc Natl Acad Sci USA. 2005;102:2555–60.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
50.Vignais PM, Billoud B. Occurrence, classification, and biological function of hydrogenases: an overview. Chem Rev. 2007;107:4206–72.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.Marreiros BC, Batista AP, Duarte AM, Pereira MM. A missing link between complex I and group 4 membrane-bound [NiFe] hydrogenases. BBA-Bioenerg. 2013;1827:198–209.CAS 
Article 

Google Scholar 
52.Martin WF. Older than genes: the acetyl CoA pathway and origins. Front Microbiol. 2020;11:817.PubMed 
PubMed Central 
Article 

Google Scholar 
53.Daniel SL, Hsu T, Dean SI, Drake HL. Characterization of the H2-and CO-dependent chemolithotrophic potentials of the acetogens Clostridium thermoaceticum and Acetogenium kivui. J Bacteriol. 1990;172:4464–71.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
54.Fontaine FE, Peterson WH, McCoy E, Johnson MJ, Ritter GJ. A new type of glucose fermentation by Clostridium thermoaceticum. J Bacteriol. 1942;43:701.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
55.Drake HL, Daniel SL. Physiology of the thermophilic acetogen Moorella thermoacetica. Microbiol Res. 2004;155:869–83.Article 

Google Scholar 
56.Hügler M, Hube H, Molyneaux SJ, Vetriani C, Sievert SM. Autotrophic CO2 fixation via the reductive tricarboxylic acid cycle in different lineages within the phylum Aquificae: evidence for two ways of citrate cleavage. Environ Microbiol. 2007;9:81–92.PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
57.Hou J, Sievert SM, Wang Y, Seewald JS, Natarajan VP, Xiao X, et al. Microbial succession during the transition from active to inactive stages of deep-sea hydrothermal vent sulfide chimneys. Microbiome. 2020;8:102.PubMed 
PubMed Central 
Article 

Google Scholar 
58.Müller AL, Kjeldsen KU, Rattei T, Pester M, Loy A. Phylogenetic and environmental diversity of DsrAB-type dissimilatory (bi) sulfite reductases. ISME J. 2015;9:1152–65.PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
59.Spang A, Stairs CW, Dombrowski N, Eme L, Lombard J, Caceres EF, et al. Proposal of the reverse flow model for the origin of the eukaryotic cell based on comparative analyses of Asgard archaeal metabolism. Nat Microbiol. 2019;4:1138–48.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
60.Seitz KW, Lazar CS, Hinrichs KU, Teske AP, Baker BJ. Genomic reconstruction of a novel, deeply branched sediment archaeal phylum with pathways for acetogenesis and sulfur reduction. ISME J. 2016;10:1696–705.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
61.Orsi WD, Vuillemin A, Rodriguez P, Coskun ÖK, Gomez-Saez GV, Lavik G, et al. Metabolic activity analyses demonstrate that Lokiarchaeon exhibits homoacetogenesis in sulfidic marine sediments. Nat Microbiol. 2020;5:248–55.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
62.Raymann K, Brochier-Armanet C, Gribaldo S. The two-domain tree of life is linked to a new root for the Archaea. Proc Natl Acad Sci USA. 2015;112:6670–5.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar  More

Based on extensive geographic sampling of morphological, ddRADseq and functional trait data, we have identified eight distinct evolutionary lineages within the putative hyperdominant taxon P. heptaphyllum. We are using a “traditional” species concept—a typological or morphological species concept where we expect different species to have discrete morphological differences but also represent monophyletic groups with limited genetic admixture. Morphological traits overlap at varying levels, but genetic data suggest that this widespread group represents distinct, independently adapted lineages that diverged over the past million years. We refute the hypothesis that P. heptaphyllum s.l. is a single species and instead find evidence for recognizing eight independently evolving species (or lineages, sensu20). Updates regarding the taxonomic treatment within P. heptaphyllum s.l. are in progress (e.g. Protium cordatum Huber sensu21) and detailed descriptions of new species are in preparation as part of a taxonomic revision (Table S2, Figure S2).Even though a large number of new species have been recently described in the Neotropics9,10,22, very few studies have attempted to resolve morphologically challenging species complexes (e.g.,13,15). In the sections below, we discuss the implications of our results in the context of: (i) revisiting the concept of hyperdominance for Amazonian trees, (ii) improving richness and diversity estimates, (iii) understanding diversification within dominant tropical lineages, and (iv) refining predictions for ecosystem function.Implications for the hyperdominance phenomenonThe fact that communities often harbor a small group of demographically abundant species in addition to a much larger number of rare species is not a recent discovery (e.g.,23 as cited in3). This pattern, also called a species oligarchy and hyperdominance, was first reported for western Amazonian forests in the early 2000s24,25 and again on the pan-Amazonian scale based on data from the ATDN, one of the largest tree community datasets ever compiled in the tropics1. Hyperdominant species have captured the imagination of many tropical ecologists. First, they have been thought to be more likely to be correctly identified than rare species22, allowing for people to use them as proxies for ecosystem-wide function. For example, Fauset et al.4 emphasized that hyperdominant species were responsible for half of carbon storage and productivity in the Amazon. Second, ecologists have proposed that hyperdominants have important shared demographic properties—they often have large geographic ranges but are only dominant in one or two regions and are often habitat specialists1. In contrast, our results suggest that the hyperdominant taxon P. heptaphyllum s.l. actually consists of several lineages warranting recognition as new species that have very different functional traits, that occupy distinct geographic ranges, and that can be rare or threatened. We postulate that similar conclusions could be reached with many of the other hyperdominant species, which are also thought to be members of species complexes (e.g., Iriartea deltoidea Ruiz and Pav.26; Eschweilera coriacea (DC.) S.A. Mori27). Further study is warranted prioritizing the study of these putative species complexes to refine our understanding of dominance across tropical regions28. If fewer species are found to be true hyperdominants, this will render conservation efforts and ecosystem modeling exercises more complicated than has been discussed to date.Taxonomic relevanceProtium is one of the best studied plant groups in the Neotropics (e.g.,17,18,19). Currently, the genus consists of approximately 200 species, and their systematics has been studied by a collaborative team of taxonomists and evolutionary biologists. Protium has a wide geographic range of specimen sampling, and genomic data are available for a number of species18,29,30. Also, current species descriptions in Protium are consistently founded on both morphological and molecular phylogenetic evidence19. Furthermore, the broad intraspecific sampling of key lineages/taxa within Protium heptaphyllum s.l. gives us high confidence in our results, in contrast to other plant groups that have not yet been subjected to similar intensive systematic study but often include diverse and abundant trees in the Amazon Basin. Here, we showed that a multidisciplinary effort and relatively short time investment on a hyperdominant species complex (3–5 years) has yielded the discovery of eight new lineages warranting species status.Multiple other tree families are relatively well studied in the Neotropics (e.g., Annonaceae, Sapotaceae, Lecythidaceae, Rubiaceae, Chrysobalanaceae, Fabaceae, Melastomataceae), but are also frequently cited as containing a few species complexes with incomplete genetic divergence and requiring detailed revision (e.g., Eschweilera27, Pouteria31). We acknowledge that the field of taxonomy is dynamic (ter Steege et al. 2019) in the sense that classification and species names are likely to change as new studies are performed. Besides that, the methods applied to describe, re-establish, or invalidate taxonomic entities are frequently inconsistent among taxonomists and experts, and can vary according to several species concepts and definitions32,33. Inconsistency in taxonomic methods has been the focus of recent debates about tree species richness in Amazonia9,10. While the authors of these studies disagreed about the number of Amazonian tree species, both sides agree that there is still much work to be done in order to improve taxonomy and to increase the pace of species discovery22. We propose that the establishment of a consensus prioritized list of other species complexes from putative hyperdominant taxa such as P. heptaphyllum, s.l. would represent a significant step to contribute to efficient progress of new species discovery. With the advent of cheaper and more rapid techniques for both molecular systematics and functional ecology, we expect that interdisciplinary approaches combined with an extensive populational sampling like we have done will be highly beneficial for the study of hyperdominant species complexes and to advance the estimates of Amazonian tree species diversity.Refining the understanding of trait variation and functional responseThe seven distinct lineages for the Central/West Amazon Basin and Northeast Amazonian regions showed substantial variation in leaf and wood traits (Fig. 4). Overall, more than half of the breadth of variation in these functional traits that has been observed across 2600 Amazonian species (e.g.,34,35,36) can be found within this single species complex.Moreover, the breadth of functional trait variation also varied markedly within lineages. The low within-lineage trait variation in the Amazonian TUC (orange) population (Fig. 4) is consistent with adaptation to dry savanna environments; that is, this lineage exhibits consistent high-water use efficiency (less negative 13C) and small, dense xylem vessels. In contrast, the CRS (red) population, which also exhibits low within-lineage functional trait variation, displays values of traits consistent with functional responses to flooding in Amazonian igapó habitats where it is found, and thus has low water use efficiency and large xylem vessels. Together, these lineages contribute to the broad variation observed in the species complex, while at the same time other lineages from Central Brazil and the Atlantic Coast (light and dark blue) show high within-lineage functional trait variation and broad distribution across habitats. Taken together across the phylogeny of the entire species complex, these functional traits can be considered to be highly labile, with some lineages adapting to contrasting extreme values whereas other sister lineages retain broad functional trait variation (Fig. 4).Importantly, this variation is consistent with the evolutionary histories of each distinct lineage within the species complex. For example, the P. heptaphyllum, P. “tucuruiense” and P. “aromaticum” (yellow, orange and blue) clades experienced a relatively recent population bottleneck and then expanded, likely with traits that permitted them to radiate into the drier savanna environments where they were then able to expand (Fig. 2). Some of these lineages have become adapted to distinct habitats and are responding very differently to current environmental conditions. Moreover, we suggest that these lineages will respond very differently to future changes in climatic conditions across the region, with some of them poorly equipped for the predicted increasing frequency and intensity of droughts37.Our results have important implications for ecosystem function. Since the advent of functional trait network studies (e.g., TRY Plant Trait Database38), understanding how plant species behave in terms of their physiological performance over large scales has led to important predictions of the consequences of future scenarios of climate and land-use change39,40,41. The importance of the Amazon region for the global climate and carbon cycle42,43 highlights the need to devote substantial effort to investigating the taxonomy of hyperdominant plant taxa. For instance, if some or most of the hyperdominant taxa actually represent multiple hidden evolutionary entities, as in P. heptaphyllum s.l., a larger fraction of Amazonian tree species would contribute proportionally more to carbon storage and cycling than described by Fauset et al.4, rendering modeling exercises and management much more nuanced. In addition, a comprehensive taxonomic review of dominant tree lineages with concomitant screening of functional traits as we conducted for this species complex, would greatly improve the understanding of climate-induced functional shifts, such as described by Esquivel‐Muelbert et al.41, and have potentially important consequences for the conservation of putative rare taxa that are nested within hyperdominant species complexes.Understanding diversification in dominant lineagesThe scenario of some hyperdominant or oligarchic taxa representing multiple diverged lineages is intriguing and relevant for understanding the processes of diversification in the Amazonian flora. Based on the demographic history of P. heptaphyllum s.l., older lineages tend to be habitat specialists and less morphologically and functionally variable. In contrast, recently evolved lineages have colonized new areas and frequently experience gene flow across very large geographic distances, and they are morphologically and functionally more variable. We hypothesize that large population sizes may be associated with higher diversification rates due to the process of population expansion followed by radiation into different habitats. Therefore, dominant lineages would have better chances to speciate via habitat specialization and generate large clades than non-dominant lineages. Habitat specialization is considered to have evolved in many tropical plant groups44,45 and many different tree genera have become specialized in contrasting environments46.According to our results, P. heptaphyllum s.l. diverged from its common ancestors around five million years ago and diversified first in the Amazon region followed by an increase in population size. Many lineages subsequently became specialized to white-sand habitats, seasonally flooded forests in the Rio Negro Basin (Igapó forests), and floodplain forests (várzea forests). However, more recently, lineages dispersed into neighboring floristic domains (e.g. Cerrado and Atlantic Forest) and ecotone areas, colonizing regions where the annual precipitation is currently lower and seasonal. These populations found in Central and Coastal Brazil are genetically very similar to each other despite the large geographical distances among them. Moreover, these populations are functionally more variable than the early-diverging lineages in the Amazon. In contrast, older lineages from the Amazon basin were found to be less morphologically plastic and more isolated in terms of gene flow. White-sand ecosystems in Amazonia are characterized by a patchy and geographically disjunct distribution47 that could have inhibited dispersal of habitat-specialist populations, resulting in repeated speciation in white-sand forests in different regions.Conservation implicationsWe show that at least four newly discovered lineages, including two “resurrected” species, are geographically restricted, demographically rare and endemic to white-sand vegetation in the Amazon (e.g. P. cordatum sensu Damasco et al. 2019a, P. “tucuruiense,” P. angustifolium Swart, and P. “reticuliflorum”). Future studies aiming to review potential species complexes among hyperdominant species could indicate other threatened lineages warranting taxonomic recognition as well as new geographic areas for conservation priority. While many studies have relied on datasets compiled by plot-inventory networks, more effort should be focused on increasing the pace of taxonomic research in the Neotropics22,28. Deforestation rates in Amazonia are likely to increase in the next few years; and yet one third of the total number of tree species are likely still undescribed or undiscovered1.We believe that reports that roughly half of all biomass and carbon storage in one of the most diverse areas on Earth is stored by a very small group of “hyperdominant” species may be misinforming policy makers and stakeholders. Based on our integrative review of P. heptaphyllum s.l., we showed that genetic diversity and functional responses to environmental gradients are much greater than expected by the hyperdominance principle. Our results revealed that a single “hyperdominant” taxon contains several lineages warranting recognition as distinct species that include great functional diversity, including at least three that are relatively rare and potentially threatened. Metapopulations of a hyperdominant taxon may be interpreted as much more resilient to future global changes than relatively rare lineages within a species complex. Thus, our findings suggest not only critical issues with species level conservation for lineages that should now be considered threatened taxa and not hyper- nor dominant, but also that some of these putatively new lineages might be at risk of extinction due to future climate change and increasing deforestation. Although the distinction between species and population delimitation may not be of consequence for calculating current carbon storage, we argue that multiple distinct lineages with limited gene flow such as we describe here will likely respond very differently to current and future global changes than a larger interbreeding metapopulation of a single lineage. We therefore caution against simplistic assumptions about biogeochemical processes and ecosystem services reliant on the attractive simplicity of a putative hyperdominance phenomenon. More

1.Bear A, Rand DG. Intuition, deliberation, and the evolution of cooperation. Proc Natl Acad Sci USA. 2016;113:936–41.CAS 
PubMed 
Article 

Google Scholar 
2.Athena Aktipis C, Boddy AM, Jansen G, Hibner U, Hochberg ME, Maley CC, et al. Cancer across the tree of life: Cooperation and cheating in multicellularity. Philos Trans R Soc B Biol Sci. 2015;370:20140219.Article 

Google Scholar 
3.D’Souza G, Shitut S, Preussger D, Yousif G, Waschina S, Kost C. Ecology and evolution of metabolic cross-feeding interactions in bacteria. Nat Prod Rep. 2018;35:455–88.PubMed 
Article 

Google Scholar 
4.West SA, Diggle SP, Buckling A, Gardner A, Griffin AS. The social lives of microbes. Annu Rev Ecol Evol Syst. 2007;38:53–77.Article 

Google Scholar 
5.Cremer J, Melbinger A, Wienand K, Henriquez T, Jung H, Frey E. Cooperation in microbial populations: theory and experimental model systems. J Mol Biol. 2019;431:4599–644.CAS 
PubMed 
Article 

Google Scholar 
6.West SA, Cooper GA. Division of labour in microorganisms: an evolutionary perspective. Nat Rev Microbiol. 2016;14:716–23.CAS 
PubMed 
Article 

Google Scholar 
7.Travisano M, Velicer GJ. Strategies of microbial cheater control. Trends Microbiol. 2004;12:72–8.CAS 
PubMed 
Article 

Google Scholar 
8.Nadell CD, Drescher K, Foster KR. Spatial structure, cooperation and competition in biofilms. Nat Rev Microbiol. 2016;14:589–600.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
9.Birch J. Are kin and group selection rivals or friends? Curr Biol. 2019;29:R433–R438.CAS 
PubMed 
Article 

Google Scholar 
10.Harcombe W. Novel cooperation experimentally evolved between species. Evolution (N. Y). 2010;64:2166–72.
Google Scholar 
11.Dobay A, Bagheri HC, Messina A, Kümmerli R, Rankin DJ. Interaction effects of cell diffusion, cell density and public goods properties on the evolution of cooperation in digital microbes. J Evol Biol. 2014;27:1869–77.CAS 
PubMed 
Article 

Google Scholar 
12.MacLean RC, Gudelj I. Resource competition and social conflict in experimental populations of yeast. Nature. 2006;441:498–501.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
13.Pande S, Kaftan F, Lang S, Svatoš A, Germerodt S, Kost C. Privatization of cooperative benefits stabilizes mutualistic cross-feeding interactions in spatially structured environments. ISME J. 2016;10:1413–23.PubMed 
Article 

Google Scholar 
14.Kreft J-U. Biofilms promote altruism. Microbiology. 2004;150:2751–60.CAS 
PubMed 
Article 

Google Scholar 
15.Chuang JS, Rivoire O, Leibler S. Simpson’s paradox in a synthetic microbial system. Science. 2009;323:272–6.CAS 
PubMed 
Article 

Google Scholar 
16.Hsu RH, Clark RL, Tan JW, Ahn JC, Gupta S, Romero PA, et al. Microbial interaction network inference in microfluidic droplets. Cell Syst. 2019;9:229–42.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
17.Park J, Kerner A, Burns MA, Lin XN Microdroplet-enabled highly parallel co-cultivation of microbial communities. PLoS One. 2011;6:e17019.18.Wilson CE, Lopatkin AJ, Craddock TJA, Driscoll WW, Eldakar OT, Lopez JV, et al. Cooperation and competition shape ecological resistance during periodic spatial disturbance of engineered bacteria. Sci Rep. 2017;7:1–13.19.Hauert C, Doebeli M. Spatial structure often inhibits the evolution of cooperation in the snowdrift game. Nature. 2004;428:643–6.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
20.Rebolleda-Gómez M, Travisano M The cost of being big: local competition, importance of dispersal, and experimental evolution of reversal to unicellularity. Am Nat. 2018;192:731–44.21.Sieuwerts S, Molenaar D, SAFT VanHijum, Beerthuyzen M, MJA Stevens, PWM Janssen, et al. Mixed-culture transcriptome analysis reveals the molecular basis of mixed-culture growth in Streptococcus thermophilus and Lactobacillus bulgaricus. Appl Environ Microbiol. 2010;76:7775–84.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
22.Blasche S, Kim Y, Mars R, Kafkia E, Maansson M, Machado D, et al. Emergence of stable coexistence in a complex microbial community through metabolic cooperation and spatio-temporal niche partitioning. bioRxiv 2019;541870.23.Gobbetti M, Corsetti A, Rossi J. The sourdough microflora. Interactions of lactic acid bacteria and yeasts: metabolism of amino acids. World J Microbiol Biotechnol. 1994;10:275–9.CAS 
PubMed 
Article 

Google Scholar 
24.Müller MJI, Neugeboren BI, Nelson DR, Murray AW. Genetic drift opposes mutualism during spatial population expansion. Proc Natl Acad Sci USA. 2014;111:1037–42.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
25.Inglis RF, Biernaskie JM, Gardner A, Kümmerli R. Presence of a loner strain maintains cooperation and diversity in well-mixed bacterial communities. Proc R Soc B Biol Sci. 2016;283:20152682.Article 
CAS 

Google Scholar 
26.García J, Traulsen A. Leaving the loners alone: Evolution of cooperation in the presence of antisocial punishment. J Theor Biol. 2012;307:168–73.PubMed 
Article 

Google Scholar 
27.Bachmann H, Fischlechner M, Rabbers I, Barfa N, Branco dos Santos F, Molenaar D, et al. Availability of public goods shapes the evolution of competing metabolic strategies. Proc Natl Acad Sci. 2013;110:14302–7.CAS 
PubMed 
Article 

Google Scholar 
28.Bull JJ, Harcombe WR Population dynamics constrain the cooperative evolution of cross-feeding. PLoS One. 2009;4:e4115.29.Smith JM. Group selection and kin selection. Nature. 1964;201:1145–7.Article 

Google Scholar 
30.Otto R, ten Brink B, Veldkamp H, Konings WN. The relation between growth rate and electrochemical proton gradient of Streptococcus cremoris. FEMS Microbiol Lett. 1983;16:69–74.Article 

Google Scholar 
31.Jarvis B, Wilrich C, Wilrich P-T. Reconsideration of the derivation of most probable numbers, their standard deviations, confidence bounds and rarity values. J Appl Microbiol. 2010;109:1660–7.CAS 
PubMed 

Google Scholar 
32.Pool WA, Neves AR, Kok J, Santos H, Kuipers OP. Natural sweetening of food products by engineering Lactococcus lactis for glucose production. Metab Eng. 2006;8:456–64.CAS 
PubMed 
Article 

Google Scholar 
33.Leenhouts KJ, Gietema J, Kok J, Venema G. Chromosomal stabilization of the proteinase genes in Lactococcus lactis. Appl Environ Microbiol. 1991;57:2568–75.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Bachmann H, Kleerebezem M, Van Hylckama Vlieg JET. High-throughput identification and validation of in situ-expressed genes of Lactococcus lactis. Appl Environ Microbiol. 2008;74:4727–36.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Nordkvist M, Jensen NBS, Villadsen J. Glucose metabolism in Lactococcus lactis MG1363 under different aeration conditions: Requirement of acetate to sustain growth under microaerobic conditions. Appl Environ Microbiol. 2003;69:3462–8.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Price CE, Branco Dos Santos F, Hesseling A, Uusitalo JJ, Bachmann H, Benavente V, et al. Adaption to glucose limitation is modulated by the pleotropic regulator CcpA, independent of selection pressure strength. BMC Evol Biol. 2019;19:1–15.Article 

Google Scholar 
37.Wilson DS. A theory of group selection. Proc Natl Acad Sci USA. 1975;72:143–6.CAS 
PubMed 
Article 

Google Scholar 
38.Rumbaugh KP, Sauer K Biofilm dispersion. Nat Rev Microbiol. 2020;18:571–86.39.Morris JJ. Black queen evolution: the role of leakiness in structuring microbial communities. Trends Genet. 2015;31:475–82.CAS 
PubMed 
Article 

Google Scholar 
40.Bachmann H, Molenaar D, Kleerebezem M, van Hylckama Vlieg JET. High local substrate availability stabilizes a cooperative trait. ISME J. 2011;5:929–32.CAS 
PubMed 
Article 

Google Scholar 
41.Gore J, Youk H, van Oudenaarden A. Snowdrift game dynamics and facultative cheating in yeast. Nature. 2009;459:253–6.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
42.Pfeiffer T, Schuster S, Bonhoeffer S. Cooperation and competition in the evolution of ATP-producing pathways. Science (80-). 2001;292:504–7.CAS 
Article 

Google Scholar 
43.Marchal M, Goldschmidt F, Derksen-Müller SN, Panke S, Ackermann M, Johnson DR A passive mutualistic interaction promotes the evolution of spatial structure within microbial populations. BMC Evol Biol. 2017;17:1–14.44.Preussger D, Giri S, Muhsal LK, Oña L, Kost C Reciprocal fitness feedbacks promote the evolution of mutualistic cooperation. Curr Biol. 2020;30:3580–90.45.Harcombe WR, Chacon J, Adamowicz E, Chubiz L, Marx C Evolution of bidirectional costly mutualism from byproduct consumption. PNAS 2018;115.46.Pacheco AR, Moel M, Segrè D Costless metabolic secretions as drivers of interspecies interactions in microbial ecosystems. Nat Commun 2019;10:1–12.47.Ponomarova O, Gabrielli N, Sévin DC, Mülleder M, Zirngibl K, Bulyha K, et al. Yeast creates a niche for symbiotic lactic acid bacteria through nitrogen overflow. Cell Syst. 2017;5:345–57.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
48.Pillai P, Gouhier TC, Vollmer SV. The cryptic role of biodiversity in the emergence of host-microbial mutualisms. Ecol Lett. 2014;17:1437–46.PubMed 
Article 

Google Scholar 
49.van Tatenhove-Pel RJ, Rijavec T, Lapanje A, van Swam I, Zwering E, Hernandez-Valdes JA, et al. Microbial competition reduces metabolic interaction distances to the low µm-range. ISME J. 2021;15:688–701.50.Co AD, van Vliet S, Kiviet DJ, Schlegel S, Ackermann M. Short-range interactions govern the dynamics and functions of microbial communities. Nat Ecol Evol. 2020;4:366–75.Article 

Google Scholar 
51.Saleski TE, Kerner AR, Chung MT, Jackman CM, Khasbaatar A, Kurabayashi K, et al. Synthrophic co-culture amplification of production phenotype for high-throughput screening of microbial strain libraries. Metab Eng. 2019;54:232–43.CAS 
PubMed 
Article 

Google Scholar 
52.Sung YJ, Young JHK, Choi H Il, Kwak HS, Sim SJ Magnetophoretic sorting of microdroplets with different microalgal cell densities for rapid isolation of fast growing strains. Sci Rep. 2017;7:1–11.53.Zengler K, Toledo G, Rappe M, Elkins J, Mathur EJ, Short JM, et al. Cultivating the uncultured. Proc Natl Acad Sci USA. 2002;99:15681–6.CAS 
PubMed 
Article 

Google Scholar 
54.Kehe J, Kulesa A, Ortiz A, Ackerman CM, Thakku SG, Sellers D, et al. Massively parallel screening of synthetic microbial communities. Proc Natl Acad Sci USA. 2019;116:12804–9.CAS 
PubMed 
Article 

Google Scholar 
55.Riehl C, Frederickson ME Cheating and punishment in cooperative animal societies. Philos Trans R Soc B Biol Sci. 2016;371:20150090.56.Nepi M, Grasso DA, Mancuso S. Nectar in plant–insect mutualistic relationships: From food reward to partner manipulation. Front Plant Sci. 2018;9:1–14.Article 

Google Scholar  More

Containing over 500 billion tonnes (Pg) of carbon (C)1,2,3, peatlands are an important part of the global C cycle, and there is considerable interest in how C accumulation in these systems has varied in the past and how it might respond to future changes in climate and land management4,5,6,7,8. Carbon accumulates in a peatland because more plant material is added to it than is lost via decay. Although rapid in the near surface (often called the acrotelm), decay rates in waterlogged deeper peat (the catotelm) are much lower, allowing peat to build up. However, the reverse can happen at times; more material can be lost than is added, resulting in a decrease in a peatland’s C store and a net release of C to the atmosphere (e.g.6,9).To reconstruct the C accumulation history of a peatland10, scientists calculate what is often called ‘aCAR’11: the apparent Carbon Accumulation Rate. To estimate aCAR, it is necessary first to establish the age of the peat down a peat core; this is done by dating samples of peat from a number of depths and fitting a curve through the data12,13. The C content of contiguous or regularly-spaced layers of peat down the core also needs to be measured. aCAR is then the amount of C (per unit area) in a layer divided by the difference in age between the top and bottom of the layer. aCAR is usually plotted against time to infer how rates of C accumulation have varied during a peatland’s developmental history. Although aCAR is a widely-used metric, problems with its interpretation have been discussed over many years (e.g.10,11,14,15,16,17,18,19).aCAR has a number of problems as an indicator of the net C accumulation rate of a peatland. First, because it is a measure of the amount of C found within a peat layer at the time of coring it is dependent on the overall age of the layer; this is because decay of the layer will continue as the layer gets older—the layer loses mass and C over time since its formation. For example, the aCAR for peat that is 3,000 years old may be greater than for peat in the same profile that is 4,000 years old. However, this difference in aCAR does not mean the peatland was necessarily accumulating C more rapidly 3,000 years ago than it was 4,000 years ago: when the layer of 4,000-year old peat was only 3,000 years old, its aCAR may have been the same as that of the layer of peat that is currently 3,000 years old. This ‘ageing’ problem makes it impossible to use aCAR to determine the true net rate of C accumulation of a peatland over time. A second problem is that aCAR does not account for what else may be happening to other layers in the peat profile10,11. When a new peat layer is being formed at the peatland surface, new C is being added to the profile, but a greater amount of older C may be being lost via decay from the rest of the peat profile below the new layer. Therefore, despite new C being added, the peatland as a whole may be losing C. This net loss of C is not part of the calculation of aCAR, which is a measure solely of how much C has been added to a peatland over a period of time (notwithstanding ongoing losses from the layer because of decay). Therefore, aCAR can only ever give positive values of C accumulation within the time period spanning the layer of interest11.Further problems in the interpretation of aCAR arise when it is calculated for near-surface peat (i.e., peat from recent decades), and these problems have been termed the ‘acrotelm effect’. We recently explained10 why aCAR from near-surface peat cannot be reliably compared to the long-term rate of C accumulation obtained from a deeper layer within a peat core, which may be several centuries, or more, old. Greater values of aCAR found in the near surface of a peatland are an artefact that arises because recently-added plant litter has decomposed much less than older, deeper, peat. The artefact is another example of the peat-ageing problem noted above, but it is exacerbated in the acrotelm because of the higher decay rates in this aerated part of the peat profile. As a result, many peatland scientists choose to ignore aCAR in the upper parts of the peat profile; they are aware of the problem and do not use or attempt to interpret the increase in aCAR in progressively younger peat (e.g.20). However, misinterpretation of the effect is common in the recent literature (e.g.21,22,23,24). These studies also mistakenly assume that aCAR, when calculated for the acrotelm as a whole (by treating the acrotelm as a single layer), gives a multi-decadal average net rate of peat C accumulation for the entire peat profile, when in fact it merely describes the C content of the acrotelm. This erroneous assumption is perhaps most prominently made by Rydin and Jeglum25, who suggest that:
“… it may be useful to have a measure of peat accumulation over the last few decades …This measure is the recent rate of carbon accumulation (RERCA), which is obtained from the bulk density [which gives the mass of peat and C] down to a dated level not far from the surface. Given the recent developments in precise and accurate dating of young peat…, this is now quite possible” [text in brackets added].
As defined by Rydin and Jeglum25 RERCA is simply aCAR calculated for a single layer of peat at the peatland surface, and, as noted above, this layer may comprise all of the acrotelm. Rydin and Jeglum25 note that RERCA may be used instead of direct measurement of the C fluxes to and from a peatland (see below); in doing so, it is clear they assume RERCA is a measure of the C budget of the peatland in its entirety.Data from peat cores are, of course, essential for understanding peatlands’ C accumulation histories, but it is their use to calculate aCAR or RERCA that we wish to challenge. In discussions with peatland scientists and policy makers we have realised that the problems with aCAR and RERCA are still not fully appreciated; in particular there is confusion over why RERCA cannot be used to give an average net C accumulation rate for a peatland as a whole over recent decades (e.g.26). To address these misunderstandings, and to expand on the explanations given in10, we present and discuss here the results from a simple numerical model based on Clymo’s14 work and a more detailed computer model of peatland development. The simple peatland model is used to illustrate, from first principles, how the acrotelm effect arises. We show how an increase in aCAR in the uppermost part of the peat profile arises even when: (1) actual rates of net C accumulation for the peatland as a whole decline over time, (2) net C accumulation rates for the peatland are steady (constant) over time, and (3) net C accumulation rates are negative (there is a net loss of C from the peatland as a whole). We also show why, except in one unusual case, RERCA is not equal to the net rate of C accumulation of the peatland. Specifically, we show how the method used to calculate RERCA is based on a misapplication of the mass balance equation. In the second part of the paper, we use the DigiBog peatland development model10,27 to show how the mismatch between aCAR and actual rates of net C accumulation, explained by our simple peatland model, applies over millennial timescales. DigiBog’s outputs enable us to calculate aCAR and actual (‘true’) rates of net C accumulation for the thousands of years over which our simulated peatlands develop. Because climate over such timescales is rarely constant, we also explore the effect of changes to temperature and net rainfall (precipitation minus evapotranspiration) on net C accumulation and aCAR.In the remainder of the paper we use a third acronym to describe the rate of C accumulation in peatlands: NCB or net carbon balance11,28. aCAR and RERCA (which, as noted above, is a special case of aCAR) are both calculated for layers of peat. NCB, in contrast, includes all C additions to, and all C losses from, a peatland and may be thought of as the true rate of net C accumulation for the whole peatland (or the whole peat column) at a particular time. NCB may be obtained directly by measuring atmosphere-peatland C exchanges using flux towers and by measuring C losses in water discharging from a peatland (e.g.29,30).Conceptualising the acrotelm effectTo illustrate how the acrotelm effect arises we first use a very simple numerical model where litter produced by peatland plants is added to the peatland surface as cohorts or layers. The rate of litter production has dimensions of mass (addition) per unit area of peatland per time (M L−2 T−1). Past cohorts (M L−2) decay at a specified proportionate rate (proportion per time—T−1) in accordance with Clymo14. We consider three scenarios. In Scenario 1, which we term the ‘establishment phase’ of a peatland, new peat forms on, for example, a bare mineral surface. This phase is illustrated for a column of peat in Fig. 1. In the model shown in the figure, the peatland grows over a period of five notional timesteps (Δt1–Δt5). Throughout this period, litter production is constant at a rate of 1.0 per timestep (arbitrary mass units per unit area), while decay occurs at a fixed proportionate rate of 0.33 per timestep. The peat column is shown for each timestep and its height is proportional to the total mass of peat (M L−2) (Fig. 1). During Δt1, 1.0 mass unit of litter (per unit area) is added and no peat is lost to decay. Therefore, the net mass accumulation rate is 1.0. During Δt2 1.0 mass unit of new litter is again added, but 0.33 mass units of old litter or peat are lost from the pre-existing cohort (formed during Δt1), so that its mass is reduced to 0.67. Therefore, by the end of Δt2 the peatland contains 1.67 mass units compared to 1.0 at the end of Δt1. In other words, the net gain or accumulation rate has reduced to 0.67 in Δt2 from 1.0 in Δt1. This pattern continues and the peatland continues to grow, but at a decreasing rate, during the remaining timesteps (Δt3–Δt5) as shown in the figure. The actual C balance (NCB) of the peatland is also shown in the figure, where it is assumed that C comprises half of the mass of the litter/peat (see below).Figure 1Scenario 1 showing the establishment phase of a peatland. Changes to a single column of peat of unit area are shown for five separate timesteps (Δt). Newly-added litter/peat is shown in pale green. Older peat is shown in pale brown/orange. The numbers in each layer of peat show the mass of peat and, in brackets, the mass of carbon, both per unit area (arbitrary units). NCB denotes the actual or real rate of net C accumulation and is given by the gains of C (in new litter) minus the loss of C from the decay of the older peat cohorts for each Δt, as shown in the boxes above the columns. The arrows below the boxes represent gains (down-pointing) and losses (up-pointing) of C. aCAR is calculated for each cohort or layer at Δt5. The graph shows NCB and aCAR for each Δt. In Scenario 1, the water table always resides at the base of the peat column; there is no catotelm, and all of the peat is aerated—it is acrotelm peat.Full size imageA scientist wishing to know rates of net peat and C accumulation (NCB), and how these change over time, during the establishment phase (i.e., between Δt1 and Δt5), might measure, at each time step, C fluxes to and from the peatland directly (see “Introduction” and reference to29,30). Alternatively, they could core the whole peat profile at each of the time steps and measure the peatland’s total C content on each occasion and see how it changes over time. In practice, neither approach may be practicable because of the time span involved. Most research projects last a few years (typically three to five years), and even long-term monitoring programmes rarely exceed one or two decades. Therefore, if Δt1–Δt5 spanned more than a few decades, the period would be much too long for most studies. Because of these practical difficulties, an alternative used by some peatland scientists is to take a peat core from the peatland in its current state (Δt5 in Fig. 1) and to use the core to reconstruct the apparent C accumulation history of the peatland. This reconstruction is done by dating the peat in the profile and by measuring the mass of layers of peat for which the ages of the upper and lower boundaries of the layers are known (i.e., the duration of the interval represented by the layer is known). For the simple case in Fig. 1, we may assume that the layers for which aCAR is calculated are coincident with the cohorts of litter added every timestep (Δt).If we assume that the proportion of the peat mass that is C is 50%17, NCB is simply half the value of the mass accumulation rate discussed above. Therefore, for Δt1, 0.5 mass units of C are added and none are lost. For Δt2, 0.5 mass units of C are again added, but 0.17 C units are lost from the existing layer formed during Δt1 (loss being the mass in the layer times the decay rate: 0.5 × 0.33 = 0.17), giving a net rate of C accumulation of 0.34 (Fig. 1). As noted above, the peatland continues to gain mass and C but the rate of gain decreases with time. This decrease in rate of net peat and C accumulation is shown by the grey dashed line connecting the top of peat columns in Fig. 1 and is what would be indicated by direct measurement. In contrast, aCAR erroneously suggests that the rate of growth is increasing over time because it does not consider decay other than in the dated layer of interest. The uppermost layer of peat representing the last timestep (Δt5) appears to accumulate at a greater rate than the older cohorts that have undergone progressively more decay with age. For example, the first or deepest cohort, formed in Δt1 (layer 1 in the figure), has an initial C content of 0.5, which, through decay, becomes 0.34, 0.22, 0.15, and finally 0.1 by the end of Δt5. Ascending from this deepest layer, aCAR increases to the surface, giving a pattern that is the opposite of the real rate (NCB), as shown by the graph in the lower right of Fig. 1.Two fundamental differences between aCAR and NCB are revealed here. First, NCB is measured in ‘real time’ (the fluxes are estimated for the time period in which they occur), whereas aCAR is measured between dated layers in a peat core that may have been taken many decades or centuries after the layer was first formed (at the end of Δt5 in our example). This ‘delay’ means that aCAR does not take account of changes to a cohort of peat after its initial formation. For example, if the peatland in Fig. 1 had been cored at the end of Δt4, rather than at the end of Δt5, the aCAR calculated for the lowermost cohort of peat in the profile (layer 1) would be 0.15 and not 0.1. Therefore, aCAR is dependent on the time at which the peatland is cored. This dependency is the ageing effect noted in the Introduction.Secondly, and more importantly, unlike NCB, aCAR does not consider what happens in the whole profile, which is necessary when constructing a whole-peatland C budget. For example, NCB during Δt4 comprises litter addition but also decay losses from the litter laid down in the previous three timesteps, giving a value of 0.15 (inputs of 0.5 minus decomposition losses of 0.35—see Fig. 1). When aCAR is calculated for the same time period—i.e., Δt4—it considers only the remaining mass of new peat laid down during Δt4, giving a value of 0.34 (for a core taken at the end of Δt5).If aCAR is calculated for the column of peat as a whole at Δt5 to give RERCA (i.e., the C mass of the whole core, given by 0.5 + 0.34 + 0.22 + 0.15 + 0.10 = 1.31 units of C, divided by 5 [timesteps]), it will give the correct average NCB (0.26) for the period between Δt1 to Δt5. This correspondence in values may be regarded as unusual because the acrotelm comprises the whole peatland in Scenario 1; using the whole peat profile means that all additions and losses are accounted for, making it impossible for RERCA to give anything other than the right value. However, in most situations the acrotelm sits atop a catotelm and there won’t be a correspondence of values between RERCA and NCB, as we show below in the next two scenarios and also in the next section (‘RERCA and net C balance are not comparable in near surface peat: the peatland mass balance equation and its misuse’).It may be argued that the situation in Fig. 1 is too simple because there is no lower zone of waterlogged peat; there is no catotelm, as seen in most peatlands. In Scenario 1, the water table resides at the bottom of the peat profile—all of the peat is in the acrotelm—and all litter/peat decay occurs at the same (oxic) proportionate rate. We can, however, extend the model by assuming that older, more decayed peat, is less permeable14 and that water drains less readily through it, causing the water table to rise. Figure 2 shows one realisation of this possibility (Scenario 2). In the figure, the acrotelm is in a dynamic equilibrium: its mass remains constant, but new litter continues to be added to it at a constant rate, while the oldest cohorts at the bottom become part of the waterlogged catotelm as the water table rises above them. In Scenario 2, the cohorts that have reached a specified degree of decay (80%, or cohorts with a remaining mass of 0.2) are transferred to, or become part of, the catotelm. Therefore, cohorts added to the top of the acrotelm are buried under new litter and continue to decay until they become part of the catotelm. This situation is similar to that modelled by Clymo14, with the main difference being that peat below the water table in Scenario 2 is assumed not to decay at all (Clymo14 allowed for a low rate of decay). The peatland accumulates mass at a constant rate, as shown by the straight dashed lines fitted to the top of the peat profile in Fig. 2.Figure 2Scenario 2 showing a dynamic acrotelm of constant mass, and a steadily-thickening catotelm (blue-grey shading representing waterlogged peat). Layers 1–4 become submerged by the rising water table, and by Δt9 the acrotelm comprises layers 5–9.Full size imageWhat happens if we core the peatland in Scenario 2 at the end of Δt9 and calculate aCAR for each of layers 1–9? We see that aCAR for layers 1–5 has now changed to 0.1 compared to the ascending values obtained when the peat was cored at Δt5 (Fig. 1). This difference is because these layers have now decomposed further before becoming part of the catotelm when decay ceased. An apparent increase in rates of C accumulation is still evident, however, but now in the layers of peat formed between Δt5 and Δt9 (layers 5–9) that lie above the water table and form the acrotelm.Both aCAR and NCB are plotted against time in the graph in the lower right of Fig. 2. aCAR shows a pattern similar to that from many real peat cores: a low and relatively stable aCAR in the older parts of the peat core, with an increase as one approaches the peatland surface10. It is instructive to compare these values with NCB. NCB was high initially when the peatland first formed (Scenario 1: Δt1 to Δt5) and declined to a steady value (Scenario 2: Δt6 to Δt9). As with Scenario 1, aCAR erroneously suggests that the net rate of C accumulation has increased to the present, and only one of the aCAR values corresponds to NCB (0.1), now for Δt5 (layer 5). In contrast to Scenario 1, RERCA (again applied to the acrotelm as a whole) is now wrong, giving a value of 0.26 ((0.5 + 0.34 + 0.22 + 0.15 + 0.1)/5 [timesteps]), instead of the correct value of 0.1.Finally, in Scenario 3 we may consider what happens if the peatland experiences a drought that causes the water table to fall so that layers that were in the catotelm and below the water table are now exposed above it and undergo renewed or ‘secondary’ oxic decay9. A realisation of this situation is shown in Fig. 3, where the peatland, overall, loses mass during Δt10. During the timestep, the peatland gains 1.0 mass units of litter but loses 1.06 mass units via decay of existing layers of peat above the drought water table (the layers laid down during Δt2–Δt9), giving a net rate of accumulation of − 0.06. For C the figures are a gain of 0.5 and a loss of 0.53, giving a net loss of 0.03 mass units of C per unit area (Fig. 3). If the peatland is cored at the end of Δt10 and aCAR calculated, the same problems as identified before are evident. aCAR suggests that C accumulation is increasing over time to the present. In addition, in this scenario not only does RERCA, when applied to the acrotelm as a whole, give the wrong value of net C accumulation, it also gives the wrong sign. In Scenario 3, RERCA suggests a net C accumulation rate of 0.18 (when calculated for the now deeper acrotelm incorporating the cohorts formed between Δt2 and Δt10), when in fact the peatland as a whole has become a net source of C. Here, we repeat an important point made by11: aCAR cannot be negative.Figure 3Scenario 3 showing a net loss of peat mass caused by the secondary decay of previously waterlogged layers of peat. During the drought in Δt10 the water table falls to the top of layer 1, exposing previously ‘protected’ peat in layers 2–5 to oxic decay (secondary decay).Full size imageOther scenarios in addition to the three discussed here are possible, such as ones that include changes in rates of litter production as well as changes in decay in response to drought (a modification of Scenario 3), and these may even lead to a decrease in aCAR towards the top of a core. However, the three scenarios have, between them, sufficient generality for revealing why aCAR is an unsatisfactory measure of C accumulation in peatlands. All peatlands are expressions of the balance equation: organic matter is added via litter production and is lost via decay (and sometimes erosion—not considered here). Therefore, regardless of differences in specific production and decay rates, scenarios akin to those considered above may arise in all types of peatland. In other words, the problems identified with the use of aCAR in each scenario apply regardless of the values of litter production and the decay coefficient that are used. It is clear that aCAR is misnamed; it is not a measure of net C accumulation rate—it never can be because of the way it is calculated. In the next section we extend our analysis to show that the calculation of aCAR is based on an erroneous version of the peatland mass balance equation. For simplicity, we confine our analysis to RERCA, which, as we note severally above, is aCAR applied to the uppermost part of the peat profile spanning the most recent decades in a peatland’s history; often, the acrotelm as a whole.RERCA and NCB are not comparable in near surface peat: the peatland mass balance equation and its misuseAn advantage of the simple peatland model is that the problems associated with calculating aCAR for near-surface peat become readily apparent. While it is not difficult to grasp intuitively how the artefact of an apparent increase in rates of net C accumulation arises, the exact cause of the apparent increase can be easily identified when cohorts of litter or peat are tracked over time as in Scenarios 1 and 2 (Figs. 1, 2). The simple model is, however, even more useful when considering the problem of RERCA. Without recourse to the simple model, it seems reasonable to suggest that the mass of the acrotelm divided by its overall age (obtained by dating the peat at its base), gives a reliable ‘bulk’, or time-averaged, estimate of the net rate of C accumulation for the peatland as a whole. In Fig. 1 we show that this suggestion is correct for a newly-formed acrotelm (because, in this case, all additions and all losses are considered), and it is tempting to think that it also applies to other situations. After all, the acrotelm contains new peat added in the years since the date of the acrotelm-catotelm boundary, so this would appear to be a net gain to the peatland, especially if there is little or no decay in the underlying catotelm.Our simple peatland model shows why this apparently reasonable view is mistaken in more typical situations (more typical than Scenario 1) where the acrotelm is already in existence and does not develop from scratch, and where a catotelm is present. Scenario 2 is one such situation. Here, the extant acrotelm has a fixed mass, but is dynamic: mass is added to it via litter production and mass is lost from it via decay and transfer to the catotelm (the latter caused by water-table rise). Conceptually, the acrotelm can be thought of as a simple store. To obtain an estimate of the rate of net mass addition or loss, it is necessary to look at the change in the store’s mass over time. In equation form, where Ia is litter input rate to the acrotelm (M L−2 T−1), Oa is output rate from the acrotelm (decay as well as transfer of peat to the catotelm) (M L−2 T−1), Sa is the amount of mass in the acrotelm store (M L−2), t is time (T), and i is time level, we can write the balance equation thus:$${I}_{a}-{O}_{a}=Delta {S}_{a}/Delta t=left({S}_{a,i}-{S}_{a, i-1}right)/left({t}_{i}-{t}_{i-1}right)$$
(1)
What this equation shows is that, if we measure ΔSa/Δt, we can obtain the rate of net mass addition (Ia − Oa) in the acrotelm. In Scenario 2 (Fig. 2), we see that ΔSa between any of the time steps is zero, meaning that Ia − Oa is also zero; there is no net accumulation of peat or C in the acrotelm. For example, at the end of Δt7 Sa,7 is 2.62 (1.0 + 0.67 + 0.45 + 0.30 + 0.20) (for C, the figure is half of this). At the end of Δt8 Sa,8 has the same value (although some different cohorts are now involved because the acrotelm has migrated upwards). Therefore, the right hand side of Eq. (1) gives (2.62–2.62)/(8–7) = 0. ΔSa/Δt is zero, as is Ia − Oa.Equation 1 may be rendered wrongly as follows:$${I}_{a}-{O}_{a}={S}_{a}/Delta t={S}_{a, i}/left({t}_{i}-{t}_{i-1}right)$$
ΔSa/Δt has been replaced by Sa/Δt. Here, net peat and C accumulation is being estimated from the mass in the acrotelm at one time only. This erroneous version of Eq. (1) is what is used when calculating RERCA, where Sa,i is the current mass of the acrotelm (i.e., at ti) and ti-1 now represents the age at the base of the acrotelm. If we apply this version of Eq. (1) to Scenario 2, we obtain 2.62 [the mass of peat per unit area held in the acrotelm]/5 [the difference in age between the peat at the top and bottom of the acrotelm] = 0.524 mass units per unit area per timestep, instead of the correct ΔSa/Δt value of zero. In C terms, the value is 0.262 C units per unit area per timestep (again, instead of the correct value of zero). This erroneous version of the equation can generally only produce the right result in the specific and unusual case where the mass in the acrotelm at ti−1 (Sa,t−1) is 0 (i.e., Scenario 1).However, there is a further problem here; the change in mass of the catotelm has been ignored. As noted above, the acrotelm loss term (Oa) includes two components: the loss of peat to decay and the transfer of peat from the acrotelm to the catotelm. Only the former represents a loss from the peatland; the latter remains part of the peatland and should not, therefore, be included in the loss term when calculating the net C balance of the peatland. In other words, it is not enough to look at the acrotelm alone when estimating the C budget of the peatland as a whole, even when decay in the catotelm is zero. When estimating the net rate of C accumulation for the whole peatland, a balance equation that includes both the acrotelm and the catotelm is needed:$${I}_{a}+{I}_{c}-{O}_{a}-{O}_{c}=left(Delta {S}_{a}+{Delta S}_{c}right)/Delta t=left({S}_{a,i}-{S}_{a, i-1}+{S}_{c,i}-{S}_{c,i-1}right)/left({t}_{i}-{t}_{i-1}right)$$
(2)
where the subscript c denotes the catotelm.Calculated correctly, the net C balance of the peatland in Scenario 2 between t7 and t8, for example (see above), is therefore (1.31 − 1.31 + 0.3 − 0.2)/(8 − 7) = 0.1 as shown in Fig. 2.In Scenario 2 we could have allowed the catotelm to decay slowly at an anoxic rate, which would have meant that the rate of peat accumulation would decrease very slightly over time, but this would not alter our main finding that aCAR wrongly suggests a rapid increase in rates of accumulation. In fact, the discrepancy between aCAR and NCB would be even larger in such a situation from t5 onwards; therefore, our assumption is conservative. What this simple analysis shows is that measurements in the acrotelm alone cannot, except in special cases, be used to provide information on the overall C balance of a peatland. In other words, RERCA is based on a misuse of the balance equation: to estimate the mass balance of the peatland as a whole, it is necessary to measure all of its components.Equation (2) allows for situations where the catotelm gains mass and C and where it is a net loser. However, it can sometimes be unclear where the boundary of the acrotelm and catotelm should be drawn. For example, in Scenario 3, should the acrotelm include layers 2–5 or not? It may be preferable to think of ‘acrotelm’ and ‘catotelm’ as somewhat contrived entities31, in which case the peatland should be considered a single store, giving:$${I}_{p}-{O}_{p}={Delta S}_{p}/Delta t=left({S}_{p,i}-{S}_{p, i-1}right)/left({t}_{i}-{t}_{i-1}right)$$
(3)
where the subscript p denotes ‘peatland’.Our simple peatland model and mass balance equations demonstrate why aCAR, and the special case of RERCA (aCAR for recent peat accumulation), cannot be used to understand changes to peatland C accumulation. However, peatlands develop over millennia and include a wide range of processes including feedbacks32 that can mediate their response to climate and land use, which are not represented in the simple model. We therefore used a more detailed process-based model (DigiBog) to simulate peatland development over thousands of years and to explore the dynamics of aCAR and NCB in response to perturbations to our model’s driving data.Simulating the effect on aCAR and NCB of changes in climateWe used the DigiBog peatland development model10,27 to ‘grow’ a sloping blanket peatland from the north of England over six millennia (see “Methods”). Our model simulates the peatland as a series of linked columns of peat. These can gain or lose mass (including C) depending on the climate inputs, simulated land uses and the autogenic mechanisms of the virtual peatland. And because the model records during the simulation the height of each peat column (based on the addition of mass to the peatland surface and the change in mass of each sub-surface peat layer), we can calculate the rate of change in the mass of C at each time step ((Delta t))—i.e. we know a peat column’s or the peatland’s NCB throughout the whole developmental history of the peatland (see “Methods”). At the end of a simulation we can also take a virtual core for a column and, as previously described, use the difference in age between the top and base of the layers within it, to calculate aCAR (see “Methods”). Because NCB must be calculated at the time the C fluxes occur, it is only possible to compare these past long-term dynamics of aCAR and NCB by using a peatland model.Here we show aCAR and NCB from four simulations of the single blanket peatland (see “Methods” for details of the model set up). The results are shown in Figs. 4, 5 and 6. We used the net rainfall (precipitation minus evapotranspiration) and temperature inputs from10 for a baseline simulation (Figs. 4 and 5) and ran three modifications to the same dataset: (1) a 0.4 m reduction in annual net rainfall to simulate a long-term drought; (2) a 1.5 °C increase in air temperature to simulate a warming climate; and (3) the inputs from 1 and 2 combined (Fig. 6). All other input parameters remained unchanged from the baseline simulation (see “Methods”). To create the perturbations in driving data we linearly increased or decreased the input(s) over 100 years and allowed the simulation to continue using the modified data for a further 200 years before reversing the increase to use the original time series for the remainder of the model run (the total time of a modification—400 years—is henceforth known as the perturbation). We implemented the temperature perturbation (Fig. 6a) earlier than the one for net rainfall (Fig. 6b) so that we could see the effect of later events on aCAR and NCB (Fig. 6c) (see “Methods” for the details and timings of the driving data perturbations).Figure 4Development of the virtual blanket peatland over six millennia. (a) Water-table depth and (b) the peat surface from the virtual core at the centre of the peatland from the baseline simulation (see the main text and “Methods”).Full size imageFigure 5aCAR and NCB (20 year moving average) for the baseline simulation (see “Methods”). The aCAR values are for a virtual core taken from halfway down the modelled hillslope at the end of the simulation. NCB is calculated during each year of the simulation (i.e. at the time peat is gained or lost) and is akin to measuring C fluxes. The inset shows the typical ‘uptick’ of aCAR in recently-accumulated peat layers seen in many peat cores.Full size imageFigure 6The effect of climate perturbations on simulated aCAR and NCB (20 year moving average). (a) Increase in temperature of 1.5 °C, (b) reduction in annual net rainfall of 0.4 m, and (c) both perturbations combined. The light grey vertical bars indicate the timing and duration of the perturbations and the dark grey dashed line is where C accumulation equals zero. The increases in NCB near to the beginning and end of the net rainfall perturbation (a and b) are due to the peatland water tables falling and later rising into the zone of maximum litter addition in DigiBog’s litter production equation.Full size imageOur more detailed model shows aCAR and NCB (Fig. 5) conform to a pattern similar to the one given by our simple peatland model in Scenario 2 (Fig. 2). These dynamics are also predicted by the models used by11,15, and clearly show that aCAR is not the same as NCB. The modelled ‘uptick’ in aCAR in recently accumulated peat (towards the right-hand side of Fig. 5) is also seen in real cores taken from peatlands in a wide range of environments10. The uptick is due to the ‘acrotelm effect’ explained earlier.The climate perturbations in Fig. 6 further illustrate why aCAR should not be used to represent NCB. Although they don’t have the same values as each other, aCAR and NCB in Fig. 6a increase and decrease similarly during the period in which the temperature perturbation occurs. This correspondence is because warming has shifted the peatland’s mass balance to be more in favour of plant litter production than the losses from decomposition. In this instance it might seem reasonable that aCAR can be used to indicate NCB. However, in Fig. 6b, the picture is more complicated and aCAR and NCB produce very different responses. The reduction in net rainfall deepens the peatland’s water tables, shifting the mass balance in favour of decomposition (i.e. all losses exceed all gains), but the changes in aCAR do not coincide with the timing of the perturbation. Whilst NCB is affected at the time of the climatic drying and becomes negative (there is an overall loss of C), the effects on aCAR are offset, but at no time is aCAR negative (it cannot be, as we explain earlier and as explained by11). aCAR suggests that C accumulation has reduced before the perturbation takes place but NCB has not actually changed at this time. This mismatch is because, as well as continuing to decompose, a peat layer can be altered by events that take place many years after it was originally formed. The reduction in aCAR shown in Fig. 6b is known as secondary decomposition or decay9,33. There follows a significant increase in aCAR ‘apparently’ indicating that C accumulation is also increasing when, in fact, the peatland is losing C as shown by NCB. This apparent increase in C is because a shift to deeper water tables can increase plant production; i.e., the mass added to the peatland increases. But because aCAR does not include the C fluxes from the whole peat column it does not take account of the increase in decomposition (the mass lost), and so the total change in C stored is not seen. Our simple peatland model in Fig. 3 also demonstrates how this difference between aCAR and NCB occurs.Finally, when the perturbations are combined (Fig. 6c), the increase in aCAR around 1,600 years ago, caused by an increase in temperature, is partially wiped away by secondary decomposition before sharply declining, but NCB remains unchanged. Whilst it is likely that an assessment of aCAR would conclude that C accumulation had reduced during the time when the temperature was perturbed, the interpretation of both the timing and the magnitude of the reduction would be wrong. And the increase in aCAR starting around 600 years ago would also be misinterpreted as an increase in C accumulation rather than an overall reduction in the peatland’s C store.Implications for assessing changes in peatland C accumulation ratesOur simple peatland model, mass balance equations and DigiBog simulations, along with evidence from previous studies10,11,14,15,17,19 show that aCAR and RERCA cannot generally be used to assess changes to the rate of peatland C accumulation (NCB). Therefore, studies that use aCAR to indicate changes in peatland carbon balance processes over time (acrotelm effect) or to estimate NCB are unreliable and should be viewed with considerable circumspection.As our simple peatland model scenarios show, in general, aCAR does not equal NCB11. Because all peatlands accumulate C according to the mass balance equation—i.e. assessing a peatland’s C balance requires that all of the peatland profile is taken into account and not just a dated section of it—our results apply to all peatlands in all circumstances. The only instance when we can be sure that aCAR equals NCB is when it is calculated for the whole of a peatland’s developmental history11. However, an average C accumulation rate for the entire history of the peatland is of limited use; land managers, researchers and policy makers are usually interested in how NCB has changed over time in response to climate and land-use. Although in some other circumstances (e.g. Fig. 6a) it appears that aCAR and NCB are sufficiently similar for aCAR to be useful (and sometimes they coincide—see Fig. 6 and11), this assessment can only be made because we can calculate NCB from our model outputs and compare the two quantities. But, unless NCB is known from C flux measurements or model simulations of peatland development, the correspondence of aCAR to NCB cannot be established.The results of our simulations, and those from other studies10,11, also show how some land uses or changes to the climate may cause further mismatches in the timing, magnitude and sign of aCAR and NCB. Although acknowledging that aCAR gives an erroneous estimate of past rates of NCB, Frolking et al.11 do not advocate abandoning its use. Based on our evidence, and that provided by other studies, we suggest there is a need to go further. Given that aCAR is based on a mistaken use of the balance equation and can give the wrong sign of NCB as well as the opposite trend, we believe that it should no longer be acceptable to use aCAR to indicate changes in NCB.Our simulations produce virtual and not real peat cores, and, by necessity, all models are simplifications of reality. The perturbations to our driving data are at the high end of what might be experienced naturally but are not implausible. They allow us to see more clearly how such events might affect the timing, magnitude and sign of aCAR in comparison to NCB. If our model is configured for a different type of peatland (e.g.10 simulated a raised bog) with different driving data or changes to land use, the results for aCAR and NCB will likely be different from the ones we show here. But despite these differences we would still not be able to reliably predict NCB from aCAR.The challenge of understanding if C accumulation rates have been altered by external forcing has been discussed in the literature since14, but the implications of using aCAR to indicate NCB have recently been brought to the fore because of the imperative to assess the impact of climate change and land use on peatland C cycling. By highlighting and explaining the deficiencies of aCAR, our aim is to encourage the use of more robust and reliable approaches for calculating past actual C accumulation rates (NCB). Ideally, direct measurements of C fluxes would be used10, but such observations are not available for many sites and, where they do exist, they will cover only the last few decades at most (see above). Therefore, for C accumulation histories extending to centuries and millennia, we propose that C balance models fitted to peatland age-depth (or age-mass) curves are used to estimate if NCB has changed over time. Simple models—for example, that of14—are already used in this regard and are worthy of further investigation19,28,34,35,36. For example, several studies have derived peatland NCB at the global28,37 regional35, and local17,19 scales. The authors back-calculate NCB from the net C pool using empirical models that consider autogenic long-term peat decomposition14. In a further step, with the aim of understanding if contemporary C accumulation rates were different from past rates17 and19 compared the calculated NCB from the catotelm to predictions of peat C mass transfer at the acrotelm-catotelm boundary, using a forward model of acrotelm peat decay. That being said, these approaches cannot differentiate the effects of long-term autogenic decay on peat versus that of secondary decomposition, which could be brought about by land-use or climate change.Given the limitations of such approaches, we encourage exploration of the potential of fitting more complete ecosystem models like the Holocene Peat Model38, MILLENNIA39 and DigiBog to data from peat cores to help estimate changes in peatland function over time. Observations of peat depth and downcore humification along with the inclusion of proxy data from the peatland in question—often shown in palaeoecological studies—are also important for contextualising model outputs10.In conclusion, aCAR is an unsuitable proxy for the actual C accumulation rates of peatlands. Approaches that conceptualise peatlands as dynamic C stores—the balance of all mass additions and losses – are needed. And, as we have noted, some studies, recognising the problems of aCAR, have provided potential alternatives. However, to be useful, it is likely that existing models will need to be modified, tested and their suitability assessed, or new ones developed so that credible comparisons of the effect of climate change or land uses on peatland C accumulation rates can be made. More

In the following, model transports and projected changes are expressed as ensemble interquartile ranges with individual model details provided in Tables S1-S10, reanalysis transport estimates are provided as the range across the three products examined (Table S11), and observational transports and associated references are provided in Table S12. Projections represent differences between the 1900–2000 historical mean and 2050–2100 means from the business-as-usual SSP5-8.5 (RCP8.5) scenarios for CMIP6 (CMIP5).Indian OceanIn the Indian Ocean, the South Equatorial Current (SEC) forms the northern limb of the subtropical gyre, carrying fresh ITF water to the western basin. The SEC bifurcates east of Madagascar, forming the Northeastern and Southeastern Madagascar Currents (NMC and SMC, Fig. 1). Along the African shelf, the NMC further splits southward through the Mozambique Channel (MZC) and northwards as the East African Coastal Current (EACC). Further south, the SMC and MZC transport combine into the Agulhas Current (AC). The AC extension continues westwards beyond the African cape where it retroflects, returning most water eastwards to the Indian basin33, while a part of this water (~ 21Sv34) escapes into the South Atlantic as Agulhas Leakage.Figure 1Schematic showing projected changes in WBC transport. Background colours show the multi-model mean projected change in sea surface temperature divided by the global mean change, e.g. 150% implies a warming rate 1.5 × the global average.Full size imageThe EACC transport across the CMIP6 models (interquartile range: 16.5–19.9 Sv Fig. 2, Fig. 3a) is consistent with the observed 19 Sv peak near 5°S and lies within the broad range of reanalysis estimates (7.4–23 Sv). The simulated NMC (19.4–22.7 Sv) and SMC (−10.2 to −15.7 Sv) are generally weaker than the range of observations (27–48 Sv and 20– 30 Sv, respectively), based on multiple short-term estimates (Table S12), but span similar ranges to the reanalysis (Fig. 3). Conversely, the simulated transports through the MZC (17–24.6 Sv) are slightly stronger than observations and reanalysis (15-19 Sv and 11.8-21 Sv, respectively). The simulated MZC transport seasonality, which is maximum around austral autumn, agrees well with observations and reanalysis (Figure S1). Further south, the CMIP6 AC transport increases to 50.8–61.6 Sv near Africa’s southern tip, somewhat weaker than the observational (70–77 Sv) but overlapping the weaker reanalysis estimates (47.2–53.7 Sv). A recent 3-year campaign34 found AC transport at ~ 27°E to be strongest in austral summer and weakest in winter, although large interannual variability was evident. This seasonality is qualitatively consistent with the models and reanalysis, although the observed seasonal range ~ 15 Sv is considerably larger than in the models ~ 3 Sv (Figure S1).Figure 2Historical meridional transport (left panels) and projected meridional transport change (right panels) by latitude along western boundaries shown in the map. Red/blue/green lines are multi-model median transport or transport change for CMIP6(SSP5-8.5)/CMIP5(RCP8.5)/CMIP6(SSP1-2.6) scenarios, associated shading indicates interquartile range (for high emission scenarios only). For projection panels lines are thickened where the multi-model median change is significant at the 95% level based on a two-sided Wilcoxon signed rank test. Black vertical lines and black polygons in the central map (along the WBC paths) show the location for the zonal and meridional transports presented in Fig. 3.Full size imageFigure 3Upper panel: mean transport for selected currents averaged over the twentieth century for 25 CMIP6 models (see legend), with the horizontal black line indicating the multi-model median (MMM). The bar-whisker with black dots is the associated MMM and interquartile range for 28 CMIP5 models. Grey bars indicate the range in transports from three reanalysis products (ORAS5, GODAS and C-COR). Lower panel: associated change in transport between 2050–2100 and the twentieth century means based on SSP5-8.5 (symbols and horizontal black line), SSP1-2.6 (green bar and whisker) and RCP8.5 (black bar and whisker). Positive transports indicate northward or eastwards direction in both panels. */ + indicate transports for which the CMIP5/CMIP6 MMM projected change is significant at the 95% level based on a two-sided Wilcoxon signed rank test.Full size imagePrevious work28 showed a broad-scale projected slowdown of the south Indian Ocean circulation by the end of the twenty-first century in CMIP5 models. Their reported weakening of both the western boundary Agulhas system and eastern boundary Leeuwin Current system is consistent with our CMIP6 and CMIP5 results. There is near-unanimous agreement across CMIP6 for reduced transport for the MZC (3.3 to 5.3 Sv), SMC (0.9 to 1.8 Sv), NMC (−2.1 to −3.9 Sv) and AC (3.4–7.6 Sv) (Fig. 2, Fig. 3b). However, neither the CMIP5 nor CMIP6 models show a consistent change in the EACC. In contrast to the reduction in transport along much of the southern African coast, the westward flowing AC extension south of Africa intensifies in all models (−3.6 to −7.7 Sv at 25°E) – a ~ 15% strengthening.Atlantic BasinAt the northern extent of the South Atlantic subtropical gyre, the westward SEC bifurcates with most of its water entering the equatorward North Brazil Current (NBC)—responsible for large upper-ocean cross-equatorial heat transport35. The remainder flows southward from ~ 10°S forming the relatively weak Brazil Current (BC). In the North Atlantic, the poleward flow, partly fed by the NBC, follows the western boundary of Central America as the Caribbean and Yucatan Currents ultimately emerging via the Florida Straits to form the GS. The GS breaks away from the coast at ~ 40°N, feeding the north-eastward North Atlantic Current.BC transport estimates from observations range from −19 and −23 Sv between 36 and 38°S. CMIP6 models generally simulate the maximum BC transport between about 35-40°S with values ranging from −13.9 to −25.8 Sv, which lies in the very broad range of reanalysis transports (−8.7 to −32.5 Sv). The observed NBC transport (23–26 Sv) is slightly underestimated by the ensemble (19–22 Sv, 5-10°S) with even weaker estimates from reanalysis (11.5–18 Sv). The models simulate maximum (minimum) transport in July-Aug (April–May, Fig. 4) (observed NBC seasonality estimates are not available at the latitudes examined). The BC and NBC forms at ~ 10°S37 just north of the basin-averaged zero wind-stress curl latitude. This bifurcation typically sits about 10° too far south in the CMIP models, in part related to a systematic southward bias in the model Atlantic wind field (Fig. 2, Figure S2).Figure 4Seasonal cycle of mean transport (upper panels) and projected change (lower panels) for selected currents, where the annual mean transports have been removed. Red line/shading indicate multi-model median/interquartile range for CMIP6 models; blue line/shading/dashed line indicate multi-model median/interquartile range/interdecile range for CMIP5 models. Grey shading in upper panels indicates the range of three ocean reanalysis.Full size imageIn the Northern Hemisphere, the complex circulation of the Caribbean Sea and Gulf of Mexico is represented very differently across the coarse resolution models. Compared to observations, most models (and reanalysis) underestimate the LLWBC transport of the Yucatan Current (30 Sv) with a model range of 13.5–23.3 Sv (reanalysis: 8.5–25.6 Sv). The GS transport intensifies moving northwards ( > 90 Sv) where it diverges from the coast, with the strongest transport occurring in boreal fall38. This northward intensification is absent in the models and reanalysis: northward transports peaks at 38 to 42.6 Sv between about 28-33°N (reanalysis: 37.1–46.7 Sv). The simulated winter intensification of the GS is consistent from the western margin to the extension region (Figure S1). However, the models are generally out of phase with the observations that indicate maximum transports during summer at 26.5°N39 and in the extension region38. While the reanalysis seasonality matches the models along the coast, there is poor agreement in the extension region.In the Southern Hemisphere, the BC is projected to intensify (4.2 to 6.0 Sv), especially south of 30ºS, associated with an increased northward basin interior transport27. This intensification is consistent with intensified westerlies across the Indian Ocean basin (Figure S3, Figure S4), which can increase northward Ekman transport and intensify the Indian Ocean input to the Atlantic via Agulhas Leakage31.Conversely, WBC transports weaken northwards of ~ 15°S. The cross-equatorial NBC flow is projected to weaken (−1.7 to −4.7 Sv). Similarly, the GS reduces at all latitudes with a −4.9 to −10.8 Sv (~ 15%) decrease around the GS maximum. These changes are poorly explained by surface wind changes and are likely associated with a weakened Atlantic Meridional Overturning Circulation (see below).Pacific BasinIn the South Pacific, the broad westward SEC bifurcates at the Australian margin forming the poleward EAC and equatorward Gulf of Papua Current (GPC). The EAC partially separates from the coast near 30ºS forming the Tasman Front, which continues southward to the east of New Zealand as the East Auckland Current and East Cape Currents (ECC). The remaining EAC water feeds a series of eddies that move southwards, forming the EAC extension and Tasman Leakage that provides a high-latitude pathway of water to the Indian Basin. The northward flowing GPC feeds the NGCU that exits the Northern Solomon Sea via multiple straits providing water to the subsurface Equatorial Undercurrent (EUC)17. In the Northern Hemisphere, the MC also feeds the EUC and forms the primary source of the ITF that transports warm water into the Indian Ocean17. Further north, the Kuroshio Current (KC) extends northwards from ~ 15ºN along eastern Japan, where it eventually separates and continues eastward.The observed EAC transport reaches about −22 Sv at 27°S40, with maximum/minimum transport in austral winter/summer41. CMIP6 transports are generally similar in strength (−20.3 to −23.4 Sv) and seasonality (Figure S1) to observations. While the seasonality is similar for the reanalysis products, they tend to underestimate the transport (−7 to −17 Sv). The EAC extension transport (~ 7 Sv) and Tasman Leakage (~ 8 Sv) are, however, systematically underestimated in the models (−1.4 to −7 Sv and −0.9 to −4 Sv, respectively), with some models simulating an EAC extension with northward mean flow, related to a poor representation of regional winds25.In the Northern Hemisphere, the KC intensifies from about 15 Sv at 18°N to over 20 Sv between 25 and 30°N. The CMIP6 models systematically overestimate the transport with an interquartile range of 30.1 to 44.2 Sv, which encompasses the reanalysis estimates of 38.6 to 39 Sv. Observations suggest that KC strength is weakest during winter, to the east of Taiwan42 while models and reanalysis display minimum transports earlier in autumn (Figure S1). Further north (28°N) observed transport is minimum in autumn43, while the models show no distinct seasonality. In the extension region, surface transport is weakest in winter/spring and strongest in summer/autumn44; while the model transports tend to peak in spring.At low latitudes, observed MC transport estimates varies considerably (15 to 35 Sv, Table S12). Model and reanalysis transports lie within these estimates (−18 to −25.7 Sv and −14.6 to −21.8 Sv, respectively). Observational estimates of NGCU transport decrease from 29 Sv at 12°S to ~ 20 Sv at 1–2°N, with a large seasonality that is strongest (weakest) in austral winter (summer). In agreement with observations, the CMIP6 NGCU transports between 5–10°S are 17.4 to 25.5 Sv, with a large seasonality that peaks from July–October. Reanalysis transports are generally weaker (6–20.7 Sv) with seasonality matching the climate models. The inter-model spread in ITF transport is small compared to most other currents −11.9 to  −13.4 Sv (Fig. 3). This is slightly underestimated compared to the observed transport (15 Sv)45, with sub-1000 m transport accounting for ~ 0.5 Sv of this discrepancy. Flow strengths through the multiple ITF straits each have different seasonality, largely controlled by local monsoonal wind changes and remote oceanic forcing, resulting in a bimodal seasonality in the total ITF transport, peaking in January and July45. In the models, which do not simulate realistic flow through multiple straits, there is a single annual maximum around July, with a much larger (~ 10 Sv) seasonal range compared to observations, but consistent with reanalysis (Figure S1).While the EAC core shows no consistent future change, the EAC extension and Tasman Leakage project large intensifications: −4.6 to −7.0 Sv (35–40°S) and −4.3 to −7.4 Sv (at 146°E), respectively. Previous studies have shown a negative low-frequency relationship between the EAC extension and the Tasman Front46. Consistent with this, most CMIP models project a weakening of the ECC that is fed by the Tasman Front. In the Northern Hemisphere, there is a projected weakening of the KC and Kuroshio extension across most models, but the changes are small relative to the mean transport.In the tropics, both the GPC and NGCU project unanimous model intensifications: 0.6 to 2.8 Sv and 1.9 to 4.9 Sv, respectively. In contrast, the MC and the ITF (which the MC feeds) decrease in all models (2.3–5.6 Sv and 2.4–3.2 Sv, respectively). Similar LLWBC changes in the CMIP3 models were linked to projected basin-wide negative wind stress curl anomalies flanking the equator26. These curl anomalies are also evident in the CMIP5 and CMIP6 models (Figure S3). Conversely, the ITF weakening in CMIP5 models25 and in an eddy-permitting ocean projection29 could not be explained by regional wind changes. Instead, these studies found that the changes are related to a slowdown in deep ocean waters entering the South Pacific.For the majority of currents examined across the basins, there is no significant difference in ensemble mean historical transports between CMIP5 and CMIP6. Only the Tasman Leakage and SMC demonstrate significantly different MMM transports (Table S13). In the case of the Tasman Leakage the MMM flow reverses direction from weakly eastwards in CMIP5 to weakly westward in CMIP6. This constitutes an important regional improvement although 20% of models still have spurious eastward flow in CMIP6 (compared to 57% for CMIP5). Only the MMM projected change in the AC extension transport is significantly different between model ensembles, with the CMIP6 suggesting a 40% smaller intensification compared to CMIP5.Seasonal changesAs described, many currents exhibit seasonal transport changes that are consistent across models (Figure S1). While seasonal timing is realistic for many currents, some simulated currents, for example the GS and Kuroshio system, poorly simulate the observed seasonal phase or amplitude. For some currents, comparison is hampered by uncertainties in the observed seasonality due to short observational records and large internal variability37,47,48. Both basin wide and local winds are important in setting transport seasonality, although the influence of remote winds may be lagged due to the slow propagation of ocean waves. The seasonal phase of wind stress curl in the models is broadly similar to the ERA5 reanalysis, although large discrepancies are evident, particularly near the boundaries of regions with strong seasonality differences (Figure S5a-c). We note that CMIP transports that show poor agreement with the observations or reanalysis seasonality (e.g. GS and KC extensions, NBC) are often associated with substantial biases in wind stress curl seasonality in the regions extending eastwards of the WBCs (Figure S5c).A subset of currents also exhibits consistent projected changes in transport seasonality across both model generations (Fig. 4). In the South Indian Ocean, most models project a reduced seasonal cycle for the NMC and SMC (Figure S1). Likewise, in the South Atlantic, the models consistently simulate a substantial weakening of NBC seasonality. In the South Pacific, both the EAC extension and Tasman Leakage show an amplification in seasonality. Conversely, east of New Zealand, the seasonality of the ECC is projected to decrease. In the North Pacific, the seasonality in meridional transport where the KC separates from the coast shows consistent increases (decreases) in boreal winter (summer).These currents with modified seasonality generally occur at latitudes where the phase of the projected wind stress curl seasonality also show large projected changes (except for the MC where the changes are just upstream of the MC latitudes examined; Figure S5e). These projected changes in wind stress curl seasonality are zonally oriented, occurring at transition zones where the historical wind stress curl seasonality changes rapidly with latitude (Figure S5b), suggesting that the changes are associated with a poleward expansion of the wind fields, a well-established consequence of anthropogenic climate change10, and their associated seasonality (Figure S5f.).Emergent constraintsFor a subset of currents, there appears to be a significant inter-model relationship between historical and projected transports (Figure S6). These relationships may provide emergent constraints to narrow the uncertainty associated with the large spread in projections, although such constraints may be biased by common structural errors49. For example, for the EAC extension and Tasman Leakage, models that underestimate mean transport or that have flow in the wrong direction tend to project the largest increases in southward or westward flow, respectively. Given observed EAC extension transports (~ 7 Sv, Table S12), a more moderate future change (~ 5 Sv) may therefore be more credible than the more extreme changes projected by some models. Similarly, given an observed mean transport of ~ 15 Sv, it is likely that the ITF decrease would be 3-4 Sv rather than the more extreme model estimates.Connections to surface wind changesNeglecting friction, non-linear processes and interactions with deep ocean circulation, the depth-integrated meridional transport away from the western boundaries can be related to gradients in the surface wind field via Sverdrup dynamics7. In particular, a positive (negative) wind stress curl drives northward (southward) flow in the ocean interior. WBCs provide a return flow for much of this meridional transport and a significant part of the inter-model differences in the historical mean WBC transport can be related to differences in interior transport (Table 1, Figure S7). The offset between the WBC-interior regression lines and the one-to-one line in Figure S7 for certain currents relate to inter-basin leakage or flow compensation in the deep ocean as part of the overturning circulation. For example, the ~ 10 Sv offsets for the EAC, EAC extension, GPC and NGCU result from a leakage of water via the ITF. The offset is ~ 5 Sv less than the ITF transport as there is also a net upwelling into the upper Pacific from below 1000 m. The offsets for the Atlantic basin currents, including the GS and NBC, result from the deep return flow below 1000 m.Table 1 Correlation between: interior (to the east of the WBC) and WBC transport (column 2), interior and derived Sverdrup transport (column 3), WBC and Sverdrup transport (column 4). Associated correlations for projected changes shown in columns 4, 5 and 6. Outliers (values exceeding 3 × scaled median deviations) are removed prior to the calculation of correlations. +1EAC extension includes transport to the east of New Zealand. +2MAD includes the WBC to the east and west of Madagascar. Scatter plots of Interior vs WBC and interior vs Sverdrup transport for the combined CMIP5 & 6 ensemble shown in Figure S7. Bold correlations indicate significant correlations at 95% level, based on Spearman Rank correlation.Full size tableFor many currents intermodel difference in interior transports can be explained to some degree by differences in the surface wind field via Sverdrup dynamics (Table 1). As a result, up to 50% of the intermodel variance in WBC transports can be related to differences in the surface winds. Other factors, including different overturning rates, different inter-basin transports and non-linear dynamics must be invoked to explain the wide range of mean WBC transports.Similarly, a significant fraction of intermodel projected WBC differences can be related to changes in surface wind stress curl, for most WBCs investigated (Table 1). In general, WBC whose mean differences are well explained by differences in their surface winds tend to be those whose projected transport differences are also well explained by differences in surface wind changes. The particularly poor relationship noted for the BC probably relates to the fact that the Sverdrup calculation becomes poorly defined as the eastern boundary lies at the southern tip of Africa. Other weak relationships in the Atlantic likely stem from large projected changes in the Atlantic overturning circulation50. Indeed, projected NBC decreases in CMIP5 are largely compensated by a weakening of North Atlantic Deep Water transport27.Near-surface transportWBCs affect the distribution of marine species via the dispersal of early-life stages and modulation of local thermal regimes24,51. However, ecosystem impacts will be most sensitive to near-surface circulation changes within the euphotic zone where most marine life thrives. As such, we also examine WBC transport changes in the top 100 m of the water column.For most currents examined, the change in the near-surface flow is of the same sign as the 1000 m integrated transport. An exception is the KC system, where the full-depth WBC is projected to weaken slightly along most of its length, while the near-surface flow is projected to intensify weakly north of 25°N (Figure S8, Figure S9). Previous work suggested that this intensification is associated with differences in warming rates across the KC, leading to an enhanced baroclinic flow52. As a tight connection between the state of the KC and the regional marine food webs has been documented53, this surface intensification may have consequences for the ecosystem. In contrast, the full-depth MC, which is projected to weaken, typically intensifies near the surface south of 7°S.In general, when the direction of a WBC is aligned with the warming signal (e.g. in the subtropics), a poleward intensified WBC will assist species dispersal at poleward range edges. In contrast, when the WBC flow opposes climate change velocities (e.g. in the tropics), strengthening would hinder dispersal at the poleward edges with greater propagule dispersal at the warming, equatorward edges24. Weakened/strengthened WBCs are likely to directly modify larval transport and thermal regimes, affecting rates of poleward range shifts51. In addition, other more subtle changes such as WBC broadening or modified coastal retention or dispersal pathways may also impact marine life24. More

1.Azam F, Malfatti F. Microbial structuring of marine ecosystems. Nat Rev Micro. 2007;5:782–91.CAS 
Article 

Google Scholar 
2.Hansell DA, Carlson CA, Repeta DJ, Schlitzer R. Dissolved organic matter in the ocean a controversy stimulates new insights. Oceanography. 2009;22:202–11.Article 

Google Scholar 
3.Moran MA, Kujawinski EB, Stubbins A, Fatland R, Aluwihare LI, Buchan A, et al. Deciphering ocean carbon in a changing world. Proc Natl Acad Sci. 2016;113:3143–51.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
4.Benner R, Pakulski JD, Mccarthy M, Hedges JI, Hatcher PG, Pakulski JD, et al. Bulk chemical characteristics of dissolved organic matter in the ocean. Science. 1992;255:1561–4.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
5.Carpenter LJ, Liss PS. On temperate sources of bromoform and other reactive organic bromine gases. J Geophys Res. 2000;105:20539–47.CAS 
Article 

Google Scholar 
6.Mac Monagail M, Cornish L, Morrison L, Araújo R, Critchley AT. Sustainable harvesting of wild seaweed resources. Eur J Phycol. 2017;52:371–90.Article 

Google Scholar 
7.Abdullah MI, Fredriksen S. Production, respiration and exudation of dissolved organic matter by the kelp Laminaria hyperborea along the west coast of Norway. J Mar Biol Assoc UK. 2004;84:887–94.Article 

Google Scholar 
8.Weigel BL, Pfister CA. The dynamics and stoichiometry of dissolved organic carbon release by kelp. Ecology. 2021;102:e03221.PubMed 
Article 
PubMed Central 

Google Scholar 
9.Pfister CA, Altabet MA, Weigel BL. Kelp beds and their local effects on seawater chemistry, productivity, and microbial communities. Ecology. 2019;100:e02798.PubMed 
Article 
PubMed Central 

Google Scholar 
10.Schapira M, McQuaid CD, Froneman PW. Free-living and particle-associated prokaryote metabolism in giant kelp forests: Implications for carbon flux in a sub-Antarctic coastal area. Estuar Coast Shelf Sci. 2012;106:69–79.CAS 
Article 

Google Scholar 
11.Newell R, Lucas M, Velirnirov B, Seiderer L. Quantitative significance of dissolved organic losses following fragmentation of kelp (Ecklonia maxima and Laminaria pallida). Mar Ecol Prog Ser. 1980;2:45–59.CAS 
Article 

Google Scholar 
12.Lozada M, Diéguez MC, García PE, Bigatti G, Livore JP, Gil MN, et al. Undaria pinnatifida exudates trigger shifts in seawater chemistry and microbial communities from Atlantic Patagonian coasts. bioRxiv 2020; 2020.10.21.349233.13.Kloareg B, Quatrano RS. Structure of the cell walls of marine algae and ecophysiological functions of the matrix polysaccharides. Oceanogr Mar Biol An Annu Rev. 1988;26:259–315.
Google Scholar 
14.Gacesa P. Alginates. Carbohydr Polym. 1988;8:161–82.CAS 
Article 

Google Scholar 
15.Rehm BHA. Alginates: biology and applications. Microbiology Monographs. 2009. Springer.16.Martin M, Barbeyron T, Martin R, Portetelle D, Michel G, Vandenbol M. The cultivable surface microbiota of the brown alga Ascophyllum nodosum is enriched in macroalgal-polysaccharide-degrading bacteria. Front Microbiol. 2015;6:1–14.Article 

Google Scholar 
17.Lin JD, Lemay MA, Parfrey LW. Diverse bacteria utilize alginate within the microbiome of the giant kelp Macrocystis pyrifera. Front Microbiol. 2018;9:1–16.Article 

Google Scholar 
18.Sangwan P, Chen X, Hugenholtz P, Janssen PH. Chthoniobacter flavus gen. nov., sp. nov., the first pure-culture representative of subdivision two, Spartobacteria classis nov., of the phylum Verrucomicrobia. Appl Environ Microbiol. 2004;70:5875–81.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
19.Ji S, Wang B, Lu M, Li F. Defluviitalea phaphyphila sp. nov., a novel thermophilic bacterium that degrades brown algae. Appl Environ Microbiol. 2016;82:868–77.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
20.Thomas F, Barbeyron T, Tonon T, Génicot S, Czjzek M, Michel G. Characterization of the first alginolytic operons in a marine bacterium: from their emergence in marine Flavobacteriia to their independent transfers to marine Proteobacteria and human gut Bacteroides. Environ Microbiol. 2012;14:2379–94.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
21.Kabisch A, Otto A, König S, Becher D, Albrecht D, Schüler M, et al. Functional characterization of polysaccharide utilization loci in the marine Bacteroidetes ‘Gramella forsetii’ KT0803. ISME J. 2014;8:1492–502.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
22.Koch H, Freese HM, Hahnke R, Simon M, Wietz M. Adaptations of Alteromonas sp. 76-1 to polysaccharide degradation: A CAZyme plasmid for ulvan degradation and two alginolytic systems. Front Microbiol. 2019;10:504.PubMed 
PubMed Central 
Article 

Google Scholar 
23.Gobet A, Barbeyron T, Matard-Mann M, Magdelenat G, Vallenet D, Duchaud E, et al. Evolutionary evidence of algal polysaccharide degradation acquisition by Pseudoalteromonas carrageenovora 9T to adapt to macroalgal niches. Front Microbiol. 2018;9:1–16.Article 

Google Scholar 
24.Dudek M, Dieudonné A, Jouanneau D, Rochat T, Michel G, Sarels B, et al. Regulation of alginate catabolism involves a GntR family repressor in the marine flavobacterium Zobellia galactanivorans DsijT. Nucleic Acids Res. 2020;48:7786–7800.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
25.Koch H, Dürwald A, Schweder T, Noriega-Ortega B, Vidal-Melgosa S, Hehemann JH, et al. Biphasic cellular adaptations and ecological implications of Alteromonas macleodii degrading a mixture of algal polysaccharides. ISME J. 2019;13:92–103.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
26.Thomas F, Bordron P, Eveillard D, Michel G. Gene expression analysis of Zobellia galactanivorans during the degradation of algal polysaccharides reveals both substrate-specific and shared transcriptome-wide responses. Front Microbiol. 2017;8:1808.PubMed 
PubMed Central 
Article 

Google Scholar 
27.Zhu Y, Thomas F, Larocque R, Li N, Duffieux D, Cladière L, et al. Genetic analyses unravel the crucial role of a horizontally acquired alginate lyase for brown algal biomass degradation by Zobellia galactanivorans. Environ Microbiol. 2017;19:2164–81.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
28.Hehemann J-H, Arevalo P, Datta MS, Yu X, Corzett CH, Henschel A, et al. Adaptive radiation by waves of gene transfer leads to fine-scale resource partitioning in marine microbes. Nat Commun. 2016;7:12860.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
29.Enke TN, Datta MS, Schwartzman J, Cermak N, Schmitz D, Barrere J, et al. Modular assembly of polysaccharide-degrading marine microbial communities. Curr Biol. 2019;29:1528–35.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
30.Jain A, Krishnan KP, Begum N, Singh A, Thomas FA, Gopinath A. Response of bacterial communities from Kongsfjorden (Svalbard, Arctic Ocean) to macroalgal polysaccharide amendments. Mar Environ Res. 2020;155:104874.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
31.Mitulla M, Dinasquet J, Guillemette R, Simon M, Azam F, Wietz M. Response of bacterial communities from California coastal waters to alginate particles and an alginolytic Alteromonas macleodii strain. Environ Microbiol. 2016;18:4369–77.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
32.Wietz M, Wemheuer B, Simon H, Giebel H-A, Seibt MA, Daniel R, et al. Bacterial community dynamics during polysaccharide degradation at contrasting sites in the Southern and Atlantic Oceans. Environ Microbiol. 2015;17:3822–31.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
33.Musat N, Foster R, Vagner T, Adam B, Kuypers MMM. Detecting metabolic activities in single cells, with emphasis on nanoSIMS. FEMS Microbiol Rev. 2012;36:486–511.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
34.Reintjes G, Arnosti C, Fuchs BM, Amann R. An alternative polysaccharide uptake mechanism of marine bacteria. ISME J. 2017;11:1640–50.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Liu Y, Fang J, Jia Z, Chen S, Zhang L, Gao W. DNA stable-isotope probing reveals potential key players for microbial decomposition and degradation of diatom-derived marine particulate matter. Microbiologyopen. 2020;9:1–24.
Google Scholar 
36.Orsi WD, Smith JM, Liu S, Liu Z, Sakamoto CM, Wilken S, et al. Diverse, uncultivated bacteria and archaea underlying the cycling of dissolved protein in the ocean. ISME J. 2016;10:2158–73.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Cunliffe M, Hollingsworth A, Bain C, Sharma V, Taylor JD. Algal polysaccharide utilisation by saprotrophic planktonic marine fungi. Fungal Ecol. 2017;30:135–8.Article 

Google Scholar 
38.Alonso C, Musat N, Adam B, Kuypers M, Amann R. HISH-SIMS analysis of bacterial uptake of algal-derived carbon in the Río de la Plata estuary. Syst Appl Microbiol. 2012;35:541–8.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
39.Arandia-Gorostidi N, Alonso-Sáez L, Stryhanyuk H, Richnow HH, Morán XAG, Musat N. Warming the phycosphere: differential effect of temperature on the use of diatom-derived carbon by two copiotrophic bacterial taxa. Environ Microbiol. 2020;22:1381–96.CAS 
PubMed 
Article 

Google Scholar 
40.Thomas F, Le Duff N, Leroux C, Dartevelle L, Riera P. Isotopic labeling of cultured macroalgae and isolation of 13C-labeled cell wall polysaccharides for trophic investigations. Adv Bot Res. 2020;95:1–17.Article 

Google Scholar 
41.Hardouin K, Burlot AS, Umami A, Tanniou A, Stiger-Pouvreau V, Widowati I, et al. Biochemical and antiviral activities of enzymatic hydrolysates from different invasive French seaweeds. J Appl Phycol. 2014;26:1029–42.CAS 
Article 

Google Scholar 
42.Pernthaler A, Pernthaler J, Amann R. Fluorescence in situ hybridization and catalyzed reporter deposition for the identification of marine bacteria. Appl Environ Microbiol. 2002;68:3094–101.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
43.Manz W, Amann R, Ludwig W, Wagner M, Schleifer KH. Phylogenetic oligodeoxynucleotide probes for the major subclasses of Proteobacteria: problems and solutions. Syst Appl Microbiol. 1992;15:593–600.Article 

Google Scholar 
44.Manz W, Amann R, Ludwig W, Vancanneyt M, Schleifer KH. Application of a suite of 16S rRNA-specific oligonucleotide probes designed to investigate bacteria of the phylum cytophaga-flavobacter-bacteroides in the natural environment. Microbiology. 1996;142:1097–106.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
45.Wallner G, Amann R, Beisker W. Optimizing fluorescent in situ hybridization with rRNA‐targeted oligonucleotide probes for flow cytometric identification of microorganisms. Cytometry. 1993;14:136–43.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.Pernthaler A, Pernthaler J Fluorescene in situ hybridization for the identification of environmental microbes. In: Hilario E, Mackay J (eds). Methods in Molecular Biology. Totowa, NJ: Humana Press Inc.; 2004. pp 153–64.47.Guerquin-Kern JL, Wu T Di, Quintana C, Croisy A. Progress in analytical imaging of the cell by dynamic secondary ion mass spectrometry (SIMS microscopy). Biochim Biophys Acta – Gen Subj. 2005;1724:228–38.CAS 
Article 

Google Scholar 
48.Slodzian G, Daigne B, Girard F, Boust F, Hillion F. Scanning secondary ion analytical microscopy with parallel detection. Biol Cell. 1992;74:43–50.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
49.Schneider CA, Rasband WS, Eliceiri KW. NIH Image to ImageJ: 25 years of Image Analysis. Nat Methods. 2012;9:671–5.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
50.Stryhanyuk H, Calabrese F, Kümmel S, Musat F, Richnow HH, Musat N. Calculation of single cell assimilation rates from sip-nanosims-derived isotope ratios: a comprehensive approach. Front Microbiol. 2018;9:1–15.Article 

Google Scholar 
51.Woebken D, Burow LC, Behnam F, Mayali X, Schintlmeister A, Fleming ED, et al. Revisiting N2 fixation in Guerrero Negro intertidal microbial mats with a functional single-cell approach. ISME J. 2015;9:485–96.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
52.Ramond P, Sourisseau M, Simon N, Romac S, Schmitt S, Rigaut-Jalabert F, et al. Coupling between taxonomic and functional diversity in protistan coastal communities. Environ Microbiol. 2019;21:730–49.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
53.Thomas F, Corre E, Cébron A. Stable isotope probing and metagenomics highlight the effect of plants on uncultured phenanthrene-degrading bacterial consortium in polluted soil. ISME J. 2019;13:1814–30.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
54.Pepe-Ranney C, Campbell AN, Koechli CN, Berthrong S, Buckley DH. Unearthing the ecology of soil microorganisms using a high resolution DNA-SIP approach to explore cellulose and xylose metabolism in soil. Front Microbiol. 2016;7:1–17.Article 

Google Scholar 
55.Buckley DH, Huangyutitham V, Hsu SF, Nelson TA. Stable isotope probing with 15N achieved by disentangling the effects of genome G+C content and isotope enrichment on DNA density. Appl Environ Microbiol. 2007;73:3189–95.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
56.Birnie G Centrifugal separations in molecular and cell biology. Boston: Butterworth & Co Publishers Ltd.; 1978.57.Thomas F, Dittami SM, Brunet M, Le Duff N, Tanguy G, Leblanc C, et al. Evaluation of a new primer combination to minimize plastid contamination in 16S rDNA metabarcoding analyses of alga‐associated bacterial communities. Environ Microbiol Rep. 2020;12:30–37.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
58.Klindworth A, Pruesse E, Schweer T, Peplies J, Quast C, Horn M, et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 2013;41:1–11.Article 
CAS 

Google Scholar 
59.Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
60.Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37:852–7.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
61.Rognes T, Flouri T, Nichols B, Quince C, Mahé F. VSEARCH: A versatile open source tool for metagenomics. PeerJ. 2016;4:e2584.PubMed 
PubMed Central 
Article 

Google Scholar 
62.Mcmurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One. 2013;8:1–11.Article 
CAS 

Google Scholar 
63.Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:1–21.Article 
CAS 

Google Scholar 
64.Youngblut ND, Barnett SE, Buckley DH. HTSSIP: an r package for analysis of high throughput sequencing data from nucleic acid stable isotope probing (sip) experiments. PLoS One. 2018;13:1–8.Article 
CAS 

Google Scholar 
65.Youngblut ND, Barnett SE, Buckley DH. SIPSim: a modeling toolkit to predict accuracy and aid design of DNA-SIP experiments. Front Microbiol. 2018;9:570.PubMed 
PubMed Central 
Article 

Google Scholar 
66.Thomas F, Lundqvist LCE, Jam M, Jeudy A, Barbeyron T, Sandström C, et al. Comparative characterization of two marine alginate lyases from Zobellia galactanivorans reveals distinct modes of action and exquisite adaptation to their natural substrate. J Biol Chem. 2013;288:23021–37.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
67.Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, et al. SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol. 2012;19:455–77.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
68.Aziz RK, Bartels D, Best A, DeJongh M, Disz T, Edwards RA, et al. The RAST Server: Rapid annotations using subsystems technology. BMC Genomics. 2008;9:1–15.Article 
CAS 

Google Scholar 
69.Seemann T. Prokka: Rapid prokaryotic genome annotation. Bioinformatics. 2014;30:2068–9.CAS 
Article 

Google Scholar 
70.Zhang H, Yohe T, Huang L, Entwistle S, Wu P, Yang Z, et al. DbCAN2: A meta server for automated carbohydrate-active enzyme annotation. Nucleic Acids Res. 2018;46:W95–W101.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
71.Barrett K, Lange L. Peptide-based functional annotation of carbohydrate-active enzymes by conserved unique peptide patterns (CUPP). Biotechnol Biofuels. 2019;12:1–21.Article 

Google Scholar 
72.Almagro Armenteros JJ, Tsirigos KD, Sønderby CK, Petersen TN, Winther O, Brunak S, et al. SignalP 5.0 improves signal peptide predictions using deep neural networks. Nat Biotechnol. 2019;37:420–3.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
73.El-Gebali S, Mistry J, Bateman A, Eddy SR, Luciani A, Potter SC, et al. The Pfam protein families database in 2019. Nucleic Acids Res. 2019;47:D427–D432.CAS 
Article 
PubMed 

Google Scholar 
74.Ferguson RL, Buckley EN, Palumbo AV. Response of marine bacterioplankton to differential filtration and confinement. Appl Environ Microbiol. 1984;47:49–55.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
75.Lucas M, Newell R, Velimirov B. Heterotrophic utilisation of mucilage released during fragmentation of kelp (Ecklonia maxima and Laminaria pallida) II. Differential utilisation of dissolved organic components from kelp mucilage. Mar Ecol Prog Ser. 1981;4:43–55.CAS 
Article 

Google Scholar 
76.Koop K, Newell R, Lucas M. Biodegradation and carbon flow based on kelp (Ecklonia maxima) debris in a sandy beach microcosm. Mar Ecol Prog Ser. 1982;7:315–26.Article 

Google Scholar 
77.Barbeyron T, Thomas F, Barbe V, Teeling H, Schenowitz C, Dossat C, et al. Habitat and taxon as driving forces of carbohydrate catabolism in marine heterotrophic bacteria: Example of the model algae-associated bacterium Zobellia galactanivorans DsijT. Environ Microbiol. 2016;18:4610–27.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
78.Díez-Vives C, Nielsen S, Sánchez P, Palenzuela O, Ferrera I, Sebastián M, et al. Delineation of ecologically distinct units of marine Bacteroidetes in the Northwestern Mediterranean Sea. Mol Ecol. 2019;28:2846–59.PubMed 

Google Scholar 
79.Balmonte JP, Buckley A, Hoarfrost A, Ghobrial S, Ziervogel K, Teske A, et al. Community structural differences shape microbial responses to high molecular weight organic matter. Environ Microbiol. 2019;21:557–71.CAS 
PubMed 
Article 

Google Scholar 
80.Alonso-Sáez L, Díaz-Pérez L, Morán XAG. The hidden seasonality of the rare biosphere in coastal marine bacterioplankton. Environ Microbiol. 2015;17:3766–80.PubMed 
Article 
PubMed Central 

Google Scholar 
81.Giovannoni SJ. SAR11 bacteria: the most abundant plankton in the oceans. Ann Rev Mar Sci. 2017;9:231–55.PubMed 
Article 
PubMed Central 

Google Scholar 
82.Grote J, Cameron Thrash J, Huggett MJ, Landry ZC, Carini P, Giovannoni SJ, et al. Streamlining and core genome conservation among highly divergent members of the SAR11 clade. MBio. 2012;3:1–13.Article 
CAS 

Google Scholar 
83.Ngugi DK, Stingl U. High-quality draft single-cell genome sequence of the NS5 marine group from the coastal Red Sea. Genome Announc. 2018;6:5–6.
Google Scholar 
84.Woyke T, Xie G, Copeland A, González JM, Han C, Kiss H, et al. Assembling the marine metagenome, one cell at a time. PLoS One. 2009;4:e5299.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
85.Teeling H, Fuchs BM, Bennke CM, Krüger K, Chafee M, Kappelmann L, et al. Recurring patterns in bacterioplankton dynamics during coastal spring algae blooms. Elife. 2016;5:1–29.Article 

Google Scholar 
86.Krüger K, Chafee M, Ben Francis T, Glavina del Rio T, Becher D, Schweder T, et al. In marine Bacteroidetes the bulk of glycan degradation during algae blooms is mediated by few clades using a restricted set of genes. ISME J. 2019;13:2800–16.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
87.Pei X, Chang Y. Shen J. Cloning, expression and characterization of an endo-acting bifunctional alginate lyase of marine bacterium Wenyingzhuangia fucanilytica. Protein Expr Purif. 2019;154:44–51.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
88.Xing P, Hahnke RL, Unfried F, Markert S, Huang S, Barbeyron T, et al. Niches of two polysaccharide-degrading Polaribacter isolates from the North Sea during a spring diatom bloom. ISME J. 2015;9:1410–22.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
89.Tanaka R, Shibata T, Miyake H, Mori T, Tamaru Y, Ueda M, et al. Temporal fluctuation in the abundance of alginate-degrading bacteria in the gut of abalone Haliotis gigantea over 1 year. Aquac Res. 2016;47:2899–908.CAS 
Article 

Google Scholar 
90.Bunse C, Koch H, Breider S, Simon M, Wietz M. Sweet and magnetic: succession and CAZyme expression of marine bacterial communities encountering a mix of alginate and pectin particles. bioRxiv 2020; 2020.12.08.416354.91.Park HH, Kam N, Lee EY, Kim HS. Cloning and characterization of a novel oligoalginate lyase from a newly isolated bacterium Sphingomonas sp. MJ-3. Mar Biotechnol. 2012;14:189–202.CAS 
Article 

Google Scholar 
92.Sim PF, Furusawa G, Teh AH. Functional and structural studies of a multidomain alginate lyase from Persicobacter sp. CCB-QB2. Sci Rep. 2017;7:1–9.Article 
CAS 

Google Scholar 
93.Lyu Q, Zhang K, Zhu Q, Li Z, Liu Y, Fitzek E, et al. Structural and biochemical characterization of a multidomain alginate lyase reveals a novel role of CBM32 in CAZymes. Biochim Biophys Acta Gen Subj. 2018;1862:1862–9.94.Han W, Gu J, Cheng Y, Liu H, Li Y, Li F. Novel alginate lyase (Aly5) from a polysaccharide-degrading marine bacterium, Flammeovirga sp. strain MY04: Effects of module truncation on biochemical characteristics, alginate degradation patterns, and oligosaccharide-yielding properties. Appl Environ Microbiol. 2016;82:364–74.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
95.Kim DH, Wang D, Yun EJ, Kim S, Kim SR, Kim KH. Validation of the metabolic pathway of the alginate-derived monomer in Saccharophagus degradans 2-40T by gas chromatography–mass spectrometry. Process Biochem. 2016;51:1374–9.CAS 
Article 

Google Scholar 
96.Arnosti C, Wietz M, Brinkhoff T, Hehemann J-H, Probandt D, Zeugner L, et al. The biogeochemistry of marine polysaccharides: sources, inventories, and bacterial drivers of the carbohydrate cycle. Ann Rev Mar Sci. 2021;13:9.1–9.28.Article 

Google Scholar 
97.Grondin JM, Tamura K, Déjean G, Abbott DW, Brumer H. Polysaccharide Utilization Loci: Fuelling microbial communities. J Bacteriol. 2017;199:e00860–16.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
98.Reintjes G, Arnosti C, Fuchs B. Amann R. Selfish, sharing and scavenging bacteria in the Atlantic Ocean: a biogeographical study of bacterial substrate utilisation. ISME J. 2019;13:1119–32.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
99.Arahal DR, Lucena T, Macián MC, Ruvira MA, González JM, Lekumberri I, et al. Marinomonas blandensis sp. nov., a novel marine gammaproteobacterium. Int J Syst Evol Microbiol. 2016;66:5544–9.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
100.Pontiller B, Martínez-García S, Lundin D, Pinhassi J. Labile dissolved organic matter compound characteristics select for divergence in marine bacterial activity and transcription. Front Microbiol. 2020;11:1–19.Article 

Google Scholar  More

1.Levy SB, Marshall B. Antibacterial resistance worldwide: causes, challenges and responses. Nat Med. 2004;10:S122.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Davies J, Davies D. Origins and evolution of antibiotic resistance. Microbiol Mol Biol Rev. 2010;74:417–33.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
3.Blair JM, Webber MA, Baylay AJ, Ogbolu DO, Piddock LJ. Molecular mechanisms of antibiotic resistance. Nat Rev Microbiol. 2015;13:42.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
4.Meredith HR, Srimani JK, Lee AJ, Lopatkin AJ, You L. Collective antibiotic tolerance: mechanisms, dynamics, and intervention. Nat Chem Biol. 2015;11:182.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Vega NM, Gore J. Collective antibiotic resistance: mechanisms and implications. Curr Opin Microbiol. 2014;21:28–34. http://www.sciencedirect.com/science/article/pii/S1369527414001234, antimicrobials.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
6.Flemming HC, Wingender J, Szewzyk U, Steinberg P, Rice SA, Kjelleberg S. Biofilms: an emergent form of bacterial life. Nat Rev Microbiol. 2016;14:563.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
7.Liu J, Prindle A, Humphries J, Gabalda-Sagarra M, Asally M, Lee DyD, et al. Metabolic codependence gives rise to collective oscillations within biofilms. Nature. 2015;523:550.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Liu J, Martinez-Corral R, Prindle A, Dong-yeon DL, Larkin J, Gabalda- Sagarra M, et al. Coupling between distant biofilms and emergence of nutrient time-sharing. Science. 2017;356:638–42.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
9.Donlan RM. Biofilms and device-associated infections. Emerg Infect Dis. 2001;7:277.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
10.Costerton JW, Stewart PS, Greenberg EP. Bacterial biofilms: a common cause of persistent infections. Science. 1999;284:1318–22.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Hauert C, Doebeli M. Spatial structure often inhibits the evolution of cooperation in the snowdrift game. Nature. 2004;428:643–6.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
12.Hallatschek O, Hersen P, Ramanathan S, Nelson DR. Genetic drift at expanding frontiers promotes gene segregation. Proc Natl Acad Sci USA. 2007;104:19926–30.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
13.Korolev KS, Avlund M, Hallatschek O, Nelson DR. Genetic demixing and evolution in linear stepping stone models. Rev Mod Phys. 2010;82:1691.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
14.Müller MJ, Neugeboren BI, Nelson DR, Murray AW. Genetic drift opposes mutualism during spatial population expansion. Proc Natl Acad Sci USA. 2014;111:1037–42.PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
15.Momeni B, Brileya KA, Fields MW, Shou W. Strong inter-population co- operation leads to partner intermixing in microbial communities. eLife. 2013;2:e00230.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
16.Momeni B, Waite AJ, Shou W. Spatial self-organization favors heterotypic cooperation over cheating. eLife. 2013;2:e00960.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
17.Lavrentovich MO, Nelson DR. Asymmetric mutualism in two- and three-dimensional range expansions. Phys Rev Lett. 2014;112:138102.PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
18.Gandhi SR, Yurtsev EA, Korolev KS, Gore J. Range expansions transition from pulled to pushed waves as growth becomes more cooperative in an experimental microbial population. Proc Natl Acad Sci USA. 2016;113:6922–7.19.Kayser J, Schreck CF, Gralka M, Fusco D, Hallatschek O. Collective motion conceals fitness differences in crowded cellular populations. Nat Ecol Evol. 2019;3:125–34.Article 

Google Scholar 
20.Gandhi SR, Korolev KS, Gore J. Cooperation mitigates diversity loss in a spatially expanding microbial population. Proc Natl Acad Sci USA. 2019;116:23582–7.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
21.Datta MS, Korolev KS, Cvijovic I, Dudley C, Gore J. Range expansion promotes cooperation in an experimental microbial metapopulation. Proc Natl Acad Sci USA. 2013;110:7354–9.PubMed 
Article 
PubMed Central 

Google Scholar 
22.Kimmel GJ, Gerlee P, Brown JS, Altrock PM. Neighborhood size-effects shape growing population dynamics in evolutionary public goods games. Commun Biol. 2019;2:1–10.Article 

Google Scholar 
23.Kimmel GJ, Gerlee P, Altrock PM. Time scales and wave formation in non-linear spatial public goods games. PLoS Comput Biol. 2019;15:e1007361.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
24.Gerlee P, Altrock PM. Persistence of cooperation in diffusive public goods games. Phys Rev E. 2019;99:062412.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
25.Celik Ozgen V, Kong W, Blanchard AE, Liu F, Lu T. Spatial interference scale as a determinant of microbial range expansion. Sci. Adv. 2018;4:eaau0695.26.Steenackers HP, Parijs I, Foster KR, Vanderleyden J. Experimental evolu- tion in biofilm populations. FEMS Microbiol Rev. 2016;40:373–97.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
27.Kepler TB, Perelson AS. Drug concentration heterogeneity facili- tates the evolution of drug resistance. Proc Natl Acad Sci USA. 1998;95:11514–9.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
28.Hermsen R, Deris JB, Hwa T. On the rapidity of antibiotic resistance evolution facilitated by a concentration gradient. Proc Natl Acad Sci USA. 2012;109:10775–80.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
29.Zhang Q, Lambert G, Liao D, Kim H, Robin K, Tung CK, et al. Acceleration of emergence of bacterial antibiotic resistance in connected microenvironments. Science. 2011;333:1764–7.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
30.Greulich P, Waclaw B, Allen RJ. Mutational pathway determines whether drug gradients accelerate evolution of drug-resistant cells. Phys Rev Lett. 2012;109:088101.PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
31.Fu F, Nowak MA, Bonhoeffer S. Spatial heterogeneity in drug concen- trations can facilitate the emergence of resistance to cancer therapy. PLoS Comput Biol. 2015;11:e1004142.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
32.Moreno-Gamez S, Hill AL, Rosenbloom DI, Petrov DA, Nowak MA, Pen- nings PS. Imperfect drug penetration leads to spatial monotherapy and rapid evolution of multidrug resistance. Proc Natl Acad Sci USA. 2015;112:E2874–E2883.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
33.Baym M, Lieberman TD, Kelsic ED, Chait R, Gross R, Yelin I, et al. Spatiotemporal microbial evolution on antibiotic landscapes. Science. 2016;353:1147–51.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.De Jong MG, Wood KB. Tuning spatial profiles of selection pressure to modulate the evolution of drug resistance. Phys Rev Lett. 2018;120:238102.PubMed 
PubMed Central 
Article 

Google Scholar 
35.Lenski RE, Hattingh SE. Coexistence of two competitors on one re- source and one inhibitor: a chemostat model based on bacteria and antibiotics. J Theor Biol. 1986;122:83–93.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
36.Dugatkin LA, Perlin M, Lucas JS, Atlas R. Group-beneficial traits, frequency-dependent selection and genotypic diversity: an antibiotic resistance paradigm. Proc R Soc B. 2005;272:79–83.PubMed 
Article 
PubMed Central 

Google Scholar 
37.Clark DR, Alton TM, Bajorek A, Holden P, Dugatkin LA, Atlas RM, et al. Evolution of altruists and cheaters in near-isogenic populations of Escherichia coli. Front Biosci. 2009;14:4815.CAS 
Article 

Google Scholar 
38.Perlin MH, Clark DR, McKenzie C, Patel H, Jackson N, Kormanik C, et al. Protection of Salmonella by ampicillin-resistant Escherichia coli in the presence of otherwise lethal drug concentrations. Proc R Soc B 2009;276:3759–68.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
39.Yurtsev EA, Chao HX, Datta MS, Artemova T, Gore J. Bacterial cheating drives the population dynamics of cooperative antibiotic resistance plasmids. Mol Syst Biol. 2013. https://doi.org/10.1038/msb.2013.39.40.Koster DA, Mayo A, Bren A, Alon U. Surface growth of a motile bac- terial population resembles growth in a chemostat. J Mol Biol. 2012;424:180–91.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
41.Sorg RA, Lin L, van Doorn GS, Sorg M, Olson J, Nizet V. et al. Collective resistance in microbial communities by intracellular antibiotic deactivation. PLoS Biol. 2016;14:e2000631 https://doi.org/10.1371/journal.pbio.2000631.CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
42.Hallinen KM, Karslake J, Wood KB. Delayed antibiotic exposure induces population collapse in enterococcal communities with drug-resistant subpopulations. eLife. 2020;9:e52813.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
43.Bagge N, Hentzer M, Andersen JB, Ciofu O, Givskov M, Høiby N. Dynamics and spatial distribution of β-lactamase expression in Pseudomonas aeruginosa biofilms. Antimicrob Agents Chemother. 2004;48:1168–74.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
44.Allen B, Gore J, Nowak MA. Spatial dilemmas of diffusible public goods. eLife. 2013;2:e01169.PubMed 
PubMed Central 
Article 

Google Scholar 
45.Medaney F, Dimitriu T, Ellis RJ, Raymond B. Live to cheat another day: bacterial dormancy facilitates the social exploitation of β-lactamases. ISME J. 2016;10:778.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.Frost I, Smith WP, Mitri S, San Millan A, Davit Y, Osborne JM, et al. Cooperation, competition and antibiotic resistance in bacterial colonies. ISME J. 2018;12:1582–93.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
47.Estrela SBS. Community interactions and spatial structure shape selection on antibiotic resistant lineages. PLoS Comput Biol. 2018;14:e1006179.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
48.Amanatidou E, Matthews AC, Kuhlicke U, Neu TR, McEvoy JP, Raymond B. Biofilms facilitate cheating and social exploitation of β-lactam resistance in Escherichia coli. npj Biofilms Microbiomes. 2019;5:1–10.CAS 
Article 

Google Scholar 
49.Santos-Lopez A, Marshall CW, Scribner MR, Snyder DJ, Cooper VS. Evolutionary pathways to antibiotic resistance are dependent upon environmental structure and bacterial lifestyle. eLife. 2019;8:e47612.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
50.Seligman SJ, Hewitt WL. Kinetics of the action of ampicillin on Escherichia coli. J Bacteriol. 1963;85:1160–4.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
51.Klementiev AD, Jin Z, Whiteley M. Micron scale spatial measurement of the O2 gradient surrounding a bacterial biofilm in real time. mBio. 2020;11:e02536-20.52.van Tatenhove-Pel RJ, Rijavec T, Lapanje A, van Swam I, Zwering E, Hernandez-Valdes JA, et al. Microbial competition reduces metabolic interaction distances to the low µm-range. ISME J. 2021;15:688–70.53.Kumar RK, Meiller-Legrand T, Alcinesio A, Gonzalez D, Mavridou DA, Meacock OJ, et al. Droplet printing reveals the importance of micron-scale structure for bacterial ecology. Nat Commun. 2021;12:857.54.Gilmore MS, Clewell DB, Ike Y, Shankar N, editors. Enterococci: From Commensals to Leading Causes of Drug Resistant Infection [Internet]. Boston: Massachusetts Eye and Ear Infirmary; 2014.55.Huycke MM, Sahm DF, Gilmore MS. Multiple-drug resistant enterococci: the nature of the problem and an agenda for the future. Emerg Infect Dis. 1998;4:239.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
56.Mohamed JA, Huang DB. Biofilm formation by enterococci. J Med Microbiol. 2007;56:1581–8.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
57.Ch’ng JH, Chong KK, Lam LN, Wong JJ, Kline KA. Biofilm-associated infection by enterococci. Nat Rev Microbiol. 2018;1:82–94.58.Murray BE, Mederski-Samaroj B. Transferable β-lactamase. A new mechanism for in vitro penicillin resistance in Streptococcus faecalis. J Clin Investig. 1983;72:1168–71. https://doi.org/10.1172/JCI111042.CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
59.Rice L, Eliopoulos G, Wennersten C, Goldmann D, Jacoby G, Moellering R. Chromosomally mediated β-lactamase production and gentamicin resistance in Enterococcus faecalis. Antimicrob Agents Chemother. 1991;35:272–6.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
60.Murray BE. Beta-lactmase-producing enterococci. Antimicrob Agents Chemother. 1992;36:2355–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
61.Miller WR, Munita JM, Arias CA. Mechanisms of antibiotic resistance in enterococci. Expert Rev Antiinfect Ther. 2014;12:1221–36.CAS 
Article 

Google Scholar 
62.Dunny GM, Lee LN, LeBlanc DJ. Improved electroporation and cloning vector system for gram-positive bacteria. Appl Environ Microbiol. 1991;57:1194–201.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
63.Aymanns S, Mauerer S, van Zandbergen G, Wolz C, Spellerberg B. High-level fluorescence labeling of gram-positive pathogens. PLoS ONE. 2011;6:e19822.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
64.Zscheck KK, Murray BE. Nucleotide sequence of the β-lactamase gene from Enterococcus faecalis HH22 and its similarity to staphylococcal β-lactamase genes. Antimicrob Agents Chemother. 1991;35:1736–40.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
65.Hallinen KM, Guardiola-Flores KA, Wood KB. Fluorescent reporter plas- mids for single-cell and bulk-level composition assays in E. faecalis. PLoS ONE. 2020;15:e0232539.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
66.AU Levin-Reisman I, AU Fridman O, AU Balaban NQ. Scan- Lag: high- throughput quantification of colony growth and lag time. JoVE 2014;51456.67.Schindelin J, Arganda-Carreras I, Frise EEA. Fiji: an open-source plat- form for biological-image analysis. Nat Methods. 2012;9:676–82.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
68.Eden M. A two-dimensional growth process. Dyn Fractal Surf. 1961;4:223–39.
Google Scholar 
69.Smith WP, Davit Y, Osborne JM, Kim W, Foster KR, Pitt-Francis JM. Cell morphology drives spatial patterning in microbial communities. Proc Natl Acad Sci USA. 2017;114:E280–E286.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
70.Gralka M, Stiewe F, Farrell F, Möbius W, Waclaw B, Hallatschek O. Allele surfing promotes microbial adaptation from standing variation. Ecol Lett. 2016;19:889–98.PubMed 
PubMed Central 
Article 

Google Scholar 
71.Paulose J, Hallatschek O. The impact of long-range dispersal on gene surfing. Proc Natl Acad Sci USA. 2020;117:7584–93.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
72.Kaznatcheev A, Peacock J, Basanta D, Marusyk A, Scott JG. Fibroblasts and alectinib switch the evolutionary games played by non-small cell lung cancer. Nat Ecol Evol. 2019;3:450–6.PubMed 
PubMed Central 
Article 

Google Scholar 
73.Levin B. Frequency-dependent selection in bacterial populations. Philos Trans R Soc Lond B Biol Sci. 1988;319:459–72.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
74.Kaznatcheev A. Two conceptions of evolutionary games: reductive vs effective. bioRxiv. 2017; 231993; https://doi.org/10.1101/231993.75.van Gestel J, Weissing FJ, Kuipers OP, Kovács ÁT. Density of founder cells affects spatial pattern formation and cooperation in Bacillus subtilis biofilms. ISME J. 2014;8:2069–79.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
76.Bamford CH, Tipper C, Compton R. Diffusion-limited reactions, vol. 25. Elsevier; 1985.77.Berg HC, Purcell EM. Physics of chemoreception. Biophys J. 1977;20:193–219.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
78.Christensen H, Martin MT, Waley SG. Beta-lactamases as fully effcient enzymes. Determination of all the rate constants in the acyl-enzyme mechanism. Biochem J. 1990;266:853.CAS 
PubMed 
PubMed Central 

Google Scholar 
79.Hardy LW, Kirsch JF. Diffusion-limited component of reactions catalyzed by Bacillus cereus. Beta-lactamase I. Biochemistry. 1984;23:1275–82.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
80.Dubus A, Ledent P, Lamotte-Brasseur J, Frère JM. The roles of residues Tyr150, Glu272, and His314 in class C β-lactamases. Proteins. 1996;25:473–85.CAS 
PubMed 
PubMed Central 

Google Scholar 
81.Voladri R, Tummuru M, Kernodle DS. Structure-function relationships among wild-type variants of Staphylococcus aureus β-lactamase: importance of amino acids 128 and 216. J Bacteriol. 1996;178:7248–53.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
82.Livermore DM. β-Lactamases: quantity and resistance. Clin Microbiol Infect. 1997;3:4S10–4S19.CAS 
Article 

Google Scholar 
83.Nikaido H, Normark S. Sensitivity of Escherichia coli to various, β-lactams is determined by the interplay of outer membrane permeability and degradation by periplasmic β-lactamases: a quantitative predictive treatment. Mol Microbiol. 1987;1:29–36.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
84.Stewart PS. Diffusion in biofilms. J Bacteriol. 2003;185:1485–91.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
85.Mah TFC, O’Toole GA. Mechanisms of biofilm resistance to antimicrobial agents. Trends Microbiol. 2001;9:34–39.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
86.Lee HH, Molla MN, Cantor CR, Collins JJ. Bacterial charity work leads to population-wide resistance. Nature. 2010;467:82–85. https://doi.org/10.1038/nature09354.CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
87.Meredith HR, Andreani V, Ma HR, Lopatkin AJ, Lee AJ, Anderson DJ, et al. Applying ecological resistance and resilience to dissect bacterial antibiotic responses. Sci Adv. 2018;4:eaau1873.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
88.Udekwu KI, Parrish N, Ankomah P, Baquero F, Levin BR. Functional relationship between bacterial cell density and the effcacy of antibiotics. J Antimicrob Chemother. 2009;63:745–57. https://doi.org/10.1093/jac/dkn554.CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
89.Tan C, Smith RP, Srimani JK, Riccione KA, Prasada S, Kuehn M, et al. The inoculum effect and band-pass bacterial response to periodic antibiotic treatment. Mol Syst Biol. 2012;8:617.90.Karslake J, Maltas J, Brumm P, Wood KB. Population density modulates drug inhibition and gives rise to potential bistability of treatment outcomes for bacterial infections. PLoS Comput Biol. 2016;12:e1005098.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
91.Nadell CD, Foster XJ. Emergence of spatial structure in cell groups and the evolution of cooperation. PLoS Comput Biol. 2010;6:e1000716PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
92.Nadell CD, Drescher K, Foster KR. Spatial structure, cooperation and competition in biofilms. Nat Rev Microbiol. 2016;14:589.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
93.Yuste S, Acedo L, Lindenberg K. Reaction front in an A+BC reaction-subdiffusion process. Phys Rev E. 2004;69:036126.CAS 
Article 

Google Scholar 
94.Grebenkov DS. Diffusion toward non-overlapping partially reactive spherical traps: fresh insights onto classic problems. J Chem Phys. 2020;152:244108.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
95.Adamowicz EM, Flynn J, Hunter RC, Harcombe WR. Cross-feeding modulates antibiotic tolerance in bacterial communities. ISME J. 2018;12:2723–35.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
96.Tanouchi Y, Pai A, Buchler NE, You L. Programming stress-induced altruistic death in engineered bacteria. Mol Syst Biol. 2012;8:626.97.Xie H, Jiao Y, Fan Q, Hai M, Yang J, Hu Z, et al. Modeling three-dimensional invasive solid tumor growth in het- erogeneous microenvironment under chemotherapy. PLoS ONE. 2018;13:e0206292.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
98.Bowness R, Chaplain MA, Powathil GG, Gillespie SH. Modelling the effects of bacterial cell state and spatial location on tuberculosis treatment: insights from a hybrid multiscale cellular automaton model. J Theor Biol. 2018;446:87–100.PubMed 
PubMed Central 
Article 

Google Scholar 
99.Dai X, Xiang S, Li J, Gao Q, Yang K. Development of a colorimetric assay for rapid quantitative measurement of clavulanic acid in microbial samples. Sci China Life Sci. 2012;55:158–63.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
100.Kobayashi S, Arai S, Hayashi S, Sakaguchi T. Simple assay of β- lactamase with agar medium containing a chromogenic cephalosporin, pyridinium-2-azo-p-dimethylaniline chromophore (PADAC). Antimicrob Agents Chemother. 1988;32:1040–5.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar  More

Study siteThis study took place at a temperate reef called the “Chicken Rock” in waters off the coast of North Carolina, USA, between Cape Hatteras and Cape Lookout (Raleigh Bay; Fig. 1). The seafloor of the Chicken Rock is composed of low-relief hardbottom and sand. The Chicken Rock is approximately 37 m deep (Fig. 2) and is an ideal location for this study for three reasons. First, it has a relatively flat seafloor that allows for a high detection rate of acoustically tagged fish49. Second, a high-resolution bathymetric map was available for the area (C. Taylor, National Centers for Coastal Ocean Science). Third, many red snapper occupy the area, allowing us to catch and tag fish relatively easily. Recreational and commercial fishing occurs at the Chicken Rock year-round for a variety of species, but red snapper can only be retained during short open seasons that have occurred periodically since 2010.Data collectionWe quantified the fine-scale movements and distance off bottom for red snapper using VPS (Innovasea, Nova Scotia, Canada). VPS uses a time-difference-of-arrival algorithm to determine the location of coded acoustic transmitters that have been detected by at least three submersible acoustic receivers50. Highly precise fish positions (~ 1 m resolution) are possible if time is synchronized exactly across all receivers, which is accomplished by using sync tags that are either deployed independently throughout the receiver array or built into the receivers themselves. One downside of VPS is that data are not available in real time; receivers must be physically recovered to download data, and then data have to be sent to Vemco to determine fish positions. The advantages of VPS, however, are immense, especially in providing highly precise spatial positions each time acoustic signals are emitted from transmitters. VPS has been used many times to successfully quantify demersal fish movements27,28,49,50,54, and three-dimensional movements can be determined if pressure sensors are built into transmitters23,42.We deployed an array of 20 submersible VR2AR receivers at the Chicken Rock on 17 April 2019. Receivers were deployed in three rows of seven receivers, except for a single receiver in the northeast corner of the grid. Based on previously estimated detection distances of 200–400 m49,55, receivers were separated 200 m from each other, so the entire receiver grid occupied an area of approximately 400 × 1200 m (0.48 km2; Fig. 2). Each receiver was connected to a line between a 36-kg steel weight and a 28-cm diameter plastic float with 8.8 kg of buoyancy, with each receiver positioned approximately 3 m off the seafloor. Each VR2AR included its own sync tag for time synchronization and acoustic release so receivers could be retrieved at the end of the study. A TCM-1 current probe (Lowell Instruments) was attached to each of three receiver buoys spread out across our receiver array (Fig. 2) to collect minute-by-minute current speed and bottom water temperature.We also deployed a reference transmitter (Vemco V13T-1x) in the receiver array on 17 April 2019 (Fig. 2) to calculate sound speed velocity for VPS analyses and quantify positional error of transmitters in the receiver array by comparing its known location to its estimated positions over the course of the study. The reference transmitter was connected to a line with a weight at one end and a buoy at the other, had a 550–650 s random ping interval, and operated on a frequency of 69 kHz.A total of 44 red snapper were tagged in this study. Twenty-three red snapper were tagged on 7 May 2019, nineteen were tagged on 13 August 2019, one was tagged on 30 August 2019, and one was tagged on 22 September 2019 (Table 1). Most of these red snapper (N = 43) were caught via hook-and-line using either circle or J-style hooks, but one red snapper (tagged on 30 August 2019) was caught in a baited fish trap. Fish in good condition (i.e., no visible signs of barotrauma, jaw hooked, active) were tagged externally because external attachment is fast (i.e., greatly reducing surface time56) and externally attached transmitters are detected better than surgically implanted transmitters57. The downside is that transmitter retention is typically lower for externally attached transmitters compared to surgically implanted transmitters.We tagged red snapper with Vemco V13P-1 × transmitters that were 13 mm wide, 46 mm long, weighed 13 g in air, had a 130–230 s pulse interval, a 613 d battery life, and operated on a frequency of 69 kHz. Each transmitter also contained a pressure sensor, which was used to determine the depth of fish for each acoustic signal (accuracy = 1.7 m). Before field work began, stainless steel wire (0.89-mm diameter) was wrapped around the non-transmitting end of the transmitter, glued with marine adhesive (3 M 5200), and covered in heat shrink tubing. Approximately 15 cm of stainless steel wire that extended beyond the transmitter was straightened, and the end was sharpened.Upon capture, red snapper had their head and eyes covered in a wet towel and were measured for total length (mm). The sharpened transmitter wire was inserted laterally through the dorsal musculature of the fish approximately 2.5 cm posterior to, and 2.5 cm below, the insertion of the fish’s first dorsal spine. The wire was pushed laterally through the fish until the transmitter was pulled firmly against the fish’s left side, while the sharpened end emerged from the same spot on the right side of the fish. An aluminum washer was threaded onto the protruding wire, followed by a #1 double sleeve steel crimp, which was crimped onto the wire once the washer and crimp were held firmly on the right side of the fish. The wire beyond the crimp and wet towel were removed, the fish was attached to a weighted SeaQualizer fish release tool, and the fish was descended to a depth of approximately 31 m before being released by the device. The total surface time for each tagged red snapper was approximately 1.5 min.Data analysesWe first assessed whether potential error in red snapper positions could influence study results. For each reference tag position estimated by VPS, we calculated horizontal positional error as the difference between the known reference tag location and its estimated position based on VPS. We visualized daily horizontal positional error of the reference transmitter with a boxplot. Daily values were provided to determine if any changes in positional error occurred over time.Next, we used positional and depth data from fish that were monitored to determine the fate of each individual and classified them based on four events: tag loss, emigration, harvest, or predation48. Fish were assumed to have lost transmitters if the transmitter stopped moving; they were assumed to have emigrated if the transmitter moved to the edge of the receiver array before disappearing. Harvest was assumed if fish disappeared from within the receiver array. Predation (e.g., by sharks) was inferred from VPS data in one of three ways: (1) transmitters moved horizontally much faster than normal red snapper swimming speeds, (2) transmitters moved quickly across a wide range of depths, typically from the bottom to the surface and back, and (3) a reduced frequency of detections, as might be expected for transmitters in the abdominal cavity of a shark. VPS data were censored after the point at which any fish experienced tag loss, harvest, or predation, and only fish with 100 or more spatial positions were included in the analyses.We then estimated movement rates of each fish over time. Movement rate (m s−1) was quantified as the distance moved between each successive pair of spatial positions divided by the time between detections. One challenge with using movement rates is that straight-line movements are assumed between detections, when in reality fish may not move in straight lines. Red snapper were detected on average every 2–4 min, so this issue is less of a problem in our study compared to those using longer time intervals between detections51, but our movement rates can be considered minimum estimates. To further prevent negatively biased movement rate estimates, we excluded movement rate estimates for time intervals longer than 20 min; this decision had negligible effects on results (see Discussion).We also estimated the distance off the seafloor for all detections of acoustically tagged red snapper. We calculated distance off the bottom (m) for each fish position as the depth of the seafloor at that location minus the depth of the fish. We encountered an issue with some transmitters after tag loss whereby depth readings appeared to slowly drift towards shallower readings even though the transmitter was sitting on the bottom and not moving horizontally; in a few instances, this same depth drift issue was detected for transmitters attached to fish alive in the study area (i.e., distance off bottom was greater than zero for long periods of time, which never occurred for red snapper with working pressure sensors). We do not know the reason for these rare instances of depth drift by the pressure sensors, but out of caution we censored depth data for fish whose transmitters provided dubious depth data.We evaluated whether individual differences in movement rates or distance off the bottom were apparent. We created boxplots of movement rate and distance off bottom for each fish in the study, and tested for differences among individuals using a linear model where fish number was included as a categorical variable. We compared the Akaike information criterion (AIC) values of models including fish number with an intercept-only model where fish number was excluded, and models with the lowest AIC value (ΔAIC = 0) were considered the most parsimonious formulations58. Movement rate was positively skewed, so it was log-transformed to improve model fit. Model diagnostics (i.e., quantile–quantile, histogram of residuals, residuals versus linear predictions, response versus fitted values plots) were used to confirm that final models met assumptions of equal variance and normal residuals. We used R version 3.6.359 to carry out all statistical tests and to create all figures.Ideally, we would then test for the effects of environmental conditions and fish size on red snapper horizontal and vertical movements using a single, integrated analysis. However, models accounting for temporal autocorrelation and incorporating individual movement rate estimates from each fish as the response variable (i.e., including fish number as a random effect) did not converge, possibly due to large sample sizes (N = 346,363), so we used mean hourly values instead. The downside of this approach is that fish size had to be evaluated separately from the effects of environmental conditions, as described below. However, note that covariate relationships changed very little across a wide variety of model formulations.We tested for the effects of fish size on movement rate and distance off the bottom using generalized additive models60 (GAMs). GAMs are a regression modeling approach that relate a response variable to a single or multiple predictor variables using nonlinear, linear, or categorical functions. Mean log-transformed movement rate or distance off bottom were the response variables of these models and cubic-spline-smoothed fish total length (mm) was included as the predictor variable. As above, we compared the AIC values of models including fish size with an intercept-only model where fish size was excluded, and the model with the lowest AIC value was selected as the best model.We then assessed the influence of various environmental factors (see below) on red snapper movement rate and distance off bottom using GAMs. For these analyses, choosing the appropriate time scale for binning response and predictor data was critical. Longer time steps (i.e., day) were problematic because response and predictor variables frequently varied over much shorter time frames, while extremely short time steps (i.e., minute) were often lacking response and predictor variable information. Therefore, we used an hourly time step for this procedure. The main concern of using an hourly step is that any particular hourly time bin is likely to be more similar to the time bin nearest in time compared to a randomly selected time bin; in other words, time bins are not truly independent of one another61 (i.e., data are temporally autocorrelated). Not accounting for temporal autocorrelation that is present often leads to a negative bias in estimated regression coefficients and confidence intervals. To account for temporal autocorrelation, we used generalized additive mixed models (GAMMs) that included an autoregressive term for model errors. We used a likelihood ratio test to compare our GAMM to a GAM that did not include autoregressive errors, and in both cases GAMMs were selected over GAMs so they were used for movement and distance off bottom models.We limited our GAMMs to five predictor variables based on previous work. The first predictor variable was time of day, which we included because red snapper movements have been shown to vary over diel periods29. We included time of day (tod) as a categorical variable with three levels: day, crepuscular period, and night. Because sunrise and sunset times varied over the course of our 8-mo study, we defined crepuscular periods as a one hour period of time spanning 30 min before sunrise or sunset to 30 min after sunrise or sunset for each day of the study. Day was defined as 30 min after sunrise to 30 min before sunset, and night was defined as 30 min after sunset to 30 min before sunrise.Bottom water temperature has been shown to be strongly correlated with red snapper movements and home range size28,29, so it was included as our second predictor variable. We calculated bottom water temperature (temp; °C) as the mean bottom temperature measured across the three current probes deployed in the receiver array. Cold bottom water temperatures were observed near the conclusion of our study (December 2019) due to declining air temperatures and water column mixing, but also during periodic upwelling events that occurred from late May through early August. Upwelling is a common oceanographic feature of the region, occurring when upwelling-favorable winds are observed concurrent with the Gulf Stream being in a relatively inshore position62,63. Upwelled water that is cold and nutrient-rich is generally only found near the bottom, which tends to cause phytoplankton blooms near the bottom that decrease water clarity. From preliminary analyses of red snapper VPS data, we observed differing behaviors of fish during periods of upwelling than periods lacking upwelling. Therefore, we developed an upwelling index as our third predictor variable, which was calculated as the difference between the surface water temperature and mean bottom water temperature (upwel; °C). Surface water temperature was not available at the study site, so we obtained hourly surface temperature data from the nearest NOAA buoy (#41159), which was located ~ 85 km southwest of the study site in a similar water depth (Fig. 1). We assume that surface water temperature at the study site could be approximated with data from this buoy, which is a reasonable assumption given surface water temperature and wave heights from this buoy were strongly correlated with values from another buoy (NOAA buoy #41025) ~ 70 km northeast of the study site.The last two predictor variables involved properties of water movement at the seafloor in the study area. The fourth predictor variable was wave orbital velocity (wov; m s−1), which is a measure of the wave-generated oscillatory flow (“sloshing”) of water at the seabed. Wave orbital velocity was included because it was much more strongly correlated with gray triggerfish (Balistes capriscus) movement rates at the Chicken Rock area than either barometric pressure or bottom water temperature43, the latter of which have been shown to be more important for organisms in shallow water64,65. Wave orbital velocity was calculated following Bacheler et al.43 using the properties of surface wave period and height, which were also obtained from NOAA buoy 41159. The last predictor variable included in models was current speed (cur; cm s−1), which was calculated as the mean horizontal current speed from the three current probes deployed on the bottom in the receiver array.The GAMMs were formulated as:$$y = upalpha + f(tod) + s_{1} (temp) + s_{2} (upwel) + s_{3}(wov) + s_{4} (cur) + varepsilon ,$$
(1)
where y is either acoustically tagged red snapper log-transformed movement rate (m s−1) or distance off the bottom (m), α is the intercept, f is a categorical function, s1-4 are cubic spline smoothing functions, and (varepsilon) is the autoregressive error term accounting for temporal autocorrelation in the data.We employed model selection techniques to assess the importance of predictor variables. Specifically, we compared full models that included all five predictor variables to reduced models that included fewer predictor variables. Model comparisons were made using AIC, and models with the lowest AIC value (ΔAIC = 0) were again considered the most parsimonious. Various diagnostics of final models were examined using the “gam.check” function in the mgcv library to ensure model fit was suitable.Given the importance of upwelling to the vertical movements of red snapper (see Results section), we last include results from a conductivity-temperature-depth (CTD) cast taken in the study area from the NOAA Ship Pisces on 29 June 2019 (07:40 EDT), which occurred during a time when bottom upwelling was present. This CTD cast was conducted using a Sea-Bird SBE 9 deployed from the surface to within 1.5 m of the bottom, and depth-specific water temperature and beam transmission data were provided to highlight the vertical extent of upwelling on this particular day. Beam transmission is the fraction of a light source reaching a light detector set a distance away and is a quantitative measure of water clarity; a common feature of upwelling in the region (in addition to cold water) is declining clarity due to increased production within nutrient-rich, upwelled water near the bottom. We combine these water temperature and beam transmission data with a boxplot of red snapper distances off the bottom by hour throughout the same day the CTD cast was taken (29 June 2019).Ethics approvalThe tagging protocol was approved by the Institutional Animal Care and Use Committee (# NCA19-002) of the North Carolina Aquariums on 20 March 2019. All research activities were carried out under a Scientific Research Permit issued to Nathan Bacheler on 10 April 2017 by the Southeast Regional Office of the U.S. National Marine Fisheries Service, in accordance with the relevant guidelines and regulations on the ethical use of animals as experimental subjects. More