Ecology

1.
World Health Organisation. World Malaria Report 2018 (World Health Organization, Geneva, 2018).
2.
Caminade, C. et al. Impact of climate change on global malaria distribution. Proc. Natl Acad. Sci. 111, 3286–3291 (2014).
ADS  CAS  PubMed  Google Scholar 

3.
Shapiro, L. L., Whitehead, S. A. & Thomas, M. B. Quantifying the effects of temperature on mosquito and parasite traits that determine the transmission potential of human malaria. PLoS Biol. 15, 2003489 (2017).
Google Scholar 

4.
Jepson, W. F., Moutia, A. & Courtois, C. The malaria problem in Mauritius: the binomics of Mauritian anophelines. Bull. Entomol. Res. 38, 177–208 (1947).
CAS  PubMed  Google Scholar 

5.
Waite, J. L., Suh, E., Lynch, P. A. & Thomas, M. B. Exploring the lower thermal limits for development of the human malaria parasite, Plasmodium falciparum. Biol. Lett. 15, 20190275 (2019).
PubMed  PubMed Central  Google Scholar 

6.
Bayoh, M. N. & Lindsay, S. W. Effect of temperature on the development of the aquatic stages of Anopheles gambiae sensu stricto (Diptera: Culicidae). Bull. Entomol. Res. 93, 375–381 (2003).
CAS  PubMed  Google Scholar 

7.
Depinay, J. M. et al. A simulation model of African Anopheles ecology and population dynamics for the analysis of malaria transmission. Malar. J. 3, 29 (2004).
PubMed  PubMed Central  Google Scholar 

8.
Mordecai, E. A. et al. Optimal temperature for malaria transmission is dramatically lower than previously predicted. Ecol. Lett. 16, 22–30 (2012).
PubMed  Google Scholar 

9.
Craig, M. H., Snow, R. W. & Le Sueur, D. Climate-based distribution model of malaria transmission in sub-Saharan Africa. Parasitol. Today 15, 105–111 (1999).
CAS  PubMed  Google Scholar 

10.
Parham, P. E. & Michael, E. Modeling the effects of weather and climate change on malaria transmission. Environ. Health Perspect. 118, 620–626 (2010).
PubMed  Google Scholar 

11.
Ryan, S. J. et al. Mapping physiological suitability limits for malaria in Africa under climate change. Vector-Borne Zoonotic Dis. 15, 718–725 (2015).
PubMed  PubMed Central  Google Scholar 

12.
Smith, M. W., Macklin, M. G. & Thomas, C. J. Hydrological and geomorphological controls of malaria transmission. Earth Sci. Rev. 116, 109–127 (2013).
ADS  Google Scholar 

13.
Bomblies, A., Duchemin, J. B. & Eltahir, E. A. Hydrology of malaria: model development and application to a Sahelian village. Water Resour. Res. 44, W12445 (2008).
ADS  Google Scholar 

14.
Yamana, T. K., Bomblies, A. & Eltahir, E. A. Climate change unlikely to increase malaria burden in West Africa. Nat. Clim. Change 6, 1009–1013 (2016).
ADS  Google Scholar 

15.
Small, J., Goetz, S. J. & Hay, S. I. Climatic suitability for malaria transmission in Africa, 1911–1995. Proc. Natl Acad. Sci. USA 100, 15341–15345 (2003).
ADS  CAS  PubMed  Google Scholar 

16.
Thomas, C. J., Davies, G. & Dunn, C. E. Mixed picture for changes in stable malaria distribution with future climate in Africa. Trends Parasitol. 20, 216–220 (2004).
PubMed  Google Scholar 

17.
Ebi, K. L. et al. Climate suitability for stable malaria transmission in Zimbabwe under different climate change scenarios. Clim. Change 73, 375–393 (2005).
ADS  Google Scholar 

18.
Van Lieshout, M., Kovats, R. S., Livermore, M. T. J. & Martens, P. Climate change and malaria: analysis of the SRES climate and socio-economic scenarios. Glob. Environ. Change 14, 87–99 (2004).
Google Scholar 

19.
Kiszewski, A. et al. A global index representing the stability of malaria transmission. Am. J. Trop. Med. Hyg. 70, 486–498 (2004).
PubMed  Google Scholar 

20.
Ermert, V., Fink, A. H., Jones, A. E. & Morse, A. P. Development of a new version of the Liverpool Malaria Model. II. Calibration and validation for West Africa. Malar. J. 10, 62 (2011).
PubMed  PubMed Central  Google Scholar 

21.
Martens, W. J. M., Niessen, L. W., Rotmans, J. & McMichael, A. J. Potential impacts of global climate change on malaria risk. Environ. Health Perspect. 103, 458–464 (1995).
CAS  PubMed  PubMed Central  Google Scholar 

22.
Tanser, F., Sharp, B. L. & Le Sueur, D. Potential effect of climate change on malaria transmission in Africa. Lancet 362, 1792–9178 (2003).
PubMed  Google Scholar 

23.
Garnham, P. C. C. The incidence of malaria at high altitudes. J. Malar. Soc. 7, 275–284 (1948).
CAS  Google Scholar 

24.
Lindsay, S. W., Parson, L. & Thomas, C. J. Mapping the ranges and relative abundance of the two principal African malaria vectors, Anopheles gambiae sensu stricto and An. arabiensis, using climate data. Proc. R. Soc. B Biol. Sci. 265, 847–854 (1998).
CAS  Google Scholar 

25.
Mabaso, M. L., Craig, M., Ross, A. & Smith, T. Environmental predictors of the seasonality of malaria transmission in Africa: the challenge. Am. J. Trop. Med. Hyg. 76, 33–38 (2007).
PubMed  Google Scholar 

26.
Kibret, S. et al. Malaria impact of large dams in sub-Saharan Africa: maps, estimates and predictions. Malar. J. 14, 339 (2015).
PubMed  PubMed Central  Google Scholar 

27.
Lysenko, A. J. & Semashko, I. N. in Itogi Nauki: Medicinskaja Geografija (ed. Lebedew, A. W.) 25–146 (Academy of Sciences, Moscow, 1968) (in Russian).

28.
Ageep, T. B. et al. Spatial and temporal distribution of the malaria mosquito Anopheles arabiensis in northern Sudan: influence of environmental factors and implications for vector control. Malar. J. 8, 123 (2009).
PubMed  PubMed Central  Google Scholar 

29.
Fuller, D. O., Parenti, M. S., Hassan, A. N. & Beier, J. C. Linking land cover and species distribution models to project potential ranges of malaria vectors: an example using Anopheles arabiensis in Sudan and Upper Egypt. Malar. J. 11, 264 (2012).
PubMed  PubMed Central  Google Scholar 

30.
Gemperli, A., Vounatsou, P., Sogoba, N. & Smith, T. Malaria mapping using transmission models: application to survey data from Mali. Am. J. Epidemiol. 163, 289–297 (2005).
PubMed  Google Scholar 

31.
Snow, R. W., Noor, A. M. & Hay, S. I. Malaria in Somalia: Assembling the Evidence and Modeling Risks (University of Oxford, UK, 2006).

32.
Hulme, M., Doherty, R., Ngara, T., New, M. & Lister, D. African climate change: 1900–2100. Clim. Res. 17, 145–168 (2001).
Google Scholar 

33.
Tierney, J. E., Ummenhofer, C. C. & de Menocal, P. B. Past and future rainfall in the Horn of Africa. Sci. Adv. 1, e1500682 (2015).
ADS  PubMed  PubMed Central  Google Scholar 

34.
Reiter, P. Global warming and malaria: knowing the horse before hitching the cart. Malar. J. 7, S3 (2008).
PubMed  PubMed Central  Google Scholar 

35.
Gething, P. et al. Climate change and the global malaria recession. Nature 465, 342–345 (2010).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

36.
Yamazaki, D., Kanae, S., Kim, H. & Oki, T. A physically based description of floodplain inundation dynamics in a global river routing model. Water Resour. Res. 47, W04501 (2011).
ADS  Google Scholar 

37.
Hazeleger, W. et al. EC-Earth V2.2: description and validation of a new seamless earth system prediction model. Clim. Dyn. 39, 2611–2629 (2012).
Google Scholar 

38.
Alfieri, L. et al. Global projections of river flood risk in a warmer world. Earth’s Future 5, 171–182 (2017).
ADS  Google Scholar 

39.
Muerth, M. J. et al. On the need for bias correction in regional climate scenarios to assess climate change impacts on river runoff. Hydrol. Earth Syst. Sci. 17, 1189–1204 (2013).
ADS  Google Scholar 

