Ecology

1.
Sherr BF, Sherr EB, Caron D, Vaulot D, Worden A. Oceanic protists. Oceanography. 2007;20:130–34.
Google Scholar 
2.
Worden AZ, Follows MJ, Giovannoni SJ, Wilken S, Zimmerman AE, Keeling PJ. Rethinking the marine carbon cycle: factoring in the multifarious lifestyles of microbes. Science. 2015;347:1257594.
PubMed  Google Scholar 

3.
Jürgens K, Massana R. Protistan Grazing on Marine Bacterioplankton. In: D.L. Kirchman [ed.], Microbial ecology of the oceans. John Wiley & Sons, Inc; Hoboken, New Jersey, 2008. p. 383–441.

4.
Pernthaler J. Predation on prokaryotes in the water column and its ecological implications. Nat Rev Microbiol. 2005;3:537–46.
CAS  PubMed  Google Scholar 

5.
Boenigk J, Arndt H. Bacterivory by heterotrophic flagellates: community structure and feeding strategies. Ant van Leeuw. 2002;81:465–80.
Google Scholar 

6.
Vørs N, Buck KR, Chavez FP, Eikrem W, Hansen LE, Østergaard JB, et al. Nanoplankton of the equatorial Pacific with emphasis on the heterotrophic protists. Deep-Sea Res II. 1995;42:585–602.
Google Scholar 

7.
Massana R, Guillou L, Díez B, Pedrós-Alió C. Unveiling the organisms behind novel eukaryotic ribosomal DNA sequences from the ocean. Appl Environ Microbiol. 2002;68:4554–58.
CAS  PubMed  PubMed Central  Google Scholar 

8.
Rodríguez-Martínez R, Rocap G, Logares R, Romac S, Massana R. Low evolutionary diversification in a widespread and abundant uncultured protist (MAST-4). Mol Biol Evol. 2012;29:1393–406.
PubMed  Google Scholar 

9.
del Campo J, Balagué V, Forn I, Lekunberri I, Massana R. Culturing bias in marine heterotrophic flagellates analyzed through seawater enrichment incubations. Micro Ecol. 2013;66:489–99.
CAS  Google Scholar 

10.
Caron DA, Alexander H, Allen AE, Archibald JM, Armbrust EV, Bachy C, et al. Probing the evolution, ecology and physiology of marine protists using transcriptomics. Nat Rev Micro. 2017;15:6–20.
CAS  Google Scholar 

11.
Yutin N, Wolf MY, Wolf YI, Koonin EV. The origins of phagocytosis and eukaryogenesis. Biol Direct. 2009;4:9.
PubMed  PubMed Central  Google Scholar 

12.
Keeling PJ. The number, speed, and impact of plastid endosymbioses in eukaryotic evolution. Annu Rev Plant Biol. 2013;64:583–607.
CAS  Google Scholar 

13.
Martin WF, Tielens AGM, Mentel M, Garg SG, Gould SB. The physiology of phagocytosis in the context of mitochondrial origin. Micro Mol Biol Rev. 2017;81:e00008–17.
CAS  Google Scholar 

14.
Rosales C, Uribe-Querol E. Phagocytosis: a fundamental process in immunity. BioMed Res Int. 2017;2017:9042851.

15.
Gotthardt D, Warnatz HJ, Henschel O, Brückert F, Schleicher M, Soldati T. (2002). High-resolution dissection of phagosome maturation reveals distinct membrane trafficking phases. Mol Biol Cell. 2002;13:3508–20.
CAS  PubMed  PubMed Central  Google Scholar 

16.
Niedergang F, Grinstein S. How to build a phagosome: new concepts for an old process. Curr Opin Cell Biol. 2018;50:57–63.
CAS  PubMed  Google Scholar 

17.
Kanehisa M, Sato Y, Furumichi M, Morishima K, Tanabe M. New approach for understanding genome variations in KEGG. Nucleic Acids Res. 2019;47:D590–5.
CAS  PubMed  Google Scholar 

18.
Bozzaro S, Bucci C, Steinert M. Phagocytosis and host-pathogen interactions in Dictyostelium with a look at macrophages. Int Rev Cell Mol Biol. 2008;271:253–300.
CAS  PubMed  Google Scholar 

19.
Jacobs ME, DeSouza LV, Samaranayake H, Pearlman RE, Siu KWM, Klobutcher LA. The Tetrahymena thermophila phagosome proteome. Eukaryot Cell. 2006;5:1990–2000.
CAS  PubMed  PubMed Central  Google Scholar 

20.
Boulais J, Trost M, Landry CR, Dieckmann R, Levy ED, Soldati T, et al. Molecular characterization of the evolution of phagosomes. Mol Syst Biol. 2010;6:423.
PubMed  PubMed Central  Google Scholar 

21.
Lie AAY, Liu Z, Terrado R, Tatters AO, Heidelberg KB, Caron DA. Effect of light and prey availability on gene expression of the mixotrophic chrysophyte Ochromonas sp. BMC Genomics. 2017;18:163.
PubMed  PubMed Central  Google Scholar 

22.
Rubin ET, Cheng S, Montalbano AL, Menden-Deuen S, Rynearson TA. Transcriptomic response to feeding and starvation in a herbivorous dinoflagellate. Front Mar Sci. 2019;6:246.
Google Scholar 

23.
Fenchel T, Patterson DJ. Cafeteria roenbergensis nov. gen., nov. sp., a heterotrophic microflagellate from marine plankton. Mar Micro Food Webs. 1988;3:9–19.
Google Scholar 

24.
Schoenle A, Hohlfeld M, Rosse M, Filz P, Wylezich C, Nitsche F, et al. Global comparison of bicosoecid Cafeteria-like flagellates from the deep ocean and surface waters, with reorganization of the family Cafeteriaceae. Eur J Protistol. 2020;73:125665.
PubMed  Google Scholar 

25.
Keeling PJ, Burki F, Wilcox HM, Allam B, Allen EE, Amaral- Zettler LA, et al. The marine microbial eukaryote transcriptome Sequencing Project (MMETSP): illuminating the functional diversity of eukaryotic life in the oceans through transcriptome sequencing. PLoS Biol. 2014;12:e1001889.
PubMed  PubMed Central  Google Scholar 

26.
Hackl T, Martin R, Barenhoff K, Duponchel S, Heider D, Fischer MG. Four high-quality draft genome assemblies of the marine heterotrophic nanoflagellate Cafeteria roenbergensis. Sci Data. 2020;7:29.
CAS  PubMed  PubMed Central  Google Scholar 

27.
Anderson R, Kjelleberg S, Mcdougald D, Jürgens K. Species-specific patterns in the vulnerability of carbon-starved bacteria to protist grazing. Aquat Micro Ecol. 2011;64:105–16.
Google Scholar 

28.
de Corte D, Paredes G, Yokokawa T, Sintes E, Herndl GJ. Differential response of Cafeteria roenbergensis to different bacterial and archaeal characteristics. Micro Ecol. 2019;78:1–5.
Google Scholar 

29.
Massana R, del Campo J, Dinter C, Sommaruga R. Crash of a population of the marine heterotrophic flagellate Cafeteria roenbergensis by viral infection. Environ Microbiol. 2007;9:2660–69.
CAS  PubMed  Google Scholar 

30.
Logares R, Deutschmann IM, Junger PC, Giner CR, Krabberød AK, Schmidt TSB, et al. Disentangling the mechanisms shaping the surface ocean microbiota. Microbiome 2020;8:55.
PubMed  PubMed Central  Google Scholar 

31.
Giner CR, Pernice MC, Balagué V, Duarte CM, Gasol JM, Logares R, et al. Marked changes in diversity and relative activity of picoeukaryotes with depth in the world ocean. ISME J. 2020;14:437–49.
PubMed  Google Scholar 

32.
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Meth. 2016;13:581–83.
CAS  Google Scholar 

33.
Obiol A, Giner CR, Sánchez P, Duarte CM, Acinas SG, Massana R. A metagenomic assessment of microbial eukaryotic diversity in the global ocean. Mol Ecol Res. 2020;20:718–31.
CAS  Google Scholar 

34.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.
CAS  PubMed  PubMed Central  Google Scholar 

35.
Mangot J-F, Forn I, Obiol A, Massana R. Constant abundances of ubiquitous uncultured protists in the open sea assessed by automated microscopy. Environ Microbiol. 2018;20:3876–89.
CAS  PubMed  Google Scholar 

36.
Lekunberri I, Gasol JM, Acinas SG, Gómez-Consarnau L, Crespo BG, Casamayor EO, et al. The phylogenetic and ecological context of cultured and whole genome-sequenced planktonic bacteria from the coastal NW Mediterranean Sea. Syst Appl Microbiol. 2014;37:216–28.
PubMed  Google Scholar 

37.
Porter KG, Feig YS. The use of DAPI for identifying aquatic microfloral. Limnol Oceanogr. 1980;25:943–48.
Google Scholar 

38.
González JM, Suttle CA. Grazing by marine nanoflagellates on viruses and virus-sized particles: ingestion and digestion. Mar Ecol Prog Ser. 1993;94:1–10.
Google Scholar 

39.
Frost BW. Effects of size and concentration of food particles on the feeding behavior of the marine planktonic copepod Calanus pacificus. Limnol Oceanogr. 1972;17:805–15.
Google Scholar 

40.
Heinbokel JF. Studies on the functional role of tintinnids in the Southern California Bight. I. Grazing and growth rates in laboratory cultures. Mar Biol. 1978;47:177–89.
Google Scholar 

41.
Menden-Deuer S, Lessard EJ. 2000. Carbon to volume relationships for dinoflagellates, diatoms, and other protist plankton. Limnol Oceanogr. 2000;45:569–79.
CAS  Google Scholar 

42.
Picelli S, Faridani OR, Björklund Å, Winberg G, Sagasser S, Sandberg R. Full-length RNA-seq from single cells using Smart-seq2. Nat Protoc. 2014;9:171–81.
CAS  PubMed  Google Scholar 

43.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 2014;30:2114–20.
CAS  PubMed  PubMed Central  Google Scholar 

44.
Langmead B, Salzberg S. Fast gapped-read alignment with Bowtie 2. Nat Meth. 2012;9:357–59.
CAS  Google Scholar 

45.
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.
CAS  PubMed  Google Scholar 

46.
The UniProt Consortium. UniProt: a worldwide hub of protein knowledge. Nucleic Acids Res. 2019;47:D506–15.
Google Scholar 

47.
Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, et al. The Pfam protein families database. Nucleic Acids Res. 2012;40:D290–301.
CAS  PubMed  Google Scholar 

48.
Powell S, Szklarczyk D, Trachana K, Roth A, Kuhn M, Muller J, et al. eggNOG v3.0: orthologous groups covering 1133 organisms at 41 different taxonomic ranges. Nucleic Acids Res. 2012;40:D284–9.
CAS  PubMed  Google Scholar 

49.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinforma. 2011;12:323.
CAS  Google Scholar 

50.
Waterhouse RM, Seppey M, Simão FA, Manni M, Ioannidis P, Klioutchnikov G, et al. BUSCO applications from quality assessments to gene prediction and phylogenomics. Mol Biol Evol. 2018;35:543–48.
CAS  PubMed  Google Scholar 

51.
van Bel M, Proost S, van Neste C, Deforce D, van de Peer Y, Vandepoele K. TRAPID, an efficient online tool for the functional and comparative analysis of de novo RNA-Seq transcriptomes. Genome Biol. 2013;14:R134.
PubMed  PubMed Central  Google Scholar 

52.
Finn RD, Attwood TK, Babbitt PC, Bateman A, Bork P, Bridge AJ, et al. InterPro in 2017-beyond protein family and domain annotations. Nucleic Acids Res. 2017;45:D190–9.
CAS  PubMed  Google Scholar 

53.
Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Meth. 2015;12:59–60.
CAS  Google Scholar 

54.
Van Bel M, Diels T, Vancaester E, Kreft L, Botzki A, Van de Peer Y, et al. PLAZA 4.0: an integrative resource for functional, evolutionary and comparative plant genomics. Nucleic Acids Res. 2018;46:D1190–6.
PubMed  Google Scholar 

