Philippa Kaur

Study areas: biodiversity hotspots
We focused on eight biodiversity hotspots21: those with at least 25% of their extent within the “tropical and subtropical moist broadleaf forests” biome44 and for which we obtained at least 1000 checklists from eBird (after applying the data selection procedure described below): Atlantic Forest, Tropical Andes, Tumbes-Chocó-Magdalena, and Mesoamerica (Americas); Eastern Afromontane (Africa); Western Ghats and Sri Lanka, Indo-Burma and Sundaland (Asia). Within each hotspot, we analysed only areas overlapping the “tropical and subtropical moist broadleaf forests” biome44 (Fig. 1, Supplementary Figs. 1, 4, and 5), assumed to have been originally forested (see Supplementary Methods 4d).
Data selection: eBird checklists
We obtained bird sightings from the eBird citizen science database23. The reporting system is based on checklists, whereby the observer provides: list of birds detected; GPS location; sampling effort (whether or not all detected species are reported; sampling duration; sampling protocol, e.g., stationary point, travel, and banding; and distance travelled in case of travelling protocol); starting time of the sampling event; and number of observers.
We used the eBird dataset released in December 201845, focusing on records from 2005 to 2018, as data collected prior to 2005 were too scarce for analysis. We filtered this dataset to obtain high-quality checklists comparable in protocol and effort: we selected complete checklists only (i.e. in which observers explicitly declare having reported all bird species detected and identified); following either the “stationary points” or the “travelling counts” protocol; with durations of continuous observation of 0.5–10 h; with observers travelling distances during the checklist 60% of the 1-km buffer around the point is forested47) versus non-forest ( More

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.
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

1.
Mitchard, E. T. A. The tropical forest carbon cycle and climate change. Nature 559, 527–534 (2018).
ADS  CAS  PubMed  Google Scholar 
2.
Saatchi, S. S. et al. Benchmark map of forest carbon stocks in tropical regions across three continents. Proc. Natl Acad. Sci. USA 108, 9899–9904 (2011).
ADS  CAS  PubMed  Google Scholar 

3.
Baccini, A. et al. Estimated carbon dioxide emissions from tropical deforestation improved by carbon-density maps. Nat. Clim. Change 2, 182–185 (2012).
ADS  CAS  Google Scholar 

4.
Liu, Y. Y. et al. Recent reversal in loss of global terrestrial biomass. Nat. Clim. Change 5, 470–474 (2015).
ADS  Google Scholar 

5.
Harris, N. L. et al. Baseline map of carbon emissions from deforestation in tropical regions. Science 336, 1573–1576 (2012).
ADS  CAS  PubMed  Google Scholar 

6.
Marco, M. D., Watson, J. E. M., Currie, D. J., Possingham, H. P. & Venter, O. The extent and predictability of the biodiversity–carbon correlation. Ecol. Lett. 21, 365–375 (2018).
PubMed  Google Scholar 

7.
Giardina, F. et al. Tall Amazonian forests are less sensitive to precipitation variability. Nat. Geosci. 11, 405–409 (2018).
ADS  CAS  Google Scholar 

8.
Erb, K.-H. et al. Unexpectedly large impact of forest management and grazing on global vegetation biomass. Nature 553, 73–76 (2018).
ADS  CAS  PubMed  Google Scholar 

9.
Zarin, D. J. et al. Can carbon emissions from tropical deforestation drop by 50% in 5 years? Glob. Change Biol. 22, 1336–1347 (2016).
ADS  Google Scholar 

10.
Chaplin-Kramer, R. et al. Degradation in carbon stocks near tropical forest edges. Nat. Commun. 6, 10158 (2015).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

11.
Brandt, M. et al. Satellite passive microwaves reveal recent climate-induced carbon losses in African drylands. Nat. Ecol. Evol. 2, 827 (2018).
PubMed  Google Scholar 

12.
Fan, L. et al. Satellite-observed pantropical carbon dynamics. Nat. Plants 5, 944–951 (2019).
CAS  PubMed  Google Scholar 

13.
Mitchard, E. T. et al. Markedly divergent estimates of Amazon forest carbon density from ground plots and satellites. Glob. Ecol. Biogeogr. 23, 935–946 (2014).
PubMed  PubMed Central  Google Scholar 

14.
Mitchard, E. T. et al. Uncertainty in the spatial distribution of tropical forest biomass: a comparison of pan-tropical maps. Carbon Balance Manag. 8, 10 (2013).
PubMed  PubMed Central  Google Scholar 

15.
Avitabile, V. et al. An integrated pan-tropical biomass map using multiple reference datasets. Glob. Change Biol. 22, 1406–1420 (2016).
ADS  Google Scholar 

16.
Réjou-Méchain, M. et al. Upscaling forest biomass from field to satellite measurements: Sources of errors and ways to reduce them. Surv. Geophys. 40, 881–911 (2019).
ADS  Google Scholar 

17.
Saatchi, S. Mapping tropical forest biomass: synthesis of ground and remote sensing inventory. Consult. Rep. 2 High Carbon Stock Sci. Study (2015).

18.
Ploton, P. et al. A map of African humid tropical forest aboveground biomass derived from management inventories. Sci. Data 7, 221 (2020).
PubMed  PubMed Central  Google Scholar 

19.
Philippon, N. et al. The light-deficient climates of Western Central African evergreen forests. Environ. Res. Lett. 14, 034007 (2018).

20.
Saatchi, S. et al. Seeing the forest beyond the trees. Glob. Ecol. Biogeogr. 24, 606–610 (2015).
Google Scholar 

21.
Mermoz, S., Le Toan, T., Villard, L., Réjou-Méchain, M. & Seifert-Granzin, J. Biomass assessment in the Cameroon savanna using ALOS PALSAR data. Remote Sens. Environ. 155, 109–119 (2014).
ADS  Google Scholar 

22.
Lewis, S. L. et al. Above-ground biomass and structure of 260 African tropical forests. Philos. Trans. R. Soc. B Biol. Sci. 368, 20120295 (2013).
Google Scholar 

23.
Hansen, M. C., Potapov, P. & Tyukavina, A. Comment on “Tropical forests are a net carbon source based on aboveground measurements of gain and loss”. Science 363, eaar3629 (2019).
CAS  PubMed  Google Scholar 

24.
Baccini, A. et al. Tropical forests are a net carbon source based on aboveground measurements of gain and loss. Science 358, 230–234 (2017).
ADS  MathSciNet  CAS  PubMed  PubMed Central  MATH  Google Scholar 

25.
Breiman, L. Random forests. Mach. Learn 45, 5–32 (2001).
MATH  Google Scholar 

26.
Lyapustin, A., Wang, Y., Korkin, S. & Huang, D. MODIS Collection 6 MAIAC algorithm. Atmos. Meas. Tech. 11, 5741–5765 (2018).

27.
Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1‐km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).
Google Scholar 

