Ecology

Kwiatkowski, L. et al. Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model projections. Biogeosciences 17, 3439–3470 (2020).ADS 
CAS 

Google Scholar 
Intergovernmental Panel on Climate Change. Climate Change 2013: 5th Assessment Report of the Intergovernmental Panel on Climate Change (Cambridge University Press, 2013).
Google Scholar 
Torres, O., Kwiatkowski, L., Sutton, A. J., Dorey, N. & Orr, J. C. Characterizing mean and extreme diurnal variability of ocean CO2 system variables across marine environments. Geophys. Res. Lett. 48, 2 (2021).
Google Scholar 
Dorey, N., Lançon, P., Thorndyke, M. & Dupont, S. Assessing physiological tipping point of sea urchin larvae exposed to a broad range of pH. Glob. Change Biol. 19, 3355–3367 (2013).
Google Scholar 
Hauri, C. et al. Spatiotemporal variability and long-term trends of ocean acidification in the California current system. Biogeosci. Discuss. 9, 10371–10428 (2012).ADS 

Google Scholar 
Dupont, S. & Pörtner, H.-O. A snapshot into ocean acidification research. Mar. Biol. 160, 1765–1771 (2013).CAS 

Google Scholar 
Dupont, S. & Thorndyke, M. Chapter: Direct impacts of near-future ocean acidification on sea urchins. in Climate Change Perspective from the Atlantic: Past, Present and Future (eds. Fernández-Palacios, J. et al.) 461–485 (2013).Byrne, M. & Hernández, J. C. Chapter 16: Sea urchins in a high CO2 world: Impacts of climate warming and ocean acidification across life history stages. in Developments in Aquaculture and Fisheries Science vol. 43 281–297 (Elsevier, 2020).Kroeker, K. J., Kordas, R. L., Crim, R. N. & Singh, G. G. Meta-analysis reveals negative yet variable effects of ocean acidification on marine organisms. Ecol. Lett. 13, 1419–1434 (2010).PubMed 

Google Scholar 
L. Kelley, A., J. Lunden, J., 1 Ocean Acidification Research Center, College of Fisheries and Ocean Sciences, University of Alaska, Fairbanks, Fairbanks, AK, 99775, USA, & 2 Haverford College, Haverford, PA, 19041, USA. Meta-analysis identifies metabolic sensitivities to ocean acidification. AIMS Environ. Sci. 4, 709–729 (2017).Stumpp, M. et al. Acidified seawater impacts sea urchin larvae pH regulatory systems relevant for calcification. Proc. Natl. Acad. Sci. U. S. A. 109, 18192–18197 (2012).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Stumpp, M. et al. Digestion in sea urchin larvae impaired under ocean acidification. Nat. Clim. Change 3, 1044–1049 (2013).ADS 
CAS 

Google Scholar 
Runcie, D. E. et al. Genomic characterization of the evolutionary potential of the sea urchin Strongylocentrotus droebachiensis facing ocean acidification. Genome Biol. Evol. 8, 272 (2017).
Google Scholar 
Sewell, M. Utilization of lipids during early development of the sea urchin Evechinus chloroticus. Mar. Ecol. Prog. Ser. 304, 133–142 (2005).ADS 
CAS 

Google Scholar 
Lucas, M. I., Walker, G., Holland, D. L. & Crisp, D. J. An energy budget for the free-swimming and metamorphosing larvae of Balanus balanoides (Crustacea: Cirripedia). Mar. Biol. 55, 221–229 (1979).
Google Scholar 
Shilling, F. M., Hoegh-Guldberg, O. & Manahan, D. T. Sources of energy for increased metabolic demand during metamorphosis of the abalone Haliotis rufescens (Mollusca). Biol. Bull. 191, 402–412 (1996).CAS 
PubMed 

Google Scholar 
Meidel, S. K. & Scheibling, R. E. Effects of food type and ration on reproductive maturation and growth of the sea urchin Strongylocentrotus droebachiensis. Mar. Biol. 134, 155–166 (1999).
Google Scholar 
Pearce, C. M. & Scheibling, R. E. Induction of metamorphosis of larvae of the green sea urchin, Strongylocentrotus droebachiensis by coralline red algae. Biol. Bull. 179, 304–311 (1990).CAS 
PubMed 

Google Scholar 
Gosselin, P. & Jangoux, M. From competent larva to exotrophic juvenile: a morphofunctional study of the perimetamorphic period of Paracentrotus lividus (Echinodermata, Echinoida). Zoomorphology 118, 31–43 (1998).
Google Scholar 
Hinegardner, R. T. Growth and development of the laboratory cultured sea urchin. Biol. Bull. 137, 465–475 (1969).CAS 
PubMed 

Google Scholar 
Strathmann, R. R. Length of pelagic period in echinoderms with feeding larvae from the Northeast Pacific. J. Exp. Biol. Ecol. 34, 23–27 (1978).
Google Scholar 
Byrne, M. et al. Unshelled abalone and corrupted urchins: Development of marine calcifiers in a changing ocean. Proc. Biol. Sci. 278, 2376–2383 (2011).PubMed 

Google Scholar 
Dupont, S., Dorey, N., Stumpp, M., Melzner, F. & Thorndyke, M. Long-term and trans-life-cycle effects of exposure to ocean acidification in the green sea urchin Strongylocentrotus droebachiensis. Mar. Biol. 160, 1835–1843 (2013).CAS 

Google Scholar 
Uthicke, S. et al. Impacts of ocean acidification on early life-history stages and settlement of the coral-eating sea star Acanthaster planci. PLoS ONE 8, e82938 (2013).ADS 
PubMed 
PubMed Central 

Google Scholar 
Dupont, S., Lundve, B. & Thorndyke, M. Near future ocean acidification increases growth rate of the lecithotrophic larvae and juveniles of the sea star Crossaster papposus. J. Exp. Zool. Mol. Dev. Evol. 314, 382–389 (2010).
Google Scholar 
Lim, Y.-K., Dang, X. & Thiyagarajan, V. Transgenerational responses to seawater pH in the edible oyster, with implications for the mariculture of the species under future ocean acidification. Sci. Total Environ. 782, 146704 (2021).ADS 
CAS 
PubMed 

Google Scholar 
Hettinger, A. et al. Persistent carry-over effects of planktonic exposure to ocean acidification in the Olympia oyster. Ecology 93, 2758–2768 (2012).PubMed 

Google Scholar 
Hettinger, A. et al. Larval carry-over effects from ocean acidification persist in the natural environment. Glob. Change Biol. https://doi.org/10.1111/gcb.12307 (2013).Article 

Google Scholar 
Albright, R. & Langdon, C. Ocean acidification impacts multiple early life history processes of the Caribbean coral Porites astreoides. Glob. Change Biol. 17, 2478–2487 (2011).ADS 

Google Scholar 
Yuan, X. et al. Elevated CO2 delays the early development of scleractinian coral Acropora gemmifera. Sci. Rep. 8, 2787 (2018).ADS 
PubMed 
PubMed Central 

Google Scholar 
Maboloc, E. A. & Chan, K. Y. K. Parental whole life cycle exposure modulates progeny responses to ocean acidification in slipper limpets. Glob. Change Biol. 2, 15647. https://doi.org/10.1111/gcb.15647 (2021).Article 

Google Scholar 
Mos, B., Byrne, M. & Dworjanyn, S. A. Effects of low and high pH on sea urchin settlement, implications for the use of alkali to counter the impacts of acidification. Aquaculture 528, 735618 (2020).CAS 

Google Scholar 
Harianto, J., Aldridge, J., Torres Gabarda, S. A., Grainger, R. J. & Byrne, M. Impacts of acclimation in warm-low pH conditions on the physiology of the sea urchin Heliocidaris erythrogramma and carryover effects for juvenile offspring. Front. Mar. Sci. 7, 588938 (2021).
Google Scholar 
Houlihan, E. P., Espinel-Velasco, N., Cornwall, C. E., Pilditch, C. A. & Lamare, M. D. Diffusive boundary layers and ocean acidification: Implications for sea urchin settlement and growth. Front. Mar. Sci. 7, 577562 (2020).
Google Scholar 
Norderhaug, K. M. & Christie, H. C. Sea urchin grazing and kelp re-vegetation in the NE Atlantic. Mar. Biol. Res. 5, 515–528 (2009).
Google Scholar 
Dickson, A., Sabine, C. L. & Christian, J. R. Guide to best practices for ocean CO2 measurements. (PICES Special Publication 3;191 pp, 2007).Lavigne, H. & Gattuso, J.-P. seacarb: seawater carbonate chemistry with R. R package version 2.4. http://CRAN.R-project.org/package=seacarb. (2011).R Core Team. R: A language and environment for statistical computing. R: A language and environment for statistical computing (2017).Guillard, R. R. L. & Ryther, J. H. Studies of marine planktonic diatoms. I. Cyclotella nana Hustedt and Detonula confervacea (Cleve) Gran. Can. J. Microbiol. 8, 229–239 (1962).CAS 
PubMed 

Google Scholar 
Stumpp, M., Wren, J., Melzner, F., Thorndyke, M. & Dupont, S. CO2 induced seawater acidification impacts sea urchin larval development I: Elevated metabolic rates decrease scope for growth and induce developmental delay. Comp. Biochem. Physiol. Mol. Integr. Physiol. 160, 331–340 (2011).CAS 

Google Scholar 
His, E., Heyvang, I., Geffard, O. & De Montaudouin, X. A comparison between oyster (Crassostrea gigas) and sea urchin (Paracentrotus lividus) larval bioassays for toxicological studies. Water Res. 33, 1706–1718 (1999).CAS 

Google Scholar 
U. S. National Institutes of Health, Bethesda, Maryland, U. ImageJ, Rasband, W.S., http://imagej.nih.gov/ij/.Smith, M. M., Cruz Smith, L., Cameron, R. A. & Urry, L. The larval stages of the sea urchin, Strongylocentrotus purpuratus. J. Morphol. 269, 713–733 (2008).PubMed 

Google Scholar 
Kahm, M., Hasenbrink, G., Lichtenberg-Frate, H., Ludwig, J. & Kschischo, M. grofit: Fitting Biological Growth Curves with R. J. Stat. Softw., 33(7), 1–21. URL http://www.jstatsoft.org/v33/i07/. (2010).Pinheiro, J., Bates, D., & R-core. Package ‘nlme’: Linear and Nonlinear Mixed Effects Models. Cran-R (2018).Pan, T.-C.F., Applebaum, S. L. & Manahan, D. T. Experimental ocean acidification alters the allocation of metabolic energy. Proc. Natl. Acad. Sci. U. S. A. 112, 4696–4701 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Jager, T., Ravagnan, E. & Dupont, S. Near-future ocean acidification impacts maintenance costs in sea-urchin larvae: Identification of stress factors and tipping points using a DEB modelling approach. J. Exp. Mar. Biol. Ecol. 474, 11–17 (2016).
Google Scholar 
Hoegh-Guldberg, O. & Emlet, R. B. Energy use during the development of a lecithotrophic and a planktotrophic echinoid. Biol. Bull. 192, 27–40 (1997).CAS 
PubMed 