55.
McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012;40:4288–97.
CAS  PubMed  PubMed Central  Google Scholar 

56.
Zinger L, Gobet A, Pommier T. Two decades of describing the unseen majority of aquatic microbial diversity. Mol Ecol. 2012;21:1878–96.
PubMed  Google Scholar 

57.
Pernice MC, Forn I, Gomes A, Lara E, Alonso-Sáez L, Arrieta JM, et al. Global abundance of planktonic heterotrophic protists in the deep ocean. ISME J. 2015;9:782–92.
CAS  PubMed  Google Scholar 

58.
Eccleston-Parry JD, Leadbeater BSC. A comparison of the growth-kinetics of 6 marine heterotrophic nanoflagellates fed with one bacterial species. Mar Ecol Prog Ser. 1994;105:167–77.
Google Scholar 

59.
Arndt H, Hausmann K, Wolf M. Deep-sea heterotrophic nanoflagellates of the Eastern Mediterranean Sea: qualitative and quantitative aspects of their pelagic and benthic occurrence. Mar Ecol Prog Ser. 2003;256:45–56.
Google Scholar 

60.
Azam F, Long RA. Sea snow microcosms. Nature 2001;414:495–98.
CAS  PubMed  Google Scholar 

61.
Fenchel T. Ecology of protozoa: The biology of free-living phagotrophic protists. Science Tech Publishers, Madison and Springer-Verlag; Madison, Wisconsin, 1987.

62.
Mestre M, Ruiz-González C, Logares R, Duarte CM, Gasol JM, Sala MM. Sinking particles promote vertical connectivity in the ocean microbiome. Proc Natl Acad Sci USA. 2018;115:E6799–807.
CAS  PubMed  Google Scholar 

63.
Beisser D, Graupner N, Bock C, Wodniok S, Grosmann L, Vos M, et al. Comprehensive transcriptome analysis provides new insights into nutritional strategies and phylogenetic relationships of chrysophytes. PeerJ. 2017;5:e2832.
PubMed  PubMed Central  Google Scholar 

64.
Liu Z, Campbell V, Heidelberg KB, Caron DA. Gene expression characterizes different nutritional strategies among three mixotrophic protists. FEMS Micro Ecol. 2016;92:fiw106.
Google Scholar 

65.
Garba L, Ali MSM, Oslan SN, RNZRB AbdulRahman. Review on fatty acid desaturases and their roles in temperature acclimatisation. J Appl Sci. 2017;17:282–95.
CAS  Google Scholar 

66.
Cheng W, Lin M, Qiu M, Kong L, Xu Y, Li Y, et al. Chitin synthase is involved in vegetative growth, asexual reproduction and pathogenesis of Phytophthora capsici and Phytophthora. Environ Microbiol. 2019;21:4537–47.
CAS  PubMed  Google Scholar 

67.
Rawlings ND, Barrett AJ. Families of cysteine peptidases. Methods Enzymol. 1994;244:461–86.
CAS  PubMed  PubMed Central  Google Scholar 

68.
Rawlings ND, Barrett AJ, Finn R. Twenty years of the MEROPS database of proteolytic enzymes, their substrates and inhibitors. Nucleic Acids Res. 2015;44:D343–50.
PubMed  PubMed Central  Google Scholar 

69.
Baltscheffsky M, Schultz A, Baltscheffsky H. H+-proton-pumping inorganic pyrophosphatase: a tightly membrane-bound family. FEBS Lett. 1999;457:527–33.
CAS  PubMed  Google Scholar 

70.
Labarre A, Obiol A, Wilken S, Forn I, Massana R. Expression of genes involved in phagocytosis in uncultured heterotrophic flagellates. Limnol Oceanogr. 2020;65:S149–60.
CAS  Google Scholar 

71.
Minakami R, Sumimotoa H. Phagocytosis-coupled activation of the superoxide-producing phagocyte oxidase, a member of the NADPH oxidase (nox) family. Int J Hematol. 2006;84:193–98.
CAS  PubMed  Google Scholar  More

1.
IPCC. 2013 Supplement to the 2006 IPCC Guidelines for National Greenhouse Gas Inventories (IPCC, Wetlands, 2014).
2.
Smith, P. et al. In Climate Change 2014: Mitigation of Climate Change. Contribution of Working Group III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (eds. Edenhofer, O. et al.) (Cambridge University Press, Cambridge, UK, 2014).

3.
IPCC. Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC, Geneva, Switzerland, 2014).

4.
Ravishankara, A. R., Daniel, J. S. & Portmann, R. W. Nitrous oxide (N2O): the dominant ozone-depleting substance emitted in the 21st century. Science 326, 123–125 (2009).
CAS  Article  Google Scholar 

5.
Baumert, K., Herzog, T. & Pershing, J. Navigating the Numbers: Greenhouse Gas Data and International Climate Policy (World Resources Institute, Washington, DC, 2005).

6.
Oktarita, S., Hergoualch, K., Anwar, S. & Verchot, L. V. Substantial N2O emissions from peat decomposition and N fertilization in an oil palm plantation exacerbated by hotspots. Environ. Res. Lett. 12, 104007 (2017).
Article  Google Scholar 

7.
Yu, Z., Loisel, J., Brosseau, D. P., Beilman, D. W. & Hunt, S. J. Global peatland dynamics since the last glacial maximum. Geophys. Res. Lett. 37, L13402 (2010).
Google Scholar 

8.
Leifeld, J. & Menichetti, L. The underappreciated potential of peatlands in global climate change mitigation strategies. Nat. Commun. 9, 1071 (2018).
CAS  Article  Google Scholar 

9.
Limpens, J., Heijmans, M. M. P. D. & Berendse, F. In Boreal Peatland Ecosystems (eds Wieder, R. K. & Vitt, D. H.) 195 (Springer, Berlin, 2006).

10.
Joosten, H., Tapio-Biström, M. L. & Tol, S. Peatlands—Guidance for Climate Change Mitigation Through Conservation, Rehabilitation and Sustainable use 2nd edn (FAO & Wetlands International Landscape, ICIMOD, Kathmandu, 2012).

11.
Green, S. M. & Page, S. Tropical peatlands: current plight and the need for responsible management. Geol. Today 33, 174–179 (2017).
Article  Google Scholar 

12.
Moore, T. R. & Clarkson, B. R. Dissolved organic carbon in New Zealand peatlands. NZ J. Mar. Freshw. Res. 41, 137–141 (2007).

13.
Renou-Wilson, F., Barry, C., Müller, C. & Wilson, D. The impacts of drainage, nutrient status and management practice on the full carbon balance of grasslands on organic soils in a maritime temperate zone. Biogeosciences 11, 4361–4379 (2014).
Article  Google Scholar 

14.
Martikainen, P. J., Nykänen, H., Crill, P. & Silvola, J. Effect of a lowered water table on nitrous oxide fluxes from northern peatlands. Nature 366, 51–53 (1993).
CAS  Article  Google Scholar 

15.
Wrage-Mönnig, N. et al. Role of nitrifier denitrification in the production of nitrous oxide revisited. Soil Biol. Biochem. 123, A3–A16 (2018).
Article  Google Scholar 

16.
Pärn, J. et al. Nitrogen-rich organic soils under warm well-drained conditions are global nitrous oxide emission hotspots. Nat. Commun. 9, 1135 (2018).
Article  Google Scholar 

17.
Repo, M. E. et al. Large N2O emissions from cryoturbated peat soil in tundra. Nat. Geosci. 2, 189–192 (2009).
CAS  Article  Google Scholar 

18.
Klemedtsson, L., Von Arnold, K., Weslien, P. & Gundersen, P. Soil C/N ratio as a scalar parameter to predict nitrous oxide emissions. Global Change Biol. 11, 1142–1147 (2005).
Article  Google Scholar 

19.
Leppelt, T. et al. Nitrous oxide emission budgets and land-use-driven hotspots for organic soils in Europe. Biogeosciences 11, 6595–6612 (2014).
Article  Google Scholar 

20.
Petersen, S. O. et al. Annual emissions of CH4 and N2O, and ecosystem respiration, from eight organic soils in western Denmark managed by agriculture. Biogeosciences 9, 403–422 (2012).
CAS  Article  Google Scholar 

21.
Leifeld, J. Distribution of nitrous oxide emissions from managed organic soils under different land uses estimated by the peat C/N ratio to improve national GHG inventories. Sci. Total Environ. 631–632, 23–26 (2018).
Article  Google Scholar 

22.
van Beek, C. L. et al. Emissions of N2O from fertilized and grazed grassland on organic soil in relation to groundwater level. Nutr. Cycling Agroecosyst. 86, 331–340 (2010).
Article  Google Scholar 

23.
Maljanen, M. et al. Afforestation does not necessarily reduce nitrous oxide emissions from managed boreal peat soils. Biogeochemistry 108, 199–218 (2012).
CAS  Article  Google Scholar 

24.
Liimatainen, M. et al. Factors controlling nitrous oxide emissions from managed northern peat soils with low carbon to nitrogen ratio. Soil Biol. Biochem. 122, 186–195 (2018).
CAS  Article  Google Scholar 

25.
Säurich, A., Tiemeyer, B., Dettmann, U. & Don, A. How do sand addition, soil moisture and nutrient status influence greenhouse gas fluxes from drained organic soils? Soil Biol. Biochem. 135, 71–84 (2019).
Article  Google Scholar 

26.
Laine, J. et al. Effect of water-level drawdown on global climatic warming: northern peatlands. Ambio 25, 179–184 (1996).
Google Scholar 

27.
Laudone, G. M. et al. A model to predict the effects of soil structure on denitrification and N2O emission. J. Hydrol. 409, 283–290 (2011).
CAS  Article  Google Scholar 

28.
Wu, L. et al. Simulation of nitrous oxide emissions at field scale using the SPACSYS model. Sci. Total Environ. 530–531, 76–86 (2015).
Article  Google Scholar 

29.
Couwenberg, J. et al. Assessing greenhouse gas emissions from peatlands using vegetation as a proxy. Hydrobiologia 674, 67–89 (2011).
CAS  Article  Google Scholar 

30.
Liu, H., Zak, D., Rezanezhad, F. & Lennartz, B. Soil degradation determines release of nitrous oxide and dissolved organic carbon from peatlands. Environ. Res. Lett. 14, 094009 (2019).
CAS  Article  Google Scholar 

31.
Tanneberger, F., Joosten, H., Moen, A. & Whinam, J. In Mires and Peatlands of Europe—Status, Distribution and Conservation (eds Joosten, H., Tanneberger, F. & Moen, A.) 173–196 (Schweizerbart Science Publishers, Stuttgart, 2017).

32.
Renou-Wilson, F. et al. Rewetting degraded peatlands for climate and biodiversity benefits: results from two raised bogs. Ecol. Eng. 127, 547–560 (2019).
Article  Google Scholar 

33.
Lamers, L. P. M. et al. Ecological restoration of rich fens in Europe and North America: From trial and error to an evidence-based approach. Biol. Rev. Camb. Philos. Soc. 90, 182–203 (2015).
Article  Google Scholar 

34.
Tiemeyer et al. A new methodology for organic soils in national greenhouse gas inventories: Data synthesis, derivation and application. Ecol. Indic. 109, 105838 (2020).
CAS  Article  Google Scholar 

35.
Wilson et al. Multiyear greenhouse gas balances at a rewetted temperate peatland. Glob. Change Biol. 22, 4080–4095 (2016).
Article  Google Scholar 

36.
Evans, C. et al. Implementation of an Emission Inventory for UK Peatlands. Report to the Department for Business, Energy and Industrial Strategy, Centre for Ecology and Hydrology, Bangor. 88 (2017).

37.
Günther, A. et al. Prompt rewetting of drained peatlands reduces climate warming despite methane emissions. Nat. Commun. 11, 1644 (2020).
Article  Google Scholar 

38.
Nykänen, H., Alm, J., Lång, K., Silvola, J. & Martikainen, P. J. Emissions of CH4, N2O and CO2 from a virgin fen and a fen drained for grassland in Finland. J. Biogeogr. 22, 351–357 (1995).
Article  Google Scholar 