28.
Kühn, I. Incorporating spatial autocorrelation may invert observed patterns. Divers. Distrib. 13, 66–69 (2007).
Google Scholar 

29.
Dormann, C. F. Effects of incorporating spatial autocorrelation into the analysis of species distribution data. Glob. Ecol. Biogeogr. 16, 129–138 (2007).
Google Scholar 

30.
Roberts, D. R. et al. Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. Ecography 40, 913–929 (2017).
Google Scholar 

31.
Valavi, R., Elith, J., Lahoz‐Monfort, J. J. & Guillera‐Arroita, G. blockCV: an r package for generating spatially or environmentally separated folds for k-fold cross-validation of species distribution models. Methods Ecol. Evol. 10, 225–232 (2019).
Google Scholar 

32.
Parmentier, I. et al. Predicting alpha diversity of African rain forests: models based on climate and satellite-derived data do not perform better than a purely spatial model. J. Biogeogr. 38, 1164–1176 (2011).
Google Scholar 

33.
Baccini, A., Walker, W., Carvalho, L., Farina, M. & Houghton, R. A. Response to Comment on “Tropical forests are a net carbon source based on aboveground measurements of gain and loss”. Science 363, eaat1205 (2019).
CAS  PubMed  MATH  Google Scholar 

34.
Xu, L., Saatchi, S. S., Yang, Y., Yu, Y. & White, L. Performance of non-parametric algorithms for spatial mapping of tropical forest structure. Carbon Balance Manag. 11, 18 (2016).
PubMed  PubMed Central  Google Scholar 

35.
Dormann, C. F. et al. Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. Ecography 30, 609–628 (2007).
Google Scholar 

36.
Irwin, A. The ecologist who wants to map everything. Nature 573, 478–481 (2019).
ADS  PubMed  Google Scholar 

37.
Legendre, P. Spatial autocorrelation: trouble or new paradigm? Ecology 74, 1659–1673 (1993).
Google Scholar 

38.
Meyer, H., Reudenbach, C., Wöllauer, S. & Nauss, T. Importance of spatial predictor variable selection in machine learning applications–Moving from data reproduction to spatial prediction. Ecol. Model. 411, 108815 (2019).
Google Scholar 

39.
Strobl, C., Boulesteix, A.-L., Kneib, T., Augustin, T. & Zeileis, A. Conditional variable importance for random forests. BMC Bioinformatics 9, 307 (2008).
PubMed  PubMed Central  Google Scholar 

40.
Lefsky, M. A. et al. Estimates of forest canopy height and aboveground biomass using ICESat. Geophys. Res. Lett. 32, L22S02 (2005).

41.
Réjou-Méchain, M. et al. Upscaling Forest biomass from field to satellite measurements: sources of errors and ways to reduce them. Surv. Geophys. 40, 881–911 (2019).

42.
Mitchard, E. T. A. et al. Comment on ‘A first map of tropical Africa’s above-ground biomass derived from satellite imagery’. Environ. Res. Lett. 6, 049001 (2011).
ADS  Google Scholar 

43.
Asner, G. P. et al. High-resolution carbon mapping on the million-hectare Island of Hawaii. Front. Ecol. Environ. 9, 434–439 (2011).
Google Scholar 

44.
Asner, G. P. et al. Human and environmental controls over aboveground carbon storage in Madagascar. Carbon Balance Manag. 7, 2 (2012).
PubMed  PubMed Central  Google Scholar 

45.
Asner, G. P. et al. High-resolution mapping of forest carbon stocks in the Colombian Amazon. Biogeosciences 9, 2683 (2012).
ADS  CAS  Google Scholar 

46.
Asner, G. P. et al. Mapped aboveground carbon stocks to advance forest conservation and recovery in Malaysian Borneo. Biol. Conserv. 217, 289–310 (2018).
Google Scholar 

47.
Xu, L. et al. Spatial distribution of carbon stored in forests of the Democratic Republic of Congo. Sci. Rep. 7, 15030 (2017).
ADS  PubMed  PubMed Central  Google Scholar 