40.
Hirpa, F. A. et al. Streamflow response to climate change in the Greater Horn of Africa. Clim. Change 156, 341–363 (2019).
ADS  Google Scholar 

41.
van der Knijff, J. M., Younis, J. & de Roo, A. P. J. LISFLOOD: a GIS-based distributed model for river basin scale water balance and flood simulation. Int. J. Geograph. Inf. Sci. 24, 189–212 (2010).
Google Scholar 

42.
Burek, P., van der Knijff, J. & de Roo, A. P. J. LISFLOOD, Distributed Water Balance and Flood Simulation Model Revised User Manual (Publ. Off., Luxembourg, 2013).

43.
Lehner, B., Verdin, K. & Jarvis, A. New global hydrography derived from spaceborne elevation data. Eos, Trans. Am. Geophys. Union 89, 93–94 (2008).
ADS  Google Scholar 

44.
Wu, H. et al. A new global river network database for macroscale hydrologic modeling. Water Resour. Res. 48, W09701 (2012).
ADS  Google Scholar 

45.
Hengl, T. et al. SoilGrids1km—global soil information based on automated mapping. PLoS ONE 9, e105992 (2014).
ADS  PubMed  PubMed Central  Google Scholar 

46.
Paaijmans, K. P., Takken, W., Githeko, A. K. & Jacobs, A. F. G. The effect of water turbidity on the near-surface water temperature of larval habitats of the malaria mosquito Anopheles gambiae. Int. J. Biometeorol. 52, 747–753 (2008).
ADS  CAS  PubMed  Google Scholar 

47.
Thomas, C. J., Cross, D. E. & Bøgh, C. Landscape movements of Anopheles gambiae malaria vector mosquitoes in rural Gambia. PLoS ONE 8, e68679 (2013).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

48.
Worldpop. Africa Continental Population Datasets (2000–2020) v2.0. https://doi.org/10.5258/SOTON/WP00004 (2016).

49.
Worldpop. Africa Continental age/sex structure Population Datasets 2000/05/10/15/20 V5.0. https://www.worldpop.org/geodata/summary?id=1276 (2016).

50.
James, W. H. et al. Gridded birth and pregnancy datasets for Africa, Latin America and the Caribbean. Sci. Data 5, 180090 (2018).
CAS  PubMed  PubMed Central  Google Scholar 

51.
United Nations, Department of Economic and Social Affairs, Population Division. World Population Prospects 2019 Online ed (United Nations, 2019).

52.
Smith, M. W. et al. LIS-MAL Estimates of Hydro-Climatic Suitability for Malaria Transmission in Africa (1971–2100). [Dataset]. https://doi.org/10.5518/786 (University of Leeds, 2020). More

According to an initial keyword search, we selected the 17 most popular microbial ecology-related journals, as these were more likely to have sequence-specific data deposition instructions or requirements. We surveyed all the articles published in these journals between January 2015 and March 2019 (n = 26,927 articles, Supplementary Table S1), as concerns over data deposition practices began to grow in 201514 and were soon followed by stricter standards for data availability12. A custom-built pattern-based text extraction algorithm followed by manual curation, we selected those studies which performed 16S rRNA gene amplicon sequencing and listed INSDC-compliant accession numbers (n = 2015, Supplementary Table S1; 145,203 samples).
To confirm that our parsing algorithm did not miss accession numbers in articles containing 16S rRNA gene amplicon sequencing, we randomly selected 150 articles which mentioned 16S rRNA, but for which no accession numbers were detected, for manual inspection. Of these, one contained a misspelled accession number, two had archived their sequences in unconventional repositories (Google Drive and GEO, a gene expression database, Supplementary Data 2), and 19 were identified as having performed 16S rRNA gene amplicon sequencing, but had not included any reference to the data. We found no cases in which accession numbers or sequence data were stored in supplementary materials. From this group, we estimate that 18% of the studies in our database (n = 469) performed 16S rRNA gene amplicon sequencing but did not provide access to the data (Fig. 1a). Four studies mentioned deposition data in dbGaP18, and we could verify the existence of three of these studies. We found that an additional 6.5% of the studies had deposited their data in the Qiita19, MG-RAST20, and figshare databases (n = 14, n = 134, and n = 24 studies, respectively). Of the estimated 2,656 studies employing 16S rRNA gene amplicon sequencing, 75.9% deposited their data to an INSDC database in the period studied (Fig. 1a).
Fig. 1: Popular locations for data storage.

Data for all studies which contained 16S rRNA amplicon sequencing (a), and the V3–V4 subset (b); n = 2656 and n = 635 studies, respectively. For the entirety of the study, studies which contained amplicon sequences but did not deposit them were inferred by manually checking 150 randomly-selected articles which did not contain INSDC accession numbers or refer to alternative databases, indicated in lighter yellow. For the V3–V4 subset, studies which contained the keywords “16S rRNA”, “515”, and “806” were selected. Studies for which INSDC-compliant accession numbers were reported but which did not exist on any INSDC database are shown in lighter blue.

Full size image

To obtain more precise estimates of the percentage of articles which deposited their data in each database, we focused on the subset of 635 studies which sequenced the V3–V4 region of the 16S rRNA gene between base pairs 515 and 806 (heretofore V3–V4 subset), a target region which has gained popularity since its development and use by the Earth Microbiome Project21,22. Of these, 74.5% (n = 474) studies listed INSDC-compliant accession numbers within the article, but of these, accession numbers from 5% of the studies (n = 33) were not findable on any INSDC database. Additionally, 19% (n = 121) did not provide an identifiable link to the data, and 6.8% of the studies deposited their data in the Qiita, MG-RAST, and figshare databases (n = 9, n = 24, n = 7, respectively, Fig. 1b). Two studies provided SRA submission IDs rather than accession numbers, and were also inaccessible.
The increasing popularity of microbial community sequencing was evident in our data. Over the period studied, the number of studies in the V3–V4 subset rose from 56 in 2015 to 214 in 2018 (Supplementary Fig. 1a). The proportion of publications which claimed to deposit data to INSDC databases increased slightly over time, from 33/56 in 2015 to 172/214 in 2018 (χ2 = 6.6, p = 0.01, Supplementary Fig. 1b), suggesting an increasing tendency towards deposition in INSDC databases. Deposition to alternative databases decreased (χ2 = 14.04, p  More

1.
Gómez-Robles, A. Dental evolutionary rates and its implications for the Neanderthal–modern human divergence. Sci. Adv. 5, eaaw1268 (2019).
ADS  PubMed  PubMed Central  Google Scholar 
2.
Higham, T. et al. The timing and spatiotemporal patterning of Neanderthal disappearance. Nature 512, 306–309 (2014).
ADS  CAS  PubMed  Google Scholar 

3.
Wenzel, S. Neanderthal presence and behaviour in central and Northwestern Europe during MIS. In Developments in Quaternary Sciences Vol. 7 (eds Sirocko, F. et al.) 173–193 (Elsevier, Amsterdam, 2007).
Google Scholar 

4.
Sørensen, B. Demography and the extinction of European Neanderthals. J. Anthropol. Archaeol. 30, 17–29 (2011).
Google Scholar 

5.
Salazar-García, D. C. et al. Neanderthal diets in central and southeastern Mediterranean Iberia. Quat. Int. 318, 3–18 (2013).
Google Scholar 

6.
Wißing, C. et al. Isotopic evidence for dietary ecology of late Neandertals in North-Western Europe. Quat. Int. 411, 327–345 (2016).
Google Scholar 

7.
Nielsen, T. K. et al. Investigating Neanderthal dispersal above 55°N in Europe during the Last Interglacial Complex. Quat. Int. 431, 88–103 (2017).
Google Scholar 

8.
Nicholson, C. M. Eemian paleoclimate zones and Neanderthal landscape-use: a GIS model of settlement patterning during the last interglacial. Quat. Int. 438, 144–157 (2017).
Google Scholar 

9.
Rhodes, S. E., Starkovich, B. M. & Conard, N. J. Did climate determine Late Pleistocene settlement dynamics in the Ach Valley, SW Germany?. PLoS ONE 14, e0215172 (2019).
PubMed  PubMed Central  Google Scholar 

10.
Defleur, A. R. & Desclaux, E. Impact of the last interglacial climate change on ecosystems and Neanderthals behavior at Baume Moula-Guercy, Ardèche, France. J. Archaeol. Sci. 104, 114–124 (2019).
Google Scholar 

11.
Hovers, E. Territorial Behavior in the Middle Paleolithic of the Southern Levant. In: Settlement dynamics of the Middle Paleolithic and Middle Stone Age, Conard, N. J., editor. Kerns Verlag T¨ubingen; 123–152 (2001).

12.
Shea, J. J. Neandertals, competition, and the origin of modern human behavior in the Levant. Evol. Anthropol. 12, 173–187 (2003).
Google Scholar 

13.
Been, E. et al. The first Neanderthal remains from an open-air Middle Palaeolithic site in the Levant. Sci. Rep. 7, 2958 (2017).
ADS  PubMed  PubMed Central  Google Scholar 

14.
Vahdati Nasab, H., Clark, G. A. & Torkamandi, S. Late Pleistocene dispersal corridors across the Iranian Plateau: a case study from Mirak, a Middle Paleolithic site on the Northern Edge of the Iranian Central Desert (Dasht-e Kavir). Quat. Int. 300, 267–281 (2013).
Google Scholar 

15.
Heydari-Guran, S. Palaeolithic Landscapes of Iran. BAR International Series, 2568 (2014).

16.
Bazgir, B. et al. Understanding the emergence of modern humans and the disappearance of Neanderthals: Insights from Kaldar Cave (Khorramabad Valley, Western Iran). Sci. Rep. 7, 43460 (2017).
ADS  PubMed  PubMed Central  Google Scholar 

17.
Smith, P. E. L. Paleolithic Archaeology in Iran (University of Pennsylvania Press, Philadelphia, 1986).
Google Scholar 

18.
Heydari-Guran, S. & Ghasidian, E. The MUP Zagros Project: tracking the Middle-Upper Palaeolithic transition in the Kermanshah region, West-Central Zagros, Iran. Antiquity 91, 1–7 (2017).
Google Scholar 

19.
Heydari-Guran, S. & Ghasidian, E. Late Pleistocene hominin settlement patterns and population dynamics in the Zagros Mountains: Kermanshah region. Archaeol. Res. Asia 21, 100161 (2020).
Google Scholar 

20.
Varela, S., Lobo, J. M. & Hortal, J. Using species distribution models in paleobiogeography: a matter of data, predictors and concepts. Palaeogeogr. Palaeoclimatol. Palaeoecol. 310, 451–463 (2011).
Google Scholar 

21.
Svenning, J. C., Fløjgaard, C., Marske, K. A., Nógues-Bravo, D. & Normand, S. Applications of species distribution modeling to paleobiology. Quat. Sci. Rev. 30, 2930–2947 (2011).
ADS  Google Scholar 

22.
Guisan, A., Thuiller, W. & Zimmermann, N. E. Habitat Suitability and Distribution Models: With Applications in R (Cambridge University Press, Cambridge, 2017).
Google Scholar 

23.
Araújo, M. B. et al. Standards for distribution models in biodiversity assessments. Sci. Adv. 5, eaat858 (2019).
Google Scholar 

24.
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).
Google Scholar 