Google Scholar 
Vaïtilingon, D. et al. Effects of delayed metamorphosis and food rations on the perimetamorphic events in the echinoid Paracentrotus lividus (Lamarck, 1816) (Echinodermata). J. Exp. Mar. Biol. Ecol. 262, 41–60 (2001).
Google Scholar 
García, E., Clemente, S. & Hernández, J. C. Ocean warming ameliorates the negative effects of ocean acidification on Paracentrotus lividus larval development and settlement. Mar. Environ. Res. 110, 61–68 (2015).PubMed 

Google Scholar 
Wangensteen, O. S., Dupont, S., Casties, I., Turon, X. & Palacín, C. Some like it hot: Temperature and pH modulate larval development and settlement of the sea urchin Arbacia lixula. J. Exp. Mar. Biol. Ecol. 449, 304–311 (2013).
Google Scholar 
García, E., Clemente, S. & Hernández, J. C. Effects of natural current pH variability on the sea urchin Paracentrotus lividus larvae development and settlement. Mar. Environ. Res. 139, 11–18 (2018).PubMed 

Google Scholar 
Marshall, D. J. & Keough, M. J. Variation in the dispersal potential of non-feeding invertebrate larvae: The desperate larva hypothesis and larval size. Mar. Ecol. Prog. Ser. 255, 145–153 (2003).ADS 

Google Scholar 
Huggett, M. J., Williamson, J. E., de Nys, R., Kjelleberg, S. & Steinberg, P. D. Larval settlement of the common Australian sea urchin Heliocidaris erythrogramma in response to bacteria from the surface of coralline algae. Oecologia 149, 604–619 (2006).ADS 
PubMed 

Google Scholar 
Espinel-Velasco, N., Agüera, A. & Lamare, M. Sea urchin larvae show resilience to ocean acidification at the time of settlement and metamorphosis. Mar. Environ. Res. 159, 104977 (2020).CAS 
PubMed 

Google Scholar 
Lamare, M. & Barker, M. Settlement and recruitment of the New Zealand sea urchin Evechinus chloroticus. Mar. Ecol. Prog. Ser. 218, 153–166 (2001).ADS 

Google Scholar 
Martin, S. et al. Early development and molecular plasticity in the Mediterranean sea urchin Paracentrotus lividus exposed to CO2-driven acidification. J. Exp. Biol. 214, 1357–1368 (2011).CAS 
PubMed 

Google Scholar 
Vargas, C. A. et al. Species-specific responses to ocean acidification should account for local adaptation and adaptive plasticity. Nat. Ecol. Evol. 1, 0084 (2017).
Google Scholar 
Espinel-Velasco, N. et al. Effects of ocean acidification on the settlement and metamorphosis of marine invertebrate and fish larvae: a review. Mar. Ecol. Prog. Ser. 606, 237–257 (2018).ADS 

Google Scholar 
Briffa, M., de la Haye, K. & Munday, P. L. High CO2 and marine animal behaviour: potential mechanisms and ecological consequences. Mar. Pollut. Bull. 64, 1519–1528 (2012).CAS 
PubMed 

Google Scholar 
Gaylord, B. et al. Ocean acidification through the lens of ecological theory. Ecology 96, 3–15 (2015).PubMed 

Google Scholar  More

NEWS
01 April 2022

Tropical forests have big climate benefits beyond carbon storage

Study finds that trees cool the planet by one-third of a degree through biophysical mechanisms such as humidifying the air.

Freda Kreier

Freda Kreier

View author publications

You can also search for this author in PubMed
 Google Scholar

Twitter

Facebook

Email

Tropical forests create cloud cover that reflects sunlight and cools the air.Credit: Thomas Marent/Minden Pictures

Tropical forests have a crucial role in cooling Earth’s surface by extracting carbon dioxide from the air. But only two-thirds of their cooling power comes from their ability to suck in CO2 and store it, according to a study1. The other one-third comes from their ability to create clouds, humidify the air and release cooling chemicals.
How much can forests fight climate change?
This is a larger contribution than expected for these ‘biophysical effects’ says Bronson Griscom, a forest climate scientist at the non-profit environmental organization Conservation International, headquartered in Arlington, Virginia. “For a while now, we’ve assumed that carbon dioxide alone is telling us essentially all we need to know about forest–climate interactions,” he says. But this study confirms that tropical forests have other significant ways of plugging into the climate system, he says.The analysis, published in Frontiers in Forests and Global Change on 24 March1, could enable scientists to improve their climate models, while helping governments to devise better conservation and climate strategies.The findings underscore growing concerns about rampant deforestation across the tropics. Scientists warn that one-third of the world’s tropical forests have been mown down in the past few centuries, and another one-third has been degraded by logging and development. This, when combined with climate change, could transform vast swathes of forest into grasslands2.“This study gives us even more reasons why tropical deforestation is bad for the climate,” says Nancy Harris, forest-research director at the World Resources Institute in Washington DC.More than a carbon spongeForests are major players in the global carbon cycle because they soak up CO2 from the atmosphere as they grow. Tropical forests, in particular, store around one-quarter of all terrestrial carbon on the planet, making them “centrepieces for climate policy” in their home countries, Griscom says.
Tropical forests may be carbon sources, not sinks
“There’s clear evidence that the tropics are producing excellent climate benefits for the entire planet,” says Deborah Lawrence, an environmental scientist at the University of Virginia in Charlottesville and a co-author of the latest study. She and her colleagues analysed the cooling capacity of forests around the globe, in particular considering biophysical effects alongside carbon storage. Tropical forests, they found, can cool Earth by a whole 1 °C — and biophysical effects contribute significantly.Although scientists knew about these effects, they hadn’t understood to what extent the various factors counter global warming.Trees in the tropics provide shade, but they also act as giant humidifiers by pulling water from the ground and emitting it from their leaves, which helps to cool the surrounding area in a way similar to sweating, Griscom says.“If you go into a forest, it immediately is a considerably cooler environment,” he says.This transpiration, in turn, creates the right conditions for clouds, which like snow and ice in the Arctic, can reflect sunlight higher into the atmosphere and further cool the surroundings. Trees also release organic compounds — for example, pine-scented terpenes — that react with other chemicals in the atmosphere to sometimes create a net cooling effect.Locally coolTo quantify these effects, Lawrence and her colleagues compared how the various effects of forests around the world feed into the climate system, breaking down their contributions in bands of ten degrees of latitude. When they considered only the biophysical effects, the researchers found that the world’s forests collectively cool the surface of the planet by around 0.5 °C.
When will the Amazon hit a tipping point?
Tropical forests are responsible for most of that cooling. But this band of trees across Latin America, Central Africa and southeast Asia is under increasing pressure from climate change and deforestation. Both of these human-caused impacts can lead rainforests to dry out, says Christopher Boulton, a geographer at the University of Exeter, UK. Last month, he and his colleagues published a review2 of nearly 30 years’ worth of satellite images of the Amazon, the largest rainforest in the world. By measuring the biomass of the vegetation in the images, the team discovered that three-quarters of the Amazon is losing resilience — the ability to recover from an extreme weather event such as a drought.Threats to tropical rainforests are dangerous not only for the global climate, but also for communities that neighbour the forests, Lawrence says. She and her colleagues found that the cooling caused by biophysical effects was especially significant locally. Having a rainforest nearby can help to protect an area’s agriculture and cities from heatwaves, Lawrence says. “Every tenth of a degree matters in limiting extreme weather. And where you have forests, the extremes are minimized.”Governments across the tropics have struggled to conserve their forests despite more than two decades of global campaigns to halt deforestation, promote sustainable development and protect the climate. Lawrence says that her team’s findings make it clear that protecting forests is a matter of self-interest, and has immediate benefits for local communities.

doi: https://doi.org/10.1038/d41586-022-00934-6

ReferencesLawrence, D., Coe, M., Walker, W., Verchot, L. & Vandecar, K. Front. For. Glob. Change https://doi.org/10.3389/ffgc.2022.756115 (2022).Article 

Google Scholar 
Boulton, C. A., Lenton, T. M. & Boers, N. Nature Clim. Change 12, 271–278 (2022).Article 

Google Scholar 
Download references

Related Articles

How much can forests fight climate change?

When will the Amazon hit a tipping point?

Tropical forests may be carbon sources, not sinks

Illegal mining in the Amazon hits record high amid Indigenous protests

Subjects

Climate sciences

Climate change

Conservation biology

Latest on:

Climate sciences

Funding battles stymie ambitious plan to protect global biodiversity
News 31 MAR 22

Trends in Europe storm surge extremes match the rate of sea-level rise
Article 30 MAR 22

From the archive: fishy business in 1972 and 1922
News & Views 29 MAR 22

Climate change

Funding battles stymie ambitious plan to protect global biodiversity
News 31 MAR 22

Trends in Europe storm surge extremes match the rate of sea-level rise
Article 30 MAR 22

The race to upcycle CO2 into fuels, concrete and more
News Feature 29 MAR 22

Jobs

Postdoctoral Fellow in Electrochemical CO2 reduction

The University of British Columbia (UBC)
Kelowna, Canada

Staff Member in Project Coordination (for 19 h/week)

Jülich Research Centre (FZJ)
Erlangen-Nürnberg, Germany

Student assistant IT administration (m/f/d)

Alfred Wegener Institute – Helmholtz Centre for Polar and Marine Research (AWI)
Potsdam, Germany

Physicist (postdoctoral researcher) (all genders) in the field of Theoretical Astrophysics

Helmholtz Centre for Heavy Ion Research GmbH (GSI)
Darmstadt, Germany More

Flores, C. O. et al. Proc. Natl Acad. Sci. USA 108, 288–297 (2011).Article 

Google Scholar 
Carlson, M. C. G. et al. Nat. Microbiol. https://doi.org/10.1038/s41564-022-01088-x (2022).Article 

Google Scholar 
Flombaum, P. et al. Proc. Natl Acad. Sci. USA 110, 9824–9829 (2013).CAS 
Article 

Google Scholar 
Coleman, M. L. & Chisholm, S. W. Trends Microbiol. 15, 398–407 (2007).CAS 
Article 

Google Scholar 
Johnson, Z. I. et al. Science 311, 1737–1740 (2006).CAS 
Article 

Google Scholar 
Martiny, A. C. et al. PLoS ONE 11, e0168291 (2016).Article 

Google Scholar 
Wilhelm, S. W. & Suttle, C. A. Bioscience 49, 781–788 (1999).Article 

Google Scholar 
Follett, C. L. et al. Proc. Natl Acad. Sci. USA 119, e2110993118 (2022).CAS 
Article 

Google Scholar 
Mojica, K. D. A. et al. ISME J. 10, 500–513 (2016).CAS 
Article 

Google Scholar 
Mruwat, N. et al. ISME J. 15, 41–54 (2021).CAS 
Article 

Google Scholar  More