48.
Schepaschenko, D. et al. The Forest Observation System, building a global reference dataset for remote sensing of forest biomass. Sci. Data 6, 1–11 (2019).
Google Scholar 

49.
Chave, J. et al. Ground data are essential for biomass remote sensing missions. Surv. Geophys. 40, 863–880 (2019).
ADS  Google Scholar 

50.
van den Hoogen, J. et al. Soil nematode abundance and functional group composition at a global scale. Nature 572, 194–198 (2019).
ADS  CAS  PubMed  Google Scholar 

51.
Steidinger, B. S. et al. Climatic controls of decomposition drive the global biogeography of forest-tree symbioses. Nature 569, 404–408 (2019).
ADS  CAS  PubMed  Google Scholar 

52.
Bastin, J.-F. et al. The global tree restoration potential. Science 365, 76–79 (2019).
ADS  CAS  PubMed  Google Scholar 

53.
Trabucco, A. & Zomer, R. J. Global aridity index (global-aridity) and global potential evapo-transpiration (global-PET) geospatial database. CGIAR Consort Spat Information (2009).

54.
Wilson, A. M. & Jetz, W. Remotely sensed high-resolution global cloud dynamics for predicting ecosystem and biodiversity distributions. PLoS Biol. 14, e1002415 (2016).
PubMed  PubMed Central  Google Scholar 

55.
Bowman, D. M., Williamson, G. J., Keenan, R. J. & Prior, L. D. A warmer world will reduce tree growth in evergreen broadleaf forests: evidence from A ustralian temperate and subtropical eucalypt forests. Glob. Ecol. Biogeogr. 23, 925–934 (2014).
Google Scholar 

56.
Malhi, Y. et al. Exploring the likelihood and mechanism of a climate-change-induced dieback of the Amazon rainforest. Proc. Natl Acad. Sci. USA 106, 20610–20615 (2009).
ADS  CAS  PubMed  Google Scholar 

57.
Rennó, C. D. et al. HAND, a new terrain descriptor using SRTM-DEM: Mapping terra-firme rainforest environments in Amazonia. Remote Sens. Environ. 112, 3469–3481 (2008).
ADS  Google Scholar 

58.
Nachtergaele, F., Velthuizen, H. V., Verelst, L. & Wiberg, D. Harmonized World Soil Database (HWSD) (Food and Agriculture Organization, U. N. Rome, 2009).

59.
Defourny, P. et al. Algorithm Theoretical Basis Document for Land Cover Climate Change Initiative. Technical report (European Space Agency, 2014).

60.
Segal, M. & Xiao, Y. Multivariate random forests. WIREs Data Min. Knowl. Discov. 1, 80–87 (2011).
Google Scholar 

61.
CCI, ESA. New Release of 300 m Global Land Cover and 150 m Water Products (v.1.6.1) and new version of the User Tool (3.10) for Download (ESA CCI Land cover website, 2016). More

Global P losses from soils and soil P balances
All continents result in negative P balances (e.g., net P losses from agricultural systems, Table 1; Fig. 2) except Asia, Oceania and Australia, with Asia having a slightly positive but near zero P balance. This is in spite of high to very high chemical fertilizer inputs (with a range of 1.7 to 13 kg ha−1 yr−1 between the different continents, with national values reaching up to 14 and 19 kg ha−1 yr−1 for the European Union 15 old Member States (EU15, member states joining before 2004 mainly in western and northern Europe) and China, respectively). Most negative P balances are indicated for Africa due to very low chemical fertilizer input of 1.7 kg ha−1 yr−1 paired with high losses due to soil erosion of 9.6 kg ha−1 yr−1. South America as well as Central and Eastern Europe (NEU11; which are the new member states joining the EU after 2004 with the exception of Cyprus and Malta) also exhibit high P losses but for different reasons. South America has a very high chemical fertilizer input but also high losses due to soil erosion paired with high P exports due to organic P management (calculated as the sum of manure and residue input minus plant uptake). In contrast, the eastern European Union New Member States (NEU11) have rather low erosional losses but also very low chemical fertilizer input. With the hypothetical assumption of no replenishment due to chemical fertilizer (e.g., due to economic or technical constraints), calculation of soil P balances results in negative balances globally, as well as for all the continents and regions considered (depletion between 4 and 20 kg P ha−1 yr−1; Figs. 3 and 4; Table 1). The latter demonstrates the vulnerability of today’s global land management system and its strong dependency on chemical P fertilizers from non-renewable mineable P deposits.
Table 1 Main phosphorus (P) statistics in kg P ha−1yr−1 for all continents and selected countries.
Full size table

Fig. 2: Global average phosphorus (P) losses due to soil erosion in kg ha−1 yr−1.

The chromatic scale represents the P losses estimates, while the gray color indicates the cropland areas that were excluded from the modeling due to data unavailability. Note that classes are not regularly scale ranked but are divided into six classes using the quantile classification method. Only plant available fractions were considered. For the more residual P fractions please refer to Table 1 or Figs. 3 and 4).

Full size image

Fig. 3: Global P soil pools and depletion due to erosion.

Arrows indicate fluxes (positive: net input to soils, negative: depletion of soils). *Organic P management = sum of manure and residue input minus plant uptake. Non-plant P = non-plant available P. Inorganic and organic P give plant available fractions. Total soil P: sum of P fractions lost from soil via erosion with relative errors. No/with chemical = P balance with and without chemical fertilizer.