25.
Warren, D. L., Glor, R. E. & Turelli, M. ENMTools: a toolbox for comparative studies of environmental niche models. Ecography 33, 607–611 (2010).
Google Scholar 

26.
Banks, W. E. et al. Neanderthal extinction by competitive exclusion. PLoS ONE 3, e3972 (2008).
ADS  PubMed  PubMed Central  Google Scholar 

27.
Franklin, J., Potts, A. J., Fisher, E. C., Cowling, R. M. & Marean, C. M. Paleodistribution modeling in archaeology and paleoanthropology. Quat. Sci. Rev. 110, 1–14 (2015).
Google Scholar 

28.
Benito, B. M. et al. The ecological niche and distribution of Neanderthals during the Last Interglacial. J. Biogeogr. 44, 51–61 (2017).
Google Scholar 

29.
Giampoudakis, K. et al. Niche dynamics of Palaeolithic modern humans during the settlement of the Palaearctic. Glob. Ecol. Biogeogr. 26, 359–370 (2017).
Google Scholar 

30.
Marean, C. W. A critique of the evidence for scavenging by Neandertals and early modern humans: new data from Kobeh Cave (Zagros Mountains, Iran) and Die Kelders Cave 1 layer 10 (South Africa). J. Hum. Evol. 35, 111–136 (1998).
CAS  PubMed  Google Scholar 

31.
Marean, C. W. & Kim, S. Y. Mousterian large-mammal remains from Kobeh Cave: behavioural implications for Neanderthals and early modern humans. Curr. Anthropol. 39, S79–S114 (1998).
Google Scholar 

32.
Mashkour, M. et al. Carnivores and their prey in the Wezmeh Cave (Kermanshah, Iran): a Late Pleistocene refuge in the Zagros. Int. J. Osteoarchaeol. 19, 678–694 (2009).
Google Scholar 

33.
Biglari, F. et al. Qaleh Bozi, new evidence of late middle paleolithic occupation in the Zayandeh-Rud Basin, Esfahan Province. Archaeol. Res. Iran 7, 7–26 (2015).
Google Scholar 

34.
Pearson, R. G. & Dawson, T. P. Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful?. Glob. Ecol. Biogeogr. 12, 361–371 (2003).
Google Scholar 

35.
Melchionna, M. et al. Fragmentation of Neanderthals’ pre-extinction distribution by climate change. Palaeogeogr. Palaeoclimatol. Palaeoecol. 496, 146–154 (2018).
Google Scholar 

36.
Hole, F. & Flannery, K. V. The prehistory of Southwestern Iran: a preliminary report. Proc. Prehist. Soc. 22, 147–206 (1967).
Google Scholar 

37.
Roustaei, K. et al. Recent paleolithic surveys in Lurestan. Curr. Anthropol. 45, 692–707 (2004).
Google Scholar 

38.
Diedrich, C. Ice age spotted hyenas as Neanderthal exhumers and scavengers in Europe. Chronicles Sci. 1, 1–34 (2014).
Google Scholar 

39.
Banks, W. E. The application of ecological niche modeling methods to archaeological data in order to examine culture-environment relationships and cultural trajectories. Quaternaire 28, 271–276 (2017).
Google Scholar 

40.
Banks, W. E. et al. Eco-cultural niches of the Badegoulian: unraveling links between cultural adaptation and ecology during the Last Glacial Maximum in France. J. Anthropol. Archaeol. 30, 359–374 (2011).
Google Scholar 

41.
King, T. et al. Azokh cave hominin remains. In Azokh cave and the transcaucasian corridor, vertebrate paleobiology and paleoanthropology (eds Fernndez-Jalvo, Y. et al.) 103–106 (Springer, Dordrecht, 2016).
Google Scholar 

42.
Solecki, R. S. Prehistory in Shanidar valley, Northern Iraq. Science 139, 179–193 (1963).
ADS  CAS  PubMed  Google Scholar 

43.
Pomeroy, E. et al. Newly-discovered Neanderthal remains from Shanidar Cave, Iraqi Kurdistan, and their attribution to Shanidar 5. J. Hum. Evol. 111, 102–118 (2017).
PubMed  Google Scholar 

44.
Zanolli, C. et al. A Neanderthal from the Central Western Zagros, Iran. Structural reassessment of the Wezmeh 1 maxillary premolar. J. Hum. Evol. 135, 102643 (2019).
PubMed  Google Scholar 

45.
Coon, C. S. The seven caves. In: Archaeological Exploration in the Middle East. (Alfred Knopf, 1975).

46.
Biglari, F. & Heydari, S. Do-Ashkaft: a recent discovery Mousterian cave site in Kermanshah Plain, Iran. Antiquity 75, 487–488 (2001).
Google Scholar 

47.
Conard, N. J., Ghasidian, E. & Heydari, S. The paleolithic of Iran. In Ancient Iran (ed. Potts, D. T.) 29–48 (Oxford Press, Oxford, 2013).
Google Scholar 

48.
Power, R. C. Neanderthals and their diet. In eLS (Wiley, 2019).

49.
Perkins, D. Prehistoric fauna from Shanidar, Iraq. Science 144, 1565–1566 (1964).
ADS  PubMed  Google Scholar 

50.
Turnbull, P. F. The mammalian fauna of Warwasi rock shelter, west-central Iran. Fieldiana Geol. 33, 141–155 (1975).
Google Scholar 

51.
Hesse, B. Paleolithic Faunal Remains from Ghar-I-Khar, Western Iran (University of Alabama Press, Alabama, 1989).
Google Scholar 

52.
Karami, M., Ghadirian, T. & Faizolahi, K. The Atlas of the Mammals of Iran (Jahad Daneshgahi Press, Tehran, 2016).
Google Scholar 

53.
Naderi, S. Evolutionary history of wild goat (Capra aegagrus) and the goat (C. hircus) based on the analysis of mitochondrial and nuclear DNA polymorphism: Implications for conservation and for the origin of the domestication. Ecology, Environment. Grenoble: Université Joseph-Fourier – Grenoble I (2007).