39.
Drösler, M. et al. Klimaschutz furch Moorschutz in der Praxis (Thünen-Institut fur Agrarklimaschutz, Brauschweig, Germany, 2013).

40.
Mojeremane, W., Rees, R. M. & Mencuccini, M. The effects of site preparation practices on carbon dioxide, methane and nitrous oxide fluxes from a peaty gley soil. Forestry 85, 1–15 (2012).
Article  Google Scholar 

41.
Pronger, J., Schipper, L. A., Hill, R. B., Campbell, D. I. & McLeod, M. Subsidence rates of drained agricultural peatlands in New Zealand and the relationship with time since drainage. J. Environ. Qual. 43, 1442 (2014).
CAS  Article  Google Scholar 

42.
Hume, N. P., Fleming, M. S. & Horne, A. J. Plant carbohydrate limitation on nitrate reduction in wetland microcosms. Water Res. 36, 577–584 (2002).
CAS  Article  Google Scholar 

43.
Höper, H. et al. In Peatlands and Climate Change (ed. Strack, M.) 182–210 (International Peat Society, Jyväskylä, Finland, 2008).

44.
Davidson, E. A. The contribution of manure and fertilizer nitrogen to atmospheric nitrous oxide since 1860. Nat. Geosci. 2, 659–662 (2009).
CAS  Article  Google Scholar 

45.
Tiemeyer, B. et al. High emissions of greenhouse gases from grasslands on peat and other organic soils. Global Change Biol. 22, 4134–4149 (2016).
Article  Google Scholar 

46.
Andersen, R. et al. An overview of the progress and challenges of peatland restoration in Western Europe. Restor. Ecol. 25, 271–282 (2016).
Article  Google Scholar 

47.
Lugato, E., Paniagua, L., Jones, A., de Vries, W. & Leip, A. Complementing the topsoil information of the Land Use/Land Cover Area Frame Survey (LUCAS) with modelled N2O emissions. PLoS ONE 12, e0176111 (2017).
Article  Google Scholar 

48.
Xu, J., Morrisa, P. J., Liu, J. & Holden, J. PEATMAP: Refining estimates of global peatland distribution based on a meta-analysis. Catena 160, 134–140 (2018).
Article  Google Scholar 

49.
Pflugmacher, D., Rabe, A., Peters, M. & Hostert, P. Mapping pan-European land cover using Landsat spectral-temporal metrics and the European LUCAS survey. Remote Sens. Environ. 221, 583–595 (2019).
Article  Google Scholar 

50.
Hierderer, R. EFSA Spatial Data Version 1.1, Data Properties and Processing (Publication Office of the European Union, Luxembourg, 2012).

51.
Jones, R. J., Hiederer, R., Rusco, E. & Montanarella, L. Estimating organic carbon in the soils of Europe for policy support. Eur. J. Soil Sci. 56, 655–671 (2005).
CAS  Article  Google Scholar 

52.
Joosten, H., Tannenberger, F. & Moen, A. Mires and Peatlands of Europe (Schweizerbart Science Publishers, Stuttgart, Germany, 2017).

53.
Liu, H., Price, J. S., Rezanezhad, F. & Lennartz, B. Century-scale shifts in peat hydro-physical properties as induced by drainage Water Resource Research (2020).

54.
Figueres, C. et al. Three years to safeguard our climate. Nature 546, 593–595 (2017).
CAS  Article  Google Scholar  More

We first describe the development of the method and its validation under controlled laboratory settings for both Synechococcus and Prochlorococcus host-phage systems. Then we apply the method to quantify the contribution of T4-like and T7-like cyanophage infection to the mortality of Prochlorococcus over the diel cycle in the upper mixed layer of the North Pacific Subtropical Gyre during the summer of 2015.
Development of the iPolony method for the detection of viral infection
Our goal was to develop a sensitive, high-throughput method amendable to analysis of field samples that can assess infection by virus families of interest. The polony method (named for PCR colonies that form during the reaction) is a direct PCR-based method that detects virus DNA inside infected cells. The original polony method was developed by Mitra and Church [45] and was subsequently adapted for the enumeration of free virus particles from the two major cyanophage lineages, the T4-like cyanomyoviruses [34] and the T7-like cyanopodoviruses [33]. Here we describe further development of the method to simultaneously detect viral infection in thousands of individual cells embedded in a solid-phase gel in a high-throughput manner. We term this method for the quantification of infected cells ‘iPolony’.
The iPolony method has two steps (Fig. 1). In the first step, infection is arrested by fixation with glutaraldehyde, and target cells are isolated and concentrated (Fig. 1a). In its use here, Prochlorococcus and Synechococcus are sorted by flow cytometry based on forward angle light scattering, a proxy for cell size, and the autofluorescence of chlorophyll a and phycoerythrin, respectively. The concentration of cells is then quantified to calculate the number of cells analyzed in the downstream molecular analysis. In the second step, the cells are screened for the presence of intracellular virus DNA using the polony method (Fig. 1b). Polyacrylamide gels are cast with embedded sorted cells and a 5′-acrydite-modified primer to anchor the primer to the gel. PCR reagents are diffused into the gels. Degenerate PCR primers target a signature gene shared by the virus group of interest, in this case the DNA polymerase gene for the T7-like cyanopodoviruses [33] and the portal protein gene (g20) for the T4-like cyanomyoviruses [34]. Prior to thermal cycling, cells are permeabilized to allow amplification of virus DNA with an in-gel heat lysis step (Supplementary text, Supplementary Fig. 3). During thermal cycling, a cell that contained virus DNA results in an anchored sphere of amplification or polony. After amplification, polonies in gels are hybridized with fluorescently labeled degenerate cyanophage group-specific probes, and gels are scanned with a microarray scanner. Each individual cell is counted as infected whether there is a single or multiple virus genome copies, which enables quantification of infection in a presence–absence manner per cell. Percent infection is calculated by dividing the number of polonies by the number of cells in the gel.
Fig. 1: iPolony: a polony method for quantifying virally infected cyanobacteria.

a First, Prochlorococcus and Synechococcus are sorted based on size and their autofluorescence properties using a flow cytometer from fixed samples. b Then, thousands of sorted cells per slide are screened for the presence of intracellular viral DNA using a solid-phase PCR polony method. Percent infection is determined based on the fraction of input cells that resulted in polonies at the end of the analysis.

Full size image

Validation of the iPolony method in controlled virus growth experiments
The method was first developed and tested on model lab systems under known conditions. We assessed the ability of the iPolony method to detect virus DNA inside cyanobacterial cells throughout the infection cycle. Two model systems were used: Synechococcus sp. strain WH8109 infected with the T7-like cyanophage, Syn5, and Prochlorococcus sp. strain MIT9515 infected with the T4-like cyanophage, S-TIM4. Virus DNA was detected inside cells at all stages of the infection cycle (Fig. 2). The efficiency of detection was lower prior to virus genome replication and was 43% for Syn5 and 11% for S-TIM4 (Fig. 2a–d). Detection of infected cells rose with the onset of DNA replication, reaching 86% of maximal infection levels based on host gDNA degradation, midway through genome replication for Syn5 and at the end of genome replication for S-TIM4 (Fig. 2a–d). Since more than a single phage can enter cells in these experiments, we verified that single gene copies can be detected inside cells for a single copy host gene, rbcL, in Synechococcus WH8109, as well as for E. coli carrying a single plasmid with a cyanophage g20 copy (see Supplementary results and discussion). These findings indicate the method is sensitive enough to detect infection throughout the entire infection process even when a single molecule of phage DNA is present prior to genome replication and reaches maximal levels after the onset of DNA replication.
Fig. 2: The iPolony method detects viral infection throughout the infection cycle.

Cultures of Synechococcus sp. strain WH8109 infected by the T7-like cyanopodovirus, Syn5 (a, c, e) and Prochlorococcus sp. strain MIT9515 infected by the T4-like cyanomyovirus, S-TIM4 (b, d, f) at MOI = 3. a, b Percent infection was determined using the polony method over the infection cycle. c, d Virus DNA replication (solid lines) and host genomic DNA degradation (dashed lines) were assessed by qPCR in infected cultures. Host and virus DNA concentrations were normalized to initial or maximum concentrations, respectively. Shaded regions indicate the period of virus genome replication. Lysis was assessed from an increase in plaque forming units measured by the plaque assay for Syn5 (e) or the appearance of extracellular virus DNA measured by qPCR for S-TIM4 (f). Note that Syn5 infections shown in (c) and (e) were not synchronized at 5 min post infection as in (a) and are shown for comparison of the timing of different phases of infection. Average and standard deviation of biological triplicates are shown in all panels.

Full size image