Full size image

Fig. 4: Soil P pools and depletion due to erosion in Africa, Europe and North America.

AD = Atmospheric Deposition. CF = Chemical Fertilizer. OM = Organic P management = sum of manure and residue input minus plant uptake. Arrows indicate fluxes (positive: net input to soils, negative: depletion of soils). Non-plant P = non-plant available P. Inorganic and organic P give plant available fractions. Soil Plost: sum of P fractions lost from soil via erosion with relative errors. No/with chemical = P balance with and without chemical fertilizer.

Full size image

Our area related calculations result in an average P loss for arable soils, due to erosion by water, of approximately (5.9_{ – 0.79% }^{ + 1.17% }) kg ha−1 yr−1 globally (Fig. 3, Table 1). This is around 60% of the rates given by Smil16 who estimated 10 kg P ha−1 yr−1 from arable fields due to soil erosion by water. Our total P losses due to soil erosion by water from arable soils globally result in 6.3 Tg yr−1 with 1.5 Tg yr−1 for organic and 4.8 Tg yr−1 for inorganic P. With these values the results are at the lower end of the range discussed in the literature (between 1–19 Tg yr−1, Table 2).
Table 2 Global fluxes from soils/arable systems to waters as discussed in recent literature. For comparability, all values were normalized to 1 billion ha of arable land (with original values given in brackets).
Full size table

Liu34 in estimating net input to, and output from, cropland systems indicated a net loss of P from the world’s croplands of about 12.8 Tg  yr−1 (calculating input from atmosphere, weathering and chemical fertilizer versus output from organic P management, soil erosion and runoff), which would be, according to their calculations, the same order of magnitude as synthetic fertilizer input (13.8 Tg yr−1 for statistical year 2003/2004). Our calculations result in an approximate net soil loss due to erosion of 6.3 Tg yr−1 with a global average chemical fertilizer input of 9.2 Tg yr−1. The fluxes clearly show the critical dependence on chemical fertilizer globally, with a hypothetical net-average-area related depletion of 10.7 kg ha−1 yr−1 globally without the compensation due to chemical fertilizer, from which a loss of 5.2 kg ha−1 yr−1 stems from organic P management (sum of manure and residue input minus plant uptake) and 5.9 kg ha−1 yr−1 from soil erosion (Table 1). Continental and national erosional P losses are between 40 and 85% of total P losses from agricultural systems with the exception of Europe and Australia (16 and 19%, respectively). Globally, as well as for Africa, South America and Asia, P soil losses due to erosion are higher than losses due to organic P management.
A recent quantification of atmospheric P dust input, based on dust measurements in the Sierra Nevada, concluded that measured dust fluxes are greater than, or equal to, modern erosional outputs and a large fractional contribution relative to bedrock35. However, even though their measured and modeled maximum atmospheric P fluxes were in the same order of magnitude as the atmospheric flux data of Wang et al.36 used in this study (0.1 kg ha−1 yr−1 for North America), erosion rates from sediment trap measurements in their investigated forest ecosystems were considerably lower (maximum of 0.06 kg ha−1 yr−135) than might be expected in arable lands worldwide. As such, we can clearly not agree with the above conclusion.
Mitigation of the soils P status in the long term by decreasing the deficits of the current organic P management seems difficult and rather unlikely in many regions of the world (see discussion of continental balances below). As such, and considering the expected shortage of P supply from industrial fertilizers in the future, the evaluation of P fluxes clearly shows that soil erosion has to be limited to the feasible absolute minimum in the future.
Only a small fraction of the total soil P is plant available, because a large fraction is either bound in, adsorbed to, or made unavailable by occlusion in minerals (apatite and occluded P)37. Our modeling approach follows the approach of Yang et al.37 that builds on existing knowledge of soil P processes and data bases to provide spatially explicit estimates of different forms of naturally occurring soil P on the global scale. Yang et al.37 uses data acquired with the Hedley fractionation method38 which splits soil P into different fractions that are extracted sequentially with successively stronger reagents and which are merged into various inorganic and organic pools. There is great uncertainty to associate these fractions with functional plant uptake39,40 but some recent work aims at quantifying residence times of the different fractions17,41. We present the total P loss as well as the plant-available pool (Figs. 3 and 4, Table 1) as the most labile P which can participate to plant nutrition in a time scale up to months17 (i.e., the so-called labile inorganic P, inorganic P bound to secondary minerals, labile and stable organic P33,37). We contrast these more labile, short lived fractions with fractions considered very stable which would not be plant available at short time scales (i.e., inorganic P associated with minerals such as apatite or occluded P17, corresponding to 52% of total soil P in our estimates at the global scale).
Verification of soil P loss with a comparison to riverine P exports
Verification of our proposed P soil losses might be done indirectly via a comparison to riverine P loads. However, to verify our on-site P soil loss approach with off-site P river export data requires two prerequisites: (i) an estimation of the agricultural erosion and runoff contribution to total P loads in the rivers (as the latter will be the total sum of agricultural, urban and industrial runoff) and (ii) an assumption on sediment delivery rates (e.g., the percent of sediments reaching the rivers from the total eroded sediments) to estimate P loads to rivers from on-site P soil loss. Regarding the first prerequisite, we searched for published river P export data with a separation of the agricultural erosion and runoff contribution from total P loads in the rivers. Regarding the second prerequisite, there is a lack of models describing integrated sediment-delivery to rivers on a continental or global scale, but it has been discussed that sediment delivery ratios generally decrease as drainage area increases, ranging from roughly 30–100% in small catchments (≤0.1 km2) to 2–20% at large spatial scales (e.g., ≥1000 km2)42. For our comparison, we used a range of average sediment delivery rates between 11–30% as used in recent large scale studies from continental to global scale43,44. In doing so, we would like to point out, that even if 70–89% of the P lost from agricultural soils might be re-deposited within catchments, potential threats of P loss from soils are not reduced. Erosion (and thus P loss) occurs predominantly on agricultural soils while re-deposition will mostly occur in depositional hollows, wetlands, riparian zones or buffer strips. Thus, P is lost as a nutrient on food and feed production sites but re-deposited as a potential ecological threat to biodiversity and ecosystem health due to its eutrophication effect in less intensively or unmanaged ecosystems. Last but not least we would like to point out that RUSLE only considers soil displacement due to rill and inter-rill erosion neither considering tillage and gully erosion nor land sliding.
A comparison of calculated potential P export to rivers from our on-site soil P losses is well within the range of published riverine P exports (Table 3). Beusen et al.45 used total suspended sediment measurements from the GEMS-GLORI database to extrapolate spatially distributed sediment rates for the world’s largest rivers with land use, topography, lithology and precipitation as factors in a multiple linear regression approach accounting for soil erosion as well as sediment trapping. The associated nutrient exports for all continents, as well as global assessment were made by calibrating nutrient export to sediment rates with the model Global News45 using established correlations between sediment and nutrient concentrations. In comparing their P export with suspended sediments in rivers45 to our assessments, we underestimated P export globally as well as for all continents, with the exception of Africa (Table 3), which is, however, (i) strongly related to assumed sediment delivery rates and (ii) our RUSLE application only considering rill and inter-rill erosional processes with an unknown contribution of gullies, landslides and tillage erosion.
Table 3 On-site P loss from gross soil erosion (this study), calculated potential riverine loads with sediment delivery ratios between 11–30%43,44 and comparison to global and regional riverine P export studies. All values in kg P ha−1yr−1.
Full size table