54.
Yusefi, G. H., Faizolahi, K., Darvish, J., Safi, K. & Brito, J. C. The species diversity, distribution, and conservation status of the terrestrial mammals of Iran. J. Mammal. 100, 55–71 (2019).
Google Scholar 

55.
Hijmans, R. J. Raster: Geographic Data Analysis and Modeling. R package, (2015).

56.
R Core Team. R: A Language and Environment for Statistical Computing (R Core Team, Vienna, Austria, 2017).

57.
Jarvis, A., Reuter, H. I., Nelson, A. & Guevara, E. Hole-filled SRTM for the globe Version 4. Available from the CGIAR-CSI SRTM 90m Database https://srtm.csi.cgiar.org on 15 April 2015 (2008).

58.
Quinn, G. P. & Keough, M. J. Experimental Designs and Data Analysis for Biologists (Cambridge University Press, Cambridge, 2002).
Google Scholar 

59.
Naimi, B. Uncertainty Analysis for Species Distribution Models. R package version 1.1–15 (2015).

60.
Araújo, M. B. & New, M. Ensemble forecasting of species distributions. Trends Ecol. Evol. 22, 42–47 (2007).
PubMed  Google Scholar 

61.
McCullagh, P. & Nelder, J. A. Generalized Linear Models (Chapman and Hall, London, 1989).
Google Scholar 

62.
Hastie, T. J. & Tibshirani, R. Generalized Additive Models (Chapman and Hall, London, 1990).
Google Scholar 

63.
Ridgeway, G. The state of boosting. J. Comput. Sci. 31, 172–181 (1999).
Google Scholar 

64.
Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231–259 (2006).
Google Scholar 

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

66.
Fielding, A. H. & Bell, J. F. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 24, 38–49 (1997).
Google Scholar 

67.
Boyce, M. S., Vernier, P. R., Nielsen, S. E. & Schmiegelow, F. K. A. Evaluating resource selection functions. Ecol. Model. 157, 281–300 (2002).
Google Scholar 

68.
Le Hirzel, A. H., Lay, G., Helfer, V., Randin, C. & Guisan, A. Evaluating the ability of habitat suitability models to predict species presences. Ecol. Model. 199, 142–152 (2006).
Google Scholar 

69.
Swets, J. A. Measuring the accuracy of diagnostic systems. Science 240, 1285–1293 (1988).
ADS  MathSciNet  CAS  MATH  Google Scholar 

70.
Schoener, T. W. Anolis lizards of Bimini: resource partitioning in a complex fauna. Ecology 49, 704–726 (1968).
Google Scholar  More

Spatial regime shifts over time and space
To assess how the dominance of large predatory fish (Eurasian perch Perca fluviatilis, northern pike Esox lucius) and three-spined stickleback (Gasterosteus aculeatus) changed over time and space, we collated juvenile fish abundance data from 13073 samplings conducted during 39 years (1979–2017) in 486 bays along a 1200 km stretch of the Swedish Baltic Sea coast. We used juvenile fish surveys because they (i) include stickleback, (ii) were conducted along the entire coast, including the outer archipelago, and (iii) capture patterns of recruitment failure; a proposed driver of local perch and pike stock decline34,35,37,45. The samplings were conducted by various monitoring programs and research projects to quantify fish recruitment. Much of the data was extracted from the Swedish national database for coastal fish (http://www.slu.se/kul). Other fish species also occurred in the data, but we—like others25,34,45—focused on the most common and strongly interacting species.
The timing and placement of most of the fish surveys was not chosen with this study in mind. We therefore included data from as many surveys as possible, ensuring that the dataset covered (i) gradients in distance to open sea and wave exposure, (ii) the entire Swedish east coast, and (iii) the longest time period possible (Supplementary Fig. 7). Nearly all sampling (94%) was conducted during July-September. To achieve the best possible spatial coverage we also included some bays only sampled in October and (for one bay) June. Initial exploration of data from bays sampled monthly from June to October suggested there were no large differences in the October and June data. Nine of the 486 bays occurred much further into the archipelago (49–67 km) than the rest (20 km from the open sea prior to 1995, we only used data from 1995–2017. Finally, we used a general linear multiple regression model to test how sampling year, wave exposure (two levels) and their interaction explained the maximum distance.
The bimodality in relative predator dominance seen over time and space (Figs. 1–3) could in theory be caused by variability in local abiotic conditions, and not by predator–prey interactions. We therefore tested whether any of four abiotic conditions estimated locally—water surface temperature, salinity, turbidity, water depth—could explain the variability in deviance residuals from the binomial glm (Supplementary Fig. 8), using a regular linear model (after assessing model assumptions and ruling out multicollinearity; see above). Salinity and turbidity (but not depth and temperature) had statistically significant but weak influences (R2 = 0.11), and the bimodality clearly remained (Supplementary Fig. 9).
Finally, we explored how time, distance to open sea, wave exposure, latitude and their two- and three-way interactions influenced log-transformed abundances of (i) predators (perch and pike pooled) and (ii) stickleback, using general linear models. We identified the most parsimonious models as outlined above.
Temporal regime shift at Forsmark: using the 34-year Forsmark time series, we first tested for temporal breakpoints in logit-transformed relative predator dominance data using change-point detection (strucchange) for linear models47. This method estimates the optimal number and (if identified) position of breakpoints using the Bayesian Information Criterion (BIC). Second, we explored what temporal change(s) in perch and stickleback abundances that preceded the shift. Because of highly non-linear patterns we modeled the temporal changes using generalized additive models (GAM) as implemented in the mgvc package68.
Importance of predator–prey reversal for fish recruitment: to assess the relative importance of predator–prey reversal for perch, pike and stickleback recruitment, we used statistical model selection based on path analyses; a form of structural equation modeling that can be used to tease apart direct vs. indirect (mediated) relationships between multiple ( >2) variables, and thereby assess the relative importance of direct vs. indirect relationships in systems48,69. Initial data exploration using multiple regression showed that one bay was a clear outlier due to 0 juvenile perch and pike, generating (i) too high leverage (influence on statistical relationships), (ii) heteroscedasticity and (iii) non-normally distributed errors. Since we suspected that juveniles had already migrated out of this bay, the bay was excluded (resulting in N = 31). Removing this statistical outlier resulted in that the model fulfilled test assumptions and the overall fit more than doubled (adjusted R2 increased from 0.17 to 0.37).
Based on ecological knowledge of the study system, we expressed 14 multivariate hypotheses of the direct and indirect drivers of perch and stickleback recruitment, as graphical network models of interacting paths36,49. Due to the relatively low sample size we restricted the number of paths to 7. The two simplest models assumed that perch+pike and stickleback juvenile abundance in summer (i.e. recruitment) was influenced by adult abundance in spring (stock recruitment) and cumulative cover of rooted vegetation, while adult abundance was explained by spring cumulative cover of all vegetation species36,46, and distance to open sea or wave exposure46. The more complex models included combinations of known predator–prey interactions: perch and pike controlling adult stickleback in spring through predation36, stickleback feeding on juvenile perch and pike45, and stickleback competing with juvenile perch for zooplankton prey37. We then analyzed each model using piecewise path analysis as implemented in the piecewiseSEM package48. First, we tested the goodness of fit of each model to the data using Shipley’s test of directional separation (D-sep)70. If missing paths were identified, they were included in a new model. For the models that fitted the data (p  > 0.05) we used the Akaike’s Information Criterion corrected for small samples (AICc) (calculated using Shipley’s general approach to calculate AIC for path analysis71) to compare relative model fit. A summary of all candidate models and details of the best-fitting model are presented in Supplementary Tables 1 and 2, respectively. The strength of paths in the best-fitting model are presented using standardized path coefficients, which (based on the best-fitting model, see Fig. 3) mean that 1 SD increase in pooled adult perch and pike abundance reduces adult stickleback abundance by 0.53 SD. We also calculated the amount of variation (R2) in adult stickleback abundance, pooled juvenile perch and pike, and juvenile stickleback abundance, that was explained by the paths. Finally, we tested whether the relative predator dominance of adult and juvenile fish in this smaller dataset (N = 31) was also uni- or bimodal, using Hartigan’s dip test (for details, see above).
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

A convincing argument that weather influences COVID-19 can be formulated in three parts: (1) experimental data suggest SARS-CoV-2 persistence on surfaces or in the air is sensitive to temperature, humidity, and ultraviolet light; (2) other environmentally sensitive respiratory viruses are seasonal, and more common in winter; and therefore, (3) climatic effects could be protective over space (hot, dry places might have less transmission) and time (summer might see reduced transmission compared to winter). All three are plausible and are generally consistent, but in many places (including, and especially, on social media), the basic premise of each has been communicated to the public, and policymakers, in a way that obscures key nuance and creates false confidence.
Experimental evidence shows that SARS-CoV-2 is environmentally sensitive
Like other viruses with a lipid envelope, SARS-CoV-2 is probably sensitive to temperature, humidity, and solar radiation; this affects its ability to persist on surfaces and in air, and might have subtle impacts on transmission. But the finer details of microbiology are often lost, leading to false confidence in how lab studies could scale up to the real world. For example, studies showing that germicidal ultraviolet radiation in hospitals and laboratories (ultraviolet C (UV-C) wavelengths) kills the virus have been misconstrued as evidence that sunlight (a mix of UV-A and UV-B) would effectively neutralize the virus in outdoor public spaces, possibly at a scale detectable from case data1. Newer experimental evidence supports the hypothesis that sunlight might also have an effect on SARS-CoV-22, although a global study found only a 1% reduction in transmission linked to environmental UV radiation3. In the real world, these effects will be slight, and unlikely to set hard limits on transmission anywhere in the world4.
Other upper respiratory tract infections are seasonal, with declines in warmer months
Influenza, the common cold, and other respiratory infections show seasonal transmission that coincides with changes in temperature, humidity, and solar radiation. But seasonal epidemics are also a product of the transmissibility of a virus, the initial susceptibility of a population, and the degree and nature of immunity conferred by infections. In basic epidemiological models, stable “oscillations” like seasonal epidemic waves usually require some degree of immunity5,6; at the start of a pandemic, when transmissibility is high and immunity is low, even strong environmental drivers are unlikely to curb transmission. Previous influenza pandemics show the importance of this nuance7. “Seasonal” versus “pandemic” influenza refers not just to different epidemic phases, but entirely different viral strains, and population susceptibility to pandemic strains starts high enough for rapid spread regardless of the season. During the first wave of the 2009 A/H1N1 pandemic, epidemic growth was still possible in August, the most environmentally unfavorable point in the year, with immunity under 20%; months later, immunity—and consequently environmental sensitivity—may eventually have been high enough to lead to a winter-driven third wave8. Scientists anticipate a similar pattern for COVID-19: while the virus could develop seasonal oscillations if it becomes endemic9 (i.e., if pandemic control fails in the long term), current susceptibility to SARS-CoV-2 is high enough that summer weather is unlikely to be protective10.
Environmental drivers could plausibly create seasonal or geographic differences in COVID-19 outbreak intensity
But those drivers’ impacts are heavily confounded by immunity, interventions, human behavior, and other details that are usually left out of models, leading to potentially spurious conclusions. For example, most available contact tracing data indicate that the proportion of indoor transmission is high11,12, a pattern likely caused by a combination of social contact patterns (including both number, intensity, and duration of contacts), air circulation, and potential weather drivers like sunlight or humidity. However, when studies attempt to model links between temperature and transmission, they almost always use gridded climate data or local weather data that represents the outdoors, and is unrepresentative of indoor conditions a virus particle or aerosolized cloud would actually experience in most transmission events.
Perhaps the biggest confounder, social behavior is environmentally driven and seasonal, but is rarely weighed alongside environmental and immunity drivers as a hypothesis for why infectious diseases show seasonality13,14. For example, school terms are seasonal and have a marked influence on social mixing patterns relevant to influenza transmission, even in pandemics15,16,17. Without individual-level transmission data, it can be difficult to distinguish direct biological impacts of weather from behaviorally mediated seasonality, and in some cases, the two are blurred (e.g., vitamin D levels are driven by both weather and seasonal behavior). Confusing the two could easily lead to spurious predictions. If behavioral patterns become unpredictable—either because of externalities like social distancing restrictions or a feedback loop between science and public risk perceptions around seasonality—attempts at forecasting the pandemic based on environmental seasonality will only become more unreliable. More

NEWS Q&A
27 August 2020

The spill released a new type of low-sulfur fuel, and its ecological effects aren’t well studied, says environment advocate Jaqueline Sauzier.

Dyani Lewis

Search for this author in:

Around 1,000 tonnes of oil have leaked from the MV Wakashio off Mauritius since 6 August.Credit: Pierre Dalais/EPA-EFE/Shutterstock

When the cargo ship MV Wakashio ran aground on a coral reef on the southeast tip of Mauritius, in the Indian Ocean almost exactly a month ago, it unleashed a vast oil spill. The Japanese-owned vessel held 200 tonnes of diesel and 3,900 tonnes of fuel oil, an estimated 1,000 tonnes of which leaked into the sea when the ship’s hull cracked on 6 August. It is the first reported spill of a new type of low-sulfur fuel that has been introduced to reduce air pollution. The spill has left a 15-kilometre stretch of the coastline — an internationally recognized biodiversity hotspot — smeared with oil.
Jacqueline Sauzier, president of the non-profit Mauritius Marine Conservation Society in Phoenix, has been helping with volunteer efforts to contain the spill. She spoke to Nature about how the clean-up is progressing.
What has been the response to the spill?
Mauritius is not geared up to deal with a catastrophe of this size, so other countries have sent experts to help. A French team arrived first, from the nearby island Réunion, to erect ocean booms — floating structures that contain the spill. The United Nations sent a team including experts in oil spills and crisis management. They’ve been working with communities, the private sector and the government to coordinate clean-up efforts. Marine ecologists and others have arrived from Japan and the United Kingdom.
Mauritians were also very proactive. In one weekend, we made nearly 80 kilometres of make-shift ocean booms out of cane trash — the leftover leaves and waste from sugar-cane processing — to contain the oil. Empty bottles were put in the middle of the booms to make them float, and anchors were attached to keep them from drifting away with the current.
For ten days, people worked night and day to contain as much oil as possible so that it wouldn’t reach the shoreline, where it is more difficult to clean. We managed to contain and remove nearly 75% of the spilled oil. Only a small amount reached the shore. But there’s still the issue of water-soluble chemicals that come from the oil, but dissolve into the water and therefore aren’t scooped out with the oil that sits on the water’s surface.
What ecosystems have been affected?
When you look at images in the media, it can feel like the whole of Mauritius is under oil. But the oil reached only 15 kilometres of the 350-kilometre shoreline, so it could have been much worse.
Unfortunately, there are a lot of environmentally sensitive areas in the region affected. The ship ran aground off Pointe d’Esny and just to the north of Blue Bay Marine Park. These sites are listed under the Ramsar Convention on Wetlands of International Importance as biodiversity hotspots. Ocean currents carried the oil northwards, so fortunately there’s none in the Blue Bay Marine Park, but the mangroves on the shoreline north of Pointe d’Esny have been covered. This will definitely have an impact, because mangroves are the nursery of the marine environment.
The Île aux Aigrettes, a small island near the wreck, has also been affected. The island is home to vulnerable pink pigeons (Nesoenas mayeri) and other native birds, and Telfair’s skink (Leiolopisma telfairii). The Mauritian Wildlife Foundation in Port Louis was already working to restore the island’s unique plants and remove invasive species. The oil didn’t go onto the island itself, but chemicals might have seeped into the corals and fumes from the spill could also have an impact.
Two rivers open into the bay where the oil spill is. The brackish water at the mouths of the rivers is an important ecosystem, and the oil has managed to go up parts of the rivers. The oil slick also floated above a large and rare area of seagrass, which is home to seahorses. Although the oil didn’t touch the seagrass, we fear that chemicals in the water could reach them.

Jacqueline Sauzier is president of the Mauritius Marine Conservation Society in Phoenix.Credit: Jacqueline Sauzier

Are there particular species affected?
It is not one species that could be at risk. It is the whole ecosystem, because of the dispersal of water-soluble chemicals in the water. Filter feeders, such as corals and crustaceans and molluscs, are probably the first to be impacted. We haven’t seen lots of animals dying, but we will need to monitor for signs.
Bad weather over the past two weeks has also forced the ship against the coral reef. That pushed a lot of sand and broken coral over the reef into the lagoon, creating a sand bar just inside the reef. That could change the currents in the lagoon and will have an impact on coral growth.
The social impacts are also a big concern for us. Fishing communities living in the region cannot fish anymore, because the fish that have been caught contain high levels of arsenic.
Something that is also concerning is that we don’t know the possible long-term effects. The oil is a new low-sulfur fuel oil that is being introduced to reduce air pollution. This is the first time that type of oil has spilled, so there have been no long-term studies on the impacts.
What steps are being taken now?
As soon as the ship grounded, people started monitoring the quality of the water. So we have this baseline from before the spill and we know the target that we have to reach for remediation of the water.
Oil-spill experts are formulating a plan to clean the shoreline properly. The impact on the mangroves could be worse if the cleaning is done badly. It could also push chemicals into the sand, which could be released in warm weather a year or two from now.
The front part of the vessel has been tugged away to be sunk along the shipping route. This was the least bad option. The rear is still on the reef. It has been cleaned of fuel, but rust and paint could still cause damage. It’s also falling apart, which can break the coral
Can future spills be avoided?
We were lucky this time that the spill was small and the boat grounded when it did. First, if it had happened in April, we wouldn’t have been able to go out, because we were locked down because of COVID. Second, the cane harvest started at the end of June. So if the wreck had happened earlier, we wouldn’t have had the cane trash readily available to contain the oil.
The spill has opened people’s eyes. Mauritius lies near a shipping highway. Around 2,500 large vessels pass close to Mauritius every month. It is difficult to describe how very big the Wakashio is. From the shore, it’s as if something foreign is sitting on your veranda. But this is the third stranding that we have had in ten years. Each time, there is uproar from the community, saying large vessels are much too close to the island.
It is not the first time, but I really hope it is the last.

doi: 10.1038/d41586-020-02446-7

This interview has been edited for length and clarity.

Latest on:

Biodiversity

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

Assembling and monitoring the synthetic bacterial communities
We selected three bacterial strains, a Chryseobacterium sp., a Staphylococcus sp., and a Bacillus sp., from a large collection of soil isolates, and we screened them with fluorescence-independent flow cytometry as we did previously [17]. We measured at an acquisition speed of 14 μl min−1 for 2 min per sample and setting a threshold of 10,000 regarding the height signal of the front scatter (FSC-H). Differently than in our previous work, here we acquired all scattering profiles based on growth assays at 30 °C that was the temperature at which we performed all the experiments. In addition, we screened the growing cultures with a temporal resolution of 20, rather than 30, min. We recorded significant interactions among the strains by comparing their single and mixed growth profiles at 30 °C (Dataset 1). These interactions were mainly positive, similarly to what we found previously at other temperatures [17]. We performed all the related growth assays in biological triplicates. All flow cytometry data are available in .fcs format online (http://flowrepository.org) under the “FR-FCM-Z25Q” identifier. Henceforth, when referring to the experiments we use the term “population” to describe the cells of a given strain within a flask at a given assay and the term “community” to describe the total bacterial cells within a flask at a given assay.
Quantification of “background noise”
Our flow cytometry method for screening the mixed bacterial cultures has an accuracy of 97% for sample densities above 105 cells ml−1 [17]. However, at lower densities sampling errors and instrument inconsistencies become increasingly important because the signal-to-noise ratio drops. This can result in substantially different counts among identical samples and can artificially inflate the observed variability. Thus, it was essential to quantify this “noise” before performing the main experiments and subtract it from the observed variability when quantifying drift.
To that end, we made a series of separate experiments to quantify “noise”. In these experiments, we mixed overnight cultures of the three strains in all the seven possible combinations and in final cell densities ranging from 1.6 × 104 to 6.3 × 107 strain−1 ml−1 (corresponding to the expected range of cell densities in the main experiment, see below), and we measured repeated aliquots from the same flask to determine the coefficient of variation (CV—Fig. 1a). We treated the samples in exactly the same way as in the main experiments to include the effect of sampling errors in our calculations. We acquired in total of 99 triplicate measurements of 1–3 populations for a total of 148 observations (Fig. 1a, Dataset 2). We hypothesized that the level of “noise” should be inversely related to the cell density of the sample, because the signal-to-noise ratio decreases at low cell densities in the flow cytometer. Accordingly, we fit different functions for the dependency of CV to cell density (Supplementary Table S1). Finally, we calculated the 99.5% confidence intervals of the best-fitting function (i.e., Michaelis-Menten) using the confint function of the MASS package [18] in R [19], and we defined the false discovery rate based on the number of observations that were above the upper 99.5% confidence interval (Supplementary Fig. S1). Finally, we verified that the levels of noise determined in this study are similar to the variability recorded from technical replicate samples taken during our previous experiment where we used the same bacterial system and instrument with identical settings [17] (Supplementary Fig. S2).
Main experiments
To quantify drift, we monitored the changes in population densities across identical starting communities incubated under the same environmental conditions (Fig. 1b). To that end, we mixed the three strains in all seven possible combinations, i.e., three monocultures, three mixed cultures of two strains and a mixed culture of all three strains, and in three different starting total cell densities (5 × 104 cells ml−1, 105 cells ml−1, and 106 cells ml−1). To perform each growth assay, we first inoculated the respective strains from overnight pure cultures in a single flask containing 300 ml of Luria–Bertani medium (Sigma). To reach the desired starting total cell density, we estimated the cell density of the overnight pure cultures with flow cytometry [17] and we inoculated the respective volume. We then mixed the culture thoroughly and we immediately split the volume equally into three flasks. We next sampled 500 μl from each flask and we compared the variability in the bacterial populations across the three flasks to the expected “background noise” for the same cell density. In specific, we examined the CVs of the bacterial populations and their z-scores compared to the “background noise”, i.e., how many standard deviations an observed CV differs from the expected “background noise” CV at a given cell density. If the observed z-scores were larger than 2 (95% CI), we aborted the given experiment because it indicated that we introduced variability when we mixed and split the cultures and thus the starting cultures could not be considered identical; this happened in ~50% of the cases. If the observed z-scores were lower than 2, indicating that the recorded variability was not statistically different or was less than the expected variability based on the “background noise”, we proceeded with the experiment, incubating the three flasks in the same chamber (New Brunswick Innova 42R, Eppendorf) at 30 °C and with shaking at 80 rpm.
We recorded the starting densities (Dataset 3, z-scores between −6.68 and 1.17) and the densities every 20 min until the end of the fourth hour of incubation starting from the 60th minute. To detect and quantify drift, our main assumption was that any larger-than-expected deviations in the population densities of identical starting communities incubated under the same environmental conditions could only be because of drift. Thus, we compared the observed CVs to the expected CVs based on the “background noise” by deploying the z-score. We quantified drift using two different thresholds:
1.
The “upper threshold” that focused on excluding false-positive observations. In this quantification, we used a cutoff significance level of z  > 3 (99.5% CI) and we ignored the lowest 17.57% of positive observations (i.e., 15 observations, corresponding to the FDR level of the “background noise”) to minimize the detection of false positives.

2.
The “mean threshold” that focused on excluding false negatives and increasing detectability. In this quantification, we used a cutoff significance level of z  > 0, meaning that we scored any observation greater than the mean noise function as positive.

The “upper threshold” quantification probably overestimates drift by taking into account only the highest among the recorded CV values while the “mean threshold” quantification underestimates drift by taking into account some low CV values that are very close to the noise levels. Thus, the “upper threshold” and “mean threshold” quantifications do not represent the true levels of drift (which are hard to define whatsoever in the presence of noise) but they rather represent the upper and lower boundaries within which the true levels of drift lie.
Estimating potential growth variability due to temperature differences within the incubation chamber
To ensure that the recorded variability in the population counts was not due to slight differences in temperature within the incubation chamber, we estimated the potential variability that could have resulted if each strain grew within the extremes of the recorded temperatures in the chamber. For that, we first measured the temperature within each flask at each experimental time point, five times per flask, using a digital immersion thermometer with an accuracy of 0.1 °C. The temperature varied by 0.15 °C ± 0.08 °C on average and by 0.28 °C at maximum. We then calculated the growth rates of each of the strains under the recorded temperature extremes at each time point by interpolating from previously recorded growth rates [17]. We interpolated both with respect to time and with respect to temperature because the previous data were recorded at intervals of 30 min and 0.5 °C (the latter at a range of 25–42 °C). Finally, for each strain, we calculated how the CV in the hypothetical population densities would increase if the strains were constantly growing within the recorded temperature extremes for the duration of the experiment and if the CV was calculated from three observations (like in the real experiments) of the resulting population density distributions. We note that with this analysis, we probably overestimated the hypothetical increase in population densities because we used growth rates from mixed cultures that were generally higher than those in monocultures because of the positive interactions among the strains (Dataset 1).
Simulations
To simulate drift in complex bacterial communities, we used in silico communities with diversity and abundance distributions similar to nature [20] where drift acts with magnitude according to our experimental data. A conceptual flowchart of the simulations can be found in Supplementary Information (Supplementary Fig. S3). Each simulation involved a metacommunity of 100 communities that were connected with dispersal and that initially contained 2000 species each.
We simulated dispersal occurring in a unidirectional way within a closed system; individuals from community n disperse to the community (n + 1) and individuals from community 100 disperse back to community 1. The strength of dispersal equaled to the percentage of individuals that disperse to the respective community and it varied between 2 and 20%. Our aim in modeling dispersal in this way was to create a setting where habitat fragmentation was high and therefore drift’s importance is expected to be more pronounced [21, 22], and where there was no gain or loss of individuals from outside the metacommunity.
We simulated selection as differences in the growth rates among species within a community. The growth rates were distributed normally with a mean of 1 (resembling systems at their carrying capacity) and with a standard deviation between 0.071 and 0.167. Growth rates changed at every generation by being re-drawn from the same distribution in an effort to represent fluctuating habitats where a given species is not always favored or disfavored. Therefore, in our simulations, the standard deviation in the growth rates represents the strength of selection, because the higher it is the bigger are the differences in the growth rates in a community and the changes in the growth of a species from generation to generation. The distribution of the abundances in a community at time zero was log-normal (mean = 4, sd = 1.1) and the distribution of the abundances of a given species across all communities was normal with a standard deviation equal to the strength of selection.
The metacommunity grew for 1000 generations under given dispersal and selection conditions with drift, where drift changed the assigned growth rates at every generation according to a distribution based on the defined threshold from the experimental data (“upper” or “mean” threshold). In parallel, an initially identical metacommunity grew under the same dispersal and selection conditions but without drift, meaning that the assigned growth rates at every generation did not change further. More details and an example on how we modeled changes in growth rates due to drift are presented in the Supplementary Text in Supplementary Information.
For a given generation, we calculated the effect of drift by comparing a given community in the metacommunity that grew under drift to the same community in the metacommunity that grew without drift. In specific, we examined the β-diversity by means of the Bray–Curtis (BC) community similarity and the differences in species richness and in Pielou’s evenness among drift-impacted and drift-free communities, calculating the metacommunity-wise mean and standard deviation on all these properties. Moreover, we kept track of the extinct species at the end of each simulation and we mapped their initial relative abundances, but here we report the metacommunity-wise median because the distribution of the relative abundances is skewed (Supplementary Fig. S4). We ran simulations under 50 different scenarios resulting from five levels of selection strength over ten levels of dispersal rate. To estimate the effect of drift on Bray–Curtis similarity in metacommunities with increasing number of species, we ran the same simulation at the highest selection and lowest dispersal levels, at intermediate selection and dispersal and at the lowest selection and highest dispersal, but we changed the number of species; we ran the simulation three times in metacommunities of 500, 1000, 2000, 4000, 6000, 8000, and 10,000 species. We performed all simulations in R. All code is available on GitHub (https://github.com/sfodel/Drift).
Reported β-diversity in stochastically assembled communities in nature
To compare our simulation results with the results from environmental surveys regarding the β-diversity in stochastically assembled communities, we searched for related studies using the following two criteria: (1) the study cites the works of Stegen and colleagues [12, 13], where the term “undominated” community assembly is presented formally for microbial ecology, (2) the study reports data on the range of the observed β-diversity in terms of Bray–Curtis dissimilarity (or similarity) in stochastically assembled communities, or this range can be inferred from the data presented in that study. More

Amos W (2009) Heterozygosity and mutation rate: evidence for an interaction and its implications. Bioessays 32:82–90
Google Scholar 

Amos W, Kosanović D, Eriksson A (2015) Inter-allelic interactions play a major role in microsatellite evolution. Proc R Soc B 282:20152125
PubMed  Google Scholar 

Amos W (2016) Heterozygosity increases microsatellite mutation rate. Biol Lett 12:20150929
PubMed  PubMed Central  Google Scholar 

Azevedo L, Serrano C, Amorim A, Cooper DN (2015) Trans-species polymorphism in humans and the great apes is generally maintained by balancing selection that modulates the host immune response. Hum Genomics 9:21
PubMed  PubMed Central  Google Scholar 

Beye M, Hasselmann M, Fondrk MK, Page RE, Omholt SW (2003) The gene csd is the primary signal for sexual development in the honeybee and encodes an SR-type protein. Cell 114:419–429
CAS  PubMed  Google Scholar 

Beye M, Seelmann C, Gempe T, Hasselmann M, Vekemans X, Fondrk MK et al. (2013) Gradual molecular evolution of a sex determination switch through incomplete penetrance of femaleness. Curr Biol 23:2559–2564
CAS  PubMed  Google Scholar 

Bezabih G, Adgaba N, Hepburn HR, Pirk CWW (2014) The territorial invasion of Apis florea in Africa. Afr Entomol 22:888–890
Google Scholar 

Biewer M, Lechner S, Hasselmann M (2016) Similar but not the same: insights into the evolutionary history of paralogous sex-determining genes of the dwarf honey bee Apis florea. Heredity 116:12–22
CAS  PubMed  Google Scholar 

Chapuis MP, Plantamp C, Streiff R, Blondin L, Piou C (2015) Microsatellite evolutionary rate and pattern in Schistocerca gregaria inferred from direct observation of germline mutations. Mol Ecol 24:6107–6119
CAS  PubMed  Google Scholar 

Chookajorn T, Kachroo A, Ripoll DR, Clark AG, Nasrallah JB (2004) Specificity determinants and diversification of the Brassica self-incompatibility pollen ligand. Proc Natl Acad Sci USA 101:911–917
CAS  PubMed  Google Scholar 

Cho S, Huang ZY, Green DR, Smith DR, Zhang J (2006) Evolution of the complementary sex-determination gene of honey bees: balancing selection and trans-species polymorphisms. Genome Res 16:1366–1375
CAS  PubMed  PubMed Central  Google Scholar 

Clarke BC (1979) The evolution of genetic diversity. Proc R Soc B 205:453–474
CAS  Google Scholar 

Ding G, Xu H, Oldroyd BP, Gloag RS (2017) Extreme polyandry aids the establishment of invasive populations of a social insect. Heredity 119:381–387
CAS  PubMed  PubMed Central  Google Scholar 

Ellegren H (2004) Microsatellites: simple sequences with complex evolution. Nat Rev Genet 5:435–445
CAS  PubMed  Google Scholar 

Fay JC, Wu CI (2000) Hitchhiking under positive Darwinian selection. Genetics 155:1405–1413
CAS  PubMed  PubMed Central  Google Scholar 

Fu YX (1997) Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147:915–925
CAS  PubMed  PubMed Central  Google Scholar 

Glémin S, Gaude T, Guillemin ML, Lourmas M, Olivieri I, Mignot A (2005) Balancing selection in the wild: testing population genetics theory of self-incompatibility in the rare species Brassica insularis. Genetics 171:279–289
PubMed  PubMed Central  Google Scholar 

Gervais CE, Castric V, Ressayre A, Billiard S (2011) Origin and diversification dynamics of self-incompatibility haplotypes. Genetics 188:625–636
CAS  PubMed  PubMed Central  Google Scholar 

Gibbs MJ, Armstrong JS, Gibbs AJ (2000) Sister-scanning: a Monte Carlo procedure for assessing signals in recombinant sequences. Bioinformatics 16:573–582
CAS  PubMed  Google Scholar 

Gloag R, Christie JR, Ding G, Stephens R, Buchmann G, Oldroyd BP (2019) Workers’ sons rescue genetic diversity at the sex locus in an invasive honey bee population. Mol Ecol 28:1585–1592
PubMed  Google Scholar 

Gloag R, Ding G, Christie JR, Buchmann G, Beekman M, Oldroyd BP (2017) An invasive social insect overcomes genetic load at the sex locus. Nat Ecol Evol 1:1–6
Google Scholar 

Hasselmann M, Beye M (2004) Signatures of selection among sex-determining alleles of the honey bee. Proc Natl Acad Sci USA 101:4888–4893
CAS  PubMed  Google Scholar 

Hasselmann M, Vekemans X, Pflugfelder J, Koeniger N, Koeniger G, Tingek S et al. (2008) Evidence for convergent nucleotide evolution and high allelic turnover rates at the complementary sex determiner gene of Western and Asian honeybees. Mol Biol Evol 25:696–708
CAS  PubMed  Google Scholar 

Hedrick PW (2005) A standardized genetic differentiation measure. Evolution 59:1633–1638
CAS  PubMed  Google Scholar 

Koch V, Nissen I, Schmitt BD, Beye M (2014) Independent evolutionary origin of fem paralogous genes and complementary sex determination in hymenopteran insects. PLoS ONE 9:e91883
PubMed  PubMed Central  Google Scholar 

Koetz AH (2013) Ecology, behaviour and control of Apis cerana with a focus on relevance to the Australian incursion. Insects 4:558–592
PubMed  PubMed Central  Google Scholar 

Kosakovsky Pond SL, Posada D, Gravenor MB, Woelk CH, Frost SD (2006) Automated phylogenetic detection of recombination using a genetic algorithm. Mol Biol Evol 23:1891–1901
PubMed  Google Scholar 

Kumar S, Stecher G, Li M, Knyaz C, Tamura K (2018) MEGA X: molecular evolutionary genetics analysis across computing platforms. Mol Biol Evol 35:1547–1549
CAS  PubMed  PubMed Central  Google Scholar 

Kurbalija Novičić Z, Sayadi A, Jelić M, Arnqvist G (2020) Negative frequency dependent selection contributes to the maintenance of a global polymorphism in mitochondrial DNA. BMC Evol Biol 20:20. https://doi.org/10.1186/s12862-020-1581-2
CAS  Article  PubMed  PubMed Central  Google Scholar 

Lechner S, Ferretti L, Schöning C, Kinuthia W, Willemsen D, Hasselmann M (2014) Nucleotide variability at its limit? Insights into the number and evolutionary dynamics of the sex-determining specificities of the honey bee Apis mellifera. Mol Biol Evol 31:272–287
CAS  PubMed  Google Scholar 

Levin BR (1988) Frequency dependent selection in bacterial populations. Philos Trans R Soc Lond B 319:459–472
CAS  Google Scholar 

Li Y-C, Korol AB, Fahima T, Nevo E (2004) Microsatellites within genes: structure, function, and evolution. Mol Biol Evol 21:991–1007
CAS  PubMed  Google Scholar 

Lynch M (2015) Feedforward loop for diversity. Nature 523:414–416
CAS  PubMed  PubMed Central  Google Scholar 

Mackensen O (1951) Viability and sex determination in the honeybee (Apis mellifera L.). Genetics 36:500–509
CAS  PubMed  PubMed Central  Google Scholar 

Martin DP, Murrell B, Golden M, Khoosal A, Muhire B (2015) RDP4: Detection and analysis of recombination patterns in virus genomes. Virus Evol 1:vev003
PubMed  PubMed Central  Google Scholar 

Martin DP, Posada D, Crandall KA, Williamson C (2005) A modified bootscan algorithm for automated identification of recombinant sequences and recombination breakpoints. AIDS Res Hum Retroviruses 21:98–102
CAS  PubMed  Google Scholar 

May G, Matzke E (1995) Recombination and variation at the A mating-type of Coprinus cinereus. Mol Biol Evol 12:794–802
CAS  Google Scholar 

Moritz RFA, Haddad N, Bataieneh A, Shalmon B, Hefetz A (2010) Invasion of the dwarf honeybee Apis florea into the near East. Biol Invasions 12:1093–1099
Google Scholar 

Moritz RFA, Härtel S, Neumann P (2005) Global invasions of the western honeybee (Apis mellifera) and the consequences for biodiversity. Écoscience 12:289–301
Google Scholar 

Muirhead CA (2001) Consequences of population structure on genes under balancing selection. Evolution 55:1532–1541
CAS  PubMed  Google Scholar 

Muirhead CA, Glass NL, Slatkin M (2002) Multilocus self-recognition systems in fungi as a cause of trans-species polymorphism. Genetics 161:633–641
CAS  PubMed  PubMed Central  Google Scholar 

Nasrallah JB (1997) Evolution of the Brassica self-incompatibility locus: a look into S-locus gene polymorphisms. Proc Natl Acad Sci USA 94:9516–9519
CAS  PubMed  Google Scholar 

Nasrallah JB, Kao TH, Chen CH, Goldberg ML, Nasrallah ME (1987) Amino-acid sequence of glycoproteins encoded by three alleles of the S locus of Brassica oleracea. Nature 326:617–619
CAS  Google Scholar 

Oldroyd BP, Wongsiri S (2006) Asian honey bees. Biology, conservation and human interactions. Harvard University Press, Cambridge
Google Scholar 

Posada D, Crandall KA (2001) Evaluation of methods for detecting recombination from DNA sequences: computer simulations. Proc Natl Acad Sci USA 98:13757–13762
CAS  PubMed  Google Scholar 

Radloff SE, Hepburn C, Hepburn HR, Fuchs S, Hadisoesilo S, Tan K et al. (2010) Population structure and classification of Apis cerana. Apidologie 41:589–601
Google Scholar 

Raper JR (1966) Sexuality in fungi. (Book reviews: genetics of sexuality in higher fungi). Science 154:758
Google Scholar 

Richman AD, Kohn JR (1999) Self-incompatibility alleles from Physalis: implications for historical inference from balanced genetic polymorphisms. Proc Natl Acad Sci USA 96:168–172
CAS  PubMed  Google Scholar 

Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE et al. (2017) DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol 34:3299–3302
CAS  PubMed  Google Scholar 

Schierup MH, Bechsgaard JS, Christiansen FB (2008) Selection at work in self-incompatible Arabidopsis lyrata. II. Spatial distribution of S haplotypes in Iceland. Genetics 180:1051–1059
PubMed  PubMed Central  Google Scholar 

Schierup MH, Vekemans X, Charlesworth D (2000) The effect of subdivision on variation at multi-allelic loci under balancing selection. Genet Res 76:51–62
CAS  PubMed  Google Scholar 

Smith DR (2011) Asian Honeybees and Mitochondrial DNA. In: Hepburn HR, Radloff SE (eds) Honeybees of Asia. Springer-Verlag Berlin, Heidelberg, p 69–93
Google Scholar 

Solignac M, Vautrin D, Loiseau A, Mougel F, Baudry E, Estoup A et al. (2003) Five hundred and fifty microsatellite markers for the study of the honeybee (Apis mellifera L.) genome. Mol Ecol Notes 3:307–311
CAS  Google Scholar 

Tajima F (1989) The effect of change in population size on DNA polymorphism. Genetics 123:597–601
CAS  PubMed  PubMed Central  Google Scholar 

Takahashi J, Shimizu S, Koyama S, Kimura K, Shimizu I, Yoshida T (2009) Variable microsatellite loci isolated from the Asian honeybee, Apis cerana (Hymenoptera; Apidae). Mol Ecol Resour 9:819–821
CAS  PubMed  Google Scholar 

Takahata N, Nei M (1990) Allelic genealogy under overdominant and frequency-dependent selection and polymorphism of major histocompatibility complex loci. Genetics 124:967–978
CAS  PubMed  PubMed Central  Google Scholar 

Wallberg A, Han F, Wellhagen G, Dahle B, Kawata M, Haddad N et al. (2014) A worldwide survey of genome sequence variation provides insight into the evolutionary history of the honeybee Apis mellifera. Nat Genet 46:1081–1088
CAS  PubMed  Google Scholar 

Walsh PS, Metzger DA, Higuchi R (1991) Chelex 100 as amedium for simple extraction of DNA for PCR-based typing from forensic material. Biotechniques 10:506–513
CAS  PubMed  Google Scholar 

Wang J (2015) Does GST underestimate genetic differentiation from marker data? Mol Ecol 24:3546–3558
CAS  PubMed  Google Scholar 

Weaver S, Shank SD, Spielman SJ, Li M, Muse SV, Kosakovsky Pond SL (2018) Datamonkey 2.0: a modern web application for characterizing selective and other evolutionary processes. Mol Biol Evol 35:773–777
CAS  PubMed  PubMed Central  Google Scholar 

Woyke J (1963) What happens to diploid drone larvae in honeybee colony? J Apic Res 2:73–75
Google Scholar 

Wright S (1939) The distribution of self-sterility alleles in populations. Genetics 24:538–552
CAS  PubMed  PubMed Central  Google Scholar 

Wu Q, Han TS, Chen X, Chen JF, Zou YP, Li ZW et al. (2017) Long-term balancing selection contributes to adaptation in Arabidopsis and its relatives. Genome Biol 18:217
PubMed  PubMed Central  Google Scholar 

Yang S, Wang L, Huang J, Zang X, Yuan Y, Chen JQ et al. (2015) Parent-progeny sequencing indicates higher mutation rates in heterozygotes. Nature 523:463–467
CAS  PubMed  Google Scholar 

Yokoyama S, Nei M (1979) Population dynamics of sex-determining alleles in honey bees and self-incompatibility alleles in plants. Genetics 91:609–626
CAS  PubMed  PubMed Central  Google Scholar 

Zareba J, Blazej P, Laszkiewicz A, Sniezewski L, Majkowski M, Janik S et al. (2017) Uneven distribution of complementary sex determiner (csd) alleles in Apis mellifera population. Sci Rep. 7:2317
PubMed  PubMed Central  Google Scholar  More