Next, we assessed whether the method could accurately detect varying levels of infection by exposing Synechococcus sp. strain WH8109 to different numbers of the Syn5 phage, such that the infective virus-to-host ratios (multiplicities of infection, MOIs) ranged from 0.1 to 3, spanning infection percentages that theoretically range from 10 to 95% based on encounter theory estimates (Supplementary Table 3). The percent infection determined by the iPolony method was significantly and positively correlated with the MOI (Fig. 3a) (F = 143.1, R2 = 0.87, p  More

We primarily relied on the “tidyverse”61, “metafor”30, “lmerTest”62, “effects”63, “maptools”64, “sf”65, and “raster”66 packages of the R software67 for organizing and analyzing data, and visualizing the results.
Data assembly
To support the use of meta-analysis to quantify the effects of species richness on plant litter decomposition (i.e., species-mixing effects), we searched the literature using the ISI Web of Science (WoS) database (https://clarivate.com/products/web-of-science/; up to and including December 2018). We used a combination of “decomposition” AND “litter”. These keywords matched 12,278 publications. We then reduced this list of the literature with the keywords of “mix* litter,” OR “litter diversity,” OR “litter mix*,” OR “litter species diversity,” OR “litter species richness,” OR “species mix*,” OR “mix* plant litter,” OR “multi* species litter,” OR “mixing effect*,” resulting in 416 publications. We also searched the literature using Scopus (https://www.scopus.com/home.uri) and Google Scholar (https://scholar.google.ca/) and the latter set of keywords, resulting in 765 publications. An anonymous expert identified 15 additional publications. We decided not to include reports in the grey literature among the papers we found, as we could not verify whether, as in all WoS-listed papers, the papers had undergone independent peer review as a quality check. We read through the papers carefully to select those that focused on quantifying the effects of litter diversity on decomposition rates based on a litter-bag experiment using both mixed- and mono-species litter. We focused on the publications that reported values of either mass loss (or mass remaining) or a decomposition rate constant k that quantified the decomposition rate. Specifically, some papers for litter-mixing effects reported only for the additive effects of mixture and had no reports on mass loss or the constant k. In some cases, no sample size nor standard deviations/errors were reported. Mass loss or standard deviations/errors values were reported only for mono-species or multi-species litter. These publications were not included in the present analyses. In cases that had no report on means but instead showed median values, we estimated means and standard errors based on individual data points or percentile values of boxplots. Note that, standard deviations/errors were often invisible because they were too small to read and thus hidden by their data points; in this case, a conservative approach was adopted by using possible maximum values of deviations/errors (i.e., the size of these symbols). Based on these criteria, we extracted information required for our meta-analysis (see below) from a total of 151 studies (Supplementary Data 1). Out of these 151 studies, 131 and 45 studies measured mass loss and the decomposition rate constant k, respectively. Also see the PRISMA work flow diagram (Supplementary Fig. 6).
For the selected studies, we identified the sets of comparisons between mixed- and mono-species litter bags in each publication. That is, a single study could have multiple treatments that focused on different litter substrates, using different mesh sizes for litter bags, having multiple experiments in different locations, changing abiotic conditions, and so on. The decomposition rate should be comparable in a given treatment under the same set of environmental conditions except for differences caused by a different number of plant species. Even for a given treatment, decomposition rates were often reported multiple times in the literature; in such cases, we only included comparisons for litter retrieved at the same time (i.e., after the same incubation period). Note that in a given treatment, the same mono-species litter could be used for multiple comparisons; for instance, mass loss from a litter bag that contained a three-species mixture (species A, B, and C simultaneously) can be compared with the mass loss in three different mono-species litter bags (only species A, B, or C). This could lead to issues related to pseudo replication, which we considered based on a multilevel random effects meta-analysis (see Data analyses). As a result, we identified 6535 comparisons (across 1949 treatments from the 131 studies) for the values of mass loss and 1423 comparisons (across 504 treatments from the 45 studies) for values of the decomposition rate constant k.
After identifying the sets of comparisons, we extracted the sample size (number of litter-bag replicates, n), and the mean and standard deviation (SD) of the decomposition rate from the main text, from any tables and figures, and from the supplemental materials of the selected studies. If standard errors or 95% confidence intervals (CI) were given, we transformed them into SD values. If only figures were given, we used version 3.4 of the Webplot-digitizer software (https://automeris.io/WebPlotDigitizer/) to extract these parameters from the graphs. We also recorded the study location (longitude, latitude, elevation), climatic region (subarctic, boreal, temperate, subtropical, tropical, or other), ecosystem type (forest, grassland, shrubland, desert, wetland, stream, seagrass, lake, or ex situ microcosm), and litter substrate (leaf, root, stem, branch, straw, or other). Microcosm studies were further divided into two categories: terrestrial or aquatic. The former and latter, respectively, placed litter bags on the soil and in water. For litter substrate, most studies used terrestrial plants, but two studies used macrophytes for their decomposition experiment; these data were included in the main analysis and excluded from the subsequent analyses that focused on leaf litter. Note that authors reported the decomposition rate for mixed-species litter bags in different ways, as the values for all species together or for individual species in the same bag. Removing data based on the latter classification had little effect on our results, so we retained that data. Because of the limited data availability, we did not consider the potential influences of species richness (here, the number of litter species; more than 74% of the comparisons used two- or three-species mixtures) and instead considered species richness as a random term (see Data analyses). Based on the geographic locations of the studies, we estimated their present bioclimatic conditions based on the WorldClim database (www.worldclim.org).
Species-mixing effects on litter decomposition
We calculated the unbiased standardized mean difference (Hedges’ d)30 of the decomposition rate between the mean values for the mixed- and mono-species litter. Hedges’ d is a bias-corrected, unit-free index that expresses the magnitude of the deviation from no difference in the response variable between comparisons. Note that many studies performed multiple comparisons for the decomposition rate between mixed- and mono-species litter, potentially causing issues of pseudo-replication and non-independence in the dataset, as is often the case for ecological syntheses68. We thus applied a multilevel random effects meta-analysis30,69 to account for this problem. In a random effects model, the effect sizes for individual comparisons are weighted by the sum of two values: the inverse of the within-study variance and the between-study variance. The multilevel model can also account for a nested structure in the dataset (in which different treatments and multiple comparisons were nested within each study) and is thus appropriate for dealing with non-independence within a dataset. To calculate the effect sizes, we defined their values to be positive for comparisons in which decomposition was faster in mixed-species litter than in mono-species litter and negative when mono-species litter decomposed faster. This is based on a species gain perspective. Note that some have mono- and two-species mixtures, and others could have for instance from 1 to 16 species in their experiments. Considering mono-species litter as a control is therefore the only way to be consistent across all studies. That is, quantification based on a species loss perspective requires us to obtain the effect sizes using mixed-species litter as a control. In this case, different studies have different levels of richness for a control, making it impossible to have a quantification in a standardized way. We first calculated the effect sizes based on mass loss and the decomposition rate constant k (Fig. 1b, c); because of the limited data availability for the constant k, we only calculated the effect sizes for the entire dataset and the subset of data for leaf litter. For mass loss, we further calculated the effect sizes for the different ecosystem types in different climatic regions; subsets of the data that originated from at least three different studies were used. We conducted a multilevel mixed effects meta-regression by considering random effects due to non-independence among some data points, resulting from a nested structure of the dataset that individual data points were nested within a treatment and treatments were also nested within a study. We used the Q statistic for the test of significance.
Note that, because publication bias is a problem in meta-analysis, we visually evaluated the possibility of such a potential bias by plotting the values of the effect sizes and their variances against sample size. If there is no publication bias, studies with small sample sizes should have an increased sampling error relative to those with larger sample sizes, and the variance should decrease with increasing sample size. In addition, the effect size should be independent of the sample size. Also, there should be large variation in effect sizes at the smallest sample sizes. Prior to the meta-analysis, we confirmed that these conditions existed (Supplementary Fig. 1), and found no publication bias large enough to invalidate our analysis.
We also visually confirmed that the datasets were normally distributed based on normal quantile plots. Furthermore, following the method of Gibson et al.70, we randomly selected only one comparison per treatment and then calculated the effect sizes for the decomposition rate (based on mass loss); we repeated this procedure 10,000 times (with replacement) and found that the species-mixing effects on decomposition were significantly positive (0.247 ± 0.045 for the mean ± 95% CI; Supplementary Fig. 2). We thus confirmed that the overall results were not affected by a publication bias or by non-independence of the dataset. We additionally focused on a subset of the entire dataset, which reported mass loss data for each component species from mixtures and thus had comparisons between mass loss of mono-species litter and that of the same species in a mixture (e.g., mass loss of mono-species litter for species A, B, and C was only compared with that for species A, B, and C within a three-species mixture, respectively). We found the results for this subset of the data (Supplementary Fig. 3) almost identical to the main results (Fig. 2a, b). Note that because of the limited data availability, we did not consider the effects of using different equation forms to estimate the rate constant k and instead considered these among-study differences based on a random term in our multilevel meta-analysis. We also calculated the influence of incubation period (days) on the effect size for the decomposition rate (based on mass loss of leaf litter) for the subsets of the dataset using a mixed effects meta-regression30 (Supplementary Fig. 4). We limited this analysis to studies that retrieved litter bags at least two different time points. To account for non-independence of data points, the study identity was included as a random term.
Impacts of global change drivers on litter decomposition
Litter decomposition is primarily controlled by three factors: the environment (e.g., climate), the community of decomposers, and litter traits46,47,48. By relying on the dataset of Makkonen et al.50, who conducted a full reciprocal litter transplant experiment with 16 plant species that varied in their traits and origins (four forest sites from subarctic to tropics), we aimed to obtain a standardized equation to estimate how decomposition rate changed along a climatic gradient (hereafter, the “standardized climate–decomposition relationship”). First, by considering litter decomposition rates at the coldest location (i.e., the subarctic site) as a reference (control), we calculated the effect sizes (i.e., standardized mean difference based on Hedges’ d) for their decomposition rate constant k for all possible comparisons (i.e., between data from the coldest site and the value, one at a time, for the three other warmer sites). As explained above, the effect sizes were calculated only for the comparisons between decomposition rates of litter that were obtained using the same protocol (litter species, origin, and decomposers). The effect sizes therefore represent the climatic effects on litter decomposition after removing the effects of the decomposer community and litter traits.
We obtained the bioclimatic variables for these four study sites from the WorldClim database, and then modelled the standardized climate–decomposition relationships using a multilevel mixed effects meta-regression30,69. Annual mean temperature, the mean temperature of the wettest quarter, or precipitation of the driest quarter were used as an explanatory variable, with the protocol as a random term. This allowed us to evaluate how decomposition rate can be altered by climate along a latitudinal gradient from tropics to subarctic. We then calculated the decomposition rate constant k for our dataset in exactly the same manner as Makkonen et al.50 and calculated the effect sizes; because of limited datasets for most biomes, we performed this analysis only for the subset of our dataset that we obtained for forest biomes, and we used only leaf litter (57 studies) so the results would be comparable with those of Makkonen et al.50. We used these effect sizes and the standardized climate–decomposition relationships to convert the species-mixing effects into climatic effects; that is, based on the temperature or precipitation required to alter the effect size to the same magnitude as the litter diversity effect, we estimated the climate-equivalency of the diversity effects. This means that we relied on the slope of the relationships to quantify the diversity effects, as has been done in previous analyses38,71. We also compared these estimates to the future projections of changes in decomposition rates in response to climate change (with estimates for the 2070s for two of the representative concentration pathway scenarios of CO2 used in the 5th Climate Model Intercomparison Project; CMIP5 RCP 2.6 and 8.5) for the locations of the original studies (the coordinates of individual studies were collected using the Google Map; www.google.com). This comparison was aimed at quantifying the impacts of different global change drivers (biodiversity and climate change) on decomposition processes. Note that we additionally conducted the above analysis after excluding data with no species identity information for mono and/or mixed-species litter bags, those that measured ash-free dry mass loss (this exclusion was to be consistent with the procedure used in Makkonen et al.50), or both of them (Table S1). As an additional confirmation, we also relied on an alternative method used by Hooper et al.45 to convert the species-mixing effects into climatic effects (mean annual temperature), and compared the impacts of different global change drivers on the decomposition rate (Supplementary Fig. 5).
We further analyzed the potential changes in decomposition resulting from litter diversity and climate change. For the dataset of Makkonen et al.50, we first modelled the relationship between climate (annual mean temperature) and the decomposition rate constant k using an LMM with the protocol as a random term. The constant k was log-transformed to improve homoscedasticity. With this equation and the values of annual mean temperature at the 57 study locations in the forest biomes included in our dataset, we then estimated the expected values of k at these study locations (the present k values). Note that because the experiment of Makkonen et al.50 had no mixed-species litter, the expected k using the above LMM was used primarily to estimate the decomposition rate of mono-species litter at a given annual mean temperature. By combining this estimate with the aforementioned estimate of the climate-equivalency of the litter diversity effect (for annual mean temperature; Fig. 3a), it was possible to estimate the potential changes in the k that resulted from increasing litter diversity from mono- to mixed-species litter at the 57 forest study locations. Specifically, the diversity effect converted into the temperature effect (Fig. 3a) was added to present estimates of annual mean temperature at these study locations, and these values of temperature change were used as an explanatory variable to project changes in k based on the above LMM (the projected k values). In other words, we quantified the projected k values under a scenario in which we brought mono-species litter to areas with a warmer annual mean temperature equivalent to the temperature increase in the climate-equivalency effect. We then converted the present and projected values of the decomposition rate constant k to a mass loss per unit time, making it possible to calculate percentage changes in the decomposition rate due to litter diversification. Furthermore, using the above LMM for the temperature–decomposition relationship with the projected changes in mean annual temperature (based on the CMIP5 RCP 2.6 and 8.5 scenarios) as an explanatory variable, we projected future changes in the decomposition rate constant k at these study locations.
Again, we converted these estimates to a litter mass loss to obtain the percentage changes in the decomposition rate that would result from climate change. The rationale for using the temperature–decomposition relationship was as follows: A similar magnitude of increase in k can, in absolute terms, have different influences on decomposition rates at different locations with a different k. For instance, a change in the value of k from 0.1 to 0.2 and a change from 3.0 to 3.1 are substantially different when considered based on litter mass loss or litter amount remaining after decomposition. We thus included annual mean temperature in our model for quantifying the effects of litter diversity on decomposition. More

1.
Katz, B. J. Controlling factors on source rock development—A review of productivity, preservation, and sedimentation rate. Depos. Org. Sediments Model. Mech. Consequences SEPM Spec. Publ. 82, 7–16 (2005).
CAS  Google Scholar 
2.
Tyson, R. V. The, “productivity versus preservation” controversy: Cause, flaws, and resolution. Depos. Org. Sediments Model. Mech. Consequences SEPM Spec. Publ. 82, 17–33 (2005).
CAS  Google Scholar 

3.
Harris, N. B., Freeman, K. H., Pancost, R. D., White, T. S. & Mitchell, G. D. The character and origin of lacustrine source rocks in the Lower Cretaceous synrift section, Congo Basin, west Africa. Am. Assoc. Pet. Geol. Bull. 88, 1163–1184 (2004).
CAS  Google Scholar 

4.
Talbot, M. R., Filippi, M. L., Jensen, N. B. & Tiercelin, J.-J. An abrupt change in the African monsoon at the end of the Younger Dryas. Geochem. Geophys. Geosyst. 8(3), Q03005. https://doi.org/10.1029/2006GC001465 (2007).
ADS  Article  Google Scholar 

5.
Scholz, C. A. et al. Scientific drilling in the Great Rift Valley: The 2005 Lake Malawi Scientific Drilling Project—An overview of the past 145,000 years of climate variability in Southern Hemisphere East Africa. Palaeogeogr. Palaeoclimatol. Palaeoecol. 303, 3–19 (2011).
Google Scholar 

6.
Ellis, G. S., Katz, B. J., Scholz, C. A. & Swart, P. K. Organic sedimentation in modern lacustrine systems: A case study from Lake Malawi, East Africa. In Paying Attention to Mudrocks: Priceless! (eds. Larsen, D. et al.) 515, (Geological Society of America, 2015).

7.
Ivory, S. J. et al. East African weathering dynamics controlled by vegetation-climate feedbacks. Geology 45, 823–826 (2017).
ADS  Google Scholar 

8.
Lambiase, J. & Morley, C. Hydrocarbons in rift basins: The role of stratigraphy. Philos. Trans. R. Soc. London. Ser. A Math. Phys. Eng. Sci. 357, 877–900 (1999).
ADS  Google Scholar 

9.
Huc, A. Y., Le Fournier, J., Vandenbroucke, M. & Bessereau, G. Northern Lake Tanganyika—An example of organic sedimentation in an anoxic rift lake. In Lacustrine Basin Exploration: Case Studies and Modern Analogs (ed. Katz B. J.) 50, 169–185 (American Association of Petroleum Geologists, 1990).

10.
Katz, B. J. A survey of rift basin source rocks. Geol. Soc. Lond. Spec. Publ. 80, 213 (1995).
ADS  Google Scholar 

11.
Lyons, R. P., Scholz, C. A., Buoniconti, M. R. & Martin, M. R. Late Quaternary stratigraphic analysis of the Lake Malawi Rift, East Africa: An integration of drill-core and seismic-reflection data. Palaeogeogr. Palaeoclimatol. Palaeoecol. 303, 20–37 (2011).
Google Scholar 

12.
Harris, N. B. & Tucker, G. E. Soils, slopes and source rocks: Application of a soil chemistry model to nutrient delivery to rift lakes. Sediment. Geol. 323, 31–42 (2015).
ADS  CAS  Google Scholar 

13.
Talbot, M. R. & Johannessen, T. A high resolution palaeoclimatic record for the last 27,500 years in tropical West Africa from the carbon and nitrogen isotopic composition of lacustrine organic matter. Earth Planet. Sci. Lett. 110, 23–37 (1992).
ADS  CAS  Google Scholar 

14.
Talbot, M. R. The origins of lacustrine oil source rocks: Evidence from the lakes of tropical Africa. Geol. Soc. Lond. Spec. Publ. 40, 29–43 (1988).
ADS  Google Scholar 

15.
Petersen, H. I. et al. World-class Paleogene oil-prone source rocks from a cored lacustrine syn-rift succession, Bach Long Vi Island, Song Hong Basin, offshore northern Vietnam. J. Pet. Geol. 37, 373–389 (2014).
CAS  Google Scholar 

16.
Nytoft, H. P. et al. Novel saturated hexacyclic C34 and C35 hopanes in lacustrine oils and source rocks. Org. Geochem. 87, 107–118 (2015).
CAS  Google Scholar 

17.
Hovikoski, J. et al. Density-flow deposition in a fresh-water Lacustrine rift basin, Paleogene Bach Long Vi Graben, Vietnam. J. Sediment. Res. 86, 982–1007 (2016).
ADS  CAS  Google Scholar 

18.
Fyhn, M. B. W. et al. Linking paleogene rifting and inversion in the Northern Song Hong and Beibuwan Basins, Vietnam, with left-lateral motion on the ailao shan-red river shear zone. Tectonics 37, 2559–2585 (2018).
Google Scholar 

19.
Tuyen, N. T. Biostratigraphic report of the ENRECA‐3 well, Bach Long Vi Island prepared for the VPI‐ODA Project, Ho Chi Minh City, Vietnam Oil and Gas Group, Vietnam Petroleum Institute (VPI), Analysis Laboratory Center. (2013).

20.
Rizzi, M. et al. Hinterland setting and composition of an Oligocene deep rift-lake sequence, Gulf of Tonkin, Vietnam: Implications for petroleum source rock deposition. Mar. Pet. Geol. 111, 496–509 (2020).
CAS  Google Scholar 

21.
Dymond, J., Suess, E. & Lyle, M. Barium in deep-sea sediment: A geochemical proxy for paleoproductivity. Paleoceanography 7, 163–181 (1992).
ADS  Google Scholar 

22.
Föllmi, K. B. The phosphorus cycle, phosphogenesis and marine phosphate-rich deposits. Earth-Sci. Rev. 40, 55–124 (1996).
ADS  Google Scholar 

23.
Kump, L. R. & Arthur, M. A. Interpreting carbon-isotope excursions: Carbonates and organic matter. Chem. Geol. 161, 181–198 (1999).
ADS  CAS  Google Scholar 

24.
Tribovillard, N., Algeo, T. J., Lyons, T. & Riboulleau, A. Trace metals as paleoredox and paleoproductivity proxies: An update. Chem. Geol. 232, 12–32 (2006).
ADS  CAS  Google Scholar 

25.
Boyle, J. F. Inorganic geochemical methods in palaeolimnology. In Tracking Environmental Change Using Lake Sediments (eds Last, W. M. & Smol, J. P.) 83–141 (Springer, Berlin, 2005).
Google Scholar 

26.
Calvert, S. & Pedersen, T. F. Elemental proxies for palaeoclimatic and palaeoceanographic variability in marine sediments: Interpretation and application. Dev. Mar. Geol. 1, 567–644 (2007).
Google Scholar 

27.
Wersin, P., Höhener, P., Giovanoli, R. & Stumm, W. Early diagenetic influences on iron transformations in a freshwater lake sediment. Chem. Geol. 90, 233–252 (1991).
ADS  CAS  Google Scholar 

28.
Melles, M. et al. 2.8 million years of arctic climate change from lake El’gygytgyn, NE Russia. Science 337, 315 (2012).
ADS  CAS  PubMed  Google Scholar 

29.
Naeher, S., Gilli, A., North, R. P., Hamann, Y. & Schubert, C. J. Tracing bottom water oxygenation with sedimentary Mn/Fe ratios in Lake Zurich, Switzerland. Chem. Geol. 352, 125–133 (2013).
ADS  CAS  Google Scholar 

30.
Hatch, J. R. & Leventhal, J. S. Relationship between inferred redox potential of the depositional environment and geochemistry of the Upper Pennsylvanian (Missourian) Stark Shale Member of the Dennis Limestone, Wabaunsee County, Kansas, U.S.A.. Chem. Geol. 99, 65–82 (1992).
ADS  CAS  Google Scholar 

31.
Rimmer, S. M. Geochemical paleoredox indicators in Devonian-Mississippian black shales, Central Appalachian Basin (USA). Chem. Geol. 206, 373–391 (2004).
ADS  CAS  Google Scholar 

32.
Cuven, S., Francus, P. & Lamoureux, S. Mid to Late Holocene hydroclimatic and geochemical records from the varved sediments of East Lake, Cape Bounty, Canadian High Arctic. Quat. Sci. Rev. 30, 2651–2665 (2011).
ADS  Google Scholar 

33.
Johnson, T. C., Brown, E. T. & Shi, J. Biogenic silica deposition in Lake Malawi, East Africa over the past 150,000 years. Palaeogeogr. Palaeoclimatol. Palaeoecol. 303, 103–109 (2011).
Google Scholar 

34.
Lerman, A. Lakes: Chemistry, Geology, Physics (Springer, Berlin , 1978). https://doi.org/10.1007/978-1-4757-1152-3.
Google Scholar 

35.
Jia, J., Liu, Z., Bechtel, A., Strobl, S. A. I. & Sun, P. Tectonic and climate control of oil shale deposition in the Upper Cretaceous Qingshankou Formation (Songliao Basin, NE China). Int. J. Earth Sci. 102, 1717–1734 (2013).
CAS  Google Scholar 

36.
Meng, Q., Liu, Z., Bruch, A. A., Liu, R. & Hu, F. Palaeoclimatic evolution during Eocene and its influence on oil shale mineralisation, Fushun basin, China. J. Asian Earth Sci. 45, 95–105 (2012).
ADS  Google Scholar 

37.
Didyk, B. M., Simoneit, B. R. T., Brassell, S. C. & Eglinton, G. Organic geochemical indicators of palaeoenvironmental conditions of sedimentation. Nature 272, 216–222 (1978).
ADS  CAS  Google Scholar 

38.
Sinninghe Damsté, J. S., Ten Haven, H. L., de Leeuw, J. W. & Rullkotter, J. Pristane/phytane ratio as environmental indicator-Reply. Nature 333, 604 (1988).
Google Scholar 

39.
Sinninghe Damsté, J. S. et al. Evidence for gammacerane as an indicator of water column stratification. Geochim. Cosmochim. Acta 59, 1895–1900 (1995).
ADS  PubMed  Google Scholar 

40.
Irwin, H., Curtis, C. & Coleman, M. Isotopic evidence for source of diagenetic carbonates formed during burial of organic-rich sediments. Nature 269, 209–213 (1977).
ADS  CAS  Google Scholar 

41.
Curtis, C. Mineralogical consequences of organic matter degradation in sediments: Inorganic/organic diagenesis. In Marine Clastic Sedimentology (eds Leggett, J. K. & Zuffa, G. G.) 108–123 (Springer, Berlin, 1987). https://doi.org/10.1007/978-94-009-3241-8_6.
Google Scholar 

42.
Mozley, P. S. & Wersin, P. Isotopic composition of siderite as an indicator of depositional environment. Geology 20, 817–820 (1992).
ADS  CAS  Google Scholar 

43.
Schieber, J., Southard, J. & Thaisen, K. Accretion of mudstone beds from migrating floccule ripples. Science 318, 1760 (2007).
ADS  CAS  PubMed  Google Scholar 

44.
Macquaker, J. H. S. & Bohacs, K. M. On the accumulation of mud. Science 318, 1734 (2007).
CAS  PubMed  Google Scholar 

45.
Föllmi, K. B. & Grimm, K. A. Doomed pioneers: Gravity-flow deposition and bioturbation in marine oxygen-deficient environments. Geology 18, 1069–1072 (1990).
ADS  Google Scholar 

46.
Ingall, E. & Jahnke, R. Evidence for enhanced phosphorus regeneration from marine sediments overlain by oxygen depleted waters. Geochim. Cosmochim. Acta 58, 2571–2575 (1994).
ADS  CAS  Google Scholar 

47.
Fyhn, M. B. W., Hoang, B. H., Anh, N. T., Hovikoski, J. & Cuong, T. D. Eocene-Oligocene syn-rift deposition in the northern Gulf of Tonkin, Vietnam. Mar. Pet. Geol. 111, 390–413 (2020).
Google Scholar 

48.
Talling, P. J. On the triggers, resulting flow types and frequencies of subaqueous sediment density flows in different settings. Mar. Geol. 352, 155–182 (2014).
ADS  Google Scholar 

49.
Morley, R. A review of the Cenozoic palaeoclimate history of Southeast Asia. In Biotic Evolution and Environmental Change in Southeast Asia (ed. Rosen, B.) 79–114 (Cambridge University Press, Cambridge, 2012).
Google Scholar 

50.
Bohacs, K. M., Carroll, A. R., Neal, J. E. & Mankiewicz, P. J. Lake-Basin Type, Source Potentional, and Hydrocarbon Character:An Integrated Sequence-Stratigraphic-Geochemichal Framework. Lake basins through space and time: AAPG Studies in Geology 46, (American Association of Petroleum Geologists, 2000).

51.
Demaison, G. J. & Moore, G. T. Anoxic environments and oil source bed genesis. Org. Geochem. 2, 9–31 (1980).
CAS  Google Scholar 

52.
Cohen, A. S. Facies relationships and sedimentation in large rift lakes and implications for hydrocarbon exploration: Examples from lakes Turkana and Tanganyika. Palaeogeogr. Palaeoclimatol. Palaeoecol. 70, 65–80 (1989).
Google Scholar 

53.
Calvert, S. E., Bustin, R. M. & Pedersen, T. F. Lack of evidence for enhanced preservation of sedimentary organic matter in the oxygen minimum of the Gulf of California. Geology 20, 757–760 (1992).
ADS  CAS  Google Scholar 

54.
Katz, B. J. Lacustrine basin hydrocarbon exploration -current thoughts. J. Palaeolimnol. 26, 161–179 (2001).
ADS  Google Scholar 

55.
Bohacs, K. M. et al. Production, destruction, and dilution—The many paths to source-rock development. Spec. Publ. SEPM 82, 61–101 (2005).
CAS  Google Scholar 

56.
Lehmann, M. F., Bernasconi, S. M., Barbieri, A. & McKenzie, J. A. Preservation of organic matter and alteration of its carbon and nitrogen isotope composition during simulated and in situ early sedimentary diagenesis. Geochim. Cosmochim. Acta 66, 3573–3584 (2002).
ADS  CAS  Google Scholar 

57.
Harris, N. B. et al. Patterns of organic-carbon enrichment in a lacustrine source rock in relation to paleo-lake level, Congo Basin, West Africa. Depos. Org. Sediments Model. Mech. Consequences SEPM Spec. Publ. 82, 103–123 (2005).
CAS  Google Scholar 

58.
Bojanowski, M. J. & Clarkson, E. N. K. Origin of siderite concretions in microenvironments of methanogenesis developed in a sulfate reduction zone: An exception or a rule?. J. Sediment. Res. 82, 585–598 (2012).
ADS  CAS  Google Scholar 

59.
Berner, R. A. Early Diagenesis—A Theoretical Approach. Princeton Series in Geochemistry (Princeton University Press, Princeton, 1980). https://doi.org/10.1016/0037-0738(81)90046-4.
Google Scholar 

60.
Raiswell, R. Chemical model for the origin of minor limestone-shale cycles by anaerobic methane oxidation. Geology 16, 641–644 (1988).
ADS  CAS  Google Scholar 

61.
Canfield, D. E. Reactive iron in marine sediments. Geochim. Cosmochim. Acta 53, 619–632 (1989).
ADS  CAS  PubMed  Google Scholar 

62.
Katsev, S., Sundby, B. & Mucci, A. Modeling vertical excursions of the redox boundary in sediments: Application to deep basins of the Arctic Ocean. Limnol. Oceanogr. 51, 1581–1593 (2006).
ADS  CAS  Google Scholar 

63.
Esbensen, K. H., Guyot, D., Westad, F. & Houmoller, L. P. Multivariate Data Analysis—In Practice: An Introduction to Multivariate Data Analysis and Experimental Design. (Camo, 2002).

64.
Bish, D. L. & Post, J. E. Quantitative mineralogical analysis using the Rietveld full-pattern fitting method. Am. Mineral. 78, 932–940 (1993).
CAS  Google Scholar 

65.
Hutton, A. C. Petrographic classification of oil shales. Int. J. Coal Geol. 8, 203–231 (1987).
Google Scholar 

66.
Taylor, G. H. et al. Organic Petrology (Gebrüder Borntraeger, Berlin, 1998).
Google Scholar 

67.
International, Committee for Coal and Organic Petrology. The new inertinite classification (ICCP System 1994). Fuel 80, 459–471 (2001).
Google Scholar 

68.
Sýkorová, I. et al. Classification of huminite—ICCP System 1994. Int. J. Coal Geol. 62, 85–106 (2005).
Google Scholar 

69.
Nytoft, H. P. et al. Biomarkers of Oligocene lacustrine source rocks, Beibuwan-Song Hong basin junction, offshore northern Vietnam. Mar. Pet. Geol. 114, 104196 (2020).
CAS  Google Scholar 

70.
Radke, M., Willsch, H. & Welte, D. H. Preparative hydrocarbon group type determination by automated medium pressure liquid chromatography. Anal. Chem. 52, 406–411 (1980).
CAS  Google Scholar 

71.
Craig, H. Isotopic standards for carbon and oxygen and correction factors for mass-spectrometric analysis of carbon dioxide. Geochim. Cosmochim. Acta 12, 133–149 (1957).
ADS  CAS  Google Scholar 

72.
Friedman, I. & O’Neil, J. R. Compilation of stable isotope fractionation factors of geochemical interest. US Gov. Print. Off. 440, 1–11 (1977).
Google Scholar  More

1.
Karl D, Michaels A, Bergman B, Capone D, Carpenter E, Letelier R, et al. Dinitrogen fixation in the world’s oceans. Biogeochemistry. 2002;57/58:47–98.
CAS  Article  Google Scholar 
2.
Bombar D, Paerl RW, Riemann L. Marine non-cyanobacterial diazotrophs: moving beyond molecular detection. Trends Microbiol. 2016;24:916–27.
CAS  Article  Google Scholar 

3.
Bird C, Martinez JM, O’Donnell AG, Wyman M. Spatial distribution and transcriptional activity of an uncultured clade of planktonic diazotrophic γ-proteobacteria in the Arabian Sea. Appl Environ Microbiol. 2005;71:2079–85.
CAS  Article  Google Scholar 

4.
Moisander PH, Beinart RA, Voss M, Zehr JP. Diversity and abundance of diazotrophic microorganisms in the South China Sea during intermonsoon. ISME J. 2008;2:954–67.
CAS  Article  Google Scholar 

5.
Moisander PH, Serros T, Paerl RW, Beinart RA, Zehr JP. Gammaproteobacterial diazotrophs and nifH gene expression in surface waters of the South Pacific Ocean. The. ISME J. 2014;8:1962–73.
CAS  Article  Google Scholar 

6.
Langlois R, Großkopf T, Mills M, Takeda S, LaRoche J. Widespread distribution and expression of Gamma A (UMB), an uncultured, diazotrophic, γ-proteobacterial nifH phylotype. PLoS ONE. 2015;10:e0128912.
Article  Google Scholar 

7.
Zehr JP, Turner PJ. Nitrogen fixation: nitrogenase genes and gene expression. Methods Microbiol. 2001;30:271–86.
CAS  Article  Google Scholar 

8.
Church MJ, Short CM, Jenkins BD, Karl DM, Zehr JP. Temporal patterns of nitrogenase gene (nifH) expression in the oligotrophic North Pacific Ocean. Appl Environ Microbiol. 2005;71:5362–70.
CAS  Article  Google Scholar 

9.
Langlois RJ, Hümmer D, LaRoche J. Abundances and distributions of the dominant nifH phylotypes in the Northern Atlantic Ocean. Appl Environ Microbiol. 2008;74:1922–31.
CAS  Article  Google Scholar 

10.
Bonnet S, Rodier M, Turk-Kubo KA, Germineaud C, Menkes C, Ganachaud A, et al. Contrasted geographical distribution of N2 fixation rates and nifH phylotypes in the Coral and Solomon Seas (southwestern Pacific) during austral winter conditions. Glob Biogeochem Cycles. 2015;29:1874–92.
CAS  Article  Google Scholar 

11.
Gradoville MR, Bombar D, Crump BC, Letelier RM, Zehr JP, White AE. Diversity and activity of nitrogen-fixing communities across ocean basins. Limnol Oceanogr. 2017;62:1895–909.
Article  Google Scholar 

12.
Chen TY, Chen YL, Sheu DS, Chen HY, Lin YH, Shiozaki T. Community and abundance of heterotrophic diazotrophs in the northern South China Sea: revealing the potential importance of a new alphaproteobacterium in N2 fixation. Deep Sea Res Part I. 2019;143:104–14.
CAS  Article  Google Scholar 

13.
Moisander PH, Benavides M, Bonnet S, Berman-Frank I, White AE, Riemann L. Chasing after non-cyanobacterial nitrogen fixation in marine pelagic environments. Front Microbiol. 2017;8:1736.
Article  Google Scholar 

14.
Farnelid H, Turk-Kubo K, Ploug H, Ossolinski JE, Collins JR, Van Mooy BA, et al. Diverse diazotrophs are present on sinking particles in the North Pacific Subtropical Gyre. ISME J. 2019;13:170–82.
Article  Google Scholar 

15.
Cornejo-Castillo FM. Diversity, ecology and evolution of marine diazotrophic microorganisms. Barcelona: Universitat Politècnica de Catalunya; 2017.

16.
Delmont TO, Quince C, Shaiber A, Esen ÖC, Lee ST, Rappé MS, et al. Nitrogen-fixing populations of planctomycetes and proteobacteria are abundant in surface ocean metagenomes. Nat Microbiol. 2018;3:804–13.
CAS  Article  Google Scholar 

17.
Benavides M, Moisander PH, Daley MC, Bode A, Aristegui J. Longitudinal variability of diazotroph abundances in the subtropical North Atlantic Ocean. J Plankton Res. 2016;38:662–72.
CAS  Article  Google Scholar 

18.
Gradoville MR, Farnelid H, White AE, Turk-Kubo KA, Stewart B, Ribalet F, et al. Latitudinal constraints on the abundance and activity of the cyanobacterium UCYN-A and other marine diazotrophs in the North Pacific. Limnol Oceanogr. 2020. https://doi.org/10.1002/lno.11423.

19.
Scavotto RE, Dziallas C, Bentzon-Tilia M, Riemann L, Moisander PH. Nitrogen-fixing bacteria associated with copepods in coastal waters of the North Atlantic Ocean. Environ Microbiol. 2015;17:3754–65.
CAS  Article  Google Scholar 

20.
Carradec Q, Pelletier E, Da Silva C, Alberti A, Seeleuthner Y, Blanc-Mathieu R, et al. A global ocean atlas of eukaryotic genes. Nat Commun. 2018;9:1–3.
CAS  Article  Google Scholar 

21.
Karl DM, Church MJ, Dore JE, Letelier RM, Mahaffey C. Predictable and efficient carbon sequestration in the North Pacific Ocean supported by symbiotic nitrogen fixation. Proc Natl Acad Sci. 2012;109:1842–9.
CAS  Article  Google Scholar 

22.
Pedersen JN, Bombar D, Paerl RW, Riemann L. Diazotrophs and N2-fixation associated with particles in coastal estuarine waters. Front Microbiol. 2018;9:2759.
Article  Google Scholar 

23.
Ploug H, Kühl M, Buchholz-Cleven B, Jørgensen BB. Anoxic aggregates-an ephemeral phenomenon in the pelagic environment? Aquat Microb Ecol. 1997;13:285–94.
Article  Google Scholar 

24.
Ploug H. Small-scale oxygen fluxes and remineralization in sinking aggregates. Limnol Oceanogr. 2001;46:1624–31.
CAS  Article  Google Scholar 

25.
Bentzon-Tilia M, Severin I, Hansen LH, Riemann L. Genomics and ecophysiology of heterotrophic nitrogen-fixing bacteria isolated from estuarine surface water. MBio. 2015;6:e00929–15.
CAS  Article  Google Scholar 

26.
Bergman B, Sandh G, Lin S, Larsson J, Carpenter EJ. Trichodesmium—a widespread marine cyanobacterium with unusual nitrogen fixation properties. FEMS Microbiol Rev. 2013;37:286–302.
CAS  Article  Google Scholar 

27.
Moraru C, Lam P, Fuchs BM, Kuypers MM, Amann R. GeneFISH—an in situ technique for linking gene presence and cell identity in environmental microorganisms. Environ Microbiol. 2010;12:3057–73.
CAS  Article  Google Scholar 

28.
Batani G, Bayer K, Böge J, Hentschel U, Thomas T. Fluorescence in situ hybridization (FISH) and cell sorting of living bacteria. Sci Rep. 2019;9:1–3.
Article  Google Scholar  More

Demographic dynamics of urban systems
We now derive the population size distribution of cities in the most general way, set by their demographic dynamics. The change of the population Ni in city i during unit time Δt necessarily results from the balance of births, deaths, and migration as

$${N}_{i}(t+Delta t)={N}_{i}(t)+Delta tleft[{upsilon }_{i}{N}_{i}(t)+sum_{j = 1,jne i}^{{N}_{{rm{c}}}}({J}_{ji}-{J}_{ij})right],$$
(2)

where Nc is the total number of cities, υi = bi − di is the city’s vital rate, the difference between its per capita average birth rate, bi, and the death rate, di. The current Jij is the number of people moving from city i to city j over the time interval Δt, and vice-versa for Jji. For simplicity of notation we will work in units where Δt = 1. Eq. (2) only accounts explicitly for migration between different specific cities, which can cross political borders. Migration to city i from non-specific places outside the system, for example from non-urban regions or other nations, can be incorporated into the vital rates υi. Note that this model is very general, and naturally accommodates the inclusion of new cities over time when their sizes become non-zero, and indeed the disappearance of others, if their size vanishes.
To proceed we need to explore the dependence of the currents Jij on population size. We start by writing these currents as population rates, Jij = mijNi(t), where mij is the probability that someone in city i moves to city j over the time period. This allows us to write Eq. (2) in matrix form

$${bf{N}}(t)={bf{A}}(t){bf{N}}(t-1),$$
(3)

where ({bf{N}}(t)={[{N}_{1}(t),ldots ,{N}_{i}(t),ldots ,{N}_{{N}_{{rm{c}}}}(t)]}^{T}) is the vector of populations in each city at time t and the matrix A projects the population at time t − 1 to time t. This matrix is known in population dynamics as the environment39, with elements

$${A}_{ij}=left{begin{array}{ll}1+{upsilon }_{i}-{m}_{i}^{{rm{out}}},&{rm{if}},i=j\ {m}_{ji},&{rm{if}},ine jend{array}right.,$$
(4)

where ({m}_{i}^{{rm{out}}}=mathop{sum }nolimits_{k = 1,kne i}^{{N}_{{rm{c}}}}{m}_{ik}) is the total probability that a person migrates out of city i per unit time. We also define the population structure vector x = [xi], with xi = Ni/NT, which is the probability of finding a person in city i out of the total population NT. Note that (mathop{sum }nolimits_{i = 1}^{{N}_{{rm{c}}}}{x}_{i}=1), so that the structure vector has Nc − 1 degrees of freedom (the magnitude of x is fixed).
Equation (3) can now be solved by repeated iteration

$${bf{N}}(t)={bf{A}}(t){bf{A}}(t-1)ldots {bf{A}}(1){bf{N}}(0).$$
(5)

City size distributions in different environments
A number of important ergodic theorems40 in population dynamics tell us about the properties of the solutions under different conditions on the environments A(t). These results rely on some constraints on the properties of A; specifically, that the product of matrices in Eq. (5) is positive for sufficiently long times39,41.
First, the weak ergodic theorem guarantees that for a dynamical sequence of environments, the difference between two different initial population structure vectors decays to zero over time. This means that there is typically an asymptotic city size “distribution”, which is a function of environmental dynamics only, independent of initial conditions. When the environment is stochastic but otherwise time independent, the strong stochastic ergodic theorem states that the structure vector becomes a random variable whose probability distribution converges to a fixed stationary distribution, regardless of initial conditions. This is the sense in which most derivations of Zipf’s law apply17,33. For stochastic environments, only probability distributions of structure vectors, not vectors themselves, can be predicted. Finally, in situations where the environment is time dependent and stochastic, the weak stochastic ergodic theorem states that the difference between the probability distributions for the structure vector resulting from any two initial populations, exposed to independent sample paths, decays to zero. Again, in cases where the environment is explicitly dynamic, besides being stochastic in a stationary sense, we cannot say much about the actual probability distribution of city sizes, only that the importance of initial conditions vanishes for sufficiently long times.
To get more intuition on these results we will now show some explicit solutions. The city size distribution can be calculated exactly when the environments A are arbitrarily complicated, but static. This is the result of the strong ergodic theorem, which says that in a constant environment A, any initial population vector converges to a fixed stable structure given by the leading eigenvector of the environment39. This is guaranteed to exist if A is a strongly connected aperiodic graph, meaning that any city can be reached from any other city in a finite and diverse number of intermediate steps along the non-zero migration flows between them. This is a reasonable assumption to make about an urban system, which may even serve as a definition. The leading eigenvector of A is the eigenvector centrality of the urban system. It describes the probabilistic location of a random walker over the graph of migration flows42. This is equal to the stationary solution of numerically integrating a network described by environment A. Figure 1 shows two numerical solutions with different initial conditions, but the same environment A. The cities in the two runs converge to similar sizes, as they experience a common environmental matrix with all entries set to be the same plus a little noise. The insets show how the population structure converges from both initial conditions to the one defined by the leading eigenvector e0. Because the structure vector is the probability of finding a person out of the total population in the urban system in a specific city, we use the Kullback-Leibler (KL) divergence43({D}_{{rm{KL}}}(P| {P}_{{{bf{e}}}_{0}})={sum }_{x}P[x]{mathrm{log}},P[x]/{P}_{{{bf{e}}}_{0}}[x]) as a measure of the ‘distance’ between the initial distribution, P, and the final, ({P}_{{{bf{e}}}_{0}}). The KL divergence is a foundational quantity in information theory. It measures the information gain in describing a probabilistic state using one distribution, ({P}_{{{bf{e}}}_{0}}), versus another, P, in units of information (bits). Note that (Pto {P}_{{{bf{e}}}_{0}}) and thus ({D}_{{rm{KL}}}(P| {P}_{{{bf{e}}}_{0}})to 0) over time, see insets in Fig. 1.
Fig. 1: Temporal evolution of the relative population structure of cities in the same environment for different initial conditions.

Lines shows the trajectory of each city in terms of the fraction its population, Ni to the total NT. a shows an initial situation where all cities start out with similar sizes, whereas in b they are initiated following Zipf’s law (shown in log-scale). In both cases, the relative city size distribution eventually becomes stationary at long times (vertical red line) with the same population structure. This is given by the eigenvector e0, corresponding to the leading eigenvalue of the environment. The insets show how the system converges in both cases to this common structure set by the environment, because the Kullback–Leibler divergence ({D}_{{rm{KL}}}(P| {P}_{{{bf{e}}}_{0}})) approaches zero at late times.

Full size image

The convergence rate of the city size distribution to the leading eigenvector is given by the ratio of secondary eigenvalues to the dominant one. Explicitly, we can now write this solution as

$${bf{N}}(t)={{bf{A}}}^{t}{bf{N}}(0)to {bf{N}}(t)=sum_{i = 1}^{{N}_{{rm{c}}}}{({lambda }_{i})}^{t}{c}_{i}{{bf{e}}}_{i}.$$
(6)

where Aei = λiei, so that ei are the eigenvectors of the matrix A and λi the corresponding eigenvalues, so that λ0  > Re(λ1)  > Re(λ2)  > . . . , where Re(λi) is the real part of the eigenvalue. The projection coefficients ci are such that ({bf{N}}(0)=mathop{sum }nolimits_{i = 1}^{{N}_{{rm{c}}}}{c}_{i}{{bf{e}}}_{i}), which can be written as c = E−1N(0), where E is a matrix whose jth column vector is ej. The strong ergodic theorem states that the largest eigenvalue λ0 is positive and real and that the associated eigenvector e0 is also positive in terms of all its entries. This means that, over time, the dynamics of the urban system approaches that of the growth of its dominant eigenvector, because the projections on all other eigenvectors decay exponentially in relative terms, as

$${left(frac{{lambda }_{i}}{{lambda }_{0}}right)}^{t}={e}^{-mathrm{ln},frac{{lambda }_{0}}{{lambda }_{i}}t},{mathop{-hskip -4.3pt-hskip -4.3ptlongrightarrow}limits_{tto infty}},0$$
(7)

and exemplified in the two insets in Fig. 1. The longest characteristic time is associated with the difference in magnitude between the two largest eigenvalues ({t}_{*}=1/{mathrm{ln}},frac{{lambda }_{0}}{| {lambda }_{1}| }). Consequently, for t ≫ ({t}_{*}) we observe an extreme dimensional reduction from an initial condition characterized in general by Nc degrees of freedom to a final one with just one, set by the environment’s leading eigenvector. Using the values of migration flows in the US urban system over the last decade, we can compute the value of ({t}_{*}). It is rather long—of the order of several centuries—when measured using census data44 and a little shorter using data on tax returns45.
These results allow us to consider very general classes of dynamics and initial conditions. We therefore conclude that in general, given a specific set of demographic conditions, the total population can grow exponentially with a common growth rate across cities. In addition, the relative population distribution at late times does not resemble anything like Zipf’s law. As expected from the statements of the ergodic theorems, the population structure is set by the nature of the intercity flows (and vital rates), or equivalently by the environment A.
Stochastic demographic dynamics and symmetry breaking
The step towards a statistical solution for the city size distribution—and obtaining Zipf’s law in particular—requires the additional consideration of statistical fluctuations, which must arise in the vital and migration rates. Just as in models of statistical physics, the introduction of stochasticity implies not just fluctuating quantities but potentially also the restoration of broken symmetries46,47. In the context of urban systems, restoring broken symmetries means removing preferences for population growth in specific cities over others, associated with imbalances in vital rates and migration flows. To make these expectations explicit we write the migration rates mij as

$${J}_{ij}={m}_{ij}{N}_{i}(t)=left(frac{{s}_{ij}+{delta }_{ij}}{2}right)frac{{N}_{j}(t)}{{N}_{{rm{T}}}(t)}{N}_{i}(t).$$
(8)

Eq. (8) splits the migration flows into two parts: sij describes the symmetric flow between two cities (sij = sji), and δij the anti-symmetric part (δij = −δji). Note that any correlations between Ni and Nj are preserved in these two quantities. Making the population size of the destination city explicit will allow to diagonalize Eq. (2) and therefore to find self-consistent solutions for the structure vector. This form of the migration currents also makes contact with the gravity law of geography, which is typically written as

$${J}_{ij}={G}_{g}frac{{N}_{i}{N}_{j}}{{d}_{ij}^{gamma }}={J}_{ji}.$$
(9)

In this form, the gravity law is symmetric, corresponding to Eq. (8) when δij = 0 and sij = GgNT(t)f(dij), where Gg is the ‘gravitational’ constant, (f({d}_{ij})=1/{d}_{ij}^{gamma }) is a decaying function of distance dij and γ is the distance exponent. Only the anti-symmetric part contributes to the dynamics of population size change, since we can write Eq. (2) as

$${N}_{i}(t+1)=left[1+{upsilon }_{i}-{overline{delta }}_{i}right]{N}_{i}(t).$$
(10)

with ({overline{delta }}_{i}=mathop{sum }nolimits_{j = 1}^{{N}_{{rm{c}}}}{delta }_{ij}{x}_{j}). We can now appreciate the importance of the bi-linearity of the migration flows on both the origin and destination population as a means to diagonalize the demographic dynamics. We have achieved this however at the cost of introducing a slow non-linearity in each city’s growth rate, via their dependence on the population structure vector x. We can establish the existence of a self-consistent solution x* such that for long times

$${upsilon }_{i}-{overline{upsilon }}^{* }=sum_{j = 1}^{{N}_{{rm{c}}}}{delta }_{ij}{x}_{j}^{* },$$
(11)

where ({overline{upsilon }}^{* }=mathop{sum }nolimits_{i}^{{N}_{{rm{c}}}}{upsilon }_{i}{x}_{i}^{* }), since (mathop{sum }nolimits_{i,j = 1}^{{N}_{{rm{c}}}}{delta }_{ij}{x}_{i}{x}_{j}=0), by the conservation of the migration currents. Note that we did not assume any explicit form for the migration flows as a function of distance between places, so that our discussion includes many specific models, including different version of the gravity law.
This equation has a clear geometric meaning: The matrix δ = [δij] is anti-symmetric and real, so that it behaves as an infinitesimal generator of rotations: R(x) = x + δx. We can write the self-consistent solution in terms of these rotation as

$$R({{bf{x}}}^{* })={{bf{x}}}^{* }+({boldsymbol{upsilon }}-overline{upsilon }),$$
(12)

where (overline{upsilon }) is still a function of x and υ = [υi]. Note that these rotations are associated with very small angles, and that the structure vector is subject to constraints, such as staying positive. The first two terms of Eq. (12) ask for a general solution for x that is invariant under rotations. However, the last term, which introduces differences in relative growth rates for different cities, breaks this symmetry explicitly. This specificity of the relative growth rates leads to particular solutions that do not coincide with Zipf’s law. Fig. 2a shows the example of a numerical solution in a non-linear, non-stochastic environment, defined by Eq. (8). The final population structure in this environment is still well defined, as it reaches a steady state. However, compared to the time independent case (see Fig. 1), the final structure takes much longer to unfold and has a stronger urban hierarchy.
Fig. 2: City size distribution emerging from stochastic demographic dynamics differ significantly from those in time independent environments.

a depicts the temporal trajectory for a set of 100 cities in a non-stochastic, non-linear environment, described by Eq. (8). Despite generalizing the demographic dynamics in Fig. 1, this results in a fixed urban hierarchy at late times (vertical red line), regardless of initial conditions. However, when migration and vital rates become stochastic, the demographic dynamics no longer has a static solution at long times, b. Instead, cities constantly change their relative sizes (rank) over time. After an initial transient, the population structure fluctuates close to Zipf’s law. The insets show the population structure (red) at the point in time when the non-stochastic evolution (a) becomes approximately stationary, with Zipf’s law (dashed blue line) shown for reference.

Full size image

Symmetry restoration and the emergence of Zipf’s law
The stage is now set to consider what it takes to counter the symmetry breaking resulting from the selective structure of the environment. To do that, we must consider the statistics of the demographic rates. We start by rewriting Eq. (10) in terms of the dynamical equations for the structure vectors x as

$${x}_{i}(t+1)=left[1+{upsilon }_{i}-overline{upsilon }-{overline{delta }}_{i}right]{x}_{i}(t)=(1+{epsilon }_{i}){x}_{i}(t).$$
(13)

We now specify the relative growth rate ({epsilon }_{i}={upsilon }_{i}-overline{upsilon }-{overline{delta }}_{i}) for each xi as a statistical quantity. It is clear by construction that its average over cities vanishes, ({overline{epsilon }}_{i}=0). The only remaining question is how ϵi varies over time. There is strong empirical evidence that on the time scale of years to decades some cities can maintain larger growth rates than others. For example, presently in the US, most cities of the Southwest and South are growing fast, while other cities show slower growth or relative decay, such as many urban areas in the Midwest and Appalachia48. Moreover, it has been true in recent decades that mid-sized cities in the US have been growing faster than either the largest cities or small micropolitan areas, so that even population aggregated growth rates over certain size ranges differ. Going further back in time one could argue that the same cities of the Midwest were once booming as they became the manufacturing centers of the US, especially between the Civil War and the decades post WWII49. Similar arguments may be made for cities in other urban systems, such as in Europe50,51. Thus, it is possible that on very long time scales, T—on the order of a century or longer—the growth rates of cities equalize to very similar numbers and therefore the temporal average (langle {epsilon }_{i}rangle =frac{1}{T}mathop{sum }nolimits_{t}^{T}{epsilon }_{i}(t)to 0). This requires that the temporal average of Eq. (11) remains zero, which means that the structure vector is both rotated and dilated ’randomly’ in ways that over sufficiently long times lead to 〈ϵi〉 = 0, for each city. Numerical solutions of the demographic equations under these conditions show a population structure that closely resembles Zipf’s law for long times, see Fig. 2b.
In the following, we show that Zipf’s law is a probability density for this derived dynamics in the long run. To do this, let us characterize the variance in the temporal growth rate fluctuations as ({sigma }_{i}^{2}=langle {epsilon }_{i}^{2}rangle ={sigma }_{{upsilon }_{i}-upsilon }^{2}+{sigma }_{{overline{delta }}_{i}}^{2}+2{rm{COV}}({upsilon }_{i}-upsilon ,{overline{delta }}_{i}) , > , 0). Over the same long time scales, statistical fluctuations in the growth rate may obey a limit theorem, so that they become approximately Gaussian, with variance, ({sigma }_{i}^{2}). Then, Eq. (13) becomes simpler, and acquires a direct correspondence to familiar stochastic growth processes. It can be written in analogy to geometric Brownian motion without drift as

$$d{x}_{i}={epsilon }_{i}{x}_{i}={x}_{i}{sigma }_{i}dW$$
(14)

where W is a standard Brownian motion. We now see how demographic dynamics can be simplified through a chain of assumptions into a form consistent with typical stochastic processes leading to Zipf’s law as proposed by Gabaix20,52.
Equation (14) can be expressed as an equation for the probability density P ≡ P[xi, t∣xi(0), 0] of observing xi at time t, having started with an initial condition where xi(0) was observed at time zero,

$$frac{d}{dt}P=frac{{d}^{2}}{d{x}_{i}^{2}}{sigma }_{i}^{2}{x}_{i}^{2}P=frac{d{J}_{{x}_{i}}}{d{x}_{i}}.$$
(15)

The last term describes a probability current, ({J}_{{x}_{i}}=frac{d}{d{x}_{i}}{sigma }_{i}^{2}{x}_{i}^{2}P({x}_{i})). This equation is exactly solvable, see Methods for technical details. To find the stationary solutions, we ask that the right hand side of Eq. (15) vanishes. The first of the two solutions corresponds to a vanishing current, ({J}_{{x}_{i}}=frac{d}{d{x}_{i}}{sigma }_{i}^{2}{x}_{i}^{2}P({x}_{i})=0) and is (P=frac{c^{prime} }{{sigma }_{i}^{2}{x}_{i}^{2}}), where (c^{prime} ) is a normalization constant. If ({sigma }_{i}^{2}={sigma }^{2}) independent of x, as proposed by Gibrat’s law, we can drop the index so this becomes Zipf’s law, (P(N)={P}_{{rm{z}}}({N}_{i}) sim frac{c}{{N}^{2}}). Note for example that if ({sigma }_{i}^{2} sim {N}^{-alpha }) (α may be negative), which violates Gibrat’s law in terms of fluctuations of the growth rate, then P ~ 1/N2−α. Thus, the city size dependence of growth rate fluctuations will change the exponent in Zipf’s law and may destroy scale-invariance altogether if it is not a power law. The second solution applies for constant current, i. e. Jx = J, which leads to (P=frac{c^{primeprime} }{{sigma }^{2}x}), where c″ is a normalization constant. This represents a flow of probability across the urban hierarchy. Both solutions are attractors of the stochastic dynamics that emerge for long times, t ≫ 1/σ2.
We find that in the limit where the relative difference in vital rates vanishes due to fluctuations, the structure vector becomes rotationally invariant on average and the angular symmetry of x is effectively restored. This leads to an effective symmetry of the relative city size distribution on the time scale tr = 1/2σ2, when growth rate fluctuations vanish. The time tr can be estimated using US data to be typically very long, on the range of many centuries or even millennia. Under these conditions all cities share statistically identical dynamics and can be interchanged, resulting in Zipf’s law as the time average of population size distributions, see Fig. 3. In the same limit, the anti-symmetric components of intercity flows will average out and the gravity law will emerge in its conventional form. However, the observations of city sizes at each time are samples of this distribution and may reflect decade-long observable positive and negative preferences for certain cities. Figure 4 shows this effect for the US urban system using decennial census data from 1790 to 1990. For time periods of a few decades, the temporal average does not visibly converge to Zipf’s law: The blue line in the inset of Fig. 4 depicts this effect. However, the cumulative temporal average of the size distributions starts to approximate Zipf’s law after about five decades. The red line (inset) shows the KL divergence between the cumulative temporal average and Zipf’s law, starting from 1790 to subsequent census.
Fig. 3: Stochastic migration rates restore the symmetry of city growth rates over time, leading to Zipf’s law.

The figures show the temporal solution of the demographic dynamics (Eq. (13)) when rotational symmetry is restored through stochasticity averaged over long times. a Shows the population structure at different times (dashed lines) and the distribution with smallest divergence from Zipf’s law over a long run (red), measured by the smallest observed DKL(P∣Pz). b shows the trajectory of the Kullback-Leibler divergence away from Zipf’s law over time. After a long time period, the population structure approximates Zipf’s law and fluctuates around it. In this regime, at any single time only samples showing some deviations from Zipf’s distribution are observable: it is only on the average over many structure vectors at long times that Zipf’s law emerges as a good characterization of the city size distribution. Fig. 4 shows that this was indeed the case for the US urban system until recently.

Full size image

Fig. 4: The dynamics of the rank-size distribution for cities in the US between 1790 and 1990.

Population structure distributions for the largest 100 cities in the US53 are shown as dashed lines, while the average over all years is shown in red. The inset depicts the DKL(P∣Pz) to Zipf’s law of the population structure for each available year (blue) and of the cumulative average over all structure vectors up the specific year (red). We see that the temporal average steadily approaches Zipf’s law as more years are added. After an initial phase, the average is always closer to Zipf’s law than the distribution in any single year. Note that in recent decades, the population structure vector began consistently deviating from Zipf’s law, leading also to a small divergence of the cumulative temporal averages. This effect is to a large extent the consequence of mid-sized and large metropolitan areas growing faster over this period than the largest cities in the US urban system.

Full size image More

Protected areas
We determined the structural connectivity of the global network of PAs by quantifying intact continuous pathways (areas largely devoid of high anthropogenic pressures) between PAs. Data on PA location and boundary were obtained from the May 2019 World Database of Protected Areas (WDPA)70. We only considered PAs that had a land area of at least 10 km2. As China removed most of its PAs from the public May 2019 WDPA version, we used the April 2018 WDPA for China only, which contained the full set of Chinese PAs at the time. It is important to note that our statistics may differ from those reported by countries and territories due to methodologies and dataset differences used to measure terrestrial area of a country or territory.
Measure of human pressure
We used the latest global terrestrial human footprint (HFP) maps—a cumulative index of eight variables measuring human pressure on the global environment—to calculate the average human pressure between PAs35. While there are other human pressure maps74,75,76, the HFP is a well-accepted dataset that provided a validation analysis using scored pressures from 3114 × 1 km2 random sample plots. The root mean squared error for the 3114 validation plots was 0.125 on the normalized 0–1 scale, indicating an average error of approximately 13%. The Kappa statistic was 0.737, also indicating high concurrence between the HFP and the validation dataset. The HFP 2013 map uses the following variables: (1) the extent of built human environments, (2) population density, (3) electric infrastructure, (4) crop lands, (5) pasture lands, (6) roads, (7) railways, and (8) navigable waterways.
Navigable waterways such as rivers and lakes are included within HFP as they can act as conduits for people to access nature35. In the latest HFP (2013), rivers and lakes are included based on size and visually identified shipping traffic and shore side settlements. Venter et al. treated the great lakes of North America, Lake Nicaragua, Lake Titicaca, Lake Onega, Lake Peipus, Lake Balkash, Lake Issyk Kul, Lake Victoria, Lake Tanganyika and Lake Malawi, as they did navigable marine coasts (i.e. only considered coasts as navigable for 80 km either direction of signs of a human settlement, which were mapped as a night lights signal with a Digital Number (DN)  > 6 within 4 km of the coast35). Rivers were included if their depth was >2 m and there were night-time lights (DN  > = 6) within 4 km of their banks, or if contiguous with a navigable coast or large inland lake, and then for a distance of 80 km or until stream depth is too shallow for boats. To map rivers and their depth, Venter et al. used the hydrosheds (hydrological data and maps based on shuttle elevation derivatives at multiple scales) dataset on stream discharge, and the following formulae: stream width = 8.1× (discharge[m3/s])0.58; and velocity = 4.0 × (discharge[m3/s])0.6/(width[m]); and cross-sectional area = discharge/velocity; and depth = 1:5× area/width35.
Each human pressure was scaled from 0–10, then weighted within that range according to estimates of their relative levels of human pressure following Sanderson et al.77. The resulting standardized pressures were then summed together to create the HFP maps for all non-Antarctic land areas35.
Within the main manuscript, we defined intact land as any 1 km2 pixel with a HFP value not higher than or equal to 4. Within this threshold, all areas with a HFP score higher than 4 are defined as nonintact. While previous analyses showed that a >4 score is a key threshold above which species extinction risk greatly increases39, we recognized that there is no one true threshold, which impacts all species equally. Some species may require no human pressure to successfully disperse, while others might successfully navigate through more intensively modified landscapes. Therefore, we conducted our analyses for two additional HFP thresholds. The first used a HFP score More