To explore how environmental gradients shape the distribution of cyanophages and picocyanobacteria, we conducted high-resolution surveys in surface waters along five oceanic transects on three cruises covering thousands of kilometres in the North Pacific Ocean in the spring or early summer of 2015, 2016 and 2017 (Fig. 1a–c). These cruises, two of which were out-and-back, passed through distinct regimes from warm, saline and nutrient-poor waters of the North Pacific Subtropical Gyre to cooler, less saline and nutrient-rich waters of higher latitudes influenced by the subpolar gyre (Fig. 1d–i)27. The shift between the two gyres was marked by abrupt changes in trophic indicators such as particulate carbon concentrations (Fig. 1g) and a chlorophyll front (defined as the 0.2 mg m−3 chlorophyll contour28; Fig. 1a–c). As such, the inter-gyre transition zone, defined by salinity and temperature thresholds29 (Fig. 1d), was distinct from both the subtropical and subpolar gyre ecosystems28.Fig. 1: Gradients in environmental conditions across the North Pacific gyres.a–c, Transects of three cruises overlaid on monthly averaged satellite-derived sea-surface chlorophyll in March 2015 (a), April 2016 (b) and June 2017 (c). d, Temperature–salinity diagram showing the boundaries of the subtropical and subpolar gyres (black dashed lines) based on the salinity thresholds reported by Roden29. e–i, Temperature (e), salinity (f) as well as the levels of particulate carbon (g), phosphate (h) and nitrate + nitrite (i) as a function of latitude. The coloured dashed lines show the position of the 0.2 mg m−3 chlorophyll contour. For environmental variables plotted against temperature, see Supplementary Fig. 3.Full size imageUnexpected Prochlorococcus declineProchlorococcus concentrations in the oligotrophic waters of the subtropical gyre were 1.5–3.0 × 105 cells ml−1, comprising an average of approximately 29% of the total bacteria (Extended Data Fig. 1) and numerically dominating the phytoplankton community in all three cruises (Extended Data Fig. 2). Prochlorococcus abundance remained high in the southern region of the transition zone in 2015 and 2016, decreasing precipitously to less than 2,000 cells ml−1 north of the chlorophyll front, generally constituting 80% of cyanophages measured, with the remainder consisting of T7-like clade A and TIM5-like cyanophages (Fig. 3 and Extended Data Fig. 4). Cyanophage abundances correlated positively with total picocyanobacteria in the subtropical gyre (Pearson’s coefficient of multiple correlation (r) = 0.54, P = 0.02, n = 26; Fig. 2d), suggesting that cyanophages were limited by the availability of susceptible hosts in this region and were not regulating picocyanobacterial populations. On average, less than 1% of the cyanobacterial populations were infected (Fig. 4), with higher infection rates by T4-like cyanophages than T7-like cyanophages (Extended Data Figs. 5 and 6). These instantaneous measurements of infection were used to estimate the daily rates of mortality39 (Methods and Supplementary Discussion), which suggests that 0.5–6% of picocyanobacterial populations were lysed by viruses each day (Extended Data Fig. 7). This implicates other factors, such as grazing45, as the major causes of cyanobacterial mortality in the North Pacific Subtropical Gyre.Fig. 3: Cyanophage community composition across the North Pacific gyres.a–c, Cyanophage abundance for the March 2015 (a), April 2016 (b) and June 2017 (c) transects. Insets: T7-like clade A and TIM5-like cyanophage abundances on an expanded scale (similar to the main images, the units for the vertical axes are ×105 viruses ml−1). The grey shaded regions show the position of the virus hotspot. See Extended Data Fig. 4 for the confidence intervals and out-and-back reproducibility and Supplementary Fig. 4 for cyanophage lineages plotted against latitude.Full size imageFig. 4: Viral infection patterns of picocyanobacteria in the North Pacific Ocean.a–f, Viral infection levels (black) of Prochlorococcus (a,c,e) and Synechococcus (b,d,f) plotted against temperature for the March 2015 (a,b), April 2016 (c,d) and June 2017 (e,f) transects. Insets: infection levels on an expanded scale. The solid lines show infection (red), Prochlorococcus (green) and Synechococcus (pink) averaged and plotted for every 0.5 °C. The dashed lines and shaded regions show the position of the chlorophyll front and the virus hotspot, respectively. For plots by latitude and the upper and lower bounds of infection, see Extended Data Figs. 5 and 6.Full size imageWithin the transition zone we observed a steep latitudinal increase in the abundance of cyanophages for every transect, which we define as a cyanophage hotspot (Fig. 2c and Extended Data Figs. 2 and 4). The cyanophage abundances in this hotspot were between three- and tenfold greater than in the subtropical gyre (Fig. 2c). Notably, cyanophages were approximately 25% more abundant (an increase of approximately 5 × 105 viruses ml−1) in the hotspot on the 2017 cruise relative to the other two cruises, reaching a maximum of 2 × 106 viruses ml−1. The hotspot peaked at temperatures of 15–16 °C on all transects, regardless of the geographical location, season or the exact pattern of the Prochlorococcus and Synechococcus distributions (Fig. 2c). Notably, the numbers of T7-like clade B cyanophages increased sharply in the transition zone to become the most abundant lineage, whereas T4-like cyanophages increased more modestly (Fig. 3 and Extended Data Fig. 4). The change in the cyanophage community structure was particularly pronounced in June 2017, when T7-like cyanophages were up to 2.3-fold more abundant than T4-like cyanophages (Fig. 3c). The switch in the relative abundance of T4-like and T7-like clade B cyanophages was diagnostic of the cyanophage hotspot compared with patterns in the subtropical and subpolar gyres.To begin assessing whether cyanophages negatively affected cyanobacterial populations in the hotspot, we tested the relationship between the abundance of cyanophages and total cyanobacteria. This showed a significant negative correlation between cyanophage and cyanobacterial abundances across all three cruises (Pearson’s r = −0.56, two-sided P = 0.0005, n = 34). This relationship was particularly distinct in 2017, when cyanobacteria were at their overall lowest abundances and cyanophages at their highest (Pearson’s r = −0.65, two-sided P = 0.004, n = 18). This suggests that viruses are one of the key regulators of picocyanobacteria in the region of the hotspot. However, no significant correlation was found across all regimes and all years (Pearson’s r = −0.008, two-sided P = 0.9, n = 87; Fig. 2d), indicating that factors other than viruses are likely to be more important in regulating the abundances of cyanobacteria in other regimes.Our single-cell infection measurements allowed us to directly evaluate active viral infection and its impact on picocyanobacteria in the transition zone. Viral infection spiked in this region each year with infection levels that were an average of two- to ninefold higher than those in the subtropical gyre (Fig. 4 and Extended Data Figs. 5,6 and 8). Infection peaked within the temperature range of 12–18 °C and was associated with a concomitant dip in Prochlorococcus abundances in all three cruises (Fig. 4 and Extended Data Fig. 5). These findings provide independent support for the strong negative correlation between cell and virus abundances (Fig. 2d) being the result of virus-induced mortality.Lineage-specific infection was also distinct in the transition zone relative to the subtropical gyre. Infection by T7-like clade B cyanophages generally increased to reach (2015 and 2016) or exceed (2017) those of T4-like cyanophages (Extended Data Figs. 5 and 6). In addition, the ratio of the abundances of T7-like clade B cyanophages to the number of cells they infected was 2.6-fold greater in the hotspot than the subtropics, whereas this ratio was similar in both regions for T4-like cyanophages. Together, these results indicate that, within the hotspot, the T4-like cyanophages displayed increased levels of infection, whereas the T7-like cyanophages displayed both increased levels of infection and produced more viruses per infection, suggesting that T7-like clade B cyanophages are better adapted to conditions in the transition zone (see below).Of the three cruises, the highest levels of viral infection were observed in June 2017, with up to 9.5% and 8.9% of Prochlorococcus and Synechococcus infected, respectively (Fig. 4e,f). This dramatic increase in infection mirrored the massive decline in Prochlorococcus abundances (Fig. 4e and Extended Data Fig. 5i). We estimate that viruses killed 10–30% of Prochlorococcus and Synechococcus cells daily at these high instantaneous levels of infection (Extended Data Fig. 6) based on the expected number of infection cycles cyanophages were able to complete at the light and temperature conditions in the transition zone (Methods and Supplementary Discussion). Given that Prochlorococcus is estimated to double every 2.8 ± 0.8 d at the low temperatures in this region12, we estimate that 21–51% of the population was infected and killed in the interval before cell division. Synechococcus is expected to have faster growth rates at these temperatures, doubling every 1.1 ± 0.2 d (refs. 12,46). Thus, we estimate that less of the Synechococcus population (9–31%) was killed before division.Under quasi-steady state conditions, abiotic controls on the growth rate of Prochlorococcus are balanced by mortality due to viral lysis, grazing and other mortality agents39,45,47. Based on the high levels of virus-mediated mortality, the parallel pattern between Prochlorococcus’ death and viral infection, and the negative correlation between cyanophage and picocyanobacterial abundances in the transition zone, we propose that enhanced viral infection in 2017 disrupted this balance, leading to the unexpected decline in Prochlorococcus populations. Grazing and other mortality agents not investigated here could also have contributed to additional mortality beyond the steady state, resulting in further losses of Prochlorococcus. In contrast to Prochlorococcus, Synechococcus maintained large populations despite high levels of infection (Fig. 4f), presumably due to their faster growth rates enabling them to maintain a positive net growth despite enhanced mortality. These findings suggest that virus-mediated mortality in 2017 was an important factor in limiting the geographic range of Prochlorococcus that resulted in a massive loss of habitat of approximately 550 km.Cyanophage abundances and infection levels dropped sharply in the higher-latitude waters north of the hotspot (Figs. 2c, 4 and Extended Data Figs. 1d,h and 2). The abundances of both T7-like clade B and T4-like cyanophages declined precipitously, yet T4-like cyanophages were the dominant cyanophage lineage (Fig. 3). T7-like clade A cyanophages generally increased locally at the northern border of the hotspot and became the dominant T7-like lineage in two samples between 38 and 39.2° N in 2017 (Fig. 3c and Extended Data Fig. 4). In contrast to all other cyanophages, the abundances of TIM5-like cyanophages increased in waters north of the hotspot (Fig. 3 and Extended Data Fig. 4d,i,m) but remained a minor component of the cyanophage community. No relationship was found between cyanophage and cyanobacterial abundances (Fig. 2d), and less than 1.5% of picocyanobacteria were infected by all cyanophage lineages in these waters (Fig. 4).The cyanophage hotspot in the transition zone is a ridge of high virus activity that separates the subtropical and subpolar gyres. The reproducibility of our observations, which were separated by days to weeks within each cruise (2016 and 2017) and by years among the three cruises (Extended Data Fig. 4), indicates that this virus hotspot is a recurrent feature at the boundary of these two major gyres in the North Pacific Ocean. This suggests that the hotspot forms due to the distinctive environment of the inter-gyre transition zone creating conditions that enhance infection of picocyanobacteria and proliferation of cyanophages. Prochlorococcus in the transition zone may be prone to stress due to being close to the limits of their temperature growth range5,6, which has the potential to increase susceptibility to viral infection. Alternatively, there may be temperature-dependent trade-offs between virus decay and production that lead to replication optima within a narrow temperature range48. Cyanophage infectivity has been observed to decay more slowly at colder temperatures49, which may allow for the accumulation of infective viruses, leading to increased infection. In addition, cyanophage infections may be more productive due to enhanced nutrient supply in the transition zone27 (Fig. 1h,i) relative to the subtropics, given that the cyanophages replicate in hosts with presumably greater intracellular nutrient quota and obtain more extracellular nutrients, both of which may increase progeny production9,10. The environmental factors influencing the production and removal of viruses probably vary in intensity at different times, leading to variability in cyanophage abundance and infection levels. Thus, the putative cyanophage replication optimum in the hotspot may reflect the combined effects of temperature and nutrient conditions that are intrinsically linked to the oceanographic forces that shape the transition zone itself.Changes in the cyanophage community structure over environmental gradients are likely to reflect differences in host range, infection properties and genomic potential to remodel host metabolism9. Our data, together with previous measurements in the North Pacific Subtropical Gyre38,39, indicate that the T4-like cyanophages are the lineage best adapted to the low-nutrient waters of the subtropics (Fig. 2d–f). As these waters are inhabited by hundreds of genomically diverse subpopulations of Prochlorococcus50, the broad host range of many T4-like cyanophages18,19,22,51 may be advantageous for finding a suitable host. T4-like cyanophages also have a large and diverse repertoire of host-derived genes21,51—such as nutrient acquisition, photosynthesis and carbon-metabolism genes—that augment host metabolism52 and may increase fitness in nutrient-poor conditions in the subtropics51. In contrast, T7-like clade A and B cyanophages seem to be better adapted to conditions in the transition zone (Fig. 3). T7-like cyanophages have narrow host ranges19,22,40, with smaller genomes and fewer genes to manipulate the host metabolism23, which may allow them to replicate and produce more progeny in regions with elevated nutrient concentrations relative to subtropical conditions. The maximal abundances of TIM5-like cyanophages were found in the most productive waters at the northern end of the transects where the cyanobacterial abundances were lowest and Synechococcus was the dominant picocyanobacterium. This may be partially due to the narrow host range of TIM5-like cyanophages and their specificity for Synechococcus40,44. Our findings of reproducible lineage-specific responses to changing ocean regimes indicate that cyanophage lineages occupy distinct ecological niches.Temperature and nutrient changes occurring in the transition zone are expected to result in shifts in picocyanobacterial diversity at the sub-genus level (Supplementary Discussion), which we speculate may affect community susceptibility to viral infection. One mechanism for this may be that the picocyanobacteria that thrive in the transition zone are intrinsically more susceptible to viral infection. Another scenario may be related to trade-offs associated with the evolution of resistance to viral infection. The horizontal advection of nutrient-rich waters to the transition zone28 may select for rapidly growing cells adapted for efficient resource utilization. Viral resistance in picocyanobacteria often incurs the cost of reduced growth rates53,54. Thus, competition for nutrients in this region may favour cells with faster growth rates but increased susceptibility to viral infection. Thus, it is probable that the cyanophage distributions do not always follow the cyanobacterial patterns (Extended Data Fig. 2) because of complex interactions between lineage-specific cyanophage traits, host community structure and environmental variables, which may vary seasonally or annually as a result of interannual variability in environmental conditions (see below).Despite consistent features in cyanophage distributions across the North Pacific Ocean, cyanophage infection was higher (Fig. 4 and Extended Data Fig. 7), whereas Prochlorococcus abundances were consistently lower (Fig. 2a), across the June 2017 transects relative to the March 2015 and April 2016 transects. Seasonality and/or climate variability could explain this interannual variability, although the data currently available to assess this are sparse. Viral infection of picocyanobacteria in the subtropical gyre increased from early spring to summer, suggesting a potential seasonal pattern that may extend across the transect (Extended Data Fig. 9a). In addition, the June 2017 transect occurred during a neutral-to-negative El Niño phase with lower sea-surface temperatures relative to the 2015 and 2016 transects, which were in years of a record marine heatwave, followed by a strong El Niño55 (Extended Data Fig. 9b). In 2015 and 2016, the Prochlorococcus abundances were found to be higher than usual in the North Pacific Ocean in this (Fig. 2a) and other studies56,57. Irrespective of the underlying drivers for the observed interannual variability, we speculate that an ecosystem tipping point was reached in the hotspot under the prevailing conditions in June 2017, aided by the higher cyanophage abundances yet smaller Prochlorococcus population sizes. In this scenario, picocyanobacterial populations were subjected to high infection levels that resulted in an accumulation of cyanophages, initiating a stronger than usual positive-feedback loop between infection and virus production, and precipitating the unexpected Prochlorococcus decline. Continued observations in the North Pacific Ocean are needed to evaluate the potential link between seasonality and/or large-scale climate forcing as ultimate drivers affecting virus–host interactions.Predicting basin-scale virus dynamicsMeasurements of cyanobacterial and cyanophage abundances rely on discrete sample collection from shipboard oceanographic expeditions, which limits the geographical and seasonal extent of available data. Therefore, we developed a multiple regression model based on high-resolution satellite data of temperature and chlorophyll to predict cyanophage abundances, a key proxy of cyanobacterial infection (Pearson’s r = 0.61, two-sided P = 1.7 × 10−8, degrees of freedom = 68, n = 70). We used the model to estimate the geographical extent of the virus hotspot. The model accurately predicted the location of the hotspot and cyanophage abundances along a fourth transect in April 2019 (Supplementary Table 1), with the majority of observations falling within the 95% confidence intervals of the model predictions (Fig. 5a–c). Application of the model to the larger region predicted that the virus hotspot formed a boundary extending across the North Pacific Ocean, with lower cyanophage abundances on both sides (Fig. 5d,e and Supplementary Fig. 1). This boundary had the hallmarks of the hotspot with a core that was dominated by T7-like cyanophages and the flanking gyre regions dominated by T4-like cyanophages. Thus, this feature may be more appropriately termed a ‘hot-zone’ due to its substantial projected aerial extent. Assuming the infection levels observed in the hotspot in June 2017 were similar throughout the hot-zone, the potential habitat loss for Prochlorococcus would be about 3.2 × 106 km2, approximately half of the cumulative area loss of the Amazonian rainforest to date58.Fig. 5: Prediction of cyanophage abundances.a–c, Model-based predictions of cyanophage abundances corresponding to the empirically measured total (a), T4-like (b) and T7-like clade B (c) cyanophage abundances along a transect in the North Pacific in April 2019. The shaded regions show the 95% confidence interval for the model predictions. d,e, Predicted total cyanophages (d) and the ratio of T4-like/T7-like clade B cyanophages (e) in June 2017 in the North Pacific Ocean. The black lines indicate the cruise track. The grey areas represent regions with no values due to cloud cover or that were beyond the limits of the predictive model. The hotspot peak corresponds to yellow regions in d and red regions in e.Full size imageVirus hotspot biogeochemistryWith the ability to predict biogeographic patterns of cyanophages, we evaluated the potential biogeochemical implications of virus-mediated picocyanobacterial lysis and release of organic material in sustaining the bacterial community6,7,8,9. The aerial extent of the hot-zone (approximately 4 × 106 km2) is only 14% of the size of the subtropical gyre (2.9 × 107 km2), and yet the total virus-mediated organic matter released from picocyanobacteria in the hot-zone in June 2017 was estimated to be on par with that for the entire North Pacific Subtropical Gyre (Methods and Supplementary Discussion). We estimate that viral lysate released from picocyanobacteria in the subtropical gyre could sustain 4.4 ± 0.8% of the calculated bacterial carbon demand there (Extended Data Fig. 10). In contrast, viral lysate released in the transition zone could sustain an average of 21 ± 12% of the bacterial carbon demand, reaching 33% in some regions (Extended Data Fig. 10), assuming that the bacterial assimilation and growth efficiencies were similar between the subtropical gyre and the hotspot. Thus, local generation of cyanobacterial viral lysate in the transition zone is likely to be an important source of carbon for the heterotrophic bacterial community that can rapidly utilize large molecular weight dissolved organic matter59 and may have contributed to the increase in their abundances south of the chlorophyll front in 2017 (Extended Data Fig. 1a,e). More