An analysis of 17 large scale European catchments (250–11,000 km2) quantifying the loss of P to surface waters from an off-site perspective (P flux measurement in waters) resulted in 0.05–1.5 kg ha−1 yr−1 of P loss due to soil erosion and runoff from agricultural lands into streams and rivers (excluding a Greek catchment with a very high P export of 6 kg ha−1 yr−1)46. Recalculating our on-site soil P losses with sediment delivery ratios between 11–30% results in rates for geographic Europe between 0.1–0.4 kg ha−1 yr−1 which are at the lower end of the range of Kronvang, et al.46 (Table 3). Potential P losses as riverine export due to agricultural runoff and erosion of 143 watersheds across the U.S.47 are in the same range as the calculated P loss assessment of our study, while P loads to Lake Erie assessed in a regional study48 seem slightly higher (Table 3). However, no partitioning between agricultural, urban or industrial fluxes was possible from the latter study. An assessment of the Yangtze River Basin (with 1.8 Million km2 near 20% of the whole Chinese territory) gives modeled soil P losses (on-site perspective) between 0–196 kg ha−1 yr−1 demonstrating the huge spatial heterogeneity of on-site soil erosion rates49. The model output was calibrated with total measured nutrient loads in rivers, while the partitioning differed between dissolved point and non-point as well as adsorbed non-point pollution which can be mostly attributed to soil erosion of agricultural fields. The average P loads due to soil erosion from agricultural fields in the Yangtze River Basin compares well with the range assessed in our study for China (2.7 versus a range of 1.4 to 3.7 kg ha−1 yr−1, respectively). The same holds true for assessment comparisons of Africa50 and South America51 (Table 3), even though it should be considered that especially for these latter studies scale differs considerably from our approach.
Regional P losses and balances
Parallel to the distribution pattern and dynamics of global soil erosion by water32, P losses from soils due to water erosion are most dramatic in countries and regions with intensive agriculture and/or extreme climates (e.g., droughts followed by significant rain events or high frequencies of heavy rain storms) due to high erosivity effects52. As such, our calculations result in extremely high P losses due to erosion ( >20 kg ha−1 yr−1) in regions such as eastern China, many regions in Indonesia, parts of east and south-eastern Africa (Ethiopia, Eritrea, Mozambique), Central America and parts of South America (South-Eastern Brazil; Southern Chile, Peru (Fig. 2)). A very high P loss (10 to 20 kg ha−1 yr−1) is estimated for parts of Southern Africa (South Africa, Madagascar, Tanzania) and South America (Bolivia) and a high loss (5–10 kg ha−1 yr−1) for most of India, as well as regions in Southern Africa (Angola, Zambia) and South America (Uruguay) (Fig. 2). Even though the underlying erosion model algorithm does not calculate the net catchment output but rather the on-site displacement of soil sediments which might then be re-located to other parts of the fields or even buried at depositional places, the considered on-site field management will clearly be confronted with substantial P losses due to soil erosion by water. Only considering agronomic P inputs and outputs without including P losses due to erosion by MacDonald et al.22 resulted in a very different global P pattern: most widespread large deficits were in South America (North-Eastern countries, e.g., Argentina and Paraguay), the northern United States and Eastern Europe while the largest surpluses covered most of East Asia, Western and Southern Europe, the coastal United States, South-Eastern Brazil and Uruguay.
With average soil depletion due to erosion of 9.6 kg ha−1 yr−1 in Africa, the overall P balance is already negative by 9.7 kg ha−1 yr−1 today (Table 1, Fig. 4). As the average P depletion in Africa due to negative fluxes in organic P management equals the input fluxes from the atmosphere plus chemical fertilizer, African farmers could decrease P losses to near zero with effective soil erosion mitigation. Even though the system’s P depletion due to organic P management is relatively low in Africa (−2.2 kg ha−1 yr−1) compared to a global average (−5.2 kg ha−1 yr−1), the high overall P losses are unlikely to be covered neither from a mitigated and more sustainable organic P management nor from increased chemical fertilizer input. P fluxes due to organic P management are calculated here as the sum of manure and residue input minus plant uptake (which results in biomass export in arable systems with the exception of residues left on the field). The overall sum of plant uptake is likely to increase with increased need for food and feed parallel to a predicted population and livestock growth in Africa in the future. Many soils in sub-Saharan Africa have already been characterized as deficient for levels of plant-available P for the last decades53. Manure and residue input is simultaneously in demand in Africa today (shortage of biomass in general, low animal production and even if there is manure available, there are no means to transport it to where it is needed), which results in the recommendations of an integrated farm management with combinations of organic and inorganic fertilizers54,55,56. With the inorganic P fertilizers becoming increasingly scarce, the depletion due to organic P management can be expected to increase in Africa in the future. Simultaneously, today’s prices for chemical fertilizer can already be 2–6 times more expensive for a farmer in Africa than in Europe due to higher transport and storage costs3, even though Africa itself has the highest geological P deposits in the world (according to today’s estimates 80% of the global geological P deposits are located in Morocco and the Western Sarah8). As such, and if the political situation does not change dramatically (e.g., that the P supplies are marketed within Africa instead of being exported to US, Europe and China), the only realistic means of reducing P depletion of African soils today, and in the future, is to drastically reduce soil erosion.
We recognize that the values calculated in Table 1 and Fig. 4 are gross estimates over large scales and that spatial context and scale, especially on the African continent is important. P deficiency is a country, district, farm and soil specific issue in Africa, for example parts of east Africa and the Sahel have substantial deficiencies57. In sub Saharan Africa ~40% of soils are considered to have low nutrient reserves ( More

Regionally-downscaled climate change projections
We used a management strategy evaluation (MSE) applied to ensemble projections of a climate-enhanced multispecies stock assessment within the integrated modeling framework of the Alaska Climate Change Integrated Modeling project (ACLIM)19. For this, six high resolution downscaled projections of oceanographic and lower trophic level conditions in the Bering Sea (using the Regional Ocean Modeling System49,50) were coupled to the BESTNPZ nutrient-phytoplankton-zooplankton model51; we refer to this model complex throughout this paper as the Bering10K ROMSNPZ, or just ROMSNPZ, model. Boundary conditions were driven by three global general circulation models (GFDL-ESM2M52, CESM153, and MIROC-ESM54) projected (2006–2099) under the high-baseline emission scenario Representative Concentration Pathway 8.5 (RCP 8.5) and midrange global carbon mitigation (RCP 4.5; note, that for CESM under RCP 4.5, projections from 2080–2100 were unavailable so conditions from 2080–2099 were held constant at 2080 conditions for that scenario only) future scenarios from the Coupled Model Intercomparison Project phase 5 (CMIP5)29,55. Hermann et al.30,56 also report on downscaled hindcasts of oceanographic and lower trophic conditions in the EBS from 1970–2012 (see refs. 30,56,57 for detailed descriptions of model evaluation and performance). For each downscaled model simulation, we replicated the National Marine Fisheries Service Alaska Fisheries Science Center annual summer bottom-trawl survey in time and space in the ROMSNPZ model (using historical mean survey date at each latitude and longitude of each gridded survey station) to derive estimates of sea surface and bottom temperatures (Fig. 1). We additionally used a polygon mask of the survey area to estimate the average zooplankton abundance in the system during spring, summer, winter, and fall months. These indices were derived for each climate projection scenario, as well as a persistence scenario where conditions were held constant at the average of those for 2006–2017 from a hindcast simulation. All index projections were bias corrected to the 2006–2017 hindcast period using the delta method assuming unequal variance in the GCM projections and hindcast58 such that:

$$T_{mathop {{{mathrm{fut}}}}limits^prime ,y} = bar T_{{mathrm{hind}},overrightarrow {scriptstyle{mathrm{ref}}} } + frac{{sigma _{{mathrm{hind}},overrightarrow {scriptstyle{mathrm{ref}}} }}}{{sigma _{{mathrm{fut}}, overrightarrow{scriptstyle{mathrm{ ref}}} }}}left( {T_{{mathrm{fut}},y} – bar T_{{mathrm{fut}},overrightarrow {scriptstyle{mathrm{ref}}} }} right)$$
(1)

where (T_{mathop {{{mathrm{fut}}}}limits^prime ,y}) is the bias-corrected projected timeseries, (T_{{mathrm{fut}},y}) is the raw projected timeseries, (bar T_{{mathrm{hind}},overrightarrow {scriptstyle{mathrm{ref}}} }) is the mean of the hindcast during the reference years (overrightarrow {{mathrm{ref}}}) (2006–2017), (bar T_{{mathrm{fut}},overrightarrow {scriptstyle{mathrm{ref}}} }) is the mean of the raw projected timeseries during the reference years (overrightarrow {{mathrm{ref}}}), (sigma _{{mathrm{hind}},overrightarrow {scriptstyle{mathrm{ref}}} }) is the standard deviation of the hindcast during the reference years (overrightarrow {{mathrm{ref}}}), (sigma _{{mathrm{fut}},overrightarrow {scriptstyle{mathrm{ref}}} }) is the standard deviation of the raw projection timeseries during the reference years (overrightarrow {{mathrm{ref}}}).
Climate-enhanced multispecies stock assessment model
Bias-corrected indices were then used as covariates in the climate-enhanced multispecies stock assessment model for the Bering Sea (hereafter CEATTLE)59 to evaluate the performance of alternative management approaches on future fish biomass and catch. CEATTLE is a climate-enhanced multispecies statistical age-structured assessment model with parameters for growth that are functions of temperature (i.e., temperature-specific average weight-at-age) and predation that are functions of temperature (via a bioenergetics-based predation sub-model)59,60,61. Since 2016, the model has been used operationally in the Bering sea as a supplement to the annual BSAI pollock stock assessment61. Various configurations of CEATTLE are possible; for this study we chose one where temperature-specific predator and prey interactions influenced natural mortality, temperature influenced weight-at-age, and the spawner-recruit relationship was a function of physical and biological future conditions as well as random variability (i.e., a climate-informed multispecies model). We fit the model using penalized maximum likelihood to survey biomass, diet, and fishery harvest data for three groundfish species pollock, Pacific cod, and arrowtooth flounder from the EBS in the EBS over the period 1979–2017. We also used the Bering10K ROMSNPZ model to produce detailed hindcasts of temperature for the period 1970–2017. We used hindcast-extracted timeseries from the ROMSNPZ model and CEATTLE model estimates of recruitment ((R_{i,y,l})) and spawning biomass ((B_{i,y – 1})) in hindcast year y for each species i to fit a climate-enhanced logistic recruitment per spawner model36, such that:

$$ln left( {hat R_{i,y}} right) = alpha _i – beta _{0,i}B_{i,y – 1} + ln left( {B_{i,y – 1}} right) + {{{bf{B}}}}_i{{{bf{X}}}} + varepsilon _{i,y}$$
(2)

where ({{{bf{B}}}}_i{{{bf{X}}_l}}) is the summed product of each covariate parameter (beta _{ij}) and the corresponding environmental covariate (X_{j,y}) for each bias-corrected environmental index (j = ( {1,2…n_j} )). We selected indices representative of ecological conditions important for groundfish recruitment in the Bering sea39; spring and fall large zooplankton abundances, survey replicated bottom temperature, and extent of the residual cold pool of extremely dense and cold sea water that persists across the EBS shelf following spring sea ice retreat. We assumed normally distributed (in log space) residual errors for each species ((varepsilon _{i,y} sim Nleft( {0,sigma _i^2} right))). The CEATTLE model was then projected forward where ROMSNPZ indices from individual projections drove growth, predation, and recruitment in each future simulation year36,62.
Evaluation of harvest management approaches
Previous authors have defined EM (i.e., the incorporation of ecosystem information into marine resource management) as a continuum between two paradigms of management and focus18. On one end is within-sector single-species management that considers ecosystem information (EAFM) and on the other is cross-sectoral whole of ecosystem management (i.e., EBM). EBFM is intermediate between these and is defined by quantitative incorporation of ecosystem interactions into assessment models and target setting (EBFM). Most fisheries management in the Bering Sea can be characterized as EBFM or EAFM, with increasing trends toward cross-sectoral coordination at the scale of EBM. Here we focus on one aspect on this scale of potential management options, operational EBFM and EAFM as captured through the CEATTLE multispecies stock assessment model and harvest policies decisions made annually under the constraint of the 2 MT cap (modeled via the ATTACH model).
MSE is a process of “assessing the consequences of a range of management strategies or options and presenting the results in a way which lays bare the tradeoffs in performance across a range of management objectives”63. MSE has been frequently used to evaluate alternative management strategies based on single-species estimation methods64. It is increasingly used to evaluate ecosystem management performance, although these evaluations are far less commonplace due to the complexity of modeling and assessing the performance of ecosystem level metrics64. Importantly, MSE “does not seek to proscribe an optimal strategy or decision”63, rather it aims to describe the uncertainty and tradeoffs inherent in alternative strategies and scenarios. In this case, through a series of workshops, we worked with managers and stakeholders to identify priority scenarios and outputs19. From this, risk, sensitivity, and uncertainty under contrasting climate scenarios were requested outputs of the analysis, as was the performance of current climate-naive EBFM policies.
A key component of MSE is identifying and quantifying uncertainty (i.e., process, observation, estimation, model, and implementation error) and representing it using an operating model. In the case of this MSE, the focus was on process error uncertainty due to variation in recruitment about the fitted stock-recruitment relationship, one major source of model error in the form of alternative climate scenarios, and implementation error. The MSE does not account for estimation error (uncertainty in the parameters of the operating model) nor observation error. This is because the estimates of recruitment and spawning biomass from CEATTLE for the BSAI are very precise (see Fig. 10 in ref. 60) and the estimation and operating models are therefore very similar. Thus, CEATTLE is the operating model for this MSE and implicitly the estimation method. In this approach we assume that while allowing for observation error would have increased overall error, the effect would have been minor compared to the investigated uncertainties. Future analyses using a full MSE (i.e., separate operating and estimation models) could evaluate the effect of observation error, but perhaps more importantly, the potential for model error, whereby the population dynamics model (on which the estimation method is based) differs from that of the operating model such that the estimates on which management decisions are made are biased relative to the true values in the operating model.
Given this we summarized the relative change in catch and biomass for the three species in the model under the following fishing scenarios (Fig. 1): (a) projections without harvest ((F_{i,y} = 0)) in each year y of scenario l for each species i, (b) projections under target harvest rate (Supplementary Fig. 7 left) and with a sloping harvest control rule (HCR) (Supplementary Fig. 7 right), (c) as in 2 but with the constraint of a 2 MT cap applied dynamically to the three focal species only.
Under the North Pacific Fishery Management Council (NPFMC) constraint of the 2 MT cap on cumulative total annual catch, realized harvest (i.e., catch) and specification of individual species harvest limits known as Total Allowable Catch (TAC; metric tons) are a function of the acceptable biological catch (ABC) for the given species, as well as ABC of other valuable species in the aggregate complex19,34 (https://github.com/amandafaig/catchfunction). TAC must be set at or below ABC for each species, therefore TAC of individual species are traded-off with one another to avoid exceeding the 2 MT cap. From 1981 to 1983, the TAC of pollock was reduced significantly below the ABC and in 1984 the 2 MT cap became part of the BSAI fishery management plan21,34,65. Pacific cod regulations have changed markedly over recent decades and it was only in the 1990s that in many years the catch and TAC approached its ABC. Thus, we used the socioeconomic ATTACH model (the R package ATTACHv1.6.0 is available with permission at https://github.com/amandafaig/catchfunction34) to model realized catch in each simulation year as a function of CEATTLE assessment estimates of ABC (tons) for pollock, Pacific cod, and arrowtooth flounder under future projections (2018–2100). This entailed three steps for each future simulation year (y):
1.
project the population forward from (y – 1) to (y) using estimated parameters from the multispecies mode of the CEATTLE model fit to data from 1979 to 2017 and recruitment based on biomass in simulation year (y) and future environmental covariates from the ROMSNPZ model downscaled projections (see “Methods” above) to determine ({mathrm{ABC}}_{i,y,l}) for each species (i) under each scenario (l) given the sloping harvest control rule for pollock, Pacific cod, and arrowtooth flounder in each simulation year (y);

2.
use ({mathrm{ABC}}_{i,y,l}) of each species from step 1 as inputs to the ATTACH model in order to determine the North Pacific Marine Fishery Council Total Allowable Catch (TACi,y,l) for the given simulation l year y;

3.
use TACi,y,l from step 2 to estimate catch (tons) in the simulation year (Fig. 1); remove catch from the population and advance the simulation forward 1 year.

Determine the annual ABC
We used end-of-century projections (2095–2099) to derive a maximum sustainable yield (MSY) proxy for future harvest recommendations (ABCi,y,l) for each scenario l. To replicate current management, we used a climate-specific harvest control rule that uses climate-naive unfished and target spawning biomass reference points ((B_{0,i}) and (B_{{mathrm{target},i}}), respectively) and corresponding harvest rates ((F_{i,y} = 0) and (F_{i,y} = F_{{mathrm{target}}}) and (B_{i,y,l}) in each simulation l year (y) for each species (i)

$${mathrm{ABC}}_{i,y,l} = mathop {sum }limits_a^{A_i} left( {frac{{S_{i,a}F_{{mathrm{ABC}},i,y,l}}}{{Z_{i,a,y,l}}}left( {1 – e^{ – Z_{i,a,y,l}}} right)N_{i,a,y,l}W_{i,a,y,l}} right)$$
(3)

where (W_{i,a,y,l}), (N_{i,a,y}), and (Z_{i,a,y,l}) is the climate-simulation specific annual weight, number, and mortality (i.e., influenced through temperature effects on recruitment, predation, and growth) at age (a) for (A_i) ages in the model, and (S_{i,a}) is the average fishery age selectivity from the estimation period 1979–201759,60. (F_{{mathrm{ABC}},i,y,l}) and was determined in each simulation timestep using an iterative approach66 whereby we: (i) first determined average (B_{0,i}) values in years 2095–2099 by projecting the model forward without harvest (i.e., (F_{i,y} = 0)) for each species under the persistence scenario. We then (ii) iteratively solved for the harvest rate that results in an average spawning biomass ((B_{i,y})) during 2095–2099, that is, 40% of (B_{0,i}) (i.e., (F_{{mathrm{target}},i})) for pollock and Pacific cod simultaneously, with arrowtooth flounder (F_{i,y}) set to the historical average (as historical F for arrowtooth flounder ≪(F_{40% })); once (F_{{mathrm{target}},i}) for pollock and Pacific cod were found, we then iteratively solved for (F_{{mathrm{target}},i}) for arrowtooth flounder (Supplementary Fig. 7 left panel)59,60. Last, (iii) to derive a climate-informed ({mathrm{ABC}}_{i,y,l}) in each simulation year, the North Pacific Marine Fisheries Council (hereafter, “Council”) Tier 3 sloping harvest control rule with an ecosystem cutoff at (B_{20% }) was applied to adjust (F_{{mathrm{ABC}},i,y,l}) lower than (F_{{mathrm{target}},i}) if the simulation specific (climate-informed) (B_{i,y,l}) was lower than 40% of the climate-naive (B_{0,i}) at the start of a given year or set to 0 if (B_{i,y,l} , More