Philippa Kaur

1.
Pinch, G. Handbook of Egyptian Mythology (ABC-CLIO, Santa Barbara, 2002).
Google Scholar 
2.
Ikram, S. Animals in ancient Egyptian religion: belief, identity, power, and economy. In The Oxford Handbook of Zooarchaeolog (ed. Viner-Daniels, S.) 452–465 (Oxford University Press, Oxford, 2017).
Google Scholar 

3.
Ikram, S. An eternal aviary: bird mummies from ancient egypt. In Between heaven and earth: Birds in Ancient Egypt 232 (ed. Rozenn, B. L.) (Rozenn Bailleul-Le Suer, London, 2012).
Google Scholar 

4.
Von den Driesch, A., Kessler, D., Steinmann, F., Berteaux, V. & Peters, J. Mummified, Deified and Buried at Hermopolis Mgna: The Sacred Birds from Tuna el-Gebel, Middle Egypt. in Ägypten und Levante vol. 15 203–44 (Manfred Bietak, Vienna, 2005).

5.
de Diodore, S. Introduction générale. in Bibliothèque historique vol. 1 74 (1993).

6.
Strabon. Le voyage en Egypte. vol. 4 (1997).

7.
Ray, J. D. The archives of Hor. Texts from Excavations 2 (Londres EES, London, 1976).
Google Scholar 

8.
Martin, G. T. The Sacred Animal Necropolis at North Saqqâra: The Southern Dependencies of the Main Temple Complex Vol. 50 (Egypt Exploration Society, London, 1981).
Google Scholar 

9.
Meeks, D. Les couveuses artificielles en Égypte. in Techniques et économie antiques et médiévales: le temps de l’innovation: colloque international (C.N.R.S.) Aix-en-Provence, 21–23 mai 1996, 1997, Travaux du Centre Camille Jullian 21 132–136 (Errance, Providence, 1997).

10.
de Davies, N. G. Ancient Egyptian Paintings (University Press, Cambridge, 1936).
Google Scholar 

11.
Wasef, S. et al. Mitogenomic diversity in Sacred Ibis Mummies sheds light on early Egyptian practices. PLoS ONE 14, e0223964 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

12.
Amiot, R. et al. Oxygen isotope fractionation between bird bone phosphate and drinking water. Sci. Nat. 104, 47 (2017).
Article  CAS  Google Scholar 

13.
Craig, H. & Gordon, L. I. Deuterium and Oxygen 18 Variations in the Ocean and the Marine Atmosphere (Laboratorio di geologia nucleare, Pisa, 1965).
Google Scholar 

14.
Dansgaard, W. Stable isotopes in precipitation. Tellus 16, 436–468 (1964).
ADS  Article  Google Scholar 

15.
IAEA/WMO. Global Network of Isotopes in Precipitation. https://www.iaea.org/water (2019).

16.
IAEA/WMO. Global Network of Isotopes in Rivers. Global Network of Isotopes in Rivers https://www.iaea.org/water (2019).

17.
Stewart, M. K. Stable isotope fractionation due to evaporation and isotopic exchange of falling waterdrops: applications to atmospheric processes and evaporation of lakes. J. Geophys. Res. 80, 1133–1146 (1975).
ADS  CAS  Article  Google Scholar 

18.
Hobson, K. A. & Clark, R. G. Assessing avian diets using stable isotopes II: factors influencing diet-tissue fractionation. Condor Ornithol. Appl. 94, 189–197 (1992).
Google Scholar 

19.
Passey, B. H. et al. Carbon isotope fractionation between diet, breath CO2, and bioapatite in different mammals. J. Archaeol. Sci. 32, 1459–1470 (2005).
Article  Google Scholar 

20.
Mizutani, H., Fukuda, M. & Kabaya, Y. 13C and 15N enrichment factors of feathers of 11 species of adult birds. Ecology 73, 1391–1395 (1992).
Article  Google Scholar 

21.
Angst, D. et al. Isotopic and anatomical evidence of an herbivorous diet in the Early Tertiary giant bird Gastornis. Implications for the structure of Paleocene terrestrial ecosystems. Naturwissenschaften 101, 313–322 (2014).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

22.
Hobson, K. A. Reconstructing Avian diets using stable-carbon and nitrogen isotope analysis of egg components: patterns of isotopic fractionation and turnover. The Condor 97, 752–762 (1995).
Article  Google Scholar 

23.
Peterson, B. J. & Fry, B. Stable isotopes in ecosystem studies. Annu. Rev. Ecol. Syst. 18, 293–320 (1987).
Article  Google Scholar 

24.
Caccamise, D. F., Reed, L. M., Castelli, P. M., Wainright, S. & Nichols, T. C. Distinguishing migratory and resident Canada geese using stable isotope analysis. J. Wildl. Manag. 64, 1084–1091 (2000).
Article  Google Scholar 

25.
Hebert, C. E., Bur, M., Sherman, D. & Shutt, J. L. Sulfur isotopes link overwinter habitat use and breeding condition in double-crested cormorants. Ecol. Appl. 18, 561–567 (2008).
PubMed  Article  PubMed Central  Google Scholar 

26.
Lott, C. A., Meehan, T. D. & Heath, J. A. Estimating the latitudinal origins of migratory birds using hydrogen and sulfur stable isotopes in feathers: influence of marine prey base. Oecologia 134, 505–510 (2003).
ADS  PubMed  Article  PubMed Central  Google Scholar 

27.
DeNiro, M. J. & Epstein, S. Influence of diet on the distribution of nitrogen isotopes in animals. Geochim. Cosmochim. Acta 45, 341–351 (1981).
ADS  CAS  Article  Google Scholar 

28.
Blum, J. D., Taliaferro, E. H. & Holmes, R. T. Changes in Sr/Ca, Ba/Ca and 87Sr/86Sr ratios between trophic levels in two forest ecosystems in the northeastern USA. Biochemistry 49, 87–101 (2000).
CAS  Google Scholar 

29.
Capo, R. C., Stewart, B. W. & Chadwick, O. A. Strontium isotopes as tracers of ecosystem processes: theory and methods. Geoderma 82, 197–225 (1998).
ADS  CAS  Article  Google Scholar 

30.
Graustein, W. C. 87Sr/86Sr ratios measure the sources and flow of strontium in terrestrial ecosystems. In Stable Isotopes in Ecological Research (eds Rundel, P. W. et al.) 491–512 (Springer, New York, 1989).
Google Scholar 

31.
Blum, J. D., Taliaferro, E. H. & Holmes, R. T. Determining the sources of calcium for migratory songbirds using stable strontium isotopes. Oecologia 126, 569–574 (2001).
ADS  PubMed  Article  PubMed Central  Google Scholar 

32.
Libby, W. F., Berger, R., Mead, J. F., Alexander, G. V. & Ross, J. F. Replacement rates for human tissue from atmospheric radiocarbon. Science 146, 1170–1172 (1964).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

33.
Thompson, C. & Ballou, E. Studies of metabolic turnover with tritium as a tracer IV. Metabolically inert lipide and protein fractions from the rat. J. Biol. Chem. 200, 731–743 (1953).
CAS  PubMed  PubMed Central  Google Scholar 

34.
Leeson, S. & Walsh, T. Feathering in commercial poultry II. Factors influencing feather growth and feather loss. Worlds Poult. Sci. J. World Poult. Sci J 60, 52–63 (2004).
Article  Google Scholar 

35.
De la Hera, I., Pérez-Tris, J. & Tellería, J. Migratory behaviour affects the trade-off between feather growth rate and feather quality in a passerine bird. Biol. J. Linn. Soc. 97, 98–105 (2009).
Article  Google Scholar 

36.
Bacon Wood, H. Growth bars in feathers. 3016 North Second St Harrisburd Pa. 67, (1949).

37.
Burke, W. H. et al. Variation of seawater 87Sr/86Sr throughout Phanerozoic time. Geology 10, 516–519 (1982).
ADS  CAS  Article  Google Scholar 

38.
Prokoph, A., Shields, G. A. & Veizer, J. Compilation and time-series analysis of a marine carbonate δ18O, δ13C, 87Sr/86Sr and δ34S database through Earth history. Earth-Sci. Rev. 87, 113–133 (2008).
ADS  CAS  Article  Google Scholar 

39.
Touzeau, A. et al. Egyptian mummies record increasing aridity in the Nile valley from 5500 to 1500yr before present. Earth Planet. Sci. Lett. 375, 92–100 (2013).
ADS  CAS  Article  Google Scholar 

40.
Touzeau, A. et al. Diet of ancient Egyptians inferred from stable isotope systematics. J. Archaeol. Sci. 46, 114–124 (2014).
CAS  Article  Google Scholar 

41.
Grimes, V. & Pellegrini, M. A comparison of pretreatment methods for the analysis of phosphate oxygen isotope ratios in bioapatite. Rapid Commun. Mass Spectrom. 27, 375–390 (2013).
CAS  PubMed  Article  PubMed Central  Google Scholar 

42.
Lécuyer, C. Oxygen isotope analysis of phosphate. In Handbook of Stable Isotope Analytical Techniques 482–496 (Elsevier, 2004).

43.
Lécuyer, C., Grandjean, P., O’Neil, J. R., Cappetta, H. & Martineau, F. Thermal excursions in the ocean at the Cretaceous—Tertiary boundary (northern Morocco): δ18O record of phosphatic fish debris. Palaeogeogr. Palaeoclimatol. Palaeoecol. 105, 235–243 (1993).
Article  Google Scholar 

44.
Fourel, F. et al. Carbon and oxygen isotope variability among foraminifera and ostracod carbonated shells. Ann. Univ. Mariae Curie-Sklodowska Sect. Phys. 70, 133 (2015).
Google Scholar 

45.
Coplen, T. B. et al. After two decades a second anchor for the VPDB δ13C scale. Rapid Commun. Mass Spectrom. 20, 3165–3166 (2006).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

46.
Friedman, I., O’neil, J. & Cebula, G. Two new carbonate stable-isotope standards. Geostand. Newsl. 6, 11–12 (1982).
CAS  Article  Google Scholar 

47.
Hut, G. Consultants’ group meeting on stable isotope reference samples for geochemical and hydrological investigations. Consult. Group Meet. Stable Isot. Ref. Samples Geochem. Hydrol. Investig. (1987).

48.
Stichler, W. Interlaboratory comparison of new materials for carbon and oxygen isotope ratio measurements. Int. At. Energy Agency Vienna 825, 67–74 (1993).
Google Scholar 

49.
Bondetti, M., Porcier, S., Ménager, M. & Vieillescazes, C. Analyse chimique de la composition de baumes provenant de momies animales égyptiennes. In Creatures of Earth, Water, and Sky Essays on Animals in Ancient Egypt and Nubia (ed. Pasquali, S.) (SideStone Press, Leiden, 2019).
Google Scholar 

50.
De Muynck, D., Huelga-Suarez, G., Van-Heghe, L., Degryse, P. & Vanhaecke, F. Systematic evaluation of a strontium-specific extraction chromatographic resin for obtaining a purified Sr fraction with quantitative recovery from complex and Ca-rich matrices. J. Anal. At. Spectrom. 24, 1498–1510 (2009).
Article  CAS  Google Scholar 

51.
Richardin, P., Coudert, M., Gandolfo, N. & Vincent, J. Radiocarbon dating of mummified human remains: application to a series of coptic mummies from the Louvre Museum. Radiocarbon 55, 345–352 (2013).
CAS  Article  Google Scholar 

52.
Richardin, P., Perraud, A., Hertzog, J., Madrigal, K. & Berthet, D. Radiocarbon dating of a series of Egyptian mummies heads from Confluences Museum. Radiocarbon 59, 609–619 (2017).
CAS  Article  Google Scholar 

53.
Richardin, P., Porcier, S., Ikram, S., Louarn, G. & Berthet, D. Cats, Crocodiles, cattle, and more: initial steps toward establishing a chronology of ancient egyptian animal mummies. Radiocarbon 59, 595–607 (2017).
Article  Google Scholar 

54.
Richardin, P. & Coudert, M. Datation par le carbone 14 et restes humains. Une étude de cas: la momie dorée de Dunkerque. Techné 44, 75–78 (2016).
Google Scholar 

55.
Moreau, C. et al. Research and development of the artemis 14C AMS facility: status report. Radiocarbon 55, 331–337 (2013).
CAS  Article  Google Scholar 

56.
Bronk-Ramsey, C. Analysis of chronological information and radiocarbon calibration : the program OxCal. Archaeol. Comput. Newsl. 41, 11–16 (1994).
Google Scholar 

57.
Dupras, T. L. & Schwarcz, H. P. Strangers in a strange land: stable isotope evidence for human migration in the Dakhleh Oasis Egypt. J. Archaeol. Sci. 28, 1199–1208 (2001).
Article  Google Scholar 

58.
Thompson, A. H., Richards, M. P., Shortland, A. & Zakrzewski, S. R. Isotopic palaeodiet studies of Ancient Egyptian fauna and humans. J. Archaeol. Sci. 32, 451–463 (2005).
Article  Google Scholar 

59.
Copley, M. S. et al. Short- and long-term foraging and foddering strategies of domesticated animals from Qasr Ibrim Egypt. J. Archaeol. Sci. 31, 1273–1286 (2004).
Article  Google Scholar 

60.
Iacumin, P., Bocherens, H., Mariotti, A. & Longinelli, A. An isotopic palaeoenvironmental study of human skeletal remains from the Nile Valley. Palaeogeogr. Palaeoclimatol. Palaeoecol. 126, 15–30 (1996).
Article  Google Scholar 

61.
Macko, S. A. et al. Documenting the diet in ancient human populations through stable isotope analysis of hair. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 354, 65–76 (1999).
CAS  PubMed  PubMed Central  Article  Google Scholar 

62.
Blake, R. E., Oneil, J. R. & Garcia, G. A. Oxygen isotope systematics of biologically mediated reactions of phosphate: I. Microbial degradation of organophosphorus compounds. Geochim. Cosmochim. Acta 61, 4411–4422 (1997).
ADS  CAS  Article  Google Scholar 

63.
Lécuyer, C. et al. Stable isotope composition and rare earth element content of vertebrate remains from the Late Cretaceous of northern Spain (Laño): did the environmental record survive?. Palaeogeogr. Palaeoclimatol. Palaeoecol. 193, 457–471 (2003).
Article  Google Scholar 

64.
Trueman, C., Chenery, C., Eberth, D. A. & Spiro, B. Diagenetic effects on the oxygen isotope composition of bones of dinosaurs and other vertebrates recovered from terrestrial and marine sediments. J. Geol. Soc. 160, 895–901 (2003).
ADS  CAS  Article  Google Scholar 

65.
Zazzo, A., Lécuyer, C., Sheppard, S. M. F., Grandjean, P. & Mariotti, A. Diagenesis and the reconstruction of paleoenvironments: a method to restore original δ18O values of carbonate and phosphate from fossil tooth enamel. Geochim. Cosmochim. Acta 68, 2245–2258 (2004).
ADS  CAS  Article  Google Scholar 

66.
Zazzo, A., Lécuyer, C. & Mariotti, A. Experimentally-controlled carbon and oxygen isotope exchange between bioapatites and water under inorganic and microbially-mediated conditions. Geochim. Cosmochim. Acta 68, 1–12 (2004).
ADS  CAS  Article  Google Scholar 

67.
Stanton Thomas, K. J. & Carlson, S. J. Microscale δ18O and δ13C isotopic analysis of an ontogenetic series of the hadrosaurid dinosaur Edmontosaurus: implications for physiology and ecology. Palaeogeogr. Palaeoclimatol. Palaeoecol. 206, 257–287 (2004).
Article  Google Scholar 

68.
Kohn, M. J., Morris, J. & Olin, P. Trace element concentrations in teeth—a modern Idaho baseline with implications for archeometry, forensics, and palaeontology. J. Archaeol. Sci. 40, 1689–1699 (2013).
CAS  Article  Google Scholar 

69.
Balter, V. & Lécuyer, C. Determination of Sr and Ba partition coefficients between apatite and water from 5 to 60°C: a potential new thermometer for aquatic paleoenvironments. Geochim. Cosmochim. Acta 68, 423–432 (2004).
ADS  CAS  Article  Google Scholar 

70.
Buzon, M. R., Simonetti, A. & Creaser, R. A. Migration in the Nile Valley during the New Kingdom period: a preliminary strontium isotope study. J. Archaeol. Sci. 34, 1391–1401 (2007).
Article  Google Scholar 

71.
Meyburg, B. U., Kirwan, G. M. & Garcia, E. F. J. Greater Spotted Eagle (Clanga clanga). In Handbook of the Birds of the World Alive (ed. de Juana, E.) (Lynx Edicions, Barcelona, 2019).
Google Scholar 

72.
Orta, J., Boesman, P., Kirwan, G. M. & Marks, J. S. Long-legged Buzzard (Buteo rufinus). In Handbook of the Birds of the World Alive (ed. de Juana, E.) (Lynx Edicions, Barcelona, 2019).
Google Scholar 

73.
Svensson, L. & Madge, S. (eds) Handbook of the Middle East, and North Africa: The Birds of the Western Palearctic (Oxford University Press, Oxford, 1977).
Google Scholar 

74.
Hancock, J., Kushlan, J. A. & Kahl, M. P. Storks (Ibises and Spoonbills of the World. Academic Press, Cambridge, 2010).
Google Scholar 

75.
Matheu, E., del Hoyo, J., Christie, D. A., Kirwan, G. M. & Garcia, E. F. J. African Sacred Ibis (Threskiornis aethiopicus). In Handbook of the Birds of the World Alive (ed. de Juana, E.) (Lynx Edicions, Barcelona, 2019).
Google Scholar 

76.
Audouin, V. Explication sommaire des planches d’oiseaux de l’Egypte et de la Syrie, publiées par J.-C. Savigny, membre de l’Institut, offrant un exposé des caractères naturels des genres avec la distinction des espèces. in Description de l’Egypte, ou Recueil des observations et des recherches qui ont été faites en Egypte pendant l’expédition de l’armée française, publié par les ordres de sa Majesté–L’Empereur Napoléon le Grand vol. 1 251–318 (1809).

77.
Savigny, J.-C. Histoire Naturelle et Mythologique de l’Ibis. (1805).

78.
Lortet, C. E. & Gaillard, C. L. faune momifiée de l’Ancienne Egypte. Archives du Muséum d’histoire naturelle de Lyon 8, 1–205 (1903).
Google Scholar 

79.
De Margerie, E. Fonction biomécanique des microstructures osseuses chez les oiseaux / Biomechanical function of bone microstructure in birds. C.R. Palevol 5, 619–628 (2006).
Article  Google Scholar 

80.
Kerley, E. R. The microscopic determination of age in human bone. Am. J. Phys. Anthropol. 23, 149–163 (1965).
CAS  PubMed  Article  PubMed Central  Google Scholar 

81.
Peek, S. & Clementz, M. T. Sr/Ca and Ba/Ca variations in environmental and biological sources: a survey of marine and terrestrial systems. Geochim. Cosmochim. Acta 95, 36–52 (2012).
ADS  CAS  Article  Google Scholar 

82.
Zuberogoitia, I., Zabala, J. & Martínez, J. E. Moult in birds of prey: a review of current knowledge and future challenges for research. Ardeola 65, 183–207 (2018).
Article  Google Scholar 

83.
Kovacs, G. Occurrence of the Long-legged buzzard (Buteo rufinus) in the Hortobagy between 1976 and 1991. Aquila 99, 41–48 (1992).
Google Scholar 

84.
Bloom, P. H. & Clark, W. S. Molt and sequence of plumages of Golden eagles and a technique for in-hand ageing. North Am. Bird Bander 26, 97–116 (2001).
Google Scholar 

85.
Lowe, K. W., Clark, A. & Clark, R. A. Body measurements, plumage and moult of the Sacred ibis in south Africa. Ostrich 56, 111–116 (1985).
Article  Google Scholar 

86.
Van Neer, W. Evolution of prehistoric fishing in the Nile Valley. J. Afr. Archaeol. 2, 251–269 (2004).
Article  Google Scholar 

87.
Porcier, S. et al. Wild crocodiles hunted to make mummies in Roman Egypt: Evidence from synchrotron imaging. J. Archaeol. Sci. 110, 105009 (2019).
Article  Google Scholar 

88.
Kjellén, N. Moult in relation to migration in birds-a review. Ornis Svec. 4, 1–24 (1994).
Google Scholar 

89.
Iacumin, P., Bocherens, H., Mariotti, A. & Longinelli, A. Oxygen isotope analyses of co-existing carbonate and phosphate in biogenic apatite: a way to monitor diagenetic alteration of bone phosphate?. Earth Planet. Sci. Lett. 142, 1–6 (1996).
ADS  CAS  Article  Google Scholar  More

This study of recreation in two geographically distinct regions of the United States indicates that social media can predict visitation at recreation sites on public land. We conclude that social media data can be applied with moderate success to estimate visitation at sites that are unmonitored or otherwise lack on-site counts, even in new regions. A basic visitation model that relies solely on generic predictors (e.g., weather, holidays, and seasonality) is only modestly successful due to regional differences in visitor behavior (Model 1). Performance is improved by including relationships between visitation and social media (Flickr, Instagram, and Twitter), even when these relationships are transferred from a different region (WWA, Model 2). These results are consistent with prior research findings that social media counts are correlated with on-site visitor counts from public lands8,9,24,32, and extend earlier findings by showing the potential for statistical models to estimate absolute numbers of visitors at unmonitored sites with parameters derived from social media. This is evidence of patterns in how visitors use and share social media. Furthermore, it suggests that measurable variables associated with social media use could support transferable models for accurately estimating visitation to public lands across large geographies.
We expected to detect regional differences in social media use given regional distinctions in climate, land management, population density, mobile phone signal coverage, and the demographics of visitors. Contrary to our expectations, we observed that the rate of posting to social media about recreation visits is similar across sites in NNM and WWA, and both regions displayed positive correlations between each of these social media data sources and observed visitation (Fig. 3). Furthermore, a model parameterized with social media use in WWA (Model 2) explains 45% of the variation in visitation across all 13 sites in NNM and 79% of the variability in visitation at the subset of sites that had social media posts. Re-parameterizing Model 2 with a portion of NNM visitation data does not improve its performance (Model 3) until these data are considered at the site level (Model 5). There is a noisy but consistent relationship between a destination’s popularity with visitors and its popularity on social media, regardless of whether the site is in NNM or WWA. We interpret this as evidence that visitors are equally likely to share their recreation experiences in NNM or WWA, despite the regional differences in the types of recreation opportunities and the people who are visiting.
Although there are consistent relationships between social media use and visitation at a regional scale, we see large site-to-site variability in how visitors use social media within both regions. Beyond simple differences in numbers of posts by site, the proportion of people who post about their visits varies by destination, ranging from 7% of visitors to Kasha–Katuwe Tent Rocks National Monument posting on Instagram to zero recorded social media user-days at several trails in the Valles Caldera National Preserve. As a result, models with random effects that allow sites to have unique relationships between social media, weather, and visitation (Models 4–5) perform substantially better than models that assume social media use is consistently related to visitation for all types of sites (Models 1–3). This is especially true at the most and least visited sites (Figs. 4, 5), where visitors may be sharing social media differently than they do at moderately visited sites, and responding differently to other conditions such as weather or holidays. These results indicate that while data on social media use are helpful for predicting visitation with moderate certainty in an otherwise unknown region (Model 1 vs. 2), their utility for estimating visitation is less clear when local data on the effects of environmental and institutional conditions such as weather and holidays are available to parameterize site-specific models.
A primary goal of this study is to test approaches for estimating visitation over relatively small areas in order to explore the limits of the data and methods. We find that six of the 13 sites in NNM—representing individual trails or groups of trails within a larger park—lack social media during the study period. Model 2, which depends on social media data to estimate visitation at unmonitored sites, consequently performs relatively poorly at these six sites (Fig. 4). Generally, sites that have sparse social media data tend to receive few visitors, but there are exceptions. Alcove House in Bandelier National Monument, for instance, is a very highly visited site that lacks Instagram images in our study because the site does not appear as a prescribed location for Instagram users. Our informal observation is that many visitors instead share untagged photographs of Alcove House or assign their images to other relevant place names such as Bandelier National Monument. Clearly, there are thresholds to where and how social media can be leveraged for visitor estimation. Our research suggests that future studies and visitation models could be improved by accounting for the popularity of sites on social media in the study design12. Visitation models that include predictors derived from social media (Models 2, 3, and 5, here) will likely out-perform alternative models for estimating use at popular sites or when longer time series are available. At locations with low or no posting activity, where social media contributes less to visitor estimates, it could be more useful to collect on-site data such as vehicle and pedestrian counts. Further research is necessary to understand what combination of on-site, social media, environmental, and other data is most valuable at different spatial and temporal scales.
These observations suggest that variability in correlations between social media and on-site visitor counts seen here and in previous studies8,12,20 is derived from local factors influencing visitors’ day-to-day decisions about whether and how to share a destination on social media. Choices about whether to post to social media are likely influenced by the characteristics of the local site—perhaps its topography, amenities, predominant activity, or unique natural features—and the ways that people relate to these features33. The characteristics of the visitors and the relative contributions of the natural versus the social experience in the motivation for the trip may play a role. For example, if a given type of site attracts visitors wishing to “unplug” and have a nature-based outing (a calm forest glade, say), there may be fewer posts per visit than for a site attracting visitors who desire a social experience within a natural setting (a famous scenic overlook, for instance). Another possibility is that the prevailing popularity of certain destinations on social media creates a positive feedback, whereby new visitors feel compelled to share content about their visit in response to the posts of others or the local hashtags that may make it easier or more enticing to post. Variability could also arise from the recent trend towards discouraging visitors from posting geolocated content and attracting attention to less popular or back-country sites that are not equipped to sustain higher use, although this is probably of minor importance, currently.
This is the first study to our knowledge that develops and tests models for estimating absolute numbers of visitors at unmonitored recreation sites or times using multiple social media data sources with differential effects. Building on earlier research exploring relationships of park visitation with numbers of posts to multiple social media platforms9,12,20,27, the present study tests whether models with a mixture of predictors to represent varying effects of three online platforms can estimate visitation in novel situations. We find that each social media data source contributes information that explains a statistically significant portion of the variability in visitation and improves the accuracy of the estimate. This is the case not only for Instagram, which captures 3–4% of visitor-days at our research sites in WWA and NNM, but also for Flickr and Twitter, with relatively small amounts of content shared ( More

1.
Arctic Monitoring and Assessment Programme. Snow, Water, Ice and Permafrost in the Arctic (SWIPA) 2017 (Arctic Monitoring and Assessment Programme (AMAP), 2017).
2.
Chapin, F. S. 3rd et al. Role of land-surface changes in arctic summer warming. Science 310, 657–660 (2005).
ADS  CAS  PubMed  Google Scholar 

3.
Tape, K. D., Christie, K., Carroll, G. & O’donnell, J. A. Novel wildlife in the Arctic: the influence of changing riparian ecosystems and shrub habitat expansion on snowshoe hares. Glob. Change Biol. 22, 208–219 (2016).
ADS  Google Scholar 

4.
Downing, A. & Cuerrier, A. A synthesis of the impacts of climate change on the First Nations and Inuit of Canada. Indian J. Tradit. Knowl. 10, 57–70 (2011).
Google Scholar 

5.
National Academies of Sciences. Understanding Northern Latitude Vegetation Greening and Browning: Proceedings of a Workshop (The National Academies Press, 2019).

6.
Bjorkman, A. D. et al. Plant functional trait change across a warming tundra biome. Nature 562, 57–62 (2018).
ADS  CAS  PubMed  Google Scholar 

7.
Elmendorf, S. C. et al. Plot-scale evidence of tundra vegetation change and links to recent summer warming. Nat. Clim. Change 2, 453–457 (2012).
ADS  Google Scholar 

8.
Gauthier, G. et al. Long-term monitoring at multiple trophic levels suggests heterogeneity in responses to climate change in the Canadian Arctic tundra. Philos. Trans. R. Soc. Ser. B 368, 20120482 (2013).
Google Scholar 

9.
Myers-Smith, I. H. et al. Eighteen years of ecological monitoring reveals multiple lines of evidence for tundra vegetation change. Ecol. Monogr. 89, e01351 (2019).
Google Scholar 

10.
Tape, K. D., Hallinger, M., Welker, J. M. & Ruess, R. W. Landscape heterogeneity of shrub expansion in Arctic Alaska. Ecosystems 15, 711–724 (2012).
CAS  Google Scholar 

11.
Pattison, R. R., Jorgenson, J. C., Raynolds, M. K. & Welker, J. M. Trends in NDVI and Tundra Community Composition in the Arctic of NE Alaska Between 1984 and 2009. Ecosystems 18, 707–719 (2015).
Google Scholar 

12.
Gamm, C. M. et al. Declining growth of deciduous shrubs in the warming climate of continental western Greenland. J. Ecol. 106, 640–654 (2018).
CAS  Google Scholar 

13.
Forchhammer M. Sea-ice induced growth decline in Arctic shrubs. Biol. Lett. 13, 20170122 (2017).

14.
Street, L., Shaver, G., Williams, M. & Van Wijk, M. What is the relationship between changes in canopy leaf area and changes in photosynthetic CO2 flux in arctic ecosystems? J. Ecol. 95, 139–150 (2007).
Google Scholar 

15.
Raynolds, M. K., Walker, D. A., Epstein, H. E., Pinzon, J. E. & Tucker, C. J. A new estimate of tundra-biome phytomass from trans-Arctic field data and AVHRR NDVI. Remote Sens. Lett. 3, 403–411 (2012).
Google Scholar 

16.
Berner, L. T., Jantz, P., Tape, K. D. & Goetz, S. J. Tundra plant aboveground biomass and shrub dominance mapped across the North Slope of Alaska. Environ. Res. Lett. 13, 035002 (2018).
ADS  Google Scholar 

17.
Bhatt, U. S. et al. Changing seasonality of panarctic tundra vegetation in relationship to climatic variables. Environ. Res. Lett. 12, 1–18 (2017).
Google Scholar 

18.
Guay, K. C. et al. Vegetation productivity patterns at high northern latitudes: a multi-sensor satellite data assessment. Glob. Change Biol. 20, 3147–3158 (2014).
ADS  Google Scholar 

19.
Pinzon, J. & Tucker, C. A non-stationary 1981–2012 AVHRR NDVI3g time series. Remote Sens. 6, 6929–6960 (2014).
ADS  Google Scholar 

20.
Ju, J. & Masek, J. G. The vegetation greenness trend in Canada and US Alaska from 1984–2012 Landsat data. Remote Sens. Environ. 176, 1–16 (2016).
ADS  Google Scholar 

21.
Karlsen, S. R., Anderson, H. B., Van der Wal, R. & Hansen, B. B. A new NDVI measure that overcomes data sparsity in cloud-covered regions predicts annual variation in ground-based estimates of high arctic plant productivity. Environ. Res. Lett. 13, 025011 (2018).
ADS  Google Scholar 

22.
McManus, kM. et al. Satellite-based evidence for shrub and graminoid tundra expansion in northern Quebec from 1986 to 2010. Glob. Change Biol. 18, 2313–2323 (2012).
ADS  Google Scholar 

23.
Frost, G. V., Epstein, H. & Walker, D. Regional and landscape-scale variability of Landsat-observed vegetation dynamics in northwest Siberian tundra. Environ. Res. Lett. 9, 025004 (2014).
ADS  Google Scholar 

24.
Raynolds, M. K. & Walker, D. A. Increased wetness confounds Landsat-derived NDVI trends in the central Alaska North Slope region, 1985–2011. Environ. Res. Lett. 11, 085004 (2016).
ADS  Google Scholar 

25.
Gorelick, N. et al. Google Earth Engine: planetary-scale geospatial analysis for everyone. Remote Sens. Environ. 202, 18–27 (2017).
ADS  Google Scholar 

26.
Zhu, Z., Wang, S. & Woodcock, C. E. Improvement and expansion of the Fmask algorithm: cloud, cloud shadow, and snow detection for Landsats 4–7, 8, and Sentinel 2 images. Remote Sens. Environ. 159, 269–277 (2015).
ADS  Google Scholar 

27.
Pastick, N. J. et al. Spatiotemporal remote sensing of ecosystem change and causation across Alaska. Glob. Change Biol. 25, 1171–1189 (2019).
ADS  Google Scholar 

28.
Walker, D. et al. Phytomass, LAI, and NDVI in northern Alaska: relationships to summer warmth, soil pH, plant functional types, and extrapolation to the circumpolar Arctic. J. Geophys. Res. 108, 8169 (2003).
Google Scholar 

29.
Lucht, W. et al. Climatic control of the high-latitude vegetation greening trend and Pinatubo effect. Science 296, 1687–1689 (2002).
ADS  CAS  PubMed  Google Scholar 

30.
Fraser, R. H., Lantz, T. C., Olthof, I., Kokelj, S. V. & Sims, R. A. Warming-induced shrub expansion and lichen decline in the Western Canadian. Arct. Ecosyst. 17, 1151–1168 (2014).
Google Scholar 

31.
Bonney, M. T., Danby, R. K. & Treitz, P. M. Landscape variability of vegetation change across the forest to tundra transition of central Canada. Remote Sens. Environ. 217, 18–29 (2018).
ADS  Google Scholar 

32.
Cuerrier, A., Brunet, N. D., Gérin-Lajoie, J., Downing, A. & Lévesque, E. The study of Inuit knowledge of climate change in Nunavik, Quebec: a mixed methods approach. Hum. Ecol. 43, 379–394 (2015).
Google Scholar 

33.
Forbes, B. C. & Stammler, F. Arctic climate change discourse: the contrasting politics of research agendas in the West and Russia. Polar Res. 28, 28–42 (2009).
Google Scholar 

34.
Forbes, B. C., Fauria, M. M. & Zetterberg, P. Russian Arctic warming and ‘greening’ are closely tracked by tundra shrub willows. Glob. Change Biol. 16, 1542–1554 (2010).
ADS  Google Scholar 

35.
Tape, K., Sturm, M. & Racine, C. The evidence for shrub expansion in Northern Alaska and the Pan-Arctic. Glob. Change Biol. 12, 686–702 (2006).
ADS  Google Scholar 

36.
Ropars, P. & Boudreau, S. Shrub expansion at the forest–tundra ecotone: spatial heterogeneity linked to local topography. Environ. Res. Lett. 7, 015501 (2012).
ADS  Google Scholar 

37.
Myers-Smith, I. H. et al. Complexity revealed in the greening of the Arctic. Nat. Clim. Change 10, 106–117 (2020).
ADS  Google Scholar 

38.
Park, T. et al. Changes in growing season duration and productivity of northern vegetation inferred from long-term remote sensing data. Environ. Res. Lett. 11, 084001 (2016).
ADS  Google Scholar 

39.
Riihimäki, H., Heiskanen, J. & Luoto, M. The effect of topography on arctic-alpine aboveground biomass and NDVI patterns. Int. J. Appl. Earth Obs. Geoinf. 56, 44–53 (2017).
ADS  Google Scholar 

40.
Fraser, R. H., Olthof, I., Lantz, T. C. & Schmitt, C. UAV photogrammetry for mapping vegetation in the low-Arctic. Arct. Sci. 2, 79–102 (2016).
Google Scholar 

41.
Berner, L. T., Beck, P. S. A., Bunn, A. G. & Goetz, S. J. Plant response to climate change along the forest-tundra ecotone in northeastern Siberia. Glob. Change Biol. 19, 3449–3462 (2013).
Google Scholar 

42.
Myers-Smith, I. H. et al. Climate sensitivity of shrub growth across the tundra biome. Nat. Clim. Change 5, 887–891 (2015).
ADS  Google Scholar 

43.
Bjorkman, A. D., Vellend, M., Frei, E. R. & Henry, G. H. Climate adaptation is not enough: warming does not facilitate success of southern tundra plant populations in the high Arctic. Glob. Change Biol. 23, 1540–1551 (2017).
ADS  Google Scholar 

44.
Post, E. & Pedersen, C. Opposing plant community responses to warming with and without herbivores. Proc. Natl Acad. Sci. USA 105, 12353–12358 (2008).
ADS  CAS  PubMed  Google Scholar 

45.
Yu, Q., Epstein, H., Engstrom, R. & Walker, D. Circumpolar arctic tundra biomass and productivity dynamics in response to projected climate change and herbivory. Glob. Change Biol. 23, 3895–3907 (2017).
ADS  Google Scholar 

46.
Liljedahl, A. K. et al. Pan-Arctic ice-wedge degradation in warming permafrost and its influence on tundra hydrology. Nat. Geosci. 9, 312–318 (2016).
ADS  CAS  Google Scholar 

47.
Perreault, N., Levesque, E., Fortier, D. & Lamarque, L. J. Thermo-erosion gullies boost the transition from wet to mesic tundra vegetation. Biogeosciences 13, 1237–1253 (2016).
ADS  Google Scholar 

48.
Grant, R. F., Mekonnen, Z. A., Riley, W. J., Arora, B. & Torn, M. S. Mathematical modelling of Arctic Polygonal Tundra with Ecosys: 2. Microtopography determines how CO2 and CH4 exchange responds to changes in temperature and precipitation. J. Geophys. Res. 122, 3174–3187 (2017).
CAS  Google Scholar 

49.
Phoenix, G. K. & Bjerke, J. W. Arctic browning: extreme events and trends reversing arctic greening. Glob. Change Biol. 22, 2960–2962 (2016).
ADS  Google Scholar 

50.
Treharne, R., Bjerke, J. W., Tømmervik, H., Stendardi, L. & Phoenix, G. K. Arctic browning: Impacts of extreme climatic events on heathland ecosystem CO2 fluxes. Glob. Change Biol. 25, 489–503 (2018).
ADS  Google Scholar 

51.
Forbes, B. C. et al. High resilience in the Yamal-Nenets social–ecological system, west Siberian Arctic, Russia. Proc. Natl Acad. Sci. USA 106, 22041–22048 (2009).
ADS  CAS  PubMed  Google Scholar 

52.
Mekonnen, Z. A., Riley, W. J. & Grant, R. F. Accelerated nutrient cycling and increased light competition will lead to 21st century shrub expansion in North American Arctic tundra. J. Geophys. Res. 123, 1683–1701 (2018).
CAS  Google Scholar 

53.
Rocha, A. V. et al. The footprint of Alaskan tundra fires during the past half-century: implications for surface properties and radiative forcing. Environ. Res. Lett. 7, 044039 (2012).
ADS  Google Scholar 

54.
Giglio, L., Boschetti, L., Roy, D. P., Humber, M. L. & Justice, C. O. The Collection 6 MODIS burned area mapping algorithm and product. Remote Sens. Environ. 217, 72–85 (2018).
ADS  PubMed  PubMed Central  Google Scholar 

55.
Hu, F. S. et al. Arctic tundra fires: natural variability and responses to climate change. Front. Ecol. Environ. 13, 369–377 (2015).
Google Scholar 

56.
Mack, M. C. et al. Carbon loss from an unprecedented Arctic tundra wildfire. Nature 475, 489–492 (2011).
ADS  CAS  PubMed  Google Scholar 

57.
Jones, B. M. et al. Identification of unrecognized tundra fire events on the north slope of Alaska. J. Geophys. Res. 118, 1334–1344 (2013).
Google Scholar 

58.
Loranty, M. M. et al. Siberian tundra ecosystem vegetation and carbon stocks four decades after wildfire. J. Geophys. Res. 119, 2144–2154 (2014).
CAS  Google Scholar 

59.
Natali, S. M. et al. Large loss of CO2 in winter observed across the northern permafrost region. Nat. Clim. Change 9, 852–857 (2019).
ADS  CAS  Google Scholar 

60.
Schuur, E. et al. Climate change and the permafrost carbon feedback. Nature 520, 171–179 (2015).
ADS  CAS  PubMed  Google Scholar 

61.
Pearson, R. G. et al. Shifts in Arctic vegetation and associated feedbacks under climate change. Nat. Clim. Change 3, 673–677 (2013).
ADS  Google Scholar 

62.
Loranty, M. M., Goetz, S. J. & Beck, P. S. A. Tundra vegetation effects on pan-Arctic albedo. Environ. Res. Lett. 6, 024014 (2011).
ADS  Google Scholar 

63.
Loranty, M. M. et al. Reviews and syntheses: changing ecosystem influences on soil thermal regimes in northern high-latitude permafrost regions. Biogeosciences 15, 5287–5313 (2018).
ADS  CAS  Google Scholar 

64.
Tape, K. D., Gustine, D. D., Ruess, R. W., Adams, L. G. & Clark, J. A. Range expansion of moose in Arctic Alaska linked to warming and increased shrub habitat. PLoS ONE 11, e0152636 (2016).
PubMed  PubMed Central  Google Scholar 

65.
Tape, K. D., Jones, B. M., Arp, C. D., Nitze, I. & Grosse, G. Tundra be dammed: beaver colonization of the Arctic. Glob. Change Biol. 24, 4478–4488 (2018).
ADS  Google Scholar 

66.
Joly, K., Jandt, R. R. & Klein, D. R. Decrease of lichens in Arctic ecosystems: the role of wildfire, caribou, reindeer, competition and climate in north‐western Alaska. Polar Res. 28, 433–442 (2009).
Google Scholar 

67.
Macias-Fauria, M., Forbes, B. C., Zetterberg, P. & Kumpula, T. Eurasian Arctic greening reveals teleconnections and the potential for structurally novel ecosystems. Nat. Clim. Change 2, 613–618 (2012).
ADS  Google Scholar 

68.
Wesche, S. D. & Chan, H. M. Adapting to the impacts of climate change on food security among Inuit in the Western Canadian Arctic. EcoHealth 7, 361–373 (2010).
PubMed  Google Scholar 

69.
Kuhnlein, H. V. & Chan, H. M. Environment and contaminants in traditional food systems of northern indigenous peoples. Annu. Rev. Nutr. 20, 595–626 (2000).
CAS  PubMed  Google Scholar 

70.
Tucker, C. J. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sens. Environ. 8, 127–150 (1979).
ADS  Google Scholar 

71.
Virtanen, R. et al. Where do the treeless tundra areas of northern highlands fit in the global biome system: toward an ecologically natural subdivision of the tundra biome. Ecol. Evol. 6, 143–158 (2016).
PubMed  Google Scholar 

72.
Masek, J. G. et al. A Landsat surface reflectance dataset for North America, 1990-2000. IEEE Geosci. Remote Sens. Lett. 3, 68–72 (2006).
ADS  Google Scholar 

73.
Vermote, E., Justice, C., Claverie, M. & Franch, B. Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product. Remote Sens. Environ. 185, 46–56 (2016).
ADS  PubMed  Google Scholar 

74.
Python Software Foundation. Python Language Software Version 3.7.3. https://www.python.org/ (2020).

75.
Foga, S. et al. Cloud detection algorithm comparison and validation for operational Landsat data products. Remote Sens. Environ. 194, 379–390 (2017).
ADS  Google Scholar 

76.
Roy, D. P. et al. Characterization of Landsat-7 to Landsat-8 reflective wavelength and normalized difference vegetation index continuity. Remote Sens. Environ. 185, 57–70 (2016).
ADS  PubMed  Google Scholar 

77.
Sulla-Menashe, D., Friedl, M. A. & Woodcock, C. E. Sources of bias and variability in long-term Landsat time series over Canadian boreal forests. Remote Sens. Environ. 177, 206–219 (2016).
ADS  Google Scholar 

78.
Liaw, A. & Wiener, M. Classification and Regression by randomForest. R News 2, 18–22 (2002).
Google Scholar 

79.
Wright, M. N. & Ziegler, A. Ranger: a fast implementation of random forests for high dimensional data in C++ and R. J. Stat. Softw. 77, 1–17 (2017).
Google Scholar 

80.
Melaas, E. K. et al. Multisite analysis of land surface phenology in North American temperate and boreal deciduous forests from Landsat. Remote Sens. Environ. 186, 452–464 (2016).
ADS  Google Scholar 

81.
Markham, B. L. & Helder, D. L. Forty-year calibrated record of earth-reflected radiance from Landsat: a review. Remote Sens. Environ. 122, 30–40 (2012).
ADS  Google Scholar 

82.
Markham, B. et al. Landsat-8 operational land imager radiometric calibration and stability. Remote Sens. 6, 12275–12308 (2014).
ADS  Google Scholar 

83.
Kendall, M. G. Rank Correlation Methods 4th edn (Charles Griffin, 1975).

84.
Sen, P. K. Estimates of the regression coefficient based on Kendall’s tau. J. Am. Stat. Assoc. 63, 1379–1389 (1968).
MathSciNet  MATH  Google Scholar 

85.
Bronaugh, D. & Werner, A. zyp: Zhang + Yue-Pilon Trends Package. R Package Version 0.10-1.1. https://CRAN.R-project.org/package=zyp (2012).

86.
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2020).

87.
Rohde, R. et al. A new estimate of the average Earth surface land temperature spanning 1753 to 2011. Geoinform. Geostat. 7, https://doi.org/10.4172/2327-4581.1000101 (2013).

88.
Hansen, J., Ruedy, R., Sato, M. & Lo, K. Global surface temperature change. Rev. Geophys. 48, RG4004 (2010).
ADS  Google Scholar 

89.
Cowtan, K. & Way, R. G. Coverage bias in the HadCRUT4 temperature series and its impact on recent temperature trends. Q. J. R. Meteorol. Soc. 140, 1935–1944 (2014).
ADS  Google Scholar 

90.
Harris, I., Jones, P. D., Osborn, T. J. & Lister, D. H. Updated high-resolution grids of monthly climatic observations—the CRU TS3.10 Dataset. Int. J. Climatol. 34, 623–642 (2014).
Google Scholar 

91.
Willmott, C. J. & Matsuura, K. Terrestrial Air Temperature and Precipitation: Monthly Time Series (1900–2017) v. 5.01. http://climate.geog.udel.edu/~climate (University of Deleware, 2018).

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

93.
Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A. & Hegewisch, K. C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Sci. Data 5, 170191 (2018).
PubMed  PubMed Central  Google Scholar 

94.
Obu, J. et al. ESA Permafrost Climate Change Initiative (Permafrost_cci): Permafrost Extent for the Northern Hemisphere, v1.0. https://doi.org/10.5285/c7590fe40d8e44169d511c70a60ccbcc (Centre for Environmental Data Analysis, 2019).

95.
Obu, J. et al. ESA Permafrost Climate Change Initiative (Permafrost_cci): Permafrost Ground Temperature for the Northern Hemisphere, v1.0. https://doi.org/10.5285/c7590fe40d8e44169d511c70a60ccbcc (Centre for Environmental Data Analysis, 2019).

96.
Obu, J. et al. ESA Permafrost Climate Change Initiative (Permafrost_cci): Permafrost Active Layer Thickness for the Northern Hemisphere, v1.0. https://doi.org/10.5285/1ee56c42cf6c4ef698693e00a63795f4 (Centre for Environmental Data Analysis, 2019).

97.
Olefeldt, D. et al. Arctic Circumpolar Distribution and Soil Carbon of Thermokarst Landscapes. https://doi.org/10.3334/ORNLDAAC/1332 (ORNL DAAC, 2015).

98.
Defourny, P. et al. Land Cover Climate Change Initiative—Product User Guide Version v2. http://maps.elie.ucl.ac.be/CCI/viewer/download/ESACCI-LC-Ph2-PUGv2_2.0.pdf (European Space Agency, 2017).

99.
Rizzoli, P. et al. Generation and performance assessment of the global TanDEM-X digital elevation model. ISPRS J. Photogramm. Remote Sens. 132, 119–139 (2017).
ADS  Google Scholar 

100.
Kuhn, M. Building predictive models in R using the caret package. J. Stat. Softw. 28, 1–26 (2008).
Google Scholar 

101.
Greenwell, B. M. pdp: an R package for constructing partial dependence plots. R. J. 9, 421–436 (2017).
Google Scholar 

102.
Le Moullec, M., Buchwal, A., Wal, R., Sandal, L. & Hansen, B. B. Annual ring growth of a widespread high arctic shrub reflects past fluctuations in community-level plant biomass. J. Ecol. 107, 436–451 (2019).
Google Scholar 

103.
Bunn, A. G. A dendrochronology program library in R (dplR). Dendrochronologia 26, 115–124 (2008).
Google Scholar 

104.
Euskirchen, E., Bret-Harte, M. S., Scott, G., Edgar, C. & Shaver G. R. Seasonal patterns of carbon dioxide and water fluxes in three representative tundra ecosystems in northern Alaska. Ecosphere 3, https://doi.org/10.1890/ES1811-00202.00201 (2012).

105.
Euskirchen, E. S. et al. Interannual and seasonal patterns of carbon dioxide, water, and energy fluxes from ecotonal and thermokarst-impacted ecosystems on carbon-rich permafrost soils in Northeastern Siberia. J. Geophys. Res. 122, 2651–2668 (2017).
CAS  Google Scholar 

106.
Baldocchi, D. et al. FLUXNET: a new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor, and energy flux densities. Bull. Am. Meteorol. Soc. 82, 2415–2434 (2001).
ADS  Google Scholar 

107.
Reichstein, M. et al. On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm. Glob. Change Biol. 11, 1424–1439 (2005).
ADS  Google Scholar 

108.
Hijmans, R. J. raster: Geographic Analysis and Modeling. R package version 3.0-12. http://CRAN.R-project.org/package=raster (2019).

109.
Bivand, R., Keitt, T. & Rowlingson B. rgdal: Bindings for the ‘Geospatial’ Data Abstraction Library. R Package Version 1.4-8. https://CRAN.R-project.org/package=rgdal (2019).

110.
Bivand, R. & Lewin-Koh, N. maptools: Tools for Handling Spatial Objects. R Package Version 0.9.9. https://CRAN.R-project.org/package=maptools (2019).

111.
Dawle, M. & Srinivasan, A. data.table: Extension of ‘data.frame’. R Package Version 1.12.8. https://CRAN.R-project.org/package=data.table (2019).

112.
Wickham, H. & Francois, R. dplyr: A Grammar of Data Manipulation. R Package Version 0.8.5. https://CRAN.R-project.org/package=dplyr (2015).

113.
Wickham, H. & Henry, L. tidyr: Tidy Messy Data. R Package Version 1.0.2. https://CRAN.R-project.org/package=tidyr (2020).

114.
Sarkar, D. Lattice: Multivariate Data Visualization with R (Springer, 2008).

115.
Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, New York, 2016).

116.
Kassambara, A. ggpubr: ‘ggplot2’ Basde Publication Ready Plots. R Package Version 0.2.5. https://CRAN.R-project.org/package=ggpubr (2020). More

Threatened species assessments
We assessed the proportion of the proposed budget allocated to RM for a total of 2328 species, independently managed subspecies, or distinct populations (hereafter species): 700 in NZ, 361 in NSW, and 1267 freshwater and terrestrial species in the U.S. In all jurisdictions this included the most threatened listed species and/or those with recovery plans: species with Threatened and Endangered status in the U.S. with active recovery plans as of January 2017, species that met a series of criteria in NSW as of 2013 (e.g., excluding less threatened species that do not require any active intervention and those with a large geographic range17), and the most threatened species in New Zealand as of 2012, which included all species in the Threatened and At-Risk categories with declining populations42. In all three jurisdictions, species are listed for legal protection if they are at risk of extinction. Once listed, recovery planning (including proposed projects, management tasks, and budgets) documents are developed with the objective of securing species from extinction and recovering populations to a point that they can be de-listed. Although our dataset examining threatened species recovery planning is the most comprehensive to date, our data do not represent all spending on species—there are other activities for both management action and RM that occur at a sub-jurisdictional level or outside of government.
Estimating resources allocated to RM vs action
We gathered information on the planned costs of management tasks necessary to achieve recovery for threatened species from previously published recovery planning databases (details provided in refs. 15,16,43,44 and Supplementary Methods). Briefly, for NZ and NSW, a suite of management tasks had been developed during structured expert elicitation workshops, as part of a systematic prioritization exercise16,17. For the U.S., management tasks and their cost had been extracted from each species’ published recovery plans (Supplementary Methods15). These data represent an evolution of the implementation of a systematic and cost-effective approach to endangered species resource allocation (i.e., the Project Prioritization Protocol), beginning with NZ in 200916, and subsequently applied to NSW in 201317 and the U.S. in 201615.
For each proposed management task we used the methods description to categorize tasks as research and monitoring or action based on the definitions in IUCN classification schemes (https://www.iucnredlist.org/resources/classification-schemes, Supplementary Table 145). For NZ and NSW, using previously published datasets we used a combination of the methods description field and 4 other columns that classified the management task methods into increasingly general categories16,17,43. We used keywords such as survey, monitor, surveillance, develop techniques, inventory, research, and develop plan to search for research and monitoring tasks. We reviewed the management tasks identified by these broad search terms to ensure only research and monitoring tasks were included. We also reviewed the management tasks that were not captured by search terms to ensure no research and monitoring tasks were excluded. For the US, the methods descriptions were too complex for keyword searches. Instead, the first author and a trained technician classified each management task manually. To ensure that management tasks were being classified similarly, the first 200 tasks were classified by both observers and any uncertainty was flagged for review together.
For all jurisdictions, any methods descriptions that were vague, lacked context, or required further assumptions were excluded (2.6% of management tasks, U.S. only). Some management tasks (3.9%) were scored as both action and RM (e.g., translocate birds, action, and monitor the success of the release, RM; weed surveillance, RM, and control, action). For some management tasks, the distinction between action and RM was unclear. These tasks were discussed among the authors and the technician to reach a consensus. For example, ‘standard surveillance to detect invasive mammals’ in NZ could be considered an action, since it is required to detect and subsequently control invasions. However, we assigned it as RM because other management tasks clearly include an action component (e.g., ‘surveillance for invasive species and control if detected’) and other authors have categorized invasive species surveillance as monitoring46. Generally, management tasks to develop conservation plans are distinct from implementing plans and were thus scored as RM (K. Martin pers. comm.). Where we were unable to distinguish between RM and action, we scored as both action and RM.
For a subset of 8050 management tasks (the first 207 species) in U.S. recovery plans, we further categorized the type of RM to explore common RM activities (Supplementary Table 1). Because we found that assigning management tasks into these 17 categories was challenging without making subjective judgement calls, we did not analyze specific tasks further.
We estimated the cost of implementing each management task for each species following similar methods to those previously published, calculating costs over 50 years15, Supplementary Methods16,17. We calculated the proportion of the proposed budget allocated to RM for each species as the total cost of all management tasks scored as research or monitoring divided by the total cost of all management tasks. For management tasks that were scored as both action and RM, we multiplied the cost of the task by the average proportional difference between action and RM for each jurisdiction.
Factors affecting the proportion allocated to RM
We compared the characteristics of each species recovery plan with the proportion of proposed spending designated as RM. Characteristics available in recovery planning databases for all three jurisdictions included taxon, the estimated benefit of implementing all management tasks, and the total budget estimated for each species (Table 1, Supplementary Methods). The most general category shared among all jurisdictions was taxon, resulting in nine categories: amphibians, birds, bryophytes, fishes, fungus, invertebrates, mammals, reptiles, and vascular plants (set as a reference category). Lichens were removed from further analysis because there were only two species. For NZ and NSW, we extracted expert-elicited estimates of the benefit of implementing all management tasks, where experts were asked to consider the probability of species being secure in 50 years with and without the suite of management tasks16,17. Thus, benefit was calculated as the difference between the probability of security with and without the management tasks. For the U.S., in the absence of expert elicitation, the benefit of completing all management tasks in a recovery plan was approximated using information embedded in Recovery Priority Numbers (RPN). RPNs are an 18-category numeric rank for each species based on three categories of threat (high, moderate, and low), high or low recovery potential, and taxonomic distinctness monotypic genus, species, and subspecies,47. The limitations of using RPN to estimate the probability of persistence with or without management are discussed by Gerber et al.15 and Avery-Gomm48. To generate the total budget for each species, we used previously published total costs, which considered actions that benefited more than one species cost as shared among species projects15,16,17. In all further analysis, we removed species with a proposed budget of 0 (23 species in the U.S.) and extinct species (Guam broadbill—Myiagra freycineti and Eastern puma—Puma concolor couguar).
We explored additional characteristics unique to U.S. recovery planning documents, using U.S. data only (Table 1). These included: the federal listing status, the number of species in the recovery plan (66% of plans include multiple species), the priority assigned to each management task (1: emergency measures needed to prevent extinction, 2: measures required to stabilize a species headed for extinction, and 3: needed to delist), the estimated management task duration in years, the fiscal year the management task was implemented, the management task status (ongoing, complete, planned, discontinued), the total estimated time to recovery, and an RPN, which we used to make a new factor called ‘recovery potential’ (one of six scores based on the RPN, where the highest had a high probability of recovery and low degree of threat and the lowest had a low probability of recovery and a high degree of threat). Federal listing status was collapsed from six into three categories: endangered, threatened, and not listed (including candidate species, species removed from ESA due to recovery, or populations considered as ‘non-essential, experimental’). Taxa were assigned to eight categories: amphibians, birds, fishes, invertebrates (set as a reference category), mammals, reptiles, and flowering and non-flowering plants.
Quantitative analysis
To examine what characteristics of recovery plans are associated with the proportion of the budget allocated to RM we used beta regression in the betareg package49 in R version 3.6.150. We fit two models—one including all data, with jurisdiction included as a covariate, and one including a wider suite of covariates only available for the U.S. (Table 1). All continuous covariates were standardized by subtracting the mean and dividing by the standard deviation to ensure the resulting parameter estimates would be comparable51. We standardized the total budget of each jurisdiction separately to account for each countries’ different currency and year the budget was estimated. To improve model fit we removed five species with total budgets over 5 million dollars (five times the median: Barton Springs salamander – Eurycea sosorum, Austin blind salamander—Eurycea waterlooensis, Indiana bat—Myotis sodalis, Bull trout – Salvelinus confluentus, Grizzly Bear—Ursus arctos horribilis). Our results are robust to the inclusion or exclusion of these species.
Categorical covariates were converted to dummy variables. To select a reference category, we ran an initial model, using the category with the lowest mean proportion of budget RM as the reference. In this initial model, we selected the dummy variable with the highest variance inflation factor VIF in the car package52; as the reference in the final model. As a result, all VIF were 3). We excluded correlated covariates in successive models and chose the final model with the lowest Akaike’s Information Criterion (AIC53). The final model excluded the total proposed budget, which was correlated with the number of species in multi-species plans and the first fiscal year of the earliest RM. We consider any covariates where 95% confidence intervals around parameter estimates exclude zero to indicate a significant effect.
Estimating species recovery outcomes
To assess the relationship between the proportion of the budget allocated to RM and species recovery outcomes, we extracted a previously published index of recovery for U.S. listed species2 and developed similar indices based on annual and semi-annual reports from NZ and NSW (Supplementary Methods).
To generate the U.S. recovery index, Gerber2 calculated sums of biennial status data from reports to Congress during 1989–2011 (total of 11 status reports30). For each species, reports included whether their status was extinct, declining (scored as −1), stable (scored as 0), improving (scored as +1), or unknown. These scores were summed, generating values from −11 to 11, indicating whether species are declining or improving more frequently.
To develop recovery indices for NZ and NSW, we used similar reports through the New Zealand Threat Classification System and New South Wales Saving our Species annual report card over 4 and 5 assessment periods respectively. Assessments were annual in NSW and in NZ the periods between reports were on average every 4 years (Supplementary Methods). For each update or report card, we used a similar scoring (−1, 0, and +1) to indicate whether species were declining, stable, or improving between assessments (further details in Supplementary Methods). Note that in this analysis we were limited to a subset of the 2328 threatened species (78.5% of U.S. species, 13.5% of NZ species, and 14.7% of NSW species). Other studies have noted the limitations of recovery assessments28.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

Overall strategy to detect the relevance of Unknown taxa
A pipeline based on network analysis was developed to detect and quantitatively measure the overall and individual impact of Unknown taxa on their environmental communities (Fig. 1). Briefly, Illumina 16S rRNA sequencing fastq files belonging to four distinctive aquatic extreme environments (i.e., hot springs, hypersaline, deep sea, and polar habitats that included both Arctic and Antarctic samples) were collected from public databases and 45 different BioProjects (Fig. 2a and Supplementary Dataset S1). We included different environment types to assess general and environment-specific patterns and chose to use extreme habitats as they comprise some of the harshest and relatively understudied habitats on Earth, and therefore, are likely to contain uncharacterized organisms.
Fig. 1: Overview of the analysis pipeline.

A minimum of 250 samples was retrieved for each of the four different extreme environments—hot springs (red), hypersaline (dark green), deep sea (turquoise), and polar (blue). Sequence reads were quality filtered, assigned to a taxonomy, and clustered to OTUs. At each classification level, any unassigned, ambiguous, or uncultured OTUs were designated as Unknowns, or “microbial dark matter” (MDM). For each environment, at each classification level, the direct co-occurrence relationship between all OTUs was mathematically modeled as a network. Networks were created for each environment, across all taxonomic classification levels, including all OTUs (Original, orange), excluding MDM (Without Unknowns, light green), and excluding an equal number of random Knowns (Bootstrap, blue). Network centrality metrics (i.e., degree, betweenness, and closeness) were calculated for each node, compared, and visualized as boxplots between these network types. Hub scores were calculated for each node in the Original network and networks were recreated, resizing by hub score, where the largest size node indicates a top hub species.

Full size image

Fig. 2: Summary of environmental 16S rRNA gene data.

a Map of the sample sites used in this study. Circles symbolize sample locations and are color-coded by environment: hot springs (HS, red), hypersaline (HY, dark green), deep sea (DS, turquoise), and polar (PO, blue). b Summary of data used in this study. OTUs counts are provided at the genus level. c Proportion of “microbial dark matter” (MDM) OTUs for each environment labeled as unassigned (dark blue), uncultured (dark red), and ambiguous (yellow) after SILVA and UCLUST-based taxonomic assignment to OTU. d Venn diagram of shared OTUs in four extreme environments. Each pie chart depicts the proportion of unique OTUs that were Known (lighter shade) and Unknown (darker shade) for each environment, with the bottom-most pie chart showing combined data for all environments. e Prevalence curves indicate the number of unique OTUs consistently present at an increasing number of samples. Dotted lines signify the prevalence of MDM OTUs and solid lines signify the prevalence of Known OTUs.

Full size image

After quality filtering, reads were mapped to OTUs and annotated against the SILVA (v128) reference database [38] by an open-reference strategy, i.e., allowing the detection of Unknown OTUs [39]. Over two million Known and novel taxa were observed with the vast majority of the taxa annotated as Bacteria. Only the bacterial taxa or OTUs unclassified at the domain-level were targeted for downstream network analyses to demonstrate the feasibility of this approach across ecosystems. As the term “microbial dark matter” can have a broad meaning, here, we define Unknown taxa as uncultured, unassigned, or ambiguous by the reference database at each taxonomic classification level (e.g., phylum to genus). For each environment, networks reflecting across-samples co-occurrence relationships between all taxa, Known and Unknown, were constructed and referred to as the “Original” networks (Fig. 1). To assess the role of the Unknown taxa on network structure and properties, the Unknown nodes were removed from the ‘Original’ network and a new network, referred to as ‘Without Unknown’ was reconstructed. To ensure that changes in network properties were not caused just by the number of nodes, a third network, referred to as the “Bootstrap” network, was created where a random set of nodes of the same number as the Unknown OTUs was removed. The relevance of the Unknown taxa was assessed by comparing changes in degree, closeness, and betweenness scores between the three network types and by evaluating the frequency of Unknowns as top hubs within each of the “Original” environmental networks.
A similar and significant fraction of Unknown taxa populates diverse environments
To assess whether there were distinctive patterns or trends of Unknown taxa within the four targeted environments, data from each habitat type were mined from several geographical locations across the globe (Fig. 2a and Supplementary Dataset S1). Reads were collected from the online repositories National Center for Biotechnology Information Sequence Read Archive (NCBI SRA) and Joint Genome Institute Genomes Online Database (JGI Gold). For each environment, between 255 and 286 16S rRNA samples and between 51 and 57 million reads were included in the analysis, resulting in a total of 219,980,340 reads from 1086 samples (Fig. 2b).
After processing, quality filtering, and OTU assignment steps were completed, there were 2102,595 unique OTUs totaling 164,896,127 amplicon read counts derived from these samples (Fig. 2b). The relative proportions of Unknown OTUs, which were designated as unassigned, ambiguous, or uncultured by the SILVA database, were compared between environments. Results indicated that all environments showed similar relative contributions of these three Unknown types, with unassigned and uncultured OTUs making up the majority of the Unknown component composition (Fig. 2c). Regardless of environment type, within each of the four habitats, between 25 and 38% of unique OTUs were cataloged as Unknown (Fig. 2b, d). Samples collected from polar habitats were significantly enriched (Fisher’s Exact Test p value < 0.05) in Unknown OTUs despite having the highest total read counts. The higher proportion of Unknown OTUs in the polar samples compared to the other habitats likely reflects the less-well-characterized biological diversity of these Arctic and Antarctic ecosystems. Next, the proportion of shared Known and Unknown OTUs between environments was evaluated. Most OTUs, regardless of assigned or unassigned taxonomic status, were environment-specific, with only 11,318 out of the 2,102,595 OTUs present in all four of the environments (Fig. 2d). The majority of shared OTUs were observed between the hypersaline and polar environments and between hypersaline and hot springs environments, possibly reflecting the widespread distribution of hypersaline habitats across diverse thermal zones. Unsurprisingly, polar and hot springs environmental samples shared the least number of OTUs (Fig. 2d). Given this low common OTU pool across environments, network analyses were applied to each environment independently. Last, the prevalence (i.e., the percentage of samples with nonzero counts in which an OTU was detected) of Known and Unknown OTUs was evaluated at the genus level to assess the consistency of OTU detection within each environment. The OTU matrix was sparse, with the majority of taxa observed in ≤50 (~20%) samples (Fig. 2e). However, prevalence curves were generally very similar for Known and Unknown OTUs in all four environments, indicating that Unknown OTUs are not necessarily rarer than already characterized species (Fig. 2e). Moreover, we confirmed that Unknown OTUs, like Known OTUs, were generally present and consistent across multiple studies within the same environment and did not tend to concentrate in any particular project (Supplementary Figs. S1–4). Consequently, these results indicated that a network analysis of these data would be a reflection of the co-occurrence structure of the community and not of potential compositional bias. Network analysis of OTU abundance at different taxonomic levels reveals the connectivity of unknown microbes Having demonstrated that the Unknown taxa comprise a substantial proportion of unique OTUs and have comparable abundances to Known taxa within a community, network metrics were used to effectively compare the ecological relevance of both Known and Unknown taxa in subsequent networks. Microbial association networks were constructed that featured only significant co-occurrence correlation relationships for OTUs with a notable prevalence in each environment, meaning that any OTUs that were not detected in a sufficient number of samples were removed. To select a suitable prevalence threshold, the proportion of Known and Unknown taxa across a range of sample percentages was evaluated (Supplementary Fig. S5). Across all taxonomic levels, the Known and Unknown taxa of hot springs and polar habitats were more prevalent than those of hypersaline and deep sea communities; therefore, a slightly more stringent prevalence threshold (40%) was chosen for the former, and a lower threshold value (30%) chosen for the latter environments. These thresholds resulted in the retention of a similar fraction of data from the initial OTU count (Table 1), with 102–297 nodes present per environment, making the networks both suitably large and comparable. Table 1 Breakdown of node and edge attributes across extreme environmental networks. Full size table The SpiecEasi Meinshausen–Buhlmann (MB) neighborhood algorithm was then used to construct networks (see “Methods”) that contained at least 100 nodes (i.e., OTUs) per environment and had edge-to-node ratios that varied from 1.9 and 2.9 (Table 1). Although most OTUs that passed the prevalence criteria became elements of the networks, no relationship between the initial data size (i.e., number of samples and taxa) and the interconnectedness (i.e., nodes/edges ratio) of the resulting network was observed (Table 1). The two environments with the highest number of initial OTUs (i.e., hypersaline and deep sea) had the lowest number of prevalent members, 193 and 102, respectively, and very different edge/node ratios, 2.7 and 1.9, respectively, indicating that high prevalence does not necessarily correlate with high co-occurrence. Similarly, the polar and hot springs networks retained a high number of prevalent OTUs but differed in edge number, yet again, indicating that the observed network structures were the result of the intrinsic properties of each environment and were not dependent on the sampling procedure. Interestingly, when evaluating the proportion of Unknown edges and Unknown-Unknown connections at the genus level, similar patterns were observed across environments. Between 45 and 62% of all connected nodes were Unknown OTUs and a higher proportion of Unknown-Unknown versus Known-Known links was present at the genus level (Table 1), for all environments but hot springs, where a higher proportion of Known-Unknown and Known-Known links was observed. The results of this global analysis of network construction and composition demonstrate that although the general community connectivity might be environment-specific, the relative contribution of Known and Unknown taxa to these networks is similar. Once again, network properties were not a direct outcome of sampling biases, but more likely, reflect the biology of their respective ecosystems. Unknown taxa play important roles in interconnectedness and connectivity of extreme environmental microbial networks Next, the position and neighborhood of Unknown nodes were examined. At the phylum level, Unknown taxa were present in the hot spring, hypersaline, and polar environments, but were not found in the deep sea network (Supplementary Fig. S6). The class level was the first taxonomic classification rank in which Unknown taxa were present in all environmental networks. To accurately assess the role of Unknowns, we evaluated class-level networks and observed that the hypersaline and polar Unknown OTUs created distinctive clusters, whereas hot springs and deep sea Unknown OTUs were intermixed with Known taxa (Fig. 3a). Hypersaline and polar Unknowns consistently appeared to be isolated and peripherally located compared to the centrally positioned hot springs and deep sea Unknown nodes across almost all classification ranks (Supplementary Fig. S6). Consequently, these results suggest that the clustering pattern is unrelated to a higher abundance of Unknowns and is more environment dependent. For example, the targeted hot spring environments had the highest number of Unknown OTUs yet showed the most dispersed connections between the these taxa (Fig. 3a). Thus, the inclusion of the Unknown taxa in our environmental networks models was, as anticipated, not solely the consequence of their level of prevalence, but rather a reflection of a particular abundance co-variation pattern. Fig. 3: Analysis of environmental network taxa interconnectedness. a Microbial networks at class taxonomic classification level. Nodes are colored by class assignment, with gray nodes representing Unknown taxa at the class level. b Bar graphs of the co-occurrence relationships (i.e., edges) of Unknown OTUs with other taxa at the class level within each environmental network. Y axis labels and colors signify the different classes with which Unknowns were found to co-occur. Unknown-Unknown relationships are represented in gray. Full size image For all four environments, Unknown OTUs had more frequently shared edges among them than with classified taxa (Fig. 3b). Unknown OTUs were found to frequently co-occur with each other within each environmental network, although the frequency of within-class interactions for Unknowns at the class level was found to be statistically no greater than all other within-class interactions for each environment (Supplementary Fig. S7). In fact, the pattern of a high frequency of shared edges among members of the same class held true for known classes as well (Supplementary Fig. S8). In accordance with other studies, these results demonstrate that OTUs of the same taxonomic classification most frequently co-occur with each other [40, 41]. Furthermore, the high frequency of shared edges between Unknown classes suggests that Unknown OTUs might be taxonomically related. To ensure that the observations found were reproducible, robust, and not biased by earlier steps of the analysis, the diversity and position of Unknown taxa in the networks were examined for several parameters. Although the number of Unknown nodes changed at each level (Supplementary Fig. S6), the environment-specific network patterns observed at the class level (Fig. 3a) were retained at other taxonomic levels. In addition, to determine whether the topology of the network was a direct consequence of our correlation metric of choice or the prevalence threshold, three other network construction approaches were used: SparCC [11], CClasso [42], and Pearson correlation. Network analyses were performed across a range of prevalence thresholds (15–35%, at 5% increments). Again, regardless of the network construction approach or sample percentage applied, network shape and Unknown OTU position remained consistent and each environment exhibited a distinctive pattern of Unknown taxa inclusion. For example, Unknown nodes continued to occupy peripheral positions for hypersaline and polar networks, whereas nodes in hot springs and deep sea environments were more centrally located when applying different correlation metrics (Supplementary Fig. S9). Although networks appeared “noisy” at more lenient prevalence thresholds (15–20%; Supplementary Figs. S10–13), the networks and positioning of Unknowns at higher percent thresholds were similar in appearance to the “Original” networks for all environments. Based on these results, we found that our general analysis strategy was robust across parameter choices and, therefore, these networks captured critical features of the relationships among taxa within each distinctive environment. Microbial dark matter acts as unifiers and frequent hubs within extreme environmental networks Although these results suggest that Unknown taxa were highly interconnected, these observations did not reveal how the presence of Unknowns affected the overall community structure. To more fully understand this role, we analyzed how network properties changed in the presence and absence of Unknown OTU nodes. We evaluated changes in degree, betweenness, and closeness, as different network metrics reveal different aspects of the relevance of nodes within their networks. This approach has the potential to ascertain whether certain Unknown OTUs were more centrally positioned (e.g., due to high closeness scores), more essential for joining other taxa (e.g., high betweenness), or simply more prevalent and likely to co-occur with others (e.g., high degree). To control for the effect of node removal and distinguish effects of Unknown taxa from network size, networks were generated that excluded several randomly picked Known nodes equal to the number of Unknown OTUs. This process was repeated 100 times to create a distribution of “Null” or “Bootstrap” networks for statistical comparisons. Differences in network parameters between networks without Unknown OTUs and the “Original” or “Bootstrap” networks were determined by the Wilcoxon test and p values were adjusted using the Holm method [43]. Strikingly, removal of the Unknown taxa caused a statistically significant impact on all measured network metrics in all four studied environments (Table 2, Fig. 4, and Supplementary Figs. S14–24). In the polar environments, for example, removal of Unknown OTUs caused a significant decrease in degree (p value  More

1.
Carson, H. S., Cook, G. S., López-Duarte, P. C. & Levin, L. A. Evaluating the importance of demographic connectivity in a marine metapopulation. Ecology 92, 1972–1984. https://doi.org/10.1890/11-0488.1 (2011).
PubMed  Article  Google Scholar 
2.
Jackson, J. B. C., Donovan, M. K., Cramer, K. L., Lam, V. & Lam, W. Status and trends of Caribbean coral reefs: 1970–2012. Glob. Coral Reef Monit. Network, IUCN, Gland. Switz. 306 (2014).

3.
Palumbi, S. R. Population genetics, demographic connectivity, and the design of marine reserves. 13, 146–158, https://doi.org/10.1890/1051-0761(2003)013[0146:PGDCAT]2.0.CO;2 (2003).

4.
Botsford, L. W. et al. Connectivity and resilience of coral reef metapopulations in marine protected areas: matching empirical efforts to predictive needs. Coral Reefs 28, 327–337. https://doi.org/10.1007/s00338-009-0466-z (2009).
ADS  PubMed  PubMed Central  CAS  Article  Google Scholar 

5.
Galindo, H. M., Olson, D. B. & Palumbi, S. R. Seascape genetics: a coupled oceanographic-genetic model predicts population structure of Caribbean corals. Curr. Biol. 16, 1622–1626. https://doi.org/10.1016/j.cub.2006.06.052 (2006).
PubMed  CAS  Article  Google Scholar 

6.
Rippe, J. P. et al. Population structure and connectivity of the mountainous star coral, Orbicella faveolata, throughout the wider Caribbean region. Ecol. Evol. 7, 9234–9246. https://doi.org/10.1002/ece3.3448 (2017).
PubMed  PubMed Central  Article  Google Scholar 

7.
Studivan, M. S. & Voss, J. D. Population connectivity among shallow and mesophotic Montastraea cavernosa corals in the Gulf of Mexico identifies potential for refugia. Coral Reefs 37, 1183–1196. https://doi.org/10.1007/s00338-018-1733-7 (2018).
ADS  Article  Google Scholar 

8.
Baums, I. B., Johnson, M. E., Devlin-Durante, M. K. & Miller, M. W. Host population genetic structure and zooxanthellae diversity of two reef-building coral species along the Florida Reef Tract and wider Caribbean. Coral Reefs 29, 835–842. https://doi.org/10.1007/s00338-010-0645-y (2010).
ADS  Article  Google Scholar 

9.
Serrano, X. M. et al. Long distance dispersal and vertical gene flow in the Caribbean brooding coral Porites astreoides. Sci. Rep. 6, 21619. https://doi.org/10.1038/srep21619 (2016).
ADS  PubMed  PubMed Central  CAS  Article  Google Scholar 

10.
Serrano, X. et al. Geographic differences in vertical connectivity in the Caribbean coral Montastraea cavernosa despite high levels of horizontal connectivity at shallow depths. Mol. Ecol. 23, 4226–4240. https://doi.org/10.1111/mec.12861 (2014).
PubMed  CAS  Article  Google Scholar 

11.
Bongaerts, P. et al. Deep reefs are not universal refuges: reseeding potential varies among coral species. Sci. Adv. 3, e1602373. https://doi.org/10.1126/sciadv.1602373 (2017).
ADS  PubMed  PubMed Central  Article  Google Scholar 

12.
Eckert, R. J., Studivan, M. S. & Voss, J. D. Populations of the coral species Montastraea cavernosa on the Belize Barrier Reef lack vertical connectivity. Sci. Rep. 9, 7200. https://doi.org/10.1038/s41598-019-43479-x (2019).
ADS  PubMed  PubMed Central  CAS  Article  Google Scholar 

13.
Goodbody-Gringley, G., Vollmer, S. V., Woollacott, R. M. & Giribet, G. Limited gene flow in the brooding coral Favia fragum (Esper, 1797). Mar. Biol. 157, 2591–2602. https://doi.org/10.1007/s00227-010-1521-6 (2010).
Article  Google Scholar 

14.
Goodbody-Gringley, G., Woollacott, R. M. & Giribet, G. Population structure and connectivity in the Atlantic scleractinian coral Montastraea cavernosa (Linnaeus, 1767). Mar. Ecol. 33, 32–48. https://doi.org/10.1111/j.1439-0485.2011.00452.x (2012).
ADS  CAS  Article  Google Scholar 

15.
Nunes, F. L. D., Norris, R. D. & Knowlton, N. Long distance dispersal and connectivity in Amphi-Atlantic corals at regional and basin scales. PLoS ONE 6, e22298. https://doi.org/10.1371/journal.pone.0022298 (2011).
ADS  PubMed  PubMed Central  CAS  Article  Google Scholar 

16.
Brazeau, D. A., Lesser, M. P. & Slattery, M. Genetic structure in the coral, Montastraea cavernosa: assessing genetic differentiation among and within mesophotic reefs. PLoS ONE 8, e65845. https://doi.org/10.1371/journal.pone.0065845 (2013).
ADS  PubMed  PubMed Central  CAS  Article  Google Scholar 

17.
Foster, N. L. et al. Connectivity of Caribbean coral populations: complementary insights from empirical and modelled gene flow. Mol. Ecol. 21, 1143–1157. https://doi.org/10.1111/j.1365-294X.2012.05455.x (2012).
PubMed  Article  Google Scholar 

18.
García-Machado, E., Ulmo-Díaz, G., Castellanos-Gell, J. & Casane, D. Patterns of population connectivity in marine organisms of Cuba. Bull. Mar. Sci. 94, 193–211. https://doi.org/10.5343/bms.2016.1117 (2018).
Article  Google Scholar 

19.
Ulmo-Díaz, G. et al. Genetic differentiation in the mountainous star coral Orbicella faveolata around Cuba. Coral Reefs 37, 1217–1227. https://doi.org/10.1007/s00338-018-1722-x (2018).
ADS  Article  Google Scholar 

20.
Creary, M. et al. Status of coral reefs in the northern Caribbean and western Atlantic GCRMN node in 2008. in Status of Coral Reefs of the World (ed. Wilkinson, C.) 239–252 (2008).

21.
Claro, R., Reshetnikov, Y. S. & Alcolado, P. M. Physical attributes of coastal Cuba. Ecol. Mar. fishes Cuba 1–20 (2001).

22.
Whittle, D. & Rey Santos, O. Protecting Cuba’s environment: efforts to design and implement effective environmental laws and policies in Cuba. Cuban Stud. 37, 73–103. https://doi.org/10.1353/cub.2007.0018 (2006).
Article  Google Scholar 

23.
Caballero, H., Alcolado, P. M. & Semidey, A. Condición de los arrecifes de coral frente a costas con asentamientos humanos y aportes terrígenos: el caso del litoral habanero. Cuba. Rev. Ciencias Mar. y Costeras 1, 49. https://doi.org/10.15359/revmar.1.3 (2009).
Article  Google Scholar 

24.
González-Díaz, P. et al. Status of Cuban coral reefs. Bull. Mar. Sci. 94, 229–247. https://doi.org/10.5343/bms.2017.1035 (2018).
Article  Google Scholar 

25.
Alcolado, P. M., Caballero, H. & Perera, S. Tendencia del cambio en el cubrimiento vivo por corales pétreos en los arrecifes coralinos de Cuba. Ser. Ocean. 5, 1–14 (2009).
Google Scholar 

26.
Toth, L. T. et al. Do no-take reserves benefit Florida’s corals? 14 years of change and stasis in the Florida Keys National Marine Sanctuary. Coral Reefs 33, 565–577. https://doi.org/10.1007/s00338-014-1158-x (2014).
ADS  Article  Google Scholar 

27.
Zlatarski, V. N. & Estalella, N. M. Los esclaractinios de Cuba. (2017).

28.
Zlatarski, V. N. Investigations on mesophotic coral ecosystems in Cuba (1970–1973) and Mexico (1983–1984). CICIMAR Oceánides 33, 27–43 (2018).
Google Scholar 

29.
Reed, J. et al. Cuba’s mesophotic coral reefs and associated fish communities. Rev. Investig. Mar. 38, 56–125 (2018).
Google Scholar 

30.
Baisre, J. A. An overview of Cuban commercial marine fisheries: the last 80 years. Bull. Mar. Sci. 94, 359–375. https://doi.org/10.5343/bms.2017.1015 (2018).
Article  Google Scholar 

31.
Gil-Agudelo, D. L. et al. Coral reefs in the Gulf of Mexico large marine ecosystem: conservation status, challenges, and opportunities. Front. Mar. Sci. 6, 807. https://doi.org/10.3389/fmars.2019.00807 (2020).
Article  Google Scholar 

32.
Perera Valderrama, S. et al. Marine protected areas in Cuba. Bull. Mar. Sci. 94, 423–442. https://doi.org/10.5343/bms.2016.1129 (2018).
Article  Google Scholar 

33.
NOAA. Sister Sanctuary: Memorandum of Understanding. (2015).

34.
Budd, A. F., Nunes, F. L. D., Weil, E. & Pandolfi, J. M. Polymorphism in a common Atlantic reef coral (Montastraea cavernosa) and its long-term evolutionary implications. Evol. Ecol. 26, 265–290. https://doi.org/10.1007/s10682-010-9460-8 (2012).
Article  Google Scholar 

35.
Bongaerts, P., Ridgway, T., Sampayo, E. M. & Hoegh-Guldberg, O. Assessing the `deep reef refugia’ hypothesis: focus on Caribbean reefs. Coral Reefs 29, 309–327. https://doi.org/10.1007/s00338-009-0581-x (2010).
Article  Google Scholar 

36.
Reed, J. K. Deepest distribution of Atlantic hermaptypic corals discovered in the Bahamas. Proc. Fifth Int. Coral Reef Congr. 6, 249–254 (1985).
Google Scholar 

37.
Szmant, A. M. Sexual reproduction by the Caribbean reef corals Montastrea annularis and M. cavernosa. Mar. Ecol. Prog. Ser. 74, 13–25. https://doi.org/10.3354/meps074013 (1991).
ADS  Article  Google Scholar 

38.
Highsmith, R. C., Lueptow, R. L. & Schonberg, S. C. Growth and bioerosion of three massive corals on the Belize barrier reef. Mar. Ecol. Prog. Ser. 13, 261–271 (1983).
ADS  Article  Google Scholar 

39.
Kitchen, S. A., Crowder, C. M., Poole, A. Z., Weis, V. M. & Meyer, E. De novo assembly and characterization of four anthozoan (Phylum Cnidaria) transcriptomes. G3 (Genes|Genomes|Genetics). 5, 2441–2452, https://doi.org/10.1534/g3.115.020164 (2015).
PubMed  PubMed Central  CAS  Article  Google Scholar 

40.
Matz Lab. Montastraea cavernosa annotated genome. (2018).

41.
Drury, C., Pérez Portela, R., Serrano, X. M., Oleksiak, M. & Baker, A. C. Fine-scale structure among mesophotic populations of the great star coral Montastraea cavernosa revealed by SNP genotyping. Ecol. Evol. 1–11, https://doi.org/10.1002/ece3.6340 (2020).

42.
Ellegren, H. Microsatellites: simple sequences with complex evolution. Nat. Rev. Genet. 5, 435–445. https://doi.org/10.1038/nrg1348 (2004).
PubMed  CAS  Article  Google Scholar 

43.
Joshi, D., Ram, R. N. & Lohani, P. Microsatellite markers and their application in fisheries. Int. J. Adv. Agric. Sci. Technol. 4, 67–104 (2017).
Google Scholar 

44.
Jarne, P. & Lagoda, P. J. L. Microsatellites, from molecules to populations and back. Trends Ecol. Evol. 11, 424–429. https://doi.org/10.1016/0169-5347(96)10049-5 (1996).
PubMed  CAS  Article  Google Scholar 

45.
Flores-Rentería, L. & Krohn, A. Scoring microsatellite loci. in Methods in molecular biology (ed. Kantartzi, S. K.) 319–336, https://doi.org/10.1007/978-1-62703-389-3_21 (Elsevier, 2013).

46.
Andrews, K. R., Good, J. M., Miller, M. R., Luikart, G. & Hohenlohe, P. A. Harnessing the power of RADseq for ecological and evolutionary genomics. Nat. Rev. Genet. 17, 81–92. https://doi.org/10.1038/nrg.2015.28 (2016).
PubMed  PubMed Central  CAS  Article  Google Scholar 

47.
Davey, J. L. & Blaxter, M. W. RADseq: next-generation population genetics. Brief. Funct. Genomics 9, 416–423. https://doi.org/10.1093/bfgp/elq031 (2010).
PubMed  CAS  Article  Google Scholar 

48.
Wang, S., Meyer, E., McKay, J. K. & Matz, M. V. 2b-RAD: a simple and flexible method for genome-wide genotyping. Nat. Methods 9, 808–810. https://doi.org/10.1038/nmeth.2023 (2012).
PubMed  CAS  Article  Google Scholar 

49.
Bradbury, I. R. et al. Transatlantic secondary contact in Atlantic Salmon, comparing microsatellites, a single nucleotide polymorphism array and restriction-site associated DNA sequencing for the resolution of complex spatial structure. Mol. Ecol. 24, 5130–5144. https://doi.org/10.1111/mec.13395 (2015).
PubMed  CAS  Article  Google Scholar 

50.
Jeffries, D. L. et al. Comparing RADseq and microsatellites to infer complex phylogeographic patterns, an empirical perspective in the Crucian carp, Carassius carassius. L. Mol. Ecol. 25, 2997–3018. https://doi.org/10.1111/mec.13613 (2016).
PubMed  Article  Google Scholar 

51.
Bohling, J., Small, M., Von Bargen, J., Louden, A. & DeHaan, P. Comparing inferences derived from microsatellite and RADseq datasets: a case study involving threatened bull trout. Conserv. Genet. 20, 329–342. https://doi.org/10.1007/s10592-018-1134-z (2019).
CAS  Article  Google Scholar 

52.
Thornhill, D. J., Xiang, Y., Fitt, W. K. & Santos, S. R. Reef endemism, host specificity and temporal stability in populations of symbiotic dinoflagellates from two ecologically dominant Caribbean corals. PLoS ONE https://doi.org/10.1371/journal.pone.0006262 (2009).
PubMed  PubMed Central  Article  Google Scholar 

53.
Eckert, R. J., Reaume, A. M., Sturm, A. B., Studivan, M. S. & Voss, J. D. Depth influences Symbiodiniaceae associations among Montastraea cavernosa corals on the Belize Barrier Reef. Front. Microbiol. 11, 1–13. https://doi.org/10.3389/fmicb.2020.00518 (2020).
Article  Google Scholar 

54.
Hume, B. C. C. et al. SymPortal: a novel analytical framework and platform for coral algal symbiont next-generation sequencing ITS2 profiling. Mol. Ecol. Resour. 19, 1063–1080. https://doi.org/10.1111/1755-0998.13004 (2019).
PubMed  PubMed Central  CAS  Article  Google Scholar 

55.
Pochon, X., Putnam, H. M., Burki, F. & Gates, R. D. Identifying and characterizing alternative molecular markers for the symbiotic and free-living dinoflagellate genus Symbiodinium. PLoS ONE https://doi.org/10.1371/journal.pone.0029816 (2012).
PubMed  PubMed Central  Article  Google Scholar 

56.
LaJeunesse, T. C. & Thornhill, D. J. Improved resolution of reef-coral endosymbiont (Symbiodinium) species diversity, ecology, and evolution through psbA non-coding region genotyping. PLoS ONE https://doi.org/10.1371/journal.pone.0029013 (2011).
PubMed  PubMed Central  Article  Google Scholar 

57.
Manzello, D. P. et al. Role of host genetics and heat-tolerant algal symbionts in sustaining populations of the endangered coral Orbicella faveolata in the Florida Keys with ocean warming. Glob. Change Biol. 25, 1016–1031. https://doi.org/10.1111/gcb.14545 (2018).
ADS  Article  Google Scholar 

58.
Warner, M. E., LaJeunesse, T. C., Robison, J. D. & Thur, R. M. The ecological distribution and comparative photobiology of symbiotic dinoflagellates from reef corals in Belize: potential implications for coral bleaching. Limnol Ocean. https://doi.org/10.4319/lo.2006.51.4.1887 (2006).
Article  Google Scholar 

59.
Finney, J. C. et al. The relative significance of host–habitat, depth, and geography on the ecology, endemism, and speciation of coral endosymbionts in the genus Symbiodinium. Microb. Ecol. 60, 250–263. https://doi.org/10.1007/s00248-010-9681-y (2010).
PubMed  Article  Google Scholar 

60.
Bongaerts, P. et al. Deep down on a Caribbean reef: lower mesophotic depths harbor a specialized coral-endosymbiont community. Sci. Rep. https://doi.org/10.1038/srep07652 (2015).
PubMed  PubMed Central  Article  Google Scholar 

61.
Reed, J. et al. Cruise report Cuba’s twilight zone reefs: Remotely Operated Vehicle surveys of deep/mesophotic coral reefs and associated fish communities of Cuba. (2017).

62.
Arriaza, L. et al. Modelación numérica de corrientes marinas alrededor del occidente de Cuba. Serie Oceanológica. 10, 11–22. (2012).
Google Scholar 

63.
Gordon, A. & Hannon, G. FASTX-Toolkit. FASTQ/A short-reads pre-processing tools. (2010). Available at: https://hannonlab.cshl.edu/fastx_toolkit/.

64.
Bayer, T. et al. Symbiodinium transcriptomes: genome insights into the dinoflagellate symbionts of reef-building corals. PLoS ONE https://doi.org/10.1371/journal.pone.0035269 (2012).
PubMed  PubMed Central  Article  Google Scholar 

65.
Davies, S. W., Ries, J. B., Marchetti, A. & Castillo, K. D. Symbiodinium functional diversity in the coral Siderastrea siderea is influenced by thermal stress and reef environment, but not ocean acidification. Front. Mar. Sci. 5, 1–14. https://doi.org/10.3389/fmars.2018.00150 (2018).
Article  Google Scholar 

66.
Ladner, J. T., Barshis, D. J. & Palumbi, S. R. Protein evolution in two co-occurring types of Symbiodinium: an exploration into the genetic basis of thermal tolerance in Symbiodinium clade D. BMC Evol. Biol. 12, 217. https://doi.org/10.1186/1471-2148-12-217 (2012).
PubMed  PubMed Central  CAS  Article  Google Scholar 

67.
Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. https://doi.org/10.1038/nmeth.1923 (2012).
PubMed  PubMed Central  CAS  Article  Google Scholar 

68.
Korneliussen, T. S., Albrechtsen, A. & Nielsen, R. ANGSD: Analysis of Next Generation Sequencing Data. BMC Bioinformatics 15, 356. https://doi.org/10.1186/s12859-014-0356-4 (2014).
PubMed  PubMed Central  Article  Google Scholar 

69.
Peakall, R. & Smouse, P. E. GenALEx 65: Genetic analysis in Excel. Population genetic software for teaching and research-an update. Bioinformatics 28, 2537–2539. https://doi.org/10.1093/bioinformatics/bts460 (2012).
PubMed  PubMed Central  CAS  Article  Google Scholar 

70.
Excoffier, L., Smouse, P. E. & Quattro, J. M. Analysis of Molecular Variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics 131, 479–491 (1992).
PubMed  PubMed Central  CAS  Google Scholar 

71.
Benjamini, Y. & Hochberg, Y. Controlling the False Discovery Rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 57, 289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x (1995).
MathSciNet  MATH  Article  Google Scholar 

72.
Prevosti, A., Ocaña, J. & Alonso, G. Distances between populations of Drosophila subobscura, based on chromosome arrangement frequencies. Theor. Appl. Genet. 45, 231–241. https://doi.org/10.1007/BF00831894 (1975).
PubMed  CAS  Article  Google Scholar 

73.
Kamvar, Z. N., Tabima, J. F. & Gr̈unwald, N. J. ,. Poppr: An R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ 2014, 1–14. https://doi.org/10.7717/peerj.281 (2014).
Article  Google Scholar 

74.
R Core Team. R: A language and environment for statistical computing. (2019).

75.
Jombart, T. & Ahmed, I. adegenet 13–1: new tools for the analysis of genome-wide SNP data. Bioinformatics 27, 3070–3071. https://doi.org/10.1093/bioinformatics/btr521 (2011).
PubMed  PubMed Central  CAS  Article  Google Scholar 

76.
Mantel, N. The detection of disease clustering and a generalized regression approach. Cancer Res. 27, 209–220 (1967).
PubMed  CAS  Google Scholar 

77.
Pritchard, J. K., Stephens, M. & Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 155, 945–959 (2000).
PubMed  PubMed Central  CAS  Google Scholar 

78.
Besnier, F. & Glover, K. A. ParallelStructure: A R package to distribute parallel runs of the population genetics program STRUCTURE on multi-core computers. PLoS ONE 8, 1–5. https://doi.org/10.1371/journal.pone.0070651 (2013).
CAS  Article  Google Scholar 

79.
Earl, D. A. & vonHoldt, B. M. STRUCTURE HARVESTER: A website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361. https://doi.org/10.1007/s12686-011-9548-7 (2012).
Article  Google Scholar 

80.
Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A. & Cresko, W. A. Stacks: an analysis tool set for population genomics. Mol. Ecol. 22, 3124–3140. https://doi.org/10.1111/mec.12354 (2013).
PubMed  PubMed Central  Article  Google Scholar 

81.
Pembleton, L. W., Cogan, N. O. I. & Forster, J. W. StAMPP: An R package for calculation of genetic differentiation and structure of mixed-ploidy level populations. Mol. Ecol. Resour. 13, 946–952. https://doi.org/10.1111/1755-0998.12129 (2013).
PubMed  CAS  Article  Google Scholar 

82.
Oksanen, J. et al. Vegan: community ecology package. (2019).

83.
Alexander, D. H. D., Novembre, J. & Lange, K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664. https://doi.org/10.1101/gr.094052.109 (2009).
PubMed  PubMed Central  CAS  Article  Google Scholar 

84.
Skotte, L., Korneliussen, T. S. & Albrechtsen, A. Estimating individual admixture proportions from next generation sequencing data. Genetics 195, 693–702. https://doi.org/10.1534/genetics.113.154138 (2013).
PubMed  PubMed Central  CAS  Article  Google Scholar 

85.
Foll, M. & Gaggiotti, O. A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics 180, 977–993. https://doi.org/10.1534/genetics.108.092221 (2008).
PubMed  PubMed Central  Article  Google Scholar 

86.
Evanno, G., Regnaut, S. & Goudet, J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611–2620. https://doi.org/10.1111/j.1365-294X.2005.02553.x (2005).
PubMed  CAS  Article  Google Scholar 

87.
Centurioni, L. R. & Niiler, P. P. On the surface currents of the Caribbean Sea. Geophys. Res. Lett. 30, 10–13. https://doi.org/10.1029/2002GL016231 (2003).
Article  Google Scholar 

88.
Candela, J. et al. The flow through the gulf of Mexico. J. Phys. Oceanogr. 49, 1381–1401. https://doi.org/10.1175/JPO-D-18-0189.1 (2019).
ADS  Article  Google Scholar 

89.
Kourafalou, V., Androulidakis, Y., Le Hénaff, M. & Kang, H. S. The Dynamics of Cuba Anticyclones (CubANs) and interaction with the Loop Current/Florida Current system. J. Geophys. Res. Ocean. 122, 7897–7923. https://doi.org/10.1002/2017JC012928 (2017).
ADS  Article  Google Scholar 

90.
Arriaza, L. et al. Marine current estimations in southeast Cuban shelf. Ser. Ocean. 4, 1–10 (2008).
Google Scholar 

91.
Carracedo-Hidalgo, D., Reyes-Perdomo, D., Calzada-estrada, A., Chang-Domínguez, D. & Rodríguez-Pupo, A. Characterization of sea currents in sea adjacent to Cuba . Main trends in the last years. Rev. Cuba. Meteorol. 25, (2019).

92.
Frys, C. et al. Fine-scale coral connectivity pathways in the Florida Reef Tract: implications for conservation and restoration. Front. Mar. Sci. 7, 1–42. https://doi.org/10.3389/fmars.2020.00312 (2020).
Article  Google Scholar 

93.
Kuba, A. Transgenerational effects of thermal stress: impacts on and beyond coral reproduction. (Nova Southeastern University, 2016).

94.
Claro, R., Lindeman, K. C., Kough, A. S. & Paris, C. B. Biophysical connectivity of snapper spawning aggregations and marine protected area management alternatives in Cuba. Fish. Oceanogr. 28, 33–42. https://doi.org/10.1111/fog.12384 (2019).
Article  Google Scholar 

95.
Holstein, D. M., Paris, C. B. & Mumby, P. J. Consistency and inconsistency in multispecies population network dynamics of coral reef ecosystems. Mar. Ecol. Prog. Ser. 499, 1–18. https://doi.org/10.3354/meps10647 (2014).
ADS  Article  Google Scholar 

96.
Szmant, A. M. Reproductive ecology of Caribbean reef corals. Coral Reefs 5, 43–53. https://doi.org/10.1007/BF00302170 (1986).
ADS  Article  Google Scholar 

97.
Drury, C. et al. Genomic variation among populations of threatened coral: Acropora cervicornis. BMC Genomics 17, 286. https://doi.org/10.1186/s12864-016-2583-8 (2016).
PubMed  PubMed Central  CAS  Article  Google Scholar 

98.
Devlin-Durante, M. K. & Baums, I. B. Genome-wide survey of single-nucleotide polymorphisms reveals fine-scale population structure and signs of selection in the threatened Caribbean elkhorn coral. Acropora palmata. PeerJ 5, e4077. https://doi.org/10.7717/peerj.4077 (2017).
PubMed  CAS  Article  Google Scholar 

99.
Wang, J., Feng, C., Jiao, T., Von Wettberg, E. B. & Kang, M. Genomic signature of adaptive divergence despite strong nonadaptive forces on Edaphic Islands: a case study of Primulina juliae. Genome Biol. Evol. 9, 3495–3508. https://doi.org/10.1093/gbe/evx263 (2017).
PubMed  PubMed Central  Article  Google Scholar 

100.
Ramos-Silva, P. et al. The skeletal proteome of the coral Acropora millepora: the evolution of calcification by co-option and domain shuffling. Mol. Biol. Evol. 30, 2099–2112. https://doi.org/10.1093/molbev/mst109 (2013).
PubMed  PubMed Central  CAS  Article  Google Scholar 

101.
Takeuchi, T., Yamada, L., Shinzato, C., Sawada, H. & Satoh, N. Stepwise evolution of coral biomineralization revealed with genome-wide proteomics and transcriptomics. PLoS ONE https://doi.org/10.1371/journal.pone.0156424 (2016).
PubMed  PubMed Central  Article  Google Scholar 

102.
Aranda, M. et al. Differential sensitivity of coral larvae to natural levels of ultraviolet radiation during the onset of larval competence. Mol. Ecol. 20, 2955–2972. https://doi.org/10.1111/j.1365-294X.2011.05153.x (2011).
PubMed  CAS  Article  Google Scholar 

103.
Reynolds, W. S., Schwarz, J. A. & Weis, V. M. Symbiosis-enhanced gene expression in cnidarian-algal associations: cloning and characterization of a cDNA, sym32, encoding a possible cell adhesion protein. Comp. Biochem. Physiol. Part A Mol. Integr. Physiol. 126, 33–44. https://doi.org/10.1016/S0742-8413(00)00099-2 (2000).
CAS  Article  Google Scholar 

104.
Iguchi, A. et al. Apparent involvement of a β1 type integrin in coral fertilization. Mar. Biotechnol. 9, 760–765. https://doi.org/10.1007/s10126-007-9026-0 (2007).
PubMed  CAS  Article  Google Scholar 

105.
Lesser, M. P. et al. Photoacclimatization by the coral Montastraea cavernosa in the mesophotic zone: light, food, and genetics. Ecology 91, 990–1003. https://doi.org/10.1890/09-0313.1 (2010).
PubMed  Article  Google Scholar 

106.
Klepac, C. et al. Seasonal stability of coral-Symbiodinium associations in the subtropical coral habitat of St. Lucie Reef, Florida. Mar. Ecol. Prog. Ser. 532, 137–151. https://doi.org/10.3354/meps11369 (2015).
ADS  Article  Google Scholar 

107.
Polinski, J. M. & Voss, J. D. Evidence of photoacclimatization at mesophotic depths in the coral-Symbiodinium symbiosis at Flower Garden Banks National Marine Sanctuary and McGrail Bank. Coral Reefs 37, 779–789. https://doi.org/10.1007/s00338-018-1701-2 (2018).
ADS  Article  Google Scholar 

108.
LaJeunesse, T. C. et al. Systematic revision of Symbiodiniaceae highlights the antiquity and diversity of coral endosymbionts. Curr. Biol. 28, 1–11. https://doi.org/10.1016/j.cub.2018.07.008 (2018).
CAS  Article  Google Scholar 

109.
Swain, T. D., Chandler, J., Backman, V. & Marcelino, L. Consensus thermotolerance ranking for 110 Symbiodinium phylotypes: an exemplar utilization of a novel iterative partial-rank aggregation tool with broad application potential. Funct. Ecol. 31, 172–183. https://doi.org/10.1111/1365-2435.12694 (2017).
Article  Google Scholar 

110.
Hodel, R. G. J. et al. Adding loci improves phylogeographic resolution in red mangroves despite increased missing data: comparing microsatellites and RAD-Seq and investigating loci filtering. Sci. Rep. 7, 1–14. https://doi.org/10.1038/s41598-017-16810-7 (2017).
CAS  Article  Google Scholar 

111.
Lemopoulos, A. et al. Comparing RADseq and microsatellites for estimating genetic diversity and relatedness — implications for brown trout conservation. Ecol. Evol. 9, 2106–2120. https://doi.org/10.1002/ece3.4905 (2019).
PubMed  PubMed Central  Article  Google Scholar 

112.
Puckett, E. E. Variability in total project and per sample genotyping costs under varying study designs including with microsatellites or SNPs to answer conservation genetic questions. Conserv. Genet. Resour. 9, 289–304. https://doi.org/10.1007/s12686-016-0643-7 (2017).
Article  Google Scholar 

113.
Hale, M. L., Burg, T. M. & Steeves, T. E. Sampling for microsatellite-based population genetic studies: 25 to 30 individuals per population is enough to accurately estimate allele frequencies. PLoS ONE 7, e45170. https://doi.org/10.1371/journal.pone.0045170 (2012).
ADS  PubMed  PubMed Central  CAS  Article  Google Scholar 

114.
Luikart, G., Sherwin, W. B., Steele, B. M. & Allendorf, F. W. Usefulness of molecular markers for detecting population bottlenecks via monitoring genetic change. Mol. Ecol. 7, 963–974. https://doi.org/10.1046/j.1365-294x.1998.00414.x (1998).
PubMed  CAS  Article  Google Scholar 

115.
Willing, E.-M., Dreyer, C. & van Oosterhout, C. Estimates of genetic differentiation measured by FST do not necessarily require large sample sizes when using many SNP markers. PLoS ONE 7, e42649. https://doi.org/10.1371/journal.pone.0042649 (2012).
ADS  PubMed  PubMed Central  CAS  Article  Google Scholar 

116.
Nazareno, A. G., Bemmels, J. B., Dick, C. W. & Lohmann, L. G. Minimum sample sizes for population genomics: an empirical study from an Amazonian plant species. Mol. Ecol. Resour. 17, 1136–1147. https://doi.org/10.1111/1755-0998.12654 (2017).
PubMed  CAS  Article  Google Scholar  More

EDITORIAL
22 September 2020

As China prepares to take on a crucial role in the governance of global biodiversity, its researchers must be at the table.

China’s researchers have shown that the country’s largest body of fresh water, Poyang Lake (pictured) in Jianxi province, has been drying out, in part because of the Three Gorges Dam.Credit: Fu Jianbin/Xinhua/ZUMA Wire

Last week, the United Nations confirmed that the world has failed, again, to achieve its goals to protect nature. This grim conclusion was delivered in the fifth edition of the United Nations Global Biodiversity Outlook report.
The report from the UN Convention on Biological Diversity reviewed progress towards 20 biodiversity targets that the convention’s participating countries set for themselves in Aichi, Japan, a decade ago (www.cbd.int/gbo).
None of the targets, which include making progress towards the sustainable harvesting of fish, controlling the spread of invasive species and preventing the extinction of threatened wildlife, will have been achieved by the deadline at the end of this year.
This is no time for regret or apology, but for urgency to act. Last year, an analysis by the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services revealed that some one million plant and animal species are at risk of extinction. And the wildlife charity WWF’s latest Living Planet Index, published earlier this month (see go.nature.com/32wzvdz), was similarly sobering, stating that vertebrate populations monitored between 1970 and 2017 have declined by an average of 68%.
All nations must do more, but some of the greatest responsibility now rests on the shoulders of China: the nation, along with the leaders of the UN biodiversity convention, will jointly host the next Conference of the Parties (COP) in Kunming next year. That summit, originally scheduled for this year, is where biodiversity targets for the next decade must be set.
As we have written before, the previous targets were destined to fail, in part because their format made progress hard to measure, and because countries did not need to report on what they were doing. This must now change. The targets, furthermore, need to be more closely aligned with the UN System of Environmental Economic Accounting, which is becoming the global standard for environmental reporting. Without these changes, the next set of biodiversity targets will almost certainly fail again.
At the same time, China’s biodiversity scientists and policy researchers should be at the table, too, as plans for Kunming start to take shape. The country has decades of experience of studying how to — and how not to — balance economic development with controlling species and ecosystems loss. The world needs to hear these stories, in all their complexity.
Learn from China
The Global Biodiversity Outlook report confirms that known species are on an accelerated path to extinction, with cycad and coral species among the groups most at risk. The report shows that, although deforestation has slowed in the past decade, forests are still being splintered by agriculture, tree-felling and urban growth. Such fragmentation will further harm biodiversity and increase carbon emissions.
Demand for food and agricultural production continue to be the main drivers of biodiversity loss. And governments are not helping. On average, they invest some US$500 billion per year in initiatives that harm the environment — eclipsing financing for biodiversity projects by a factor of 6, the report says.
China has a set of experiences that could help the world learn valuable lessons. Its rapid economic growth lifted a generation out of poverty; however, this created a cascade of environmental problems, not least elevated pollution in the air and on land. People in China rightly questioned their leaders for underestimating — if not downplaying — the environmental and social impacts of its industrialization. Partly in response, China’s authorities have been working with researchers from China and around the world to chart a greener way forward.
For example, national and local administrations have been devising and experimenting with environmental targets, and creating mechanisms for monitoring and reporting progress towards them — albeit with mixed success.
China’s national biodiversity strategy includes creating what it calls ‘redlines’ — areas where human activities are restricted to protect biodiversity — across the country.
Then there’s China’s US$6-trillion Belt and Road Initiative — a massive programme to build roads, ports and infrastructure, which will run through natural habitats across Asia, Europe and Africa. Much of this investment did not initially come with safeguards to mitigate environmental risks — but these are now being actively studied.
And last but not least, China has a large community of researchers working to quantify, in monetary terms, the value of natural capital and ecosystem services, so that people and policymakers can more clearly understand that nature’s services to people do not come for free.
On 30 September, heads of governments will meet at the UN for a day of talks on biodiversity, ahead of next year’s Kunming COP. Nature spoke to a number of representatives of national delegations who plan to attend this meeting, including researchers and non-governmental observers. All want the Kunming COP to succeed in bringing nations together and reaching an agreement on targets that are measurable and meaningful. But they expressed concern over the limited public engagement from China’s government about its goals or strategy for Kunming — and the relatively limited involvement of its researchers in the process so far.
Scientists in China have been central to their country’s conservation and economic-development journey. Their collective experience on what works, and what doesn’t, can provide important learning opportunities for countries as they look to slow down and eventually reverse bio-diversity and ecosystem loss. These researchers are in the academy of sciences; in universities; in the academy of environmental planning; and in the community of Chinese and international non-governmental organizations.
Many are also active in the China Council for Inter-national Cooperation on Environment and Development, an organization located in both Canada and China, which last week concluded a two-day conference presenting its latest research outputs. This important but little-known advisory body, now nearly three decades old, has been instrumental in connecting China’s environmental-science and environmental-policy communities with international counterparts.
Next year will be the first time that China has hosted an international environmental meeting — similar to the 2015 Paris climate accords — where the stakes are too high to fail. It must draw on its rich diversity of talent and experience. Other nations’ researchers must be equally forthcoming with their knowledge. All sides must put aside political differences to agree on ambitious targets, ways to achieve them and methods to measure that progress.
The best way to preserve and revive biodiversity is to acknowledge where we’ve all failed it before, to learn from that and to try again, together.

Nature 585, 481-482 (2020)
doi: 10.1038/d41586-020-02697-4

Latest on:

Biodiversity

Government

An essential round-up of science news, opinion and analysis, delivered to your inbox every weekday.

Related Articles More

Data
All the data sets analyzed in this work have been previously published and were obtained from EBI Metagenomics69. Previous publications (Supplementart Table 1) report the original experiments and corresponding analysis. In order to test the robustness of the macroecological laws and the modeling framework presented in this work, we considered 7 data sets that differ not only for the biome considered, but also for the sequencing techniques and the pipeline used to process the data. Data sets were selected to represent a wide set of biomes. We considered only data sets with at least 50 samples with more than 104 reads. No data set was excluded a-posteriori.
Sampling and compositional data
In order to study how (relative) abundance varies across communities and species, one needs to remove the effect of sampling noise, as it is not a biologically informative source of variation. By explicitly modeling sampling (Supplementary Note 2), one finds that the probability of observing n reads of species i in a sample with N total number of reads, is given by

$${P}_{i}(n| N)=int_{0}^{1}dx,{rho }_{i}(x),left(begin{array}{*{20}{c}} {n} \ {N} end{array}right){x}^{n}{(1-x)}^{N-n},,$$
(6)

where ρi(x) is the AFD, i.e., the probability (over communities or times) that the relative abundance of i is equal to x. Note that this equation does not assume anything about independence across species or communities. It only assumes the sampling process is carried independently across communities.
Since the random variable xi, whose distribution is ρi(x), is a relative abundance, one has that ∑ixi = 1 (i.e., the data are compositional70). As discussed in Supplementary Note 2, given the range of variation of the empirical relative abundances, one can substitute Eq. (6) with

$${P}_{i}(n| N)=int_{0}^{infty }dx,{rho }_{i}(x)frac{{(xN)}^{n}}{n!}{e}^{-xN},,$$
(7)

and the condition (sum)ixi = 1 to ({sum }_{i}{bar{x}}_{i}=1), where ({bar{x}}_{i}=mathop{int}nolimits_{0}^{infty }dx,{rho }_{i}(x)x) is the mean value of xi. Under this assumption, one can also take the limits of the integration from 0 to ∞, instead of considering them from 0 to 1, as the contribution of the integrand from 1 to ∞ is negligible.
Note that, because of sampling, the average of a function f(x) over the pdf ρ(x) differs in general from the average of f(n/N) over P(n∣N)

$$int_{0}^{1}dx,rho (x)f(x), ne, sum _{n = 0}^{N}P(n| N)fleft(frac{n}{N}right)=int_{0}^{1}dx,rho (x)sum _{n = 0}^{N}fleft(frac{n}{N}right)frac{{(xN)}^{n}}{n!}{e}^{-xN},,$$
(8)

and the inequality becomes equality only if f(x) is linear. The important difference between the right- and the left-hand side is often neglected in the literature. In fact, the right-hand side is a good approximation of the left-hand size only in the limit xN ≫ 1, which is far from being realized in the data for most of the species. Supplementary Note 2 introduces a method to reconstruct the moments of ρ(x) from the moments of P(n∣N). More generally, I show that it is possible to infer the moment generating function of ρ(x) from the data, which allow to reconstruct the shape of the empirical ρ(x).
Excluding competitive exclusion
A Gamma-distributed AFD implies that all the species present in a community of a biome are present in all the communities from that biome. Therefore, when a species is not observed is because it is undetected due to sampling errors. I test this claim in two different ways. First, it is shown that one can in fact predict the occupancy of a species from its abundance fluctuations. Secondly, I show that a model without true absences is statistically more supported than a model where species are allowed to be absent.
The first way to test this hypothesis is to directly test its immediate prediction: if the absence is a consequence of sampling, one should be able to predict occupancy of a species (the probability that a species is present) simply from its average and variance of abundance (together with the total number of reads of each sample). In particular, assuming a Gamma AFD, the occupancy of species i is given by

$$langle {o}_{i}rangle =1-frac{1}{T}sum _{s}P(0| {N}_{s})=1-frac{1}{T}sum _{s = 1}^{T}{left(1+frac{{bar{x}}_{i}{N}_{s}}{{beta }_{i}}right)}^{-{beta }_{i}},,$$
(9)

where Ns is the total number of reads in sample s, T is the total number of samples, and ({beta }_{i}={bar{x}}_{i}^{2}/{sigma }_{{x}_{i}}^{2}). As shown in Fig. 2 and in Supplementary Fig. 3, this prediction well reproduces the observed occupancy across species. The prediction of Eq. (12) also matches the occupancy of temporal (longitudinal) data Supplementary Fig. 20.
The second, more rigorous, way to test the hypothesis that (most) species are always present is to use model selection. In this context we want to compare two (or more) models that aim at describing the observed number of reads of each species starting from alternative hypothesis. In particular I compare a purely Gamma AFD with a zero-inflated Gamma, which reads

$${varrho }_{i}(x| {vartheta }_{i},{beta }_{i},{bar{x}}_{i})={vartheta }_{i}delta (x)+(1-{vartheta }_{i})frac{1}{Gamma ({beta }_{i})}{left(frac{{beta }_{i}}{{bar{x}}_{i}}right)}^{{beta }_{i}}{x}^{{beta }_{i}-1}exp left(-{beta }_{i}frac{x}{{bar{x}}_{i}}right),,$$
(10)

where ϑi is the probability that a species is truly absent in a community and δ( ⋅ ) is the Dirac delta distribution. Our goal is to test whether the ϑis are significantly different from zero. Since the two models are nested, one can compare the maximum likelihood estimator in the case ϑi = 0 with the (maximum) likelihood marginalized over ϑ (which has prior μ(ϑ)). Given the number of reads ({n}_{i}^{s}) of species i in community s, with Ns total number of reads, one can compute the ratio (Supplementary Note 4)

$${ell }_{i}=frac{mathop{max }nolimits_{bar{x},beta }{prod }_{s}int dx{varrho }_{i}(x| 0,beta ,bar{x})frac{{(x{N}_{s})}^{{n}_{i}^{s}}}{{n}_{i}^{s}!}{e}^{-x{N}_{s}}}{int dvartheta ,mu (vartheta )left(mathop{max }nolimits_{bar{x},beta }{prod }_{s}int dx{varrho }_{i}(x| vartheta ,beta ,bar{x})frac{{(x{N}_{s})}^{{n}_{i}^{s}}}{{n}_{i}^{s}!}{e}^{-x{N}_{s}}right)},,$$
(11)

where μ(ϑ) is a prior over ϑ. If ℓi  > 1, the model with ϑi = 0 is more strongly supported than the model with ϑ ≠ 0. Under Beta prior with parameters 0.25 and 8, one obtains that ℓi  > 1 in 98.8% of the cases (averaged across data sets, ranging from 94.4 to 99.7%) and ℓi  > 100 in 97.5% cases (ranging from 92.8 to 99.2%). See Supplementary Note 4 for a more detailed description of the methodology and Supplementary Fig. 6 for results obtained with other priors.
Prediction of macroecological patterns
Given laws #1, #2, and #3, the probability to observe n reads of a randomly chosen species in a sample with N total reads is

$$P(n| N)=int_{-infty }^{infty }deta ,frac{Gamma (beta +n)}{n!Gamma (beta )}{left(frac{{e}^{eta }N}{beta +{e}^{eta }N}right)}^{n}{left(frac{beta }{beta +{e}^{eta }N}right)}^{beta }frac{exp left(-frac{{(eta -mu )}^{2}}{2{sigma }^{2}}right)}{sqrt{2pi {sigma }^{2}}},,$$
(12)

where (eta ={mathrm{log}},(bar{x})). All the properties of species are fully specified by its mean abundance (bar{x}={e}^{eta }). The probability of observing n reads of species with average abundance (bar{x}) in a sample with N total number of reads is therefore

$$P(n| N,bar{x})=frac{Gamma (beta +n)}{n!Gamma (beta )}{left(frac{bar{x}N}{beta +bar{x}N}right)}^{n}{left(frac{beta }{beta +bar{x}N}right)}^{beta },.$$
(13)

The predictions for the patterns shown in Fig. 3 are reported here. The full derivation of this and other patterns is presented in Supplementary Note 8.
The total number of observed species in a sample with N total number of reads can be easily calculated using Eq. (12). The probability of not observing a species is simply P(0∣N). The expected number of distinct species 〈s(N)〉 in a sample with N reads is therefore

$$langle s(N)rangle ={s}_{tot}left(1-P(0| N)right)={s}_{tot}left(1-int_{-infty }^{infty }deta ,frac{exp left(-frac{{(eta -mu )}^{2}}{2{sigma }^{2}}right)}{sqrt{2pi {sigma }^{2}}},{left(frac{beta }{beta +{e}^{eta }N}right)}^{beta }right),,$$
(14)

where stot is the total number of species in the biome (including unobserved ones, see Supplementary Note 7). Note that stot is (substantially) larger than sobs, the number of different species observed in the union of all the communities, which can instead be written as

$$langle {s}_{obs}rangle ={s}_{tot}left(1-int_{-infty }^{infty }deta ,frac{exp left(-frac{{(eta -mu )}^{2}}{2{sigma }^{2}}right)}{sqrt{2pi {sigma }^{2}}},{left(prod _{s = 1}^{T}frac{beta }{beta +{e}^{eta }{N}_{s}}right)}^{beta }right),.$$
(15)

Figure 3a shows that the prediction of Eq. (14) correctly matches the data (Supplementary Fig. 13).
The SAD, one of the most studied patterns in ecology and directly related to the Relative Species Abundance35, is defined as the fraction of species with a given abundance. According to our model, the expected SAD is given by

$$langle {Phi }_{n}(N)rangle := frac{langle {s}_{n}(N)rangle }{langle s(N)rangle }=frac{P(n| N)}{1-P(0,N)}=frac{int_{-infty }^{infty }deta ,frac{Gamma (beta +n)}{n!Gamma (beta )}{left(frac{{e}^{eta }N}{beta +{e}^{eta }N}right)}^{n}{left(frac{beta }{beta +{e}^{eta }N}right)}^{beta }frac{exp left(-frac{{(eta -mu )}^{2}}{2{sigma }^{2}}right)}{sqrt{2pi {sigma }^{2}}}}{1-int_{-infty }^{infty }deta ,{left(frac{beta }{beta +{e}^{eta }N}right)}^{beta },frac{exp left(-frac{{(eta -mu )}^{2}}{2{sigma }^{2}}right)}{sqrt{2pi {sigma }^{2}}}},,$$
(16)

where 〈sn(N)〉 is the number of species with n reads in a sample with N total number of reads. The cumulative SAD is defined as

$$langle {Phi }_{n}^{ ,{ > },}(N)rangle := sum _{m = n}^{infty }langle {Phi }_{m}(N)rangle =frac{int mathop{int}nolimits_{-infty }^{infty },{I}_{frac{{e}^{eta }N}{beta +{e}^{eta }N}}(n,beta )frac{exp left(-frac{{(eta -mu )}^{2}}{2{sigma }^{2}}right)}{sqrt{2pi {sigma }^{2}}}}{1-mathop{int}nolimits_{-infty }^{infty }eta ,{left(frac{beta }{beta +{e}^{eta }N}right)}^{beta },frac{exp left(-frac{{(eta -mu )}^{2}}{2{sigma }^{2}}right)}{sqrt{2pi {sigma }^{2}}}},,$$
(17)

where Ip(n, β) is the regularized incomplete Beta function. Figure 3b shows that the Eq. (17) captures the empirical cumulative SAD (Supplementary Fig. 17).
The occupancy probability is defined as the probability that a species is present in a given fraction of communities. This quantity has been extensively studied in a variety of contexts (from genomics71 to Lego sets and texts72) and has been more recently considered in microbial ecology37. The three macroecological laws predict (see derivation in Supplementary Note 8)

$${p}_{obs}(o)=frac{mathop{int}nolimits_{-infty }^{infty }deta ,sum_{t = 1}^{T}delta left(o-1+frac{1}{T}mathop{sum }nolimits_{s = 1}^{T}{left(frac{beta }{beta +{e}^{eta }{N}_{s}}right)}^{beta }right)frac{exp left(-frac{{(eta -mu )}^{2}}{2{sigma }^{2}}right)}{sqrt{2pi {sigma }^{2}}}mathop{prod }nolimits_{s = 1}^{T}left(1-{left(frac{beta }{beta +{e}^{eta }{N}_{s}}right)}^{beta }right)}{mathop{int}nolimits_{-infty }^{infty }deta ,frac{exp left(-frac{{(eta -mu )}^{2}}{2{sigma }^{2}}right)}{sqrt{2pi {sigma }^{2}}},mathop{prod }nolimits_{s = 1}^{T}left(1-{left(frac{beta }{beta +{e}^{eta }{N}_{s}}right)}^{beta }right)},,$$
(18)

where δ( ⋅ ) is a Dirac delta function. Figure 3c compares the prediction of Eq. (18) with the data (Supplementary Fig. 15).
Occupancy (the fraction of communities where a species is found) and abundance are not independent properties, and their relative dependence is often referred to as occupancy-abundance relationship21 Given an average (relative) abundance (bar{x}=exp (eta )), the expected occurrence is

$${langle orangle }_{eta }=1-frac{1}{T}sum _{s = 1}^{T}P(0| {N}_{s},bar{x})=1-frac{1}{T}sum _{s = 1}^{T}{left(frac{beta }{beta +bar{x}{N}_{s}}right)}^{beta },,$$
(19)

Figure 3d shows the comparison between data and predictions (Supplementary Fig. 16). These predictions are also tested for temporal (longitudinal) data in Supplementary Figs. 22–24.
Transition probabilities in longitudinal data
For longitudinal data, in addition to the stationary AFD, one can study the probability ({rho }_{i}(x^{prime} ,t+Delta t| x,t)) that a species i has abundance (x^{prime}) at time t + Δt conditioned on having abundance x at time t. Instead of focusing on the full distribution, we study its first two (conditional) central moments, i.e. the average and variance of the abundance at t + Δt conditioned to abundance x at time t. In the analysis of the data stationarity is assumed (the distribution ({rho }_{i}(x^{prime} ,t+Delta t| x,t)) depends on Δt but not on t). I also assume that the dynamics of different species are governed by similar equations that only differ in their parameters. One would like therefore to average over species, by properly rescaling their abundances. The average over species is potentially problematic, as it could add a spurious effect to the conditional averages. For instance, only species with larger fluctuations would appear for extreme values of the initial abundance. In order to avoid these problems, instead of consider the actual abundance, its cumulative probability distribution value (calculated using the empirical AFD of each species) was used, that is referred as “quantile abundance”. This is equivalent to rank the abundances of each species over communities and use the (relative) ranking of each community instead of the abundance. A value equal to 0 corresponds to the lowest observed abundance, and a value equal to 1 to the highest. By definition, the quantile abundance is always uniformly distributed.
Ruling out demographic stochasticity
Demographic stochasticity can reproduce a Gamma AFD. A birth, death, and immigration process has a Gamma as stationary distribution35. In the limit of large populations sizes, it corresponds to the following equation35

$$frac{dx}{dt}=m-(d-b)x+sqrt{(b+d)x}xi (t),,$$
(20)

where m is the migration rate, while b and d are the per-capita birth and death rate. The Gaussian white noise term ξ(t) has mean zero and time-correlation (langle xi (t)xi (t^{prime} )rangle =delta (t-t^{prime} )). The stationary distribution of this process turns out to be

$$rho (x)=frac{1}{Gamma left(2frac{m}{b+d}right)}{left(frac{b+d}{2(d-b)}right)}^{-2frac{m}{b+d}}{x}^{2frac{m}{b+d}-1}exp left(-2frac{d-b}{b+d}xright),.$$
(21)

The average abundance is equal to (bar{x}=m/(d-b)), while the variance turns out to be ({sigma }_{x}^{2}=(m/2)(b+d)/{(b-d)}^{2}). The square of the coefficient of variation would therefore be equal to (b + d)/(2m).
More generally, one can assume that all the parameters are species dependent, and the population of species i is described by

$$frac{d{x}_{i}}{dt}={m}_{i}-({d}_{i}-{b}_{i}){x}_{i}+sqrt{({b}_{i}+{d}_{i}){x}_{i}}{xi }_{i}(t),,$$
(22)

where (langle {xi }_{i}(t){xi }_{j}(t^{prime} )rangle ={delta }_{ij}delta (t-t^{prime} )) was assumed.
Taylor’s Law and the wide variation of average abundance together imply that mi/(bi + di) is constant while mi/(di − bi) varies across species on several orders of magnitudes. This imposes a constraint on the variation of parameter values across species.
For instance, one can consider the scenario where species migrate to local communities from a common species pool (metacommunity). As abundance in the metacommunity varies across species the migration rate is a species-dependent quantity. Under neutrality, the per-capita birth and death rates in the local communities are constant and independent of the identity of the species. In this case mi depends on the species, while b and d do not. One could recover the Lognormal MAD by imposing that mi is Lognormally distributed. On the other hand, this model would fail in reproducing Taylor’s law with exponent 2, as it would predict and exponent 1.
More in general, the condition imposed on the parameters corresponds to an unnatural fine-tuned relationship between migration, birth, and death rates. Variation of the average abundance is observed across, at least, 7 orders of magnitudes. In order to reproduce this variation across species and Taylor’s law with exponent 2, the range of variability of (bi − di)/(bi +  di) should be of the same order. It is unrealistic that the relative difference between birth and death rates, which have strong and direct connection to fundamental biological processes, vary so much across bacterial species. It is important to underline however, that the model of Eq. (22) can, in fact, for a proper parameterization, explain the observed variation of the data. But the choice of parameters explaining the empirical variation require for achieving this goal requires careful and unrealistic fine-tuning of the microscopic parameters.
Stochastic logistic model
The SLM is defined as

$$frac{d{x}_{i}}{dt}=frac{{x}_{i}}{{tau }_{i}}left(1-frac{{x}_{i}}{{K}_{i}}right)+sqrt{frac{{sigma }_{i}}{{tau }_{i}}}{x}_{i}{xi }_{i}(t),,$$
(23)

where ξ(t) is a Gaussian white noise term with mean zero and correlation (langle {xi }_{i}(t){xi }_{j}(t^{prime} )rangle ={delta }_{ij}delta (t-t^{prime} )). Taylor’s Law and the observed Lognormal MAD constraints the parameter value. The parameters 1/τi, Ki and σi are the intrinsic growth rate, the carrying capacity and the coefficient of variation of the growth-rate fluctuations. Taylor’s Law requires σi = σ (independently of i). Since the average abundance of the SLM is ({bar{x}}_{i}={K}_{i}(1-{sigma }_{i}/2)), if σi = σ, the average abundance and the carrying capacity turn out to be proportional to each other. The lognormal MAD implies therefore that the Kis are lognormally distributed. The stationary distribution corresponding to Eq. (23) reads

$${rho }_{i}(x)=frac{1}{Gamma (2{sigma }_{i}^{-1}-1)}{left(frac{2}{{K}_{i}{sigma }_{i}}right)}^{2{sigma }_{i}^{-1}-1}exp left(-frac{2}{{K}_{i}{sigma }_{i}}xright){x}^{2{sigma }_{i}^{-1}-2},.$$
(24)

The parameter τi does not affect stationary properties, but determines the timescale of relaxation to the stationary distribution. For small deviation of abundance from the average and for large times, the conditional expected abundance behaves as

$${langle {x}_{i}(t+Delta t)rangle }_{{x}_{i}(t)}={bar{x}}_{i}+left({x}_{i}(t)-{bar{x}}_{i}right){e}^{-frac{Delta t}{{tau }_{i}}},.$$
(25)

From the slopes of Fig. 4g one can then determine the timescales τi, which turn out to be approximately equal to 19 h. In Fig. 4 it was assumed τi = 19 h for all species.
Equation (23) can emerge as effective description of more complicated coupled equations. For instance, it is possible to show that a Lotka-Volterra system of equation with random interactions reduces to Eq. (23) (with colored noise to be self-consistently determined)42. If the coefficient of variation of the interaction coefficient does not increase with the number of species (e.g., if it is constant) then the Lotka-Volterra equations can be effectively approximated with Eq. (23).
The noise term in Eq. (23) can be interpreted as corresponding to environmental fluctuations. These fluctuations are typically known to have a characteristic timescale and are not white40,41. Supplementary Note 13 and Supplementary Fig. 25 show that colored noise in Eq. (23) does not affect significantly the predictions obtained with the SLM with white noise.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More