The generation of this map of the world’s mountains consisted of five steps (Fig. 1): (i) the identification and hierarchisation of named mountain ranges and the recording of range-specific information; (ii) the manual digitization of the ranges’ general shape; (iii) the definition of mountainous terrain (and the inventory’s outer borders) using a DEM-based algorithm; (iv) the automatic refinement of the digitized and named ranges’ inner borders; and (v) the preparation of the final layers. The resulting products consist of a refined mountain definition (GMBA Definition v2.0), two versions of the inventory (GMBA Inventory v2.0_standard & GMBA Inventory v2.0_broad), and a set of tools to work with the inventories.Step i: Identification and hierarchisation of mountain rangesIn a first step, we identified mountain ranges worldwide. To do so we adopted the mountain ranges identified in the GMBA Inventory v1.410,14 and searched existing resources in any languages for other named ranges not yet included. The ranges added could either be adjacent to, included in (child range or subrange) or including (parent range or mountain system) mountain ranges of the GMBA Inventory v1.4. The resources used for our searches included world atlases (e.g. The Times Comprehensive Atlas of the World19, Knaurs grosser Weltatlas20, Pergamon World Atlas21); topographic maps (e.g. http://legacy.lib.utexas.edu/maps/imw/, http://legacy.lib.utexas.edu/maps/onc/, https://maps.lib.utexas.edu/maps/tpc/, www.topomap.co.nz, https://norgeskart.no, www.ign.es/iberpix/visor/); encyclopaedias (www.wikipedia.org; www.britannica.com); online gazetteers and reference sites (e.g. www.wikidata.org, www.geonames.org (GeoNames), www.mindat.org); mountain classification systems (e.g. the International Standardized Mountain Subdivision of the Alps or SOIUSA for the Alps22, Alpenvereinseinteilung der Ostalpen23, Classification of the Himalaya24, www.peakbagger.com/rangindx.aspx (PEMRACS), www.carpathian-research-network.eu/ogulist, http://www.sopsr.sk/symfony-bioregio/lkpcarporog, www.dinarskogorje.com, https://bivouac.com/, https://climbnz.org.nz/); and national or regional landscape, geomorphological, or physiographic maps and publications4,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42. The full list of the consulted sources and references is available on GitHub at https://www.github.com/GMBA-biodiversity/Inventory (GMBA Mountain Inventory v2.0 References.pdf).All identified mountain ranges were recorded in a Microsoft Access relational database (“Mountain database”, see below) and given a name, a unique 5-digit identifier (GMBA_V2_ID), and the corresponding Wikidata unique resource identifier (URI), when available. This URI gives access to a range’s name as well as to its Wikipedia page URL in all available languages and lists other identifiers for given mountain ranges in a variety of other repositories such as GeoNames or PEMRACS. The primary mountain range names were based on the resources used for range identification and were preferably recorded in English. Names used nationally, locally, as well as/or by indigenous people and local communities were extracted from Wikidata and recorded in a separate attribute field.In the process of cataloguing, we attributed a parent range to each of the mapped mountain ranges. Information about parent ranges is included in PEMRACS, often also in Wikidata as a property that can be extracted though a SPARQL query, in the corresponding Wikipedia pages description, and in regional hierarchical mountain classifications that exist for the European Alps (SOIUSA), the Carpathians, and the Dinaric Alps. When no such information was available, we relied on other sources of information that we found either using a general web search (leading to specific papers, reports, or web pages on mountain ranges) or by consulting (online) topographical maps and atlases at different scales. The information about parent ranges was used to construct a hierarchy of up to 10 levels using a recursive SQL query (see Step v). The result of this step was a relational database with a hierarchy of mountain systems and (sub-) ranges (Fig. 1, “Mountain database”).Step ii: Digitization of the mountain rangesIn a second step, we digitized all identified ‘childless’ mountain ranges (i.e. smallest mapping units, called ‘Basic’ as opposed to ‘Aggregated’ in the database) in one vector GIS layer. To do so, we used the Google Maps Terrain layers (Google, n.d.) as background and the WHYMAP named rivers layer42 as spatial reference since descriptions of mountain range areal extension is often given with reference to major rivers. The digitization, which was done in QGIS43 using the WGS 84 / Pseudo-Mercator (EPSG 3857) coordinate reference system, consisted in the drawing of shapes (polygons) that roughly followed the core area of each mountain range. In general, the approximate shape and extent of the mountain ranges we digitized could be distinguished based on the terrain structure as represented by the shaded relief background that corresponded to the placement and orientation of the range’s name label on a topographical map, atlas or other resource. As the exact placement and orientation of mountain range labels in each specific source can be influenced by cartographic considerations (e.g. avoiding overlaps with other features), the final approximation of the mountain range was obtained by consulting a variety of sources for each mountain range. Occasionally, the mountain terrain’s geomorphological characteristics strongly hampered the accuracy of our visual identification of mountain subranges within larger systems. This was particularly the case in old, eroded massifs such as the Brazilian Highlands or the highlands of Madagascar, where individual mountain ranges are not separated by deep well-defined valleys and have a very complex topography. In these cases, we referred to available topographical descriptions of range extent and to the river layer (see above). Other complex regions included Borneo and the Angolan Highlands, whereas subranges in mountain systems such as the European Alps, the Himalayas, and the North American Cordillera were comparatively easy to map. Moreover, the density of currently available mountain toponymical information varied quite strongly between regions. Accordingly, regional variation in the size of the smallest mountain range map units can be considerable. The result of this step was a (manually) digitized vector layer of named mountain ranges shapes (Fig. 1, “Manual mountain shapes”).Step iii: Definition of mountainous terrainIn a third step, we defined mountainous terrain (GMBA Definition v2.0). To distinguish mountainous from non-mountainous terrain, we developed a simple algorithm which we implemented in ArcMap 10.7.144. This algorithm is based on ruggedness (defined as highest minus lowest elevation in meter) within eight circular neighbourhood analysis windows (NAWs) of different sizes (from 1 pixel (≈ 250 m) to 20 (≈ 5 km) around each point, Fig. 2c) combined with empirically derived thresholds for each NAW (Fig. 2). The decision to use multiple NAW sizes was made because calculating ruggedness based on only a small or a large NAW comes at the risk of identifying the many local irregularities typically occurring in flat or rolling terrain as mountainous or of including extensive flat ‘skirts’ through the smoothing and generalization of large NAWs3. Accordingly, our approach ensures that any point in the landscape classified as mountainous showed some level of ruggedness not only at one but across scales. This also resulted in a smooth and homogeneous delineation of mountainous terrain, very suitable for our mapping purpose.Fig. 2Elevation range thresholds for the eight neighbourhood analysis windows (NAW) and their contribution to calculations of the GMBA Definition v2.0. (a) distribution of elevation range values (ruggedness) for NAWs (numbered I to VIII) in mountain regions as defined by the geometric intersection of K1, K2 and K3. (b): plot of the minimum elevation range versus the area of the NAW (n = 920). (c) NAWs and their corresponding threshold values. (d) percent overlap between GMBA Definition v2.0 (intersection of eight NAW-threshold pairs) and area defined by each individual NAW-threshold pair. (e) percent eliminated by each NAW-threshold pair (I to VIII) from the mountain area defined by the other 7 NAW-threshold combinations. Highlighted bars in the two graphs represent the combination of three NAW-threshold pairs that results in the highest overlap with the GMBA Definition v2.0.Full size imageWe used the median value of the 7.5 arc second GMTED2010 DEM45 as our source map. To reduce the latitudinal distortion of the raster, and thus the shape and area of the NAWs, we divided the global DEM into three raster layers corresponding to three latitudinal zones (84° N to 30° N, 30° N to 30° S and 30° S to 56° S) excluding ice-covered Antarctica and projected the two high latitude zones to Lambert Azimuthal Equal Area and the equatorial zone to WGS 1984 Cylindrical Equal Area. We used these reprojected DEM layers to produce eight ruggedness layers, each using one of the eight NAWs.To determine the threshold values of our algorithm, we selected 1000 random points within the area defined by the geometric intersection (Fig. 1b) of the three commonly applied mountain definitions, i.e. the definitions by UNEP-WCMC46, GMBA15, and USGS3. These layers (referred to as K1, K2, and K3, respectively by Sayre and co-authors12) were obtained from the Global Mountain Explorer47. We eliminated 80 clearly misclassified points (i.e., points that fell within lakes, oceans, or clearly flat areas according to the shaded relief map we used as a background) and used the remaining 920 to sample the eight ruggedness layers. For each of the 8 layers, we retained the lowest of the 920 ruggedness values as the threshold for the layer’s specific NAW (Fig. 2c). The eight threshold values were then used to reclassify each of the eight layers by attributing the value 1 to all cells with a ruggedness value higher than or equal to the corresponding threshold and the value 0 to all other cells. Finally, we performed a geometric intersection (see Fig. 1b) of the eight reclassified layers to derive the new mountain definition.After these calculations, we reprojected the three raster layers to WGS84 and combined them through mosaic to new raster. To eliminate isolated cells and jagged borders, we then generalized the resulting raster map by passing a majority filter (3 × 3 pixels, majority threshold) three times. This layer corresponds to the GMBA Definition v2.0.The resulting mountain definition (GMBA Definition v2.0) distinguishes itself from previous ones because of the empirically derived thresholds method used to develop it and the use of eight NAWs. In line with the previous GMBA definition, it relies entirely on the ruggedness values within NAWs. The GMBA Definition v2.0 was used to determine the outer delineation of this inventory’s mountainous terrain. As expected, it includes neither the wide ‘skirts’ of flat or undulating land around mountain ranges nor the topographical irregularities that are both typically included when other approaches are applied. It also successfully excludes extensive areas of rolling non-mountainous terrain such as the 52,000 km2 Badain Jaran Desert sand dunes in China. However, this mountain definition is conservative and only includes the highest, most rugged cores of low mountain systems, as for example in the Central Uplands of Germany, and therefore excludes some lower hill areas still considered by some as mountains.As a further step towards generalization, we considered that small ( More

Volk T, Hoffert MI. Ocean carbon pumps: Analysis of relative strengths and efficiencies in ocean-driven atmospheric CO2 changes. In: Sundquist ET, Broecker WS. (eds). The carbon cycle and atmospheric CO2: Natural variations archean to present. American Geophysical Union, Geophysical Monograph, Washington, DC: 1985. p. 32:99–110.Scharek R, Tupas LM, Karl DM. Diatom fluxes to the deep sea in the oligotrophic North Pacific gyre at Station ALOHA. Mar Ecol-Prog Ser. 1999;182:55–67.
Google Scholar 
Simon M, Grossart H, Schweitzer B, Ploug H. Microbial ecology of organic aggregates in aquatic ecosystems. Aquat Micro Ecol. 2002;28:175–211.
Google Scholar 
Siegenthaler U, Sarmiento JL. Atmospheric carbon dioxide and the ocean. Nature. 1993;365:119–25.CAS 

Google Scholar 
Ducklow H, Steinberg DK. Upper ocean carbon export and the biological pump. Oceanography. 2001;14:50–58.
Google Scholar 
Jiao N, Herndl GJ, Hansell DA, Benner R, Kattner G, Wilhelm SW, et al. Microbial production of recalcitrant dissolved organic matter: long-term carbon storage in the global ocean. Nat Rev Microbiol. 2010;8:593–9.CAS 
PubMed 

Google Scholar 
DeLong EF, Franks DG, Alldredge AL. Phylogenetic diversity of aggregate-attached vs. free-living marine bacterial assemblages. Limnol Oceanogr. 1993;38:924–34.
Google Scholar 
Allen AE, Allen LZ, McCrow JP. Lineage specific gene family enrichment at the microscale in marine systems. Curr Opin Microbiol. 2013;16:605–17.CAS 
PubMed 

Google Scholar 
D’Ambrosio L, Ziervogel K, MacGregor B, Teske A, Arnosti C. Composition and enzymatic function of particle-associated and free-living bacteria: a coastal/offshore comparison. ISME J. 2014;8:2167–79.PubMed 
PubMed Central 

Google Scholar 
Martin JH, Knauer GA, Karl DM, Broenkow WW. VERTEX: carbon cycling in the northeast Pacific. Deep-Sea Res Part I-Oceanogr Res Pap. 1987;34:267–85.CAS 

Google Scholar 
Buesseler KO. The decoupling of production and particulate export in the surface ocean. Glob Biogeochem Cycle. 1998;12:297–310.CAS 

Google Scholar 
Schlitzer R. Applying the adjoint method for biogeochemical modeling: export of particulate organic matter in the world ocean. In: Kasibhata P, editor. Inverse Methods in Global biogeochemical Cycles. Washington, DC: American Geophysical Union; 2000. p. 114:107–24.Steinberg DK, Van Mooy BAS, Buesseler KO, Boyd PW, Kobari T, Karl DM. Bacterial vs. zooplankton control of sinking particle flux in the ocean’s twilight zone. Limnol Oceanogr. 2008;53:1327–38.
Google Scholar 
Cho BC, Azam F. Major role of bacteria in biogeochemical fluxes in the ocean’s interior. Nature. 1988;332:441–3.CAS 

Google Scholar 
Herndl GJ, Reinthaler T. Microbial control of the dark end of the biological pump. Nat Geosci. 2013;6:718–24.CAS 
PubMed 
PubMed Central 

Google Scholar 
Bergh Ø, Borsheim KY, Bratbak G, Heldal M. High abundance of viruses found in aquatic environments. Nature. 1989;340:467–8.CAS 
PubMed 

Google Scholar 
Suttle CA. Viruses in the sea. Nature. 2005;437:356–61.CAS 
PubMed 

Google Scholar 
Zhang R, Wei W, Cai L. The fate and biogeochemical cycling of viral elements. Nat Rev Microbiol. 2014;12:850–1.CAS 
PubMed 

Google Scholar 
Middelboe M, Lyck PG. Regeneration of dissolved organic matter by viral lysis in marine microbial communities. Aquat Micro Ecol. 2002;27:187–94.
Google Scholar 
Weinbauer MG, Brettar I, Hofle MG. Lysogeny and virus-induced mortality of bacterioplankton in surface, deep, and anoxic marine waters. Limnol Oceanogr. 2003;48:1457–65.
Google Scholar 
Fuhrman JA. Marine viruses and their biogeochemical and ecological effects. Nature. 1999;399:541–8.CAS 
PubMed 

Google Scholar 
Jover LF, Effler TC, Buchan A, Wilhelm SW, Weitz JS. The elemental composition of virus particles: implications for marine biogeochemical cycles. Nat Rev Microbiol. 2014;12:519–28.CAS 
PubMed 

Google Scholar 
Bongiorni L, Magagnini M, Armeni M, Noble R, Danovaro R. Viral production, decay rates, and life strategies along a trophic gradient in the North Adriatic Sea. Appl Environ Microbiol. 2005;71:6644–50.CAS 
PubMed 
PubMed Central 

Google Scholar 
Weinbauer MG, Bettarel Y, Cattaneo R, Luef B, Maier C, Motegi C, et al. Viral ecology of organic and inorganic particles in aquatic systems: avenues for further research. Aquat Micro Ecol. 2009;57:321–41.CAS 

Google Scholar 
Tian Y, Cai L, Xu Y, Luo T, Zhao Z, Wang Q, et al. Stability and infectivity of allochthonous viruses in deep sea: A long-term high pressure simulation experiment. Deep-Sea Res Part I-Oceanogr Res Pap. 2020;161:103302.
Google Scholar 
Lara E, Vaqué D, Sà EL, Boras JA, Gomes A, Borrull E, et al. Unveiling the role and life strategies of viruses from the surface to the dark ocean. Sci Adv. 2017;3:e1602565.PubMed 
PubMed Central 

Google Scholar 
Zhang R, Li Y, Yan W, Wang Y, Cai L, Luo T, et al. Viral control of biomass and diversity of bacterioplankton in the deep sea. Commun Biol. 2020;3:256.CAS 
PubMed 
PubMed Central 

Google Scholar 
Woźniak SB, Stramski D, Stramska M, Reynolds RA, Wright VM, Miksic EY, et al. Optical variability of seawater in relation to particle concentration, composition, and size distribution in the nearshore marine environment at Imperial Beach, California. J Geophys Res. 2010;115:C08027.
Google Scholar 
White AE, Letelier RM, Whitmire AL, Barone B, Bidigare RR, Church MJ, et al. Phenology of particle size distributions and primary productivity in the North Pacific subtropical gyre (Station ALOHA). J Geophys Res-Oceans. 2015;120:7381–99.PubMed 
PubMed Central 

Google Scholar 
Vaulot D, Courties C, Partensky F. A simple method to preserve oceanic phytoplankton for flow cytometric analyses. Cytom Part A. 1989;10:629–35.CAS 

Google Scholar 
Chen X, Liu H, Weinbauer M, Chen B, Jiao N. Viral dynamics in the surface water of the western South China Sea in summer 2007. Aquat Micro Ecol. 2011;63:145–60.
Google Scholar 
Wei W, Zhang R, Peng L, Liang Y, Jiao N. Effects of temperature and photosynthetically active radiation on virioplankton decay in the western Pacific Ocean. Sci Rep. 2018;8:1525–34.PubMed 
PubMed Central 

Google Scholar 
Marie D, Partensky F, Vaulot D, Brussaard C. Numeration of phytoplankton, bacteria and viruses in marine samples. Curr Protoc Cytom. 1999;11:1–15.
Google Scholar 
Marie D, Brussaard CPD, Thyrhaug R, Bratbak G, Vaulot D. Enumeration of marine viruses in culture and natural samples by flow cytometry. Appl Environ Microbiol. 1999;65:45–52.CAS 
PubMed 
PubMed Central 

Google Scholar 
Brussaard CP. Optimization of procedures for counting viruses by flow cytometry. Appl Environ Microbiol. 2004;70:1506–13.CAS 
PubMed 
PubMed Central 

Google Scholar 
Wilhelm SW, Brigden SM, Suttle CA. A dilution technique for the direct measurement of viral production: a comparison in stratified and tidally mixed coastal waters. Micro Ecol. 2002;43:168–73.CAS 

Google Scholar 
Weinbauer MG, Rowe JM, Wilhelm SW. Determining rates of virus production in aquatic systems by the virus reduction approach. In: Wilhelm SW, Weinbauer MG, Suttle CA. (eds). Manual of Aquatic Viral Ecology. American Society of Limnology and Oceanography Inc., Waco, TX: 2010. p. 1–8.Chen X, Wei W, Wang J, Li H, Sun J, Ma R, et al. Tide driven microbial dynamics through virus-host interactions in the estuarine ecosystem. Water Res. 2019;160:118–29.CAS 
PubMed 

Google Scholar 
Luef B, Luef F, Peduzzi P. Online program ‘vipcal’ for calculating lytic viral production and lysogenic cells based on a viral reduction approach. Environ Microbiol Rep. 2009;1:78–85.CAS 
PubMed 
PubMed Central 

Google Scholar 
Winget DM, Helton RR, Williamson KE, Bench SR, Williamson SJ. Repeating patterns of virioplankton production within an estuarine ecosystem. Proc Natl Acad Sci USA. 2011;108:11506–11.CAS 
PubMed 
PubMed Central 

Google Scholar 
Wei W, Wang N, Cai L, Zhang C, Jiao N, Zhang R. Impacts of freshwater and seawater mixing on the production and decay of virioplankton in a subtropical estuary. Micro Ecol. 2019;78:843–54.CAS 

Google Scholar 
Noble RT, Fuhrman JA. Virus decay and its causes in coastal waters. Appl Environ Microbiol. 1997;63:77–83.CAS 
PubMed 
PubMed Central 

Google Scholar 
Suttle CA, Chen F. Mechanisms and rates of decay of marine viruses in seawater. Appl Environ Microbiol. 1992;58:3721–9.CAS 
PubMed 
PubMed Central 

Google Scholar 
Rowe JM, Saxton MA, Cottrell MT, DeBruyn JM, Berg GM, Kirchman DL, et al. Constraints on viral production in the Sargasso Sea and North Atlantic. Aquat Micro Ecol. 2008;52:233–44.
Google Scholar 
Evans C, Pearce I, Brussaard CP. Viral-mediated lysis of microbes and carbon release in the sub-Antarctic and Polar Frontal zones of the Australian Southern Ocean. Environ Microbiol. 2009;11:2924–34.CAS 
PubMed 

Google Scholar 
De Corte D, Sintes E, Winter C, Yokokawa T, Reinthaler T, Herndl GJ. Links between viral and prokaryotic communities throughout the water column in the (sub)tropical Atlantic Ocean. ISME J. 2010;4:1431–42.PubMed 

Google Scholar 
Li Y, Lou T, Sun J, Cai L, Liang Y, Jiao N, et al. Lytic viral infection of bacterioplankton in deep waters of the western Pacific Ocean. Biogeosciences. 2014;11:2531–42.
Google Scholar 
Liang Y, Zhang Y, Zhang Y, Luo T, Rivkin R, Jiao N. Distributions and relationships of virio- and picoplankton in the epi-, meso- and bathypelagic zones of the Western Pacific Ocean. FEMS Microbiol Ecol. 2017;93:fiw238.PubMed 

Google Scholar 
Wommack KE, Colwell RR. Virioplankton: viruses in aquatic ecosystems. Microbiol Mol Biol Rev. 2000;64:69–114.CAS 
PubMed 
PubMed Central 

Google Scholar 
Parikka KJ, Le Romancer M, Wauters N, Jacquet S. Deciphering the virus-to-prokaryote ratio (VPR): insights into virus-host relationships in a variety of ecosystems. Biol Rev. 2016;92:1081–1100.PubMed 

Google Scholar 
Parada V, Herndl GJ, Weinbauer MG. Viral burst size of heterotrophic prokaryotes in aquatic systems. J Mar Biol Assoc UK. 2006;86:613–21.
Google Scholar 
Yuan D. A numerical study of the South China Sea deep circulation and its relation to the Luzon Strait transport. Acta Oceano Sin. 2002;21:187–202.
Google Scholar 
Tian J, Yang Q, Zhao W. Enhanced diapycnal mixing in the South China Sea. J Phys Oceanogr. 2009;39:3191–203.
Google Scholar 
Alford MH, Lien R, Simmons H, Klymak J, Ramp S, Yang YJ, et al. Speed and evolution of nonlinear internal waves transiting the South China Sea. J Phys Oceanogr. 2010;40:1338–55.
Google Scholar 
Parada V, Sintes E, Van Aken HM, Weinbauer MG, Herndl GJ. Viral abundance, decay, and diversity in the meso- and bathypelagic waters of the north atlantic. Appl Environ Microbiol. 2007;73:4429–38.CAS 
PubMed 
PubMed Central 

Google Scholar 
De Corte D, Sintes E, Yokokawa T, Reinthaler T, Herndl GJ. Links between viruses and prokaryotes throughout the water column along a North Atlantic latitudinal transect. ISME J. 2012;6:1566–77.PubMed 
PubMed Central 

Google Scholar 
Zachary A. An ecological study of bacteriophages of Vibrio natriegens. Appl Environ Microbiol. 1978;24:321–4.CAS 

Google Scholar 
Motegi C, Nagata T. Enhancement of viral production by addition of nitrogen or nitrogen plus carbon in subtropical surface waters of the South Pacific. Aquat Micro Ecol. 2007;48:27.
Google Scholar 
Bratbak G, Egge JK, Heldal M. Viral mortality of the marine alga Emiliania huxleyi (Haptophyceae) and termination of algal blooms. Mar Ecol-Prog Ser. 1993;93:39–48.
Google Scholar 
Motegi C, Kaiser K, Benner R, Weinbauer MG. Effect of P-limitation on prokaryotic and viral production in surface waters of the Northwestern Mediterranean Sea. J Plankton Res. 2015;37:16–20.CAS 

Google Scholar 
Hewson I, O’Neil JM, Fuhrman JA, Dennison WC. Virus-like particle distribution and abundance in sediments and overmaying waters along eutrophication gradients in two subtropical estuaries. Limnol Oceanogr. 2001;46:1734–46.
Google Scholar 
Wilson WH, Mann NH. Lysogenic and lytic viral production in marine microbial communities. Aquat Micro Ecol. 1997;13:95–100.
Google Scholar 
Paul JH. Prophages in marine bacteria: dangerous molecular time bombs or the key to survival in the seas? ISME J. 2008;2:579–89.CAS 
PubMed 

Google Scholar 
Chibani-Chennoufi S, Bruttin A, Dillmann ML, Brussow H. Phage-host interaction: an ecological perspective. J Bacteriol. 2004;186:3677–86.CAS 
PubMed 
PubMed Central 

Google Scholar 
Weinbauer MG. Ecology of prokaryotic viruses. FEMS Microbiol Rev. 2004;28:127–81.CAS 
PubMed 

Google Scholar 
Williamson SJ, Paul JH. Nutrient stimulation of lytic phage production in bacterial populations of the Gulf of Mexico. Aquat Micro Ecol. 2004;36:9–17.
Google Scholar 
Williamson SJ, Paul JH. Environmental factors that influence the transition from lysogenic to lytic existence in the ϕHSIC/Listonella pelagia marine phage–host system. Micro Ecol. 2006;52:217–25.CAS 

Google Scholar 
Cissoko M, Desnues A, Bouvy M, Sime-Ngando T, Verling E, Bettarel Y. Effects of freshwater and seawater mixing on virio- and bacterioplankton in a tropical estuary. Freshw Biol. 2008;53:1154–62.
Google Scholar 
Bettarel Y, Bouvier T, Agis M, Bouvier C, Van Chu T, Combe M, et al. Viral distribution and life strategies in the Bach Dang Estuary, Vietnam. Micro Ecol. 2011;62:143–54.
Google Scholar 
Shkilnyj P, Koudelka GB. Effect of salt shock on stability of λimm434 lysogens. J Bacteriol. 2007;189:3115–23.CAS 
PubMed 
PubMed Central 

Google Scholar 
Tuomi P, Fagerbakke KM, Bratbak G, Heldal M. Nutritional enrichment of a microbial community: the effects on activity, elemental composition, community structure and virus production. FEMS Microbiol Ecol. 1995;16:23–134.
Google Scholar 
Dell’Anno A, Corinaldesi C, Danovaro R. Virus decomposition provides an important contribution to benthic deep-sea ecosystem functioning. Proc Natl Acad Sci USA. 2015;112:E2014–E2019.PubMed 
PubMed Central 

Google Scholar 
Mojica KD, Brussaard CP. Factors affecting virus dynamics and microbial host-virus interactions in marine environments. FEMS Microbiol Ecol. 2014;89:495–515.CAS 
PubMed 

Google Scholar 
Zweifel UL. Factors controlling accumulation of labile dissolved organic carbon in the Gulf of Riga. Estuar Coast Shelf Sci. 1999;48:357–70.CAS 

Google Scholar 
Pomeroy LR, Wiebe WJ. Temperature and substrates as interactive limiting factors for marine heterotrophic bacteria. Aquat Micro Ecol. 2001;23:187–204.
Google Scholar 
Ploug H, Grossart H, Azam F, Jørgensen BB. Photosynthesis, respiration, and carbon turnover in sinking marine snow from surface waters of Southern California Bight: implications for the carbon cycle in the ocean. Mar Ecol-Prog Ser. 1999;179:1–11.CAS 

Google Scholar 
Azam F, Malfatti F. Microbial structuring of marine ecosystems. Nature. 2007;5:782–91.CAS 

Google Scholar 
De Corte D, Sintes E, Yokokawa T, Lekunberri I, Herndl GJ. Large-scale distribution of microbial and viral populations in the South Atlantic Ocean. Environ Microbiol Rep. 2016;8:305–15.PubMed 
PubMed Central 

Google Scholar 
Yang YH, Yokokawa T, Motegi C, Nagata T. Large-scale distribution of viruses in deep waters of the Pacific and Southern Oceans. Aquat Micro Ecol. 2014;71:193–202.
Google Scholar 
Labonté JM, Swan BK, Poulos B, Luo H, Koren S, Hallam SJ, et al. Single-cell genomics-based analysis of virus-host interactions in marine surface bacterioplankton. ISME J. 2015;9:2386–99.PubMed 
PubMed Central 

Google Scholar 
Martinez-Hernandez F, Fornas Ò, Lluesma Gomez M, Garcia-Heredia I, Maestre-Carballa L, López-Pérez M, et al. Single-cell genomics uncover Pelagibacter as the putative host of the extremely abundant uncultured 37-F6 viral population in the ocean. ISME J. 2019;13:232–6.CAS 
PubMed 

Google Scholar 
Mruwat N, Carlson MCG, Goldin S, Ribalet F, Kirzner S, Hulata Y, et al. A single-cell polony method reveals low levels of infected Prochlorococcus in oligotrophic waters despite high cyanophage abundances. ISME J. 2021;15:41–54.CAS 
PubMed 

Google Scholar 
Peduzzi P, Weinbauer M. Effect of concentrating the virus-rich 2–200 nm size fraction of seawater on the formation of algal flocs (marine snow). Limnol Oceanogr. 1993;38:1562–5.
Google Scholar 
Uitz J, Stramski D, Baudoux A, Reynolds RA, Wright VM, Dubranna J, et al. Variations in the optical properties of a particle suspension associated with viral infection of marine bacteria. Limnol Oceanogr. 2010;55:2317–30.
Google Scholar 
Sullivan MB, Weitz JS, Wilhelm SW. Viral ecology comes of age. Environ Microbiol Rep. 2017;9:33–35.PubMed 

Google Scholar 
Laber CP, Hunter JE, Carvalho F, Collins JR, Hunter EJ, Schieler BM, et al. Coccolithovirus facilitation of carbon export in the North Atlantic. Nat Microbiol. 2018;3:537–47.CAS 
PubMed 

Google Scholar 
Kranzler CF, Brzezinski MA, Cohen NR, Lampe RH, Maniscalco M, Till CP, et al. Impaired viral infection and reduced mortality of diatoms in iron-limited oceanic regions. Nat Geosci. 2021;4:231–7.
Google Scholar 
Hewson I, Fuhrman JA. Viriobenthos production and virioplankton sorptive scavenging by suspended sediment particles in coastal and pelagic waters. Micro Ecol. 2003;46:337–47.CAS 

Google Scholar  More

Study origin and designThe MeadoWatch project (MW) is a project run collaboratively between the University of Washington (UW) and the United States National Park Service to monitor the phenology of alpine and subalpine wildflower species across large elevational gradients in Mount Rainier National Park (Fig. 2). MW was established in 2013 with the goal of understanding long-term effects of climate change on Mount Rainier National Park wildflower communities using community-science approaches. The first MW transect was established along Reflection Lakes, Skyline, and Paradise Glacier trail system in 2013 (9–11 plots). In 2015, MW expanded to include a second transect (15–17 plots) along the Glacier Basin trail (Fig. 1a). The MW transects span around 5 km each, over a 400 m altitudinal gradient (Reflection Lakes: 1490m–1889m a.s.l.; Glacier Basin: 1460m–1831m a.s.l.)Fig. 2Alpine meadows, plot extension, and target species. (a) Species-rich alpine meadow in Mount Rainier National Park (Mount Tahoma), showing many of the target species in the foreground. (b) MW volunteer coordinator Anna Wilson at a plot, indicating the arm span that defines the plot area (personal likeness used with confirmed consent). (c) Species composition and proportion of reports per species in each of the transects; species common to both trails are highlighted with striped shadowing. Photographs: A. John (a), L. Felker (b).Full size imagePlots are located along the side of each trail, marked with a colored survey marker. The area of each plot is defined by the arm-span of volunteers when they position themselves over the plot marker looking away from the trail (Fig. 2b). While less accurate than marking the corners of plots, this approach was used to avoid establishing permanent structures in wilderness areas within the National Park. The surveyed area in each plot is, on average, 1.25 m2. Each plot is also equipped with temperature sensors (HOBO Pendant Logger, Onset Computer Corp.) buried approximately 4 cm below the ground. Sensors are placed at the start of each fall season and removed at the beginning of each summer season for data downloading. The HOBO sensors provide an estimate for the date of snow disappearance and in-situ temperature at 3 hours intervals. Once plots are covered in snow, soil temperatures remain at 0 °C and show no diurnal variation, so that daily changes in temperatures above 1 °C can be used to determine the earliest date without snow cover20. We use these approaches to provide dates of snow appearance and disappearance, snow cover duration, and minimum soil temperatures for each year and plot. Occasionally, temperature data during the snow disappearing window were lost due to sensor failure or loss of sensors (which occurs because plots are not permanently marked and/or well-meaning visitors remove sensors). This, and the lack of temperature sensors in the first year of the project, resulted in approx. 20% of cases of missing data. In those cases, we used a data imputation method to estimate the missing values based on data from nearby plots and a parallel temperature data collection with 890 total observations. These estimates were highly reliable in filling the data gaps (see Appendix C in16 for further details).Focal speciesWe originally targeted 16 native wildflower species along each transect, which were chosen based on their abundance, ease of identification, and presence in the plot. Four of those target species were present on both transects. In 2016 we replaced one species with a different one (see further information below), making for a total of 17 species monitored (Fig. 2c). The focal species are: American bistort* (Polygonum bistortoides), avalanche lily (Erythronium montanum), bracted lousewort* (Pedicularis bracteosa), broadleaf arnica (Arnica latifolia), cascade aster (Aster ledophyllus; synonym Eucephalus ledophyllus), glacier lily (Erythronium grandiflorum), Gray’s lovage (Ligusticum grayi), magenta paintbrush (Castilleja parviflora), mountain daisy (Erigenon peregrinus; synonym Erigeron glacialis), northern microseris (Microseris alpestris; synonym Nothocalais alpestris), scarlet paintbrush (Castilleja miniata), sharptooth angelica (Angelica arguta), sitka valerian* (Valeriana sitchensis), subalpine lupine* (Lupinus arcticus; synonym Lupinus latifolius var. subalpinus), tall bluebell (Mertensia paniculata), Canby’s licorice-root (Ligusticum canbyi), and western anemone (Anemone occidentalis). Asterisks denote species monitored along both trails.Due to challenges in species identification, we dropped Canby’s licorice-root (Ligusticum canbyi) as a target species in 2016. Consequently, Ligusticum canbyi has limited replication in the database (Fig. 2c). While we included the phenological records of this species for the sake of completeness, we recommend focusing on the other 16 species, which are both better represented (in terms of data coverage) and are free of any potential misidentification issues.For additional information on the species, methods, identification cues, and image resources see: http://www.meadowatch.org, https://www.youtube.com/channel/UCGBFTKxf8FIWswMDxBavpuQ, and the appendices therein16.Data collection and volunteer trainingDuring the summer months, MW volunteers and scientists collect reproductive phenology data with a frequency between 3 and 9 trail reports per week. Each report records the presence or absence of 4 phenophases for each target species present in each of the plots. The phenophases are ‘budding’, ‘flowering’, ‘ripening fruit’, and ‘releasing seed’. Phenophases were defined as follows:BuddingThe beginning growth of the flower which has not yet opened. A plant is considered budding if buds are present, but the stamen and pistils are not yet visible and available to pollinators.FloweringThe generally “showy” part of the plant that holds the reproductive parts (stamens and pistils). A plant is considered flowering when at least one flower is open, and the stamens and pistils are visible and available for pollination and reproduction.Ripening fruitThe fruit develops from the female part of the flower following successful pollination. In the target species, fruits can take many forms, from hard, fleshy capsules, juicy berries, to a feathery tuft on the end of a seed. A plant is in the ripening fruit stage when reproductive parts on at least one reproductive stalk are non-functional and the formation of the fruit part is clearly ongoing (visible), but seeds are not yet fully mature and not yet being dispersed.Releasing seedAfter the fruit ripens, seeds are released to be dispersed by gravity, wind, or animals. A plant is considered in the releasing seed stage if seeds are actively being released on at least one reproductive stalk (but there are still seeds present).A full description, including illustrations for each species’ phenophase and identification cues is available in http://www.meadowatch.org/volunteer-resources.html, as well as in Annex 1 – Supplementary Documentation. Multiple phenophases can be present simultaneously, depending on the species, and are noted independently. Additionally, volunteers are also asked to record the presence of snow (‘snow covered plot’, ‘partially covered plot’, or ‘snow-free plot’), and, since 2017, the presence of damage by herbivory (‘presence of browsed stems’) on each plot.In years not impacted by the SARS-Cov-2 pandemic MW volunteers attend an in-person 3-hour botanical and phenological training session taught by UW scientists at the beginning of each sampling season. Volunteers were also provided with detailed species-identification guides, including an extensive description of sampling methods and location of the plots. The trainings for the 2020 and 2021 seasons were held virtually via a series of online training videos. In these years, volunteers took a quiz on wildflower phenology, plant identification and data collection methods after viewing these videos and were required to ‘pass’ a certain threshold to volunteer (unlimited attempts were allowed). During these virtual trainings, volunteers were provided with digital copies of the species identification guides, with many returning volunteers using printed guides they had kept from previous years.At the end of their phenological hike, volunteers submit their data sheets either by depositing them in lockboxes located within the park, or by scanning and emailing them directly to mwatch@uw.edu. Data are then entered manually and stored in the UW repositories after being checked for consistency at the end of each sampling season.The parallel data collection by members of UW’s Hille Ris Lambers group (including PI, postdoctoral researchers, graduate students, and trained interns) acted as the following: (i) a quality-control, i.e., allowing us to compare the consistency in phenology assessments between volunteers and scientists, and (ii) a way to increase the temporal resolution and scale of the data, e.g., by reducing early season gaps and ‘weekend bias’17. This parallel expert sampling was carried out around once a week between 2013 and 2020, showing great consistency between the two groups. For detailed comparisons between volunteers and scientists’ assessments see the data validation section (as well as Appendix E in16).Processed dataIn addition to the raw phenological data, we also provide here parameters to construct the year, species, and plot-specific flowering phenology based on the timing of snow disappearance (as in16). Models describe unimodal probability distributions that were fitted with maximum likelihood models to binomial flowering data from each species, year, and plot. These curves have been used to estimate peak flowering dates and diversity and link them to reported visitor experiences16. Here, we provide the 3 parameters defining the unimodal curve of flowering probability per species i, plot j and year k: the duration of flowering (𝛿ijk), the maximum probability of flowering (𝜇ijk), and peak flowering (in DOY – ρijk); following the equations described in16 and https://github.com/ajijohn/MeadoWatch).The parameters of these probability distribution curves are ready-to-use values that can be broadly and easily used to estimate floral compositional change over past seasons due to changing environmental conditions—for example, to inform plant-pollinator interaction networks if combined with pollinator behavioral data (e.g.21). More

Disease data sourceAll analyses were conducted at the administrative level, and the exact list of known zoonotic diseases is recorded in the GIDEON database22. GIDEON is currently the most comprehensive and frequently updated infectious disease outbreak database reporting epidemics of human infectious diseases at the global scale and has been widely used in global zoonosis studies42,43 (Last access date, November 9, 2020). The administrative designations used in our analyses were based on the Global Administrative Areas (GADM) database (www.gadm.org, downloaded on November 8, 2020), which includes very detailed boundary data for global countries and major island groups.Pattern and correlates of zoonosis events worldwideNumber of zoonosis eventsGIDEON defines human infectious disease reservoirs as any animal, plant, or substrate supporting the survival and reproduction of infectious agents and promoting transmission to potential susceptible hosts. Its host category therefore includes all human-specific, zoonotic, multihost, and environmental agents. As our main aim was to test the role of established alien animal species in the emergence of zoonotic diseases, we focused on a total of 161 diseases specified in GIDEON’s host designations and definitions as nonhuman zoonotic (n = 115) and multihost (n = 46) diseases (Supplementary Data 1) and excluded diseases with human-specific hosts that do not need animals to persist or be transmitted. The infectious agents of nonhuman zoonotic diseases complete their entire lifecycle in nonhuman hosts but may have the potential to spillover and infect human populations. Infectious agents of multihost diseases can use both human and animal hosts for their development and reproduction. We measured the number of zoonosis events for each jurisdiction according to five host taxonomic groups: mammals, birds, invertebrates, reptiles and amphibians. These zoonoses were mainly caused by bacteria, viruses, parasitic animals and fungi. We excluded zoonoses from the Algae (3 diseases) due to low sample sizes in GIDEON.Correlates of the number of zoonosis eventsClimatic variablesFollowing a previous study21, we used global environmental stratification (GEnS) as a composite bioclimatic variable generated by stratifying the Earth’s surface into zones with similar climates44. The GEnS database was constructed based on a total of 125 strata across 18 global environmental zones with a spatial resolution of 30 arc seconds (equivalent to approximately 0.86 km2 at the equator). The values in GEnS range from 1 to 18 with a higher value indicating warmer and wetter conditions.Human population densityWe used human population density as one general anthropogenic factor reflecting propagule pressure and human-assisted pathogen movements1,21,45. Human population size data and the land area of each jurisdiction were collected from World Bank Open Data from 2011 to 2020 (available at https://data.worldbank.org/indicator/SP.POP.TOTL, accessed on November 18, 2020). We then calculated the human population density using the human population size divided by the land area.Native potential host richness and biodiversity lossData on the richness of native amphibians, birds, and mammals were derived from the Biodiversity Mapping website (https://biodiversitymapping.org/wordpress/index.php/home/, accessed on August 19, 2020), which were based on studies from Jenkins et al. (2013)’s and Pimm et al. (2014)46,47. The map of reptile diversity is based on an updated database of the global spatial distribution of reptiles48. All diversity maps for each taxon were generated through the calculation of grid-based richness at a spatial resolution of 10 km × 10 km in ArcGIS46. We did not include native invertebrate richness, as global maps for most invertebrate taxa are not yet available. For the loss of native biodiversity, we followed the previous study by first extracting the list of threatened species (NT, EN and VU categories evaluated by the IUCN Red List, access on May 10th, 2021)29, and then calculated the number of threatened species for each taxon distributes in each administrative unite as a proxy of biodiversity loss.Richness of established alien zoonotic host speciesWe quantified the richness of established alien animal species from the five main taxonomic groups (mammals, birds, reptiles, amphibians and invertebrates) based on 4,522 establishment events of 795 alien animals in each of 201 jurisdictions according to various databases. Data on 262 established alien reptiles and amphibians were compiled from multiple publications, including Kraus’s compendium49 and other recent updates50. Data on 337 established alien birds after removing all migratory bird species as vagrants were collected from the Global Avian Invasions Atlas (GAVIA)51, which is a comprehensive database of the global distribution of established alien birds. Data on 119 established alien mammals were obtained from the Introduced Mammals of the World database52 and the more recent update53. Data on 77 terrestrial alien invertebrates (66 insects and 11 other groups) across 7 taxa with native and invaded range information were obtained from the Global Invasive Species Database (GISD, http://www.iucngisd.org/gisd/, accessed on July 1, 2020). We calculated the richness of both zoonotic and non-zoonotic alien host species for each order. We first conducted an intensive literature review for each established alien species of each of the four taxa to determine whether they transmit pathogens to humans (Supplementary Data 2). The identification of zoonotic or non-zoonotic host may be influenced by under-sampling in the literature. We therefore incorporated the latest synthesis of human-infecting pathogens in the ‘CLOVER’ dataset to identify zoonotic and non-zoonotic animal hosts54. The CLOVER dataset compiled GMPD255, EID256, HP323 and Shaw57 databases and is currently the most comprehensive dataset on host-pathogen associations. Based on this information, we then categorized each alien species as a ‘zoonotic host’ or ‘non-zoonotic host’. The records of the established alien species were assigned to GADM jurisdictions, and we calculated the richness of the established alien zoonotic and non-zoonotic host species for each taxonomic group within each jurisdiction. In order to increase the statistical power, we conducted subsequent modeling analyses based on four mammalian orders (i.e., Carnivora, Cetartiodactyla, Lagomorpha, and Rodentia), five avian groups (i.e., waterfowl including five orders: Anseriformes, Gruiformes, Pelecaniformes, Phoenicopteriformes and Suliformes; Columbiformes, Galliformes, Passeriformes, Psittaciformes), the order Diptera of the invertebrates, and herpetofauna as a whole, which have established alien populations in at least 50 administrative units.Climate changeWe extracted historical monthly mean temperature and precipitation data recorded between 1901 and 2009 from the University of East Anglia Climate Research Unit (CRU, https://sites.uea.ac.uk/cru/, accessed on November 30, 2020)58. This database provides historical global-scale yearly climatic data with the finest resolution of 0.5° grids. We generated the temperature and precipitation values for all grids in each jurisdiction, calculated the slope of the temperature and precipitation for the time series of the years 1901 to 2009 for each grid and generated the averages based on all grids within each jurisdiction.Anthropogenic land-use changeWe downloaded global land-use data from the Anthromes v2 Dataset (Anthropogenic Biomes version 2, accessed on October 15, 2020) in ESRI GRID format59. We used the 1900 and 2000 data to calculate the temporal changes in land use. By using the reclassify and raster function in ArcGIS, we calculated the percentage of grids in which the land-use type changed to a more anthropogenically influenced type from 1900 to 2000 for each jurisdiction, including 15 scenarios: Wildlands to Seminatural, Wildlands to Rangelands, Wildlands to Croplands, Wildlands to Villages, Wildlands to Dense Settlements, Seminatural to Rangelands, Seminatural to Croplands, Seminatural to Villages, Seminatural to Dense Settlements, Rangelands to Croplands, Rangelands to Villages, Rangelands to Dense Settlements, Croplands to Villages, Croplands to Dense Settlements, and Villages to Dense Settlements.Sampling effort, reporting bias and incomplete dataA potential issue in quantifying the effects of different predictor variables on the number of zoonosis events is the need to account for the differences in survey effort, reporting bias and incomplete disease data among regions1,21,28. There is a high probability that zoonosis discovery is spatially biased by uneven levels of surveillance across countries, as the global allocation of scientific resources has been focused on rich and developed countries. We thus included the Infectious Disease Vulnerability Index (IDVI), which is a comprehensive metric reflecting the demographic, health care, public health, socioeconomic, and political factors that may have an impact on the capacity of surveillance and detection of infectious diseases in each country60. Second, we followed the methods of a previous study21 to control for reporting biases. We incorporated PubMed citations per disease for each jurisdiction using a Python-based PubCrawler21. In addition, we added the longitude and latitude of the geographic centroid of administrative units to control for spatial autocorrelation as there would be a higher probability of having similar diseases in nearby than distant administrative units61.Statistical analysisThe number of zoonosis events, native potential host richness, established alien animal richness and human population density were log-transformed to improve linearity. A potential issue in our data analysis is that the numbers of zoonosis events and the numbers of native and alien animal species are strongly influenced by geographical area, as larger countries or regions may host more native or alien animal species and more disease events. We therefore calculated the density of native or alien species richness and the number of zoonosis events using the total number divided by the geographical area of each jurisdiction. Furthermore, the number of zoonosis events may also be influenced by the degree of local disease surveillance. We thus obtained the residuals from a regression correlating zoonosis event density and all disease event density, and used them as the dependent variable for further analyses (Fig. 1). As some of our variables may be expected to be nonlinear, we performed generalized additive mixed model (GAMM) analyses following Mollentze & Streicker 2020’s framework25 to quantify the relationships between different predictor variables and the number of zoonosis events. We started with a full model with zoonosis event density controlling for overall disease surveillance as the response variable and 13 smoothed fixed effects (Fig. 1 and Supplementary Data 4): GEnS, human population density, density of native species richness, biodiversity loss, density of alien zoonotic host richness, density of alien non-zoonotic host richness, climate (temperature and precipitation) change, land-use change, IDVI, PubMed citations, longitude and latitude of geographic centroid of administrative units. The reason why we included the density of alien non-zoonotic host richness as a covariate is because this variable can serve as a positive control for propagule pressure, allowing us to more explicitly test whether zoonotic alien hosts contribute to zoonoses beyond propagule pressure associated with non-zoonotic alien hosts, which cannot directly increase zoonotic diseases. These predictor variables were not highly collinear as their correlation coefficients based on Pearson rank correlation analyses were all More