Ecology

1.
SAUP. Sea Around Us. http://www.seaaroundus.org/data/ (2020).
2.
Coll, M., Navarro, J., Olson, R. J. & Christensen, V. Assessing the trophic position and ecological role of squids in marine ecosystems by means of food-web models. Deep Sea Res. Part II Top. Stud. Oceanogr. 95, 21–36 (2013).
ADS  Article  Google Scholar 

3.
Hastie, L. et al. Cephalopods in the north-eastern Atlantic: Species, biogeography, ecology, exploitation and conservation. In Oceanography and Marine Biology (eds. Gibson, R., Atkinson, R. & Gordon, J.) vol. 20092725, 111–190 (CRC Press, Boca Raton, 2009).

4.
Piatkowski, U. & Pierce, G. J. Impact of cephalopods in the food chain and their interaction with the environment and fisheries: An overview. Fish. Res. 6, 5–10 (2001).
Article  Google Scholar 

5.
Pierce, G. J. et al. A review of cephalopod–environment interactions in European Seas. Hydrobiologia 612, 49–70 (2008).
Article  Google Scholar 

6.
André, J., Haddon, M. & Pecl, G. T. Modelling climate-change-induced nonlinear thresholds in cephalopod population dynamics. Glob. Change Biol. 16, 2866–2875 (2010).
ADS  Article  Google Scholar 

7.
Jereb, P. et al. Cephalopod biology and fisheries in Europe: II. Species Accounts. ICES Cooper. Res. Rep. 325, 1–360 (2015).
Google Scholar 

8.
Sims, D. W., Genner, M. J., Southward, A. J. & Hawkins, S. J. Timing of squid migration reflects North Atlantic climate variability. Proc. R. Soc. Lond. Ser. B Biol. Sci. 268, 2607–2611 (2001).

9.
Dorey, N. et al. Ocean acidification and temperature rise: Effects on calcification during early development of the cuttlefish Sepia officinalis. Mar. Biol. 160, 2007–2022 (2013).
CAS  Article  Google Scholar 

10.
Rodhouse, P. G. K. et al. Environmental effects on cephalopod population dynamics. In Advances in Marine Biology vol. 67, 99–233 (Elsevier, Amsterdam, 2014).

11.
Millar, R. J. et al. Emission budgets and pathways consistent with limiting warming to 1.5 °C. Nat. Geosci. 10, 741–747 (2017).
ADS  CAS  Article  Google Scholar 

12.
Otto, F. E. L., Frame, D. J., Otto, A. & Allen, M. R. Embracing uncertainty in climate change policy. Nat. Clim. Change 5, 917–920 (2015).
ADS  Article  Google Scholar 

13.
Meinshausen, M. et al. The RCP greenhouse gas concentrations and their extensions from 1765 to 2300. Clim. Change 109, 213–241 (2011).
ADS  CAS  Article  Google Scholar 

14.
van Vuuren, D. P. et al. The representative concentration pathways: An overview. Clim. Change 109, 5–31 (2011).
ADS  Article  Google Scholar 

15.
Gissi, E. et al. A review of the combined effects of climate change and other local human stressors on the marine environment. Sci. Total Environ. 755, 142564 (2021).
ADS  CAS  PubMed  Article  Google Scholar 

16.
Beaugrand, G. et al. Prediction of unprecedented biological shifts in the global ocean. Nat. Clim. Change 9, 237–243 (2019).
ADS  Article  Google Scholar 

17.
Jorda, G. et al. Ocean warming compresses the three-dimensional habitat of marine life. Nat. Ecol. Evol. 4, 109–114 (2020).
PubMed  Article  Google Scholar 

18.
Vidal, E. A. G., DiMarco, F. P., Wormuth, J. H. & Lee, P. G. Influence of temperature and food availability on survival, growth and yolk utilization in hatchling squid. Bull. Mar. Sci. 71, 915–931 (2002).
Google Scholar 

19.
Doubleday, Z. A. et al. Global proliferation of cephalopods. Curr. Biol. 26, R406–R407 (2016).
CAS  PubMed  Article  Google Scholar 

20.
van der Kooij, J., Engelhard, G. H. & Righton, D. A. Climate change and squid range expansion in the North Sea. J. Biogeogr. 43, 2285–2298 (2016).
Article  Google Scholar 

21.
Jin, Y., Jin, X., Gorfine, H., Wu, Q. & Shan, X. Modeling the oceanographic impacts on the spatial distribution of common cephalopods during autumn in the yellow sea. Front. Mar. Sci. 7, (2020).

22.
Pang, Y. et al. Variability of coastal cephalopods in overexploited China Seas under climate change with implications on fisheries management. Fish. Res. 208, 22–33 (2018).
Article  Google Scholar 

23.
Le Marchand, M. et al. Climate change in the Bay of Biscay: Changes in spatial biodiversity patterns could be driven by the arrivals of southern species. Mar. Ecol. Prog. Ser. 647, 17–31 (2020).
ADS  Article  Google Scholar 

24.
Lima, F. D., Ángeles-González, L. E., Leite, T. S. & Lima, S. M. Q. Global climate changes over time shape the environmental niche distribution of Octopus insularis in the Atlantic Ocean. Mar. Ecol. Prog. Ser. 652, 111–121 (2020).
ADS  Article  Google Scholar 

25.
Xavier, J. C., Peck, L. S., Fretwell, P. & Turner, J. Climate change and polar range expansions: Could cuttlefish cross the Arctic?. Mar. Biol. 163, 78 (2016).
Article  Google Scholar 

26.
Selig, E. R. et al. Mapping global human dependence on marine ecosystems. Conserv. Lett. 12, e12617 (2019).
Article  Google Scholar 

27.
Blasiak, R. et al. Climate change and marine fisheries: Least developed countries top global index of vulnerability. PLoS ONE 12, e0179632 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

28.
FAO. The State of Mediterranean and Black Sea Fisheries. (General Fisheries Commission for the Mediterranean, 2016).

29.
Lam, V. W. Y., Cheung, W. W. L., Reygondeau, G. & Sumaila, U. R. Projected change in global fisheries revenues under climate change. Sci. Rep. 6, 32607 (2016).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

30.
Badjeck, M.-C., Perry, A., Renn, S., Brown, D. & Poulain, F. The vulnerability of fishing-dependent economies to disasters. FAO Fish. Aquac. Circ. 1081, 1–19 (2013).
Google Scholar 

31.
Allison, E. H. et al. Vulnerability of national economies to the impacts of climate change on fisheries. Fish Fish. 10, 173–196 (2009).
Article  Google Scholar 

32.
Adloff, F. et al. Mediterranean Sea response to climate change in an ensemble of twenty first century scenarios. Clim. Dyn. 45, 2775–2802 (2015).
Article  Google Scholar 

33.
Alexander, M. A. et al. Projected sea surface temperatures over the 21st century: Changes in the mean, variability and extremes for large marine ecosystem regions of Northern Oceans. Elementa Sci. Anthropocene 6, 9 (2018).
Article  Google Scholar 

34.
Gaines, S. D. et al. Improved fisheries management could offset many negative effects of climate change. Sci. Adv. 4, eaao1378 (2018).
ADS  PubMed  PubMed Central  Article  Google Scholar 

35.
Pierce, G. J. et al. Status and trends of European cephalopod stocks. In ASC 2019 ICES Conference, Gothenburg, Sweden 1 (2019).

36.
Hutchinson, G. E. Concluding remarks. Cold Spring Harb. Symp. Quant. Biol. 22, 415–427 (1957).
Article  Google Scholar 

37.
Hutchinson, G. E. An Introduction to Population Ecology (Yale University Press, New Haven, 1978).
Google Scholar 

38.
Peterson, A. & Soberón, J. Species distribution modeling and ecological niche modeling: Getting the concepts right. Natureza e Conservação 10, 1–6 (2012).
Article  Google Scholar 

39.
Colwell, R. K. & Rangel, T. F. Hutchinson’s duality: The once and future niche. Proc. Natl. Acad. Sci. 106, 19651–19658 (2009).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

40.
Buisson, L., Thuiller, W., Casajus, N., Lek, S. & Grenouillet, G. Uncertainty in ensemble forecasting of species distribution. Glob. Change Biol. 16, 1145–1157 (2010).
ADS  Article  Google Scholar 

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

42.
Hao, T., Elith, J., Guillera-Arroita, G. & Lahoz-Monfort, J. J. A review of evidence about use and performance of species distribution modelling ensembles like BIOMOD. Divers. Distrib. https://doi.org/10.1111/ddi.12892 (2019).
Article  Google Scholar 

43.
Goberville, E., Beaugrand, G., Hautekèete, N.-C., Piquot, Y. & Luczak, C. Uncertainties in the projection of species distributions related to general circulation models. Ecol. Evol. 5, 1100–1116 (2015).
PubMed  PubMed Central  Article  Google Scholar 

44.
Leroy, B. et al. Forecasted climate and land use changes, and protected areas: The contrasting case of spiders. Divers. Distrib. 20, 686–697 (2014).
Article  Google Scholar 

45.
Schickele, A. et al. Modelling European small pelagic fish distribution: Methodological insights. Ecol. Model. 416, 108902 (2020).
Article  Google Scholar 

46.
Araújo, M. B. et al. Standards for distribution models in biodiversity assessments. Sci. Adv. 5, eaat4858 (2019).
ADS  PubMed  PubMed Central  Article  Google Scholar 

47.
Barbet-Massin, M., Thuiller, W. & Jiguet, F. How much do we overestimate future local extinction rates when restricting the range of occurrence data in climate suitability models?. Ecography 33, 878–886 (2010).
Article  Google Scholar 

48.
Beaugrand, G., Luczak, C., Goberville, E. & Kirby, R. Marine biodiversity and the chessboard of life. PLoS ONE 13, e0194006 (2018).
PubMed  PubMed Central  Article  CAS  Google Scholar 

49.
Støa, B., Halvorsen, R., Mazzoni, S. & Gusarov, V. I. Sampling bias in presence-only data used for species distribution modelling: Theory and methods for detecting sample bias and its effects on models. Sommerfeltia 38, 1–53 (2018).
Article  Google Scholar 

50.
Dufresne, J.-L. et al. Climate change projections using the IPSL-CM5 Earth System Model: from CMIP3 to CMIP5. Clim. Dyn. 40, 2123–2165 (2013).
Article  Google Scholar 

51.
Voldoire, A. et al. The CNRM-CM5.1 global climate model: Description and basic evaluation. Clim. Dyn. 40, 2091–2121 (2013).
Article  Google Scholar 

52.
Hourdin, F. et al. Impact of the LMDZ atmospheric grid configuration on the climate and sensitivity of the IPSL-CM5A coupled model. Clim. Dyn. 40, 2167–2192 (2013).
Article  Google Scholar 

53.
Friedlingstein, P. et al. Uncertainties in CMIP5 climate projections due to carbon cycle feedbacks. J. Clim. 27, 511–526 (2013).
ADS  Article  Google Scholar 

54.
Wiens, J. A., Stralberg, D., Jongsomjit, D., Howell, C. A. & Snyder, M. A. Niches, models, and climate change: Assessing the assumptions and uncertainties. PNAS 106, 19729–19736 (2009).
ADS  CAS  PubMed  Article  Google Scholar 

55.
Martinez-Meyer, E. Climate change and biodiversity: Some considerations in forecasting shifts in species’ potential distributions. Biodivers. Inform. 2, 42–55 (2005).
Article  Google Scholar 

56.
Levitus, S. Climatological atlas of the world ocean. Eos Trans. Am. Geophys. Union 64, 962–963 (2011).
ADS  Article  Google Scholar 

57.
Cabanes, C. et al. The CORA dataset: Validation and diagnostics of in-situ ocean temperature and salinity measurements. Ocean Sci. 9, 1–18 (2013).
ADS  Article  Google Scholar 

58.
Giorgetta, M. A. et al. Climate and carbon cycle changes from 1850 to 2100 in MPI-ESM simulations for the Coupled Model Intercomparison Project phase 5: Climate Changes in MPI-ESM. J. Adv. Model. Earth Syst. 5, 572–597 (2013).
ADS  Article  Google Scholar 

59.
Stevens, B. et al. Atmospheric component of the MPI-M Earth System Model: ECHAM6. J. Adv. Model. Earth Syst. 5, 146–172 (2013).
ADS  Article  Google Scholar 

60.
Jones, C. D. et al. The HadGEM2-ES implementation of CMIP5 centennial simulations. Geosci. Model Dev. Discuss. 4, 689–763 (2011).
ADS  Google Scholar 

61.
Schmidt, G. A. et al. Configuration and assessment of the GISS ModelE2 contributions to the CMIP5 archive: GISS MODEL-E2 CMIP5 SIMULATIONS. J. Adv. Model. Earth Syst. 6, 141–184 (2014).
ADS  Article  Google Scholar 

62.
Beaugrand, G., Lenoir, S., Ibañez, F. & Manté, C. A new model to assess the probability of occurrence of a species, based on presence-only data. Mar. Ecol. Prog. Ser. 424, 175–190 (2011).
ADS  Article  Google Scholar 

63.
Raybaud, V., Bacha, M., Amara, R. & Beaugrand, G. Forecasting climate-driven changes in the geographical range of the European anchovy (Engraulis encrasicolus). ICES J. Mar. Sci. 74, 1288–1299 (2017).
Article  Google Scholar 

64.
Thuiller, W., Lafourcade, B., Engler, R. & Araújo, M. B. BIOMOD—A platform for ensemble forecasting of species distributions. Ecography 32, 369–373 (2009).
Article  Google Scholar 

65.
Thuiller, W., Georges, D., Engler, R. & Breiner, F. Ensemble Platform for Species Distribution Modelling. (2016).

66.
Dormann, C. F. et al. Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. Ecography 36, 27–46 (2013).
Article  Google Scholar 

67.
Lenoir, J. et al. Species better track climate warming in the oceans than on land. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-020-1198-2 (2020).
Article  PubMed  Google Scholar 

68.
Smith, W. H. F. & Sandwell, D. T. Global sea floor topography from satellite altimetry and ship depth soundings. Science 277, 1956–1962 (1997).
CAS  Article  Google Scholar 

69.
NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group. Distance to the nearest coast. https://oceancolor.gsfc.nasa.gov/docs/distfromcoast/ (01/03/2018) (2009).

70.
Hattab, T. et al. Towards a better understanding of potential impacts of climate change on marine species distribution: A multiscale modelling approach. Glob. Ecol. Biogeogr. 23, 1417–1429 (2014).
Article  Google Scholar 

71.
Varela, S., Anderson, R. P., García-Valdés, R. & Fernández-González, F. Environmental filters reduce the effects of sampling bias and improve predictions of ecological niche models. Ecography 37, 1084–1091 (2014).
Google Scholar 

72.
Ben Rais Lasram, F. et al. An open-source framework to model present and future marine species distributions at local scale. Ecol. Inform. 59, 101130 (2020).
Article  Google Scholar 

73.
Montgomery, D. C. Design and Analysis of Experiments (Wiley, Hoboken, 2005).
Google Scholar 

74.
Getz, W. M. & Wilmers, C. C. A local nearest-neighbor convex-hull construction of home ranges and utilization distributions. Ecography 27, 489–505 (2006).
Article  Google Scholar 

75.
Cornwell, W. K., Schwilk, D. W. & Ackerly, D. D. A trait-based test for habitat filtering: Convex hull volume. Ecology 87, 1465–1471 (2004).
Article  Google Scholar 

76.
Hirzel, A. H., Le 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).
Article  Google Scholar 

77.
Leroy, B. et al. Without quality presence–absence data, discrimination metrics such as TSS can be misleading measures of model performance. J. Biogeogr. 45, 1994–2002 (2018).
Article  Google Scholar 

78.
Faillettaz, R., Beaugrand, G., Goberville, E. & Kirby, R. R. Atlantic Multidecadal Oscillations drive the basin-scale distribution of Atlantic bluefin tuna. Sci. Adv. 5, eaar6993 (2019).
ADS  PubMed  PubMed Central  Article  Google Scholar 

79.
Elith, J., Ferrier, S., Huettmann, F. & Leathwick, J. The evaluation strip: A new and robust method for plotting predicted responses from species distribution models. Ecol. Model. 186, 280–289 (2005).
Article  Google Scholar 

80.
VanDerWal, J. et al. Focus on poleward shifts in species’ distribution underestimates the fingerprint of climate change. Nat. Clim. Change 3, 239–243 (2013).
ADS  Article  Google Scholar 

81.
Taylor, K. E. Summarizing multiple aspects of model performance in a single diagram. J. Geophys. Res. Atmos. 106, 7183–7192 (2001).
ADS  Article  Google Scholar 

82.
Cristofari, R. et al. Climate-driven range shifts of the king penguin in a fragmented ecosystem. Nat. Clim. Change 8, 245–251 (2018).
ADS  Article  Google Scholar 

83.
Péron, C., Weimerskirch, H. & Bost, C.-A. Projected poleward shift of king penguins’ (Aptenodytes patagonicus) foraging range at the Crozet Islands, southern Indian Ocean. Proc. Biol. Sci. 279, 2515–2523 (2012).
PubMed  PubMed Central  Article  Google Scholar 

84.
Bloor, I. S. M., Attrill, M. J. & Jackson, E. L. Chapter One—A Review of the Factors Influencing Spawning, Early Life Stage Survival and Recruitment Variability in the Common Cuttlefish (Sepia officinalis). In Advances in Marine Biology (ed. Lesser, M.) vol. 65, 1–65 (Academic Press, Cambridge, 2013).

85.
Vidal, E. A. G., Roberts, M. J. & Martins, R. S. Yolk utilization, metabolism and growth in reared Loligo vulgaris reynaudii paralarvae. Aquat. Living Resour. 18, 385–393 (2005).
Article  Google Scholar 

86.
Bouchaud, O. Energy consumption of the cuttlefish Sepia officinalis L. (Mollusca: Cephalopoda) during embryonic development, preliminary results. Bull. Mar. Sci. 49, 333–340 (1991).
Google Scholar 

87.
Laptikhovsky, V. Latitudinal and bathymetric trends in egg size variation: A new look at Thorson’s and Rass’s rules. Mar. Ecol. 27, 7–14 (2006).
ADS  Article  Google Scholar 

88.
Hengl, T., Sierdsema, H., Radović, A. & Dilo, A. Spatial prediction of species’ distributions from occurrence-only records: Combining point pattern analysis, ENFA and regression-kriging. Ecol. Model. 220, 3499–3511 (2009).
Article  Google Scholar 

89.
Clarke, M. R. The role of cephalopods in the world’s oceans: general conclusions and the future. Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 351, 1105–1112 (1996).
ADS  Article  Google Scholar 

90.
Marmion, M., Parviainen, M., Luoto, M., Heikkinen, R. K. & Thuiller, W. Evaluation of consensus methods in predictive species distribution modelling. Divers. Distrib. 15, 59–69 (2009).
Article  Google Scholar 

91.
Kissling, W. D. et al. Towards novel approaches to modelling biotic interactions in multispecies assemblages at large spatial extents: Modelling multispecies interactions. J. Biogeogr. 39, 2163–2178 (2012).
Article  Google Scholar 

92.
Clark, J. S., Gelfand, A. E., Woodall, C. & Zhu, K. More than the sum of the parts: Forest climate response from joint species distribution models. Ecol. Appl. 24, 990–999 (2014).
PubMed  Article  Google Scholar 

93.
Harris, D. J. Generating realistic assemblages with a joint species distribution model. Methods Ecol. Evol. 6, 465–473 (2015).
Article  Google Scholar 

94.
Nogués-Bravo, D. Predicting the past distribution of species climatic niches. Glob. Ecol. Biogeogr. 18, 521–531 (2009).
Article  Google Scholar 

95.
Lee, Q., Thorson, J. T., Gertseva, V. V. & Punt, A. E. The benefits and risks of incorporating climate-driven growth variation into stock assessment models, with application to Splitnose Rockfish (Sebastes diploproa). ICES J. Mar. Sci. 75, 245–256 (2018).
Article  Google Scholar 

96.
Colléter, M., Gascuel, D., Ecoutin, J.-M. & Tito de Morais, L. Modelling trophic flows in ecosystems to assess the efficiency of marine protected area (MPA), a case study on the coast of Sénégal. Ecol. Model. 232, 1–13 (2012).
Article  Google Scholar 

97.
Allen, K. R. Relation between production and biomass. J. Fish. Res. Board Can. 28, 1573–1581 (1971).
Article  Google Scholar 

98.
FAO. Review of the state of world marine fishery resources. FAO Fish. Aquac. Tech. Pap. 334 (2011).

99.
Cheung, W. W. L. et al. Transform high seas management to build climate resilience in marine seafood supply. Fish Fish. 18, 254–263 (2016).
Article  Google Scholar 

100.
Sumaila, U. R., Cheung, W. W. L., Lam, V. W. Y., Pauly, D. & Herrick, S. Climate change impacts on the biophysics and economics of world fisheries. Nat. Clim. Change 1, 449–456 (2011).
ADS  Article  Google Scholar 

101.
Barange, M. et al. Impacts of climate change on marine ecosystem production in societies dependent on fisheries. Nat. Clim. Change 4, 211–216 (2014).
ADS  Article  Google Scholar 

102.
Ojea, E., Lester, S. E. & Salgueiro-Otero, D. Adaptation of fishing communities to climate-driven shifts in target species. One Earth 2, 544–556 (2020).
Article  Google Scholar  More

Site description, chronology, and sampling
A detailed description of the site, coring methods, age-depth model reconstruction, and sampling strategy can be found in Alsos et al.41. Briefly, Lake Øvre Æråsvatnet is located on Andøya, Northern Norway (69.25579°N, 16.03517°E) (Fig. 1a, b). In 2013, two cores were collected from the deepest sediments, AND10 and AND11, which were stored at 4 °C prior to sampling. Macrofossil remains were dated, with those from AND10 all dating to within the LGM. For the longer core AND11, a Bayesian age-depth model was required to estimate the age of each layer41. In this study, we selected one sample of LGM sediments from each of the two cores. According to the Bayesian age-depth model, sample Andøya_LGM_B, from 1102 cm depth in AND11, was dated to a median age of 17,700 (range: 20,200–16,500) cal yr BP. The age of Andøya_LGM_A, from 938 cm depth in AND10, was estimated at 19,500 cal yr BP, based on the interpolated median date between two adjacent macrofossils (20 cm above: 19,940–18,980 cal yr BP, 30 cm below: 20,040–19,000 cal yr BP). As Andøya_LGM_A falls within the age range of Andøya_LGM_B, we consider the samples to be broadly contemporaneous.
Sampling, DNA extraction, library preparation, and sequencing
The two cores were subsampled at the selected layers under clean conditions, in a dedicated ancient DNA laboratory at The Arctic University Museum of Norway in Tromsø. We extracted DNA from 15 g of sediment following the Taberlet phosphate extraction protocol29 in the same laboratory. We shipped a 210 µL aliquot of each DNA extract to the ancient DNA dedicated laboratories at the Centre for GeoGenetics (University of Copenhagen, Denmark) for double-stranded DNA library construction. The extracts were concentrated to 80 µL using Amicon Ultra-15 30 kDa centrifugal filters (Merck Millipore, Darmstadt, Germany) and half of each extract (40 µL, totalling between 31.7 and 36.0 ng of DNA) was converted into Illumina-compatible libraries using established protocols10. Each library was dual-indexed via 12 cycles of PCR. The libraries were then purified using the AmpureBead protocol (Beckman Coulter, Indianapolis, IN, USA), adjusting the volume ratio to 1:1.8 library:AmpureBeads, and quantified using a BioAnalyzer (Agilent, Santa Clara, CA, USA). The indexed libraries were pooled equimolarly and sequenced on a lane of the Illumina HiSeq 2500 platform using 2 × 80 cycle paired-end chemistry.
Raw read filtering
For each sample, we merged and adapter-trimmed the paired-end reads with SeqPrep (https://github.com/jstjohn/SeqPrep/releases, v1.2) using default parameters. We only retained the resulting merged sequences, which were then filtered with the preprocess function of the SGA toolkit v0.10.15 (ref. 57) by the removal of those shorter than 35 bp or with a DUST complexity score > 1.
Metagenomic analysis of the sequence data
We first sought to obtain an overview of the taxonomic composition of the samples and therefore carried out a BLAST-based metagenomic analysis on the two filtered sequence datasets. To make the datasets more computationally manageable, we subsampled the first and last one-million sequences from the filtered dataset of each sample and analysed each separately. The data subsets were each identified against the NCBI nucleotide database (release 223) using the blastn function from the NCBI-BLAST+ suite v2.2.18+58 under default settings. For each sample, the results from the two subsets were checked for internal consistency, merged into one dataset, and loaded into MEGAN v6.12.3 (ref. 59). Analysis and visualization of the Last Common Ancestor (LCA) was carried out for the taxonomic profile using the following settings: min score = 35, max expected = 1.0E−5, min percent identity = 95%, top percent = 10%, min support percentage = 0.01, LCA = naive, min percent sequence to cover = 95%. We define sequences as the reads with BLAST hits assigned to taxa post-filtering, thus ignoring “unassigned” and “no hit” categories.
Alignment to reference genome panels
We mapped our filtered data against three different reference panels to help improve taxonomic identifications and provide insight into the sequence abundance of the identified taxa (Supplementary Data 2 and 3). The first reference panel consisted of 42 nuclear genomes that included genera expected from Northern Norway, exotic/implausible taxa for LGM Andøya, human, six Nannochloropsis species, and four strains of Mycobacterium. The inclusion of exotic taxa was to give an indication of the background spurious mapping rate, which can result from mappings to conserved parts of the genome and/or short and damaged ancient DNA molecules22,23. We included Nannochloropsis, Mycobacterium, and human genomes, due to their overrepresentation in the BLAST-based metagenomic analysis. The other two reference panels were based on either all 8486 mitochondrial or 2495 chloroplast genomes on NCBI GenBank (as of January 2018). The chloroplast dataset was augmented with 247 partial or complete chloroplast genomes generated by the PhyloNorway project60 for 2742 chloroplast genomes in total. The filtered data were mapped against each reference genome or organellar genome set individually using bowtie2 v2.3.4.1 (ref. 61) under default settings. The resulting bam files were processed with SAMtools v0.1.19 (ref. 62). We removed unmapped sequences with SAMtools view and collapsed PCR duplicate sequences with SAMtools rmdup.
For the nuclear reference panel, we reduced potential spurious or nonspecific sequence mappings by comparing the mapped sequences to both the aligned reference genome and the NCBI nucleotide database using NCBI-BLAST+, following the method used by Graham et al.9, as modified by Wang et al.11. The sequences were aligned using the following NCBI-BLAST+ settings: num_alignments = 100 and perc_identity = 90. Sequences were retained if they had better alignments, based on bit score, to reference genomes as compared to the NCBI nucleotide database. If a sequence had a better or equal match against the NCBI nucleotide database, it was removed, unless the LCA of the highest NCBI nucleotide bit score was from the same genus as the reference genome (based on the NCBI taxonID). To standardize the relative mapping frequencies to genomes of different size, we calculated the number of retained mapped sequences per Mb of genome sequence.
The sequences mapped against the chloroplast and mitochondrial reference panels were filtered and reported in a different manner than the nuclear genomes. First, to exclude any non-eukaryotic sequences, we used NCBI-BLAST+ to search sequence taxonomies and retained sequences if the LCA was, or was within, Eukaryota. Second, for the sequences that were retained, the LCA was calculated and reported in order to summarize the mapping results across the organelle datasets. LCAs were chosen as the reference sets are composed of multiple genera.
Within the Nannochloropsis nuclear reference alignments, the relative mapping frequency was highest for N. limnetica. In addition, the relative mapping frequency for other Nannochloropsis taxa was higher than those observed for the exotic taxa. This could represent the mapping of sequences that are conserved between Nannochloropsis genomes or suggest the presence of multiple Nannochloropsis taxa in a community sample. We therefore cross-compared mapped sequences to determine the number of uniquely mapped sequences per reference genome. First, we individually remapped the filtered data to six available Nannochloropsis nuclear genomes, the accession codes of which are provided in Supplementary Data 2. For each sample, we then calculated the number of sequences that uniquely mapped to, or overlapped, between each Nannochloropsis genome. We repeated the above analysis with six available chloroplast sequences (Supplementary Data 2) to get a comparable overlap estimation for the chloroplast genome.
Reconstruction of the Andøya Nannochloropsis community organellar palaeogenomes
To place the Andøya Nannochloropsis community taxon into a phylogenetic context, and provide suitable reference sequences for variant calling, we reconstructed environmental palaeogenomes for the Nannochloropsis mitochondria and chloroplast. First, the raw read data from both samples were combined into a single dataset and re-filtered with the SGA toolkit to remove sequences shorter than 35 bp, but retain low complexity sequences to assist in the reconstruction of low complexity regions in the organellar genomes. This re-filtered sequence dataset was used throughout the various steps for environmental palaeogenome reconstruction.
The re-filtered sequence data were mapped onto the N. limnetica reference chloroplast genome (NCBI GenBank accession: NC_022262.1) with bowtie2 using default settings. SAMtools was used to remove unmapped sequences and PCR duplicates, as above. We generated an initial consensus genome from the resulting bam file with BCFtools v1.9 (ref. 62), using the mpileup, call, filter, and consensus functions. For variable sites, we produced a majority-rule consensus using the –variants-only and –multiallelic-caller options, and for uncovered sites the reference genome base was called. The above steps were repeated until the consensus could no longer be improved. The re-filtered sequence data was then remapped onto the initial consensus genome sequence with bowtie2, using the above settings. The genomecov function from BEDtools v2.17.0 (ref. 63) was used to identify gaps and low coverage regions in the resulting alignment.
We attempted to fill the identified gaps, which likely consisted of diverged or difficult-to-assemble regions. For this, we assembled the re-filtered sequence dataset into de novo contigs with the MEGAHIT pipeline v1.1.4 (ref. 64), using a minimum k-mer length of 21, a maximum k-mer length of 63, and k-mer length increments of six. The MEGAHIT contigs were then mapped onto the initial consensus genome sequence with the blastn tool from the NCBI-BLAST+ toolkit. Contigs that covered the gaps identified by BEDtools were incorporated into the initial consensus genome sequence, unless a blast comparison against the NCBI nucleotide database suggested a closer match to non-Nannochloropsis taxa. We repeated the bowtie2 gap-filling steps iteratively, using the previous consensus sequence as reference, until a gap-free consensus was obtained. The re-filtered sequence data were again mapped, the resulting final assembly was visually inspected, and the consensus was corrected where necessary. This was to ensure the fidelity of the consensus sequence, which incorporated de novo-assembled contigs that could potentially be problematic, due to the fragmented nature and deaminated sites of ancient DNA impeding accurate assembly65.
Annotation of the chloroplast genome was carried out with GeSeq v1.77 (ref. 66), using the available annotated Nannochloropsis chloroplast genomes (accession codes provided in Supplementary Table 7). The resulting annotated chloroplast was visualized with OGDRAW v1.3.1 (ref. 67).
The same assembly and annotation methods outlined above were used to reconstruct the mitochondrial palaeogenome sequence, where the initial mapping assembly was based on the N. limnetica mitochondrial sequence (NCBI GenBank accession: NC_022256.1). The final annotation was carried out by comparison against all available annotated Nannochloropsis mitochondrial genomes (accession codes provided in Supplementary Table 7).
If the Nannochloropsis sequences derived from more than one taxon, then alignment to the N. limnetica chloroplast genome could introduce reference bias, which would underestimate the diversity of the Nannochloropsis sequences present. We therefore reconstructed Nannochloropsis chloroplast genomes, but using the six available Nannochloropsis chloroplast genome sequences, including N. limnetica, as seed genomes (accession codes for the reference genomes are provided in Supplementary Table 3). The assembly of the consensus sequences followed the same method outlined above, but with two modifications to account for the mapping rate being too low for complete genome reconstruction based on alignment to the non-N. limnetica reference sequences. First, consensus sequences were called with SAMtools, which does not incorporate reference bases into the consensus at uncovered sites. Second, neither additional gap filling nor manual curation was implemented.
Analysis of ancient DNA damage patterns
We checked for the presence of characteristic ancient DNA damage patterns for sequences aligned to three nuclear genomes: human, Nannochloropsis limnetica and Mycobacterium avium. We further analysed damage patterns for sequences aligned to both the reconstructed N. limnetica composite organellar genomes. Damage analysis was conducted with mapDamage v2.0.8 (ref. 68) using the following settings: –merge-reference-sequences and –length = 160.
Assembly of high- and low-frequency variant consensus sequences
The within-sample variants in each reconstructed organellar palaeogenome was explored by creating two consensus sequences, which included either high- or low-frequency variants at multiallelic sites. For each sample, the initial filtered sequence data were mapped onto the reconstructed Nannochloropsis chloroplast palaeogenome sequence with bowtie2 using default settings. Unmapped and duplicate sequences were removed with SAMtools, as above. We used the BCFtools mpileup, call, and normalize functions to identify the variant sites in the mapped dataset, using the –skip-indels, –variants-only, and –multiallelic-caller options. The resulting alleles were divided into two sets, based on either high- or low-frequency variants. High-frequency variants were defined as those present in the reconstructed reference genome sequence. Both sets were further filtered to only include sites with a quality score of 30 or higher and a coverage of at least half the average coverage of the mapping assembly (minimum coverage: Andøya_LGM_A = 22×, Andøya_LGM_B = 14×). We then generated the high- and low-frequency variant consensus sequences using the consensus function in BCFTools. The above method was repeated for the reconstructed Nannochloropsis mitochondrial genome sequence in order to generate comparable consensus sequences of high- and low-frequency variants (minimum coverage: Andøya_LGM_A = 16×, Andøya_LGM_B = 10×).
Phylogenetic analysis of the reconstructed organellar palaeogenomes
We determined the phylogenetic placement of our high- and low-frequency variant organellar palaeogenomes within Nannochloropsis, using either full mitochondrial and chloroplast genome sequences or three short loci (18S, ITS, rbcL). We reconstructed the 18S and ITS1-5.8S-ITS2 complex using DQ977726.1 (full length) and EU165325.1 (positions 147:1006, corresponding to the ITS complex) as seed sequences following the same approach that was used for the organellar palaeogenome reconstructions, except that the first and last 10 bp were trimmed to account for the lower coverage due to sequence tiling. We then called high and low variant consensus sequences as described above.
We created six alignments using available sequence data from NCBI GenBank (Supplementary Data 4) with the addition of: (1 + 2) the high- and low-frequency variant chloroplast or mitochondrial genome consensus sequences, (3) an ~1100 bp subset of the chloroplast genome for the rbcL alignment, (4 + 5) ~1800 and ~860 bp subsets of the nuclear multicopy complex for the 18S and ITS alignments, respectively, and (6) the reconstructed chloroplast genome consensus sequences derived from the alternative Nannochloropsis genome starting points. Full details on the coordinates of the subsets are provided in Supplementary Data 4. We generated alignments using MAFFT v7.427 (ref. 69) with the maxiterate = 1000 setting, which was used for the construction of a maximum likelihood tree in RAxML v8.1.12 (ref. 70) using the GTRGAMMA model and without outgroup specified. We assessed branch support using 1000 replicates of rapid bootstrapping.
Nannochloropsis variant proportions and haplogroup diversity estimation
To estimate major haplogroup diversity, we calculated the proportions of high and low variants in the sequences aligned to our reconstructed Nannochloropsis mitochondrial and chloroplast genomes. For each sample, we first mapped the initial filtered sequence data onto the high- and low-frequency variant consensus sequences with bowtie2. To avoid potential reference biases, and for each organellar genome, the sequence data were mapped separately against both frequency consensus sequences. The resulting bam files were then merged with SAMtools merge. We removed exact sequence duplicates, which may have been mapped to different coordinates, from the merged bam file by randomly retaining one copy. This step was replicated five times to examine its impact on the estimated variant proportions. After filtering, remaining duplicate sequences—those with identical mapping coordinates—were removed with SAMtools rmdup. We then called variable sites from the duplicate-removed bam files using BCFTools under the same settings as used in the assembly of the high- and low-frequency variant consensus sequences. We restricted our analyses to transversion-only variable positions to remove the impact of ancient DNA deamination artifacts. For each variable site, the proportion of reference and alternative alleles was calculated, based on comparison to the composite N. limnetica reconstructed organellar palaeogenomes. We removed rare alleles occurring at a proportion of More

Viruses and mammalian data
Viral genomic data
Complete sequences of coronaviruses were downloaded from Genbank44. Sequences labelled with the terms: ‘vaccine’, ‘construct’, ‘vector’, ‘recombinant’ were removed from the analyses. In addition, we removed those associated with experimental infections where possible. This resulted in a total of 3264 sequences for 411 coronavirus species or strains (i.e., viruses below species level on NCBI taxonomy tree). Of those, 88 were sequences of coronavirus species, and 307 sequences of strains (in 25 coronavirus species, with total number of species included = 92). Of our included species, six in total were unclassified Coronavirinae (unclassified coronaviruses).
Selection of potential mammalian hosts of coronaviruses
We processed meta-data accompanying all sequences (including partial sequences but excluding vaccination and experimental infections) of coronaviruses uploaded to GenBank to extract information on hosts (to species level) of these coronaviruses. We supplemented these data with species-level hosts of coronaviruses extracted from scientific publications via the ENHanCEd Infectious Diseases Database (EID2)45. This resulted in identification of 313 known terrestrial mammalian hosts of coronaviruses (regardless of whether a complete genome was available or not, n = 185 mammalian species for which an association with a coronavirus with complete genome was identified). We expanded this set of potential hosts by including terrestrial mammalian species in genera containing at least one known host of coronavirus, and which are known to host one or more other virus species (excluding coronaviruses, information of whether the host is associated with a virus were obtained from EID2). This results in total of 876 mammalian species which were selected.
Quantification of viral similarities
We computed three types of similarities between each two viral genomes as summarised below.
Biases and codon usage
We calculated proportion of each nucleotide of the total coding sequence length. We computed dinucleotide and codon biases46 and codon-pair bias, measured as the codon-pair score46,47 in each of the above sequences. This enabled us to produce for each genome sequence (n = 3264) the following feature vectors: nucleotide bias, dinucleotide bias, codon biased and codon-pair bias.
Secondary structure
Following alignment of sequences (using AlignSeqs function in R package Decipher48), we predicted the secondary structure for each sequence using PredictHEC function in the R package Decipher48. We obtained both states (final prediction), and probability of secondary structures for each sequence. We then computed for each 1% of the genome length both the coverage (number of times a structure was predicted) and mean probability of the structure (in the per cent of the genome considered). This enabled us to generate six vectors (length = 100) for each genome representing: mean probability and coverage for each of three possible structures—Helix (H), Beta-Sheet (E) or Coil (C).
Genome dissimilarity (distance)
We calculated pairwise dissimilarity (in effect a hamming distance) between each two sequences in our set using the function DistanceMatrix in the R package Decipher48. We set this function to penalise gap-to-gap and gap-to-letter mismatches.
Similarity quantification
We transformed the feature (traits) vectors described above into similarities matrices between coronaviruses (species or strains). This was achieved by computing cosine similarity between these vectors in each category (e.g., codon-pair usage, H coverage, E probability). Formally, for each genomic feature (n = 10) presented by vector as described above, this similarity was calculated as follows:

$${mathrm{sim}}_{{mathrm{genomic}}_{mathrm{l}}}left( {s_m,s_n} right) = {mathrm{sim}}_{{mathrm{genomic}}_{mathrm{l}}}left( {{mathbf{V}}_m^{f_l},{mathbf{V}}_n^{f_l}} right) = frac{{mathop {sum}nolimits_{i = 1}^d {left( {{mathbf{V}}_m^{f_l}[i] times {mathbf{V}}_n^{f_l}[i]} right)} }}{{sqrt {mathop {sum}nolimits_{i = 1}^d {{mathbf{V}}_m^{f_l}[i]^2} } times sqrt {mathop {sum}nolimits_{i = 1}^d {{mathbf{V}}_n^{f_l}[i]^2} } }}$$
(1)

where sm and sn are two sequences presented by two feature vectors ({mathbf{V}}_m^{f_l}) and ({mathbf{V}}_n^{f_l}) from the genomic feature space fl (e.g., codon-pair bias) of the dimension d (e.g., d = 100 for H coverage).
We then calculated similarity between each pair of virus strains or species (in each category) as the mean of similarities between genomic sequences of the two virus strains or species (e.g., mean nucleotide bias similarity between all sequences of SARS-CoV-2 and all sequences of MERS-CoV presented the final nucleotide bias similarity between SARS-CoV-2 and MERS-CoV). This enabled us to generate 11 genomic features similarity matrices (the above 10 features represented by vectors and genomic similarity matrix) between our input coronaviruses. Supplementary Fig. 1 illustrates the process.
Similarity network fusion (SNF)
We applied SNF49 to integrate the following similarities in order to reduce our viral genomic feature space: (1) nucleotide, dinucleotide, codon and codon-pair usage biases were combined into one similarity matrix—genome bias similarity. And (2) Helix (H), Beta-Sheet (E) or Coil (C) mean probability and coverage similarities (six in total) were combined into one similarity matrix—secondary structure similarity.
SNF applies an iterative nonlinear method that updates every similarity matrix according to the other matrices via nearest neighbour approach and is scalable and is robust to noise and data heterogeneity. The integrated matrix captures both shared and complementary information from multiple similarities.
Quantification of mammalian similarities
We calculated a comprehensive set of mammalian similarities. Table 3 summarises these similarities and provides justification for inclusion. Supplementary Note 1 provides full details.
Table 3 Mammalian phylogenetic, ecological and geospatial similarities.
Full size table

Quantification of network similarities
Network construction
We processed meta-data accompanying all sequences (including partial genome but excluding vaccination and experimental infections) of coronaviruses uploaded to Genbank44 (accessed 4 May 2020) to extract information on hosts (to species level) of these coronaviruses. We supplemented these data with virus–host associations extracted from publications via the EID2 Database45. This resulted in 1669 associations between 1108 coronaviruses and 545 hosts (including non-mammalian hosts). We transformed these associations into a bipartite network linking species and strains of coronaviruses with their hosts.
Quantification of topological features
The above constructed network summarises our knowledge to date of associations between coronaviruses and their hosts, and its topology expresses patterns of sharing these viruses between various hosts and host groups. Our analytical pipeline captures this topology, and relations between nodes in our network, by means of node embeddings. This approach encodes each node (here either a coronavirus or a host) with its own vector representation in a continuous vector space, which, in turn, enables us to calculate similarities between two nodes based on this representation.
We adopted DeepWalk23 to compute vectorised representations for our coronaviruses and hosts from the network connecting them. DeepWalk23 uses truncated random walks to get latent topological information of the network and obtains the vector representation of its nodes (in our case coronaviruses and their hosts) by maximising the probability of reaching a next node (i.e., probability of a virus–host association) given the previous nodes in these walks (Supplementary Note 2 lists further details).
Similarity calculations
Following the application of DeepWalk to compute the latent topological representation of our nodes, we calculated the similarity between two nodes in our network—n (vectorised as N) and m (vectorised as M), by using cosine similarity as follows24,25:

$${mathrm{sim}}_{{mathrm{network}}}left( {n,m} right) = {mathrm{sim}}_{{mathrm{network}}}left( {{mathbf{M}},{mathbf{N}}} right) = frac{{mathop {sum}nolimits_{i = 1}^d {left( {m_i times n_i} right)} }}{{sqrt {mathop {sum}nolimits_{i = 1}^d {m_i^2} } times sqrt {mathop {sum}nolimits_{i = 1}^d {n_i^2} } }}$$
(2)

where d is the dimension of the vectorised representation of our nodes: M, N; and mi and ni are the components of vectors M and N, respectively.
Similarity learning meta-ensemble—a multi-perspective approach
Our analytical pipeline stacks 12 similarity learners into testable meta-ensembles. The constituent learners can be categorised by the following three ‘points of view’ (see also Supplementary Fig. 4 for a visual description):
Coronaviruses—the virus point of view
We assembled three models derived from (a) genome similarity, (b) genome biases and (c) genome secondary structure. Each of these learners gave each coronavirus–mammalian association (( {v_i to m_j} )) a score, termed confidence, based on how similar the coronavirus vi is to known coronaviruses of mammalian species mj, compared to how similar vi is to all included coronaviruses. In other words, if vi is more similar (e.g., based on genome secondary structure) to coronaviruses observed in host mj than it is similar to all coronaviruses (both observed in mj and not), then the association (v_i to m_j) is given a higher confidence score. Conversely, if vi is somewhat similar to coronaviruses observed in mj, and also somewhat similar to viruses not known to infect this particular mammal, then the association (v_i to m_j) is given a medium confidence score. The association (v_i to m_j) is given a lower confidence score if vi is more similar to coronaviruses not known to infect mj than it is similar to coronaviruses observed in this host.
Formally, given an adjacency matrix A of dimensions (left| {mathbf{V}} right| times left| {mathbf{M}} right|) where (left| {mathbf{V}} right|) is number of coronaviruses included in this study (for which a complete genome could be found), and (left| {mathbf{M}} right|) is number of included mammals, such that for each (v_i in {mathbf{V}}) and (m_j in {mathbf{M}}), aij = 1 if an association is known to exist between the virus and the mammal, and 0 otherwise. Then for a similarity matrix simviral corresponding to each of the similarity matrices calculated above, a learner from the viral point of view is defined as follows24,25:

$${mathrm{confidence}}_{{mathrm{viral}}}( {v_i to m_j} ) = frac{{mathop {sum}nolimits_{l = 1,,l ne i}^{left| {mathbf{V}} right|} {( {{mathrm{sim}}_{{mathrm{viral}}}( {v_i,,v_l} ) times a_{lj}} )} }}{{mathop {sum}nolimits_{l = 1,,l ne i}^{left| {mathbf{V}} right|} {{mathrm{sim}}_{{mathrm{viral}}}left( {v_i,,v_l} right)} }}$$
(3)

Mammals—the host point of view
We constructed seven learners from the similarities summarised in Table 3. Each of these learners calculated for every coronavirus–mammalian association (( {v_i to m_j})) a confidence score based on how similar the mammalian species mj is to known hosts of the coronavirus vi, compared to how similar mj is to mammals not associated with vi. For instance, if mj is phylogenetically close to known hosts of vi, and also phylogenetically distant to mammalian species not known to be associated with this coronavirus, then the phylogenetic similarly learner will assign (v_i to m_j) a higher confidence score. However, if mj does not overlap geographically with known hosts of vi, then the geographical overlap learner will assign it a low (in effect 0) confidence score.
Formally, given the above-defined adjacency matrix A, and a similarity matrix simmammalian corresponding to each of the similarity matrices summarised in Table 3, a learner from the mammalian point of view is defined as follows24,25:

$${mathrm{confidence}}_{{mathrm{mammalian}}}( {v_i to m_j} ) = frac{{mathop {sum}nolimits_{l = 1,,l ne j}^{left| {mathbf{M}} right|} {( {{mathrm{sim}}_{{mathrm{mammalian}}}( {m_j,,m_l} ) times a_{il}} )} }}{{mathop {sum}nolimits_{l = 1,,l ne j}^{left| {mathbf{M}} right|} {{mathrm{sim}}_{{mathrm{mammalian}}}( {m_j,,m_l} )} }}$$
(4)

Network—the network point of view
We integrated two learners based on network similarities—one for mammals and one for coronaviruses. Formally, given the adjacency matrix A, our two learners from the network point of view as defined as follows24:

$${mathrm{confidence}}_{{mathrm{network}}_{mathbf{V}}}( {v_i to m_j} ) = frac{{mathop {sum}nolimits_{l = 1,,l ne i}^{left| {mathbf{V}} right|} {( {{mathrm{sim}}_{{mathrm{network}}}( {v_i,,v_l} ) times a_{lj}} )} }}{{mathop {sum}nolimits_{l = 1,,l ne i}^{left| {mathbf{V}} right|} {{mathrm{sim}}_{{mathrm{network}}}( {v_i,,v_l} )} }};;$$
(5)

$${mathrm{confidence}}_{{mathrm{network}}_{mathbf{M}}}( {v_i to m_j} ) = frac{{mathop {sum}nolimits_{l = 1,,l ne j}^{left| {mathbf{M}} right|} {( {{mathrm{sim}}_{{mathrm{network}}}( {m_j,,m_l} ) times a_{il}} )} }}{{mathop {sum}nolimits_{l = 1,,l ne j}^{left| {mathbf{M}} right|} {{mathrm{sim}}_{{mathrm{network}}}( {m_j,,m_l} )} }}$$
(6)

Ensemble construction
We combined the learners described above by stacking them into ensembles (meta-ensembles) using Stochastic Gradient Boosting (GBM). The purpose of this combination is to incorporate the three points of views, as well as varied aspects of the coronaviruses and their mammalian potential hosts, into a generalisable, robust model50. We selected GBM as our stacking algorithm following an assessment of seven machine-learning algorithms using held-out test sets (20% of known associations randomly selected, N = 5—Supplementary Fig. 14). In addition, GBM is known for its ability to handle non-linearity and high-order interactions between constituent learners51, and have been used to predict reservoirs of viruses46 and zoonotic hot-spots51.We performed the training and optimisation (tuning) of these ensembles using the caret R Package52.
Sampling
Our GBM ensembles comprised 100 replicate models. Each model was trained with balanced random samples using tenfold cross-validation (Supplementary Fig. 4). Final ensembles were generated by taking mean predictions (probability) of constituent models. Predictions were calculated form the mean probability at three cut-offs: >0.5 (standard), >0.75 and ≥0.9821. SD from mean probability was also generated and its values subtracted/added to predictions, to illustrate variation in the underlying replicate models.
Validation and performance estimation
We validated the performance of our analytical pipeline externally against 20 held-out test sets. Each test set was generated by splitting the set of observed associations between coronaviruses and their hosts into two random sets: a training set comprising 85% of all known associations and a test set comprising 15% of known associations. These test sets were held-out throughout the processes of generating similarity matrices; similarity learning, and assembling our learners, and were only used for the purposes of estimating performance metrics of our analytical pipeline. This resulted in 20 runs in which our ensemble learnt using only 85% of observed associations between our coronaviruses and their mammalian hosts. For each run, we calculated three performance metrics based on the mean probability across each set of 100 replicate models of the GBM meta-ensembles: AUC, true skill statistics (TSS) and F-score.
AUC is a threshold-independent measure of model predictive performance that is commonly used as a validation metric for host–pathogen predictive models21,46. Use of AUC has been criticised for its insensitivity to absolute predicted probability and its inclusion of a priori untenable prediction51,53, and so we also calculated the TSS (TSS = sensitivity + specificity − 1)54. F-score captures the harmonic mean of the precision and recall and is often used with uneven class distribution. Our approach is relaxed with respect to false positives (unobserved associations), hence the low F-score recorded overall.
We selected three probability cut-offs for our meta-ensemble: 0.50, 0.75 and 0.9821. One extreme of our cut-off range (0.5) maximises the ability of our ensemble to detect known associations (higher AUC, lower F-score). The other (0.9821) is calculated so that 90% of known positives are captured by our ensemble, while reducing the number of additional associations predicted (higher F-score, lower AUC).
Changes in network structure
We quantified the diversity of the mammalian hosts of each coronavirus in our input by computing mean phylogenetic distance between these hosts. Similarly, we captured the diversity of coronaviruses associated with each mammalian species by calculating mean (hamming) distance between the genomes of these coronaviruses. We termed these two metrics: mammalian diversity per virus and viral diversity per mammal, respectively. We aggregated both metrics at the network level by means of simple average. This enabled us to quantify changes in these diversity metrics, at the level of network, with addition of predicted links at three probability cut-offs: >0.5, >0.75 and ≥0.9821.
In addition, we captured changes in the structure of the bipartite network linking CoVs with their mammalian hosts, with the addition of predicted associations, by computing a comprehensive set of structural properties (Supplementary Note 3) at the probability cut-offs mentioned above, and comparing the results with our original network. Here we ignore properties that deterministically change with the addition of links (e.g., degree centrality, connectance; Supplementary Table 2 lists all computed metrics and changes in their values). Instead, we focus on non-trivial structural properties. Specifically, we capture changes in network stability, by measuring its nestedness55,56,57; and we quantify non-independence in interaction patterns by means of C-score58. Supplementary Note 3 provides full definition of these concepts as well as other metrics we computed for our networks.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

arising from: P. B. Ngor et al.; Scientific Reports https://doi.org/10.1038/s41598-018-27340-1 (2018).
As one of the richest sources of fisheries-related data in the lower Mekong basin, the Tonle Sap dai fishery has received considerable attention in the literature in recent years as concerns grow over the impacts of hydropower dams on fisheries, which are important for livelihoods and food security1,2,3.
Ngor et al.4 reported a decline since 2000 in the catch of larger species which tend to occupy higher trophic levels; compensatory increases in the catch of smaller species; and declines in the mean body weight (and length) of common species in the Tonle Sap dai fishery, as evidence of the effects of indiscriminate fishing or “fishing-down” of the multi-species fish assemblage in the lower Mekong basin. We provide evidence below that suggest that these apparent recent changes are more likely to reflect changing hydrological conditions than fishing-down effects, possibly caused by climate change and recently also by hydropower development.
The dai fishery has been reliably monitored since 1997–98. Without explanation, Ngor et al. excluded the first three seasons (1997–98 to 1999–2000) of monitoring data which include one of the driest fishing seasons on record (1998–99). The authors thereby created a time series beginning with the three wettest seasons (largest floods) since monitoring began (2000–1 to 2002–3) that were followed by 12 seasons of variable, but decreasing flows caused by hydropower dam construction, low rainfalls possibly resulting from climate change, and abstractions for agriculture5,6 (Fig. 1).
Figure 1

Source: Mekong River Commission Secretariat.

The flood index (FI) or flood pulse14 in the Tonle Sap Great Lake System (1997/08–2014/15). The FI is a measure of the flood extent and duration, calculated as the sum of the flooded area days above the mean flooded area from April to March of the following year2. Whilst highly variable, a downward decline (p-value = 0.06) in the FI is observed between 2000/01 (Year 2001) and 2014/15 (Year 2015) shown by solid circles. Adding the most recent data for 2016–2018 (not shown here), confirmed that a downward linear trend in the FI since the 2000/01 season is statistically significant (p-value  45 cm) excluding those with zero catch in any year. These 28 species formed approximately 16% of the total catch during the study period. We also found negative regression coefficients for all 28 species, supporting the findings of Ngor et al. However, the combined annual catch of these 28 species did not decline significantly through time (R2 = 0.22; p-value = 0.07).
We did however find that the combined annual catch of these 28 larger species varied significantly with the annual flood index (FI)—a measure of flood extent and duration (R2 = 0.46; p-value  45 cm) species and the flood index (R2 = 0.46; p-value  More

1.
Paterson, S. et al. Antagonistic coevolution accelerates molecular evolution. Nature 464, 275–278 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 
2.
Schulte, R. D., Makus, C., Hasert, B., Michiels, N. K. & Schulenburg, H. Multiple reciprocal adaptations and rapid genetic change upon experimental coevolution of an animal host and its microbial parasite. Proc. Natl Acad. Sci. USA 107, 7359–7364 (2010).
CAS  PubMed  Article  PubMed Central  Google Scholar 

3.
Koskella, B. & Lively, C. M. Evidence for negative frequency-dependent selection during experimental coevolution of a freshwater snail and a sterilizing trematode. Evolution 63, 2213–2221 (2009).
PubMed  Article  PubMed Central  Google Scholar 

4.
Decaestecker, E. et al. Host–parasite ‘Red Queen’ dynamics archived in pond sediment. Nature 450, 870–873 (2007).
CAS  PubMed  Article  PubMed Central  Google Scholar 

5.
Gómez, P. & Buckling, A. Bacteria–phage antagonistic coevolution in soil. Science 332, 106–109 (2011).
PubMed  Article  CAS  PubMed Central  Google Scholar 

6.
Refardt, D. & Ebert, D. Inference of parasite local adaptation using two different fitness components. J. Evol. Biol. 20, 921–929 (2007).
CAS  PubMed  Article  PubMed Central  Google Scholar 

7.
Duffy, M. A., Hall, S. R., Cáceres, C. E. & Ives, A. R. Rapid evolution, seasonality, and the termination of parasite epidemics. Ecology 90, 1441–1448 (2009).
PubMed  Article  Google Scholar 

8.
Springer, Y. P. Clinical resistance structure and pathogen local adaptation in a serpentine flax–flax rust interaction. Evolution 61, 1812–1822 (2007).
PubMed  Article  PubMed Central  Google Scholar 

9.
Tack, A. J. M., Laine, A.-L., Burdon, J. J., Bissett, A. & Thrall, P. H. Below-ground abiotic and biotic heterogeneity shapes above-ground infection outcomes and spatial divergence in a host–parasite interaction. New Phytol. 207, 1159–1169 (2015).
PubMed  PubMed Central  Article  Google Scholar 

10.
Wolinska, J. & King, K. C. Environment can alter selection in host–parasite interactions. Trends Parasitol. 25, 236–244 (2009).
PubMed  Article  PubMed Central  Google Scholar 

11.
Auld, S. K. J. R., Hall, S. R., Ochs, J. H., Sebastian, M. & Duffy, M. A. Predators and patterns of within-host growth can mediate both among-host competition and evolution of transmission potential of parasites. Am. Nat. 184, S77–S90 (2014).
PubMed  Article  PubMed Central  Google Scholar 

12.
Wright, R. C. T., Brockhurst, M. A. & Harrison, E. Ecological conditions determine extinction risk in co-evolving bacteria–phage populations. BMC Evol. Biol. 16, 227 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

13.
Duffy, M. A. et al. Ecological context influences epidemic size and parasite-driven evolution. Science 335, 1636–1638 (2012).
CAS  PubMed  Article  Google Scholar 

14.
Auld, S. K. J. R. & Brand, J. Environmental variation causes different (co) evolutionary routes to the same adaptive destination across parasite populations. Evol. Lett. 1, 245–254 (2017).
PubMed  PubMed Central  Article  Google Scholar 

15.
Su, M. & Boots, M. The impact of resource quality on the evolution of virulence in spatially heterogeneous environments. J. Theor. Biol. 416, 1–7 (2017).
PubMed  Article  PubMed Central  Google Scholar 

16.
Auld, S. K. J. R. & Tinsley, M. C. The evolutionary ecology of complex lifecycle parasites: linking phenomena with mechanisms. Heredity 114, 125–132 (2015).
CAS  PubMed  Article  PubMed Central  Google Scholar 

17.
Cardon, M., Loot, G., Grenouillet, G. & Blanchet, S. Host characteristics and environmental factors differentially drive the burden and pathogenicity of an ectoparasite: a multilevel causal analysis. J. Anim. Ecol. 80, 657–667 (2011).
PubMed  Article  PubMed Central  Google Scholar 

18.
Mahmud, M. A., Bradley, J. E. & MacColl, A. D. C. Abiotic environmental variation drives virulence evolution in a fish host–parasite geographic mosaic. Funct. Ecol. 31, 2138–2146 (2017).
Article  Google Scholar 

19.
Arruda, J. A., Marzolf, G. R. & Faulk, R. T. The role of suspended sediments in the nutrition of zooplankton in turbid reservoirs. Ecology 64, 1225–1235 (1983).
Article  Google Scholar 

20.
Mostowy, R. & Engelstädter, J. The impact of environmental change on host–parasite coevolutionary dynamics. Proc. R. Soc. B 278, 2283–2292 (2011).
PubMed  Article  Google Scholar 

21.
Thompson, J. N. The Geographic Mosaic of Coevolution (Univ. Chicago Press, 2005).

22.
Brett, M. T. Chaoborus and fish-mediated influences on Daphnia longispina population structure, dynamics and life history strategies. Oecologia 89, 69–77 (1992).
PubMed  Article  Google Scholar 

23.
Goss, L. B. & Bunting, D. L. Daphnia development and reproduction: responses to temperature. J. Therm. Biol. 8, 375–380 (1983).
Article  Google Scholar 

24.
Luijckx, P., Fienberg, H., Duneau, D. & Ebert, D. A matching-allele model explains host resistance to parasites. Curr. Biol. 23, 1085–1088 (2013).
CAS  PubMed  Article  Google Scholar 

25.
Bento, G. et al. The genetic basis of resistance and matching-allele interactions of a host–parasite system: the Daphnia magna–Pasteuria ramosa model. PLoS Genet. 13, e1006596 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

26.
Grosberg, R. K. Mate selection and the evolution of highly polymorphic self/nonself recognition genes. Science 289, 2111–2114 (2000).
CAS  PubMed  Article  Google Scholar 

27.
Hutchinson, G. E. The Ecological Theater and the Evolutionary Play (Yale Univ. Press, 1965).

28.
Stuart, Y. E. et al. Contrasting effects of environment and genetics generate a continuum of parallel evolution. Nat. Ecol. Evol. 1, 0158 (2017).
Article  Google Scholar 

29.
Klüttgen, B., Dülmer, U., Engels, M. & Ratte, H. ADaM, an artificial freshwater for the culture of zooplankton. Water Res. 28, 743–746 (1994).
Article  Google Scholar 

30.
Ebert, D., Zschokke-Rohringer, C. D. & Carius, H. J. Within- and between-population variation for resistance of Daphnia magna to the bacterial endoparasite Pasteuria ramosa. Proc. R. Soc. B 265, 2127–2134 (1998).
Article  Google Scholar 

31.
Auld, S. K. J. R. & Brand, J. Simulated climate change, epidemic size, and host evolution across host–parasite populations. Glob. Change Biol. 23, 5045–5053 (2017).
Article  Google Scholar 

32.
Holm, S. A simple sequentially rejective multiple test procedure. Scand. J. Stat. 6, 65–70 (1979).
Google Scholar 

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

34.
Brereton, R. G. & Lloyd, G. R. Re-evaluating the role of the Mahalanobis distance measure. J. Chemom. 30, 134–143 (2016).
CAS  Article  Google Scholar 

35.
D’Orazio, M. StatMatch: Statistical Matching or Data Fusion. R package version 1.4.0 (2019).

36.
Goslee, S. C. & Urban, D. L. The ecodist package for dissimilarity-based analysis of ecological data. J. Stat. Softw. 22, 1–22 (2007).
Article  Google Scholar 

37.
Lefcheck, J. S. piecewiseSEM: piecewise structural equation modelling in R for ecology, evolution and systematics. Methods Ecol. Evol. 7, 573–579 (2016).
Article  Google Scholar 

38.
Auld, S. K. J. R., Wilson, P. J. & Little, T. J. Rapid change in parasite infection traits over the course of an epidemic in a wild host–parasite population. Oikos 123, 232–238 (2014).
Article  Google Scholar 

39.
Shocket, M. S. et al. Parasite rearing and infection temperatures jointly influence disease transmission and shape seasonality of epidemics. Ecology 99, 1975–1987 (2018).
PubMed  Article  PubMed Central  Google Scholar 

40.
Duncan, A. B., Mitchell, S. E. & Little, T. J. Parasite-mediated selection and the role of sex and diapause in Daphnia. J. Evol. Biol. 19, 1183–1189 (2006).
CAS  PubMed  Article  Google Scholar 

41.
Auld, S. K. J. R. et al. Variation in costs of parasite resistance among natural host populations. J. Evol. Biol. 26, 2479–2486 (2013).
CAS  PubMed  Article  Google Scholar 

42.
Laine, A.-L. Evolution of host resistance: looking for coevolutionary hotspots at small spatial scales. Proc. R. Soc. B 273, 267–273 (2006).
PubMed  Article  Google Scholar 

43.
Lohse, K., Gutierrez, A. & Kaltz, O. Experimental evolution of resistance in Paramecium caudatum against the bacterial parasite Holospora undulata. Evolution 60, 1177–1186 (2006).
Article  Google Scholar 

44.
Duffy, M. A. & Sivars-Becker, L. Rapid evolution and ecological host–parasite dynamics. Ecol. Lett. 10, 44–53 (2007).
PubMed  Article  Google Scholar 

45.
Brewer, M. J., Butler, A. & Cooksley, S. L. The relative performance of AIC, AICC and BIC in the presence of unobserved heterogeneity. Methods Ecol. Evol. 7, 679–692 (2016).
Article  Google Scholar 

46.
Shipley, B. A new inferential test for path models based on directed acyclic graphs. Struct. Equ. Model. 7, 206–218 (2000).
Article  Google Scholar  More

1.
Mayr E. Populations, species, and evolution: an abridgment of animal species and evolution. Cambridge: Belknap Press of Harvard University Press; 1970.
2.
Bickford D, Lohman DJ, Sodhi NS, Ng PKL, Meier R, Winker K, et al. Cryptic species as a window on diversity and conservation. Trends Ecol Evol. 2007;22:148–55.
PubMed  Article  Google Scholar 

3.
Fišer C, Robinson CT, Malard F. Cryptic species as a window into the paradigm shift of the species concept. Mol Ecol. 2018;27:613–35.
PubMed  Article  Google Scholar 

4.
Struck TH, Feder JL, Bendiksby M, Birkeland S, Cerca J, Gusarov VI, et al. Finding evolutionary processes hidden in cryptic species. Trends Ecol Evol. 2018;33:153–63.
PubMed  Article  Google Scholar 

5.
Sarno D, Kooistra WHCF, Medlin LK, Percopo I, Zingone A. Diversity in the genus Skeletonema (Bacillariophyceae). II. An assessment of the taxonomy of S. costatum-like species with the description of four new species. J Phycol. 2005;41:151–76.
Article  Google Scholar 

6.
Gaonkar CC, Kooistra WHCF, Lange CB, Montresor M, Sarno D. Two new species in the Chaetoceros socialis complex (Bacillariophyta): C. sporotruncatus and C. dichatoensis, and characterization of its relatives. J Phycol. 2017;53:889–907.
CAS  PubMed  Article  Google Scholar 

7.
Li Y, Boonprakob A, Gaonkar CC, Kooistra WHCF, Lange CB, Hernández-Becerril D, et al. Diversity in the globally distributed diatom genus Chaetoceros (Bacillariophyceae): three new species from warm-temperate waters. PLoS ONE. 2017;12:e0168887.
PubMed  PubMed Central  Article  CAS  Google Scholar 

8.
Finlay BJ, Clarke KJ. Ubiquitous dispersal of microbial species. Nature. 1999;400:828.
CAS  Article  Google Scholar 

9.
Finlay BJ, Fenchel T. Divergent perspectives on protist species richness. Protist. 1999;150:229–33.
CAS  PubMed  Article  Google Scholar 

10.
Fenchel T, Finlay BJ. The ubiquity of small species: patterns of local and global diversity. Bioscience. 2004;54:777.
Article  Google Scholar 

11.
Fenchel T. Cosmopolitan microbes and their ‘cryptic’ species. Aquat Microb Ecol. 2005;41:49–54.
Article  Google Scholar 

12.
Miglietta MP, Faucci A, Santini F. Speciation in the sea: overview of the symposium and discussion of future directions. Integr Comp Biol. 2011;51:449–55.
PubMed  Article  Google Scholar 

13.
Kooistra WHCF, Sarno D, Balzano S, Gu H, Andersen RA, Zingone A. Global diversity and biogeography of Skeletonema species (Bacillariophyta). Protist. 2008;159:177–93.
CAS  PubMed  Article  Google Scholar 

14.
Nanjappa D, Audic S, Romac S, Kooistra WHCF, Zingone A. Assessment of species diversity and distribution of an ancient diatom lineage using a DNA metabarcoding approach. PLoS ONE. 2014;9:e103810.
PubMed  PubMed Central  Article  CAS  Google Scholar 

15.
Kaczmarska I, Mather L, Luddington IA, Muise F, Ehrman JM. Cryptic diversity in a cosmopolitan diatom known as Asterionellopsis glacialis (Fragilariaceae): Implications for ecology, biogeography, and taxonomy. Am J Bot. 2014;101:267–86.
PubMed  Article  Google Scholar 

16.
Zhao Y, Yi Z, Gentekaki E, Zhan A, Al-Farraj SA, Song W. Utility of combining morphological characters, nuclear and mitochondrial genes: An attempt to resolve the conflicts of species identification for ciliated protists. Mol Phylogenet Evol. 2016;94:718–29.
PubMed  Article  Google Scholar 

17.
Weiner A, Aurahs R, Kurasawa A, Kitazato H, Kucera M. Vertical niche partitioning between cryptic sibling species of a cosmopolitan marine planktonic protist. Mol Ecol. 2012;21:4063–73.
PubMed  Article  Google Scholar 

18.
Lamari N, Ruggiero MV, d’Ippolito G, Kooistra WHCF, Fontana A, Montresor M. Specificity of lipoxygenase pathways supports species delineation in the marine diatom genus Pseudo-nitzschia. PLoS ONE. 2013;8:e73281.
CAS  PubMed  PubMed Central  Article  Google Scholar 

19.
Škaloud P, Friedl T, Hallmann C, Beck A, Dal Grande F. Taxonomic revision and species delimitation of coccoid green algae currently assigned to the genus Dictyochloropsis (Trebouxiophyceae, Chlorophyta). J Phycol. 2016;52:599–617.
PubMed  Article  CAS  Google Scholar 

20.
de Jesus PB, Costa AL, de Castro Nunes JM, Manghisi A, Genovese G, Morabito M, et al. Species delimitation methods reveal cryptic diversity in the Hypnea cornuta complex (Cystocloniaceae, Rhodophyta). Eur J Phycol. 2019;54:135–53.
Article  CAS  Google Scholar 

21.
Díaz-Tapia P, Ly M, Verbruggen H. Extensive cryptic diversity in the widely distributed Polysiphonia scopulorum (Rhodomelaceae, Rhodophyta): molecular species delimitation and morphometric analyses. Mol Phylogenet Evol. 2020;152:106909.
PubMed  Article  Google Scholar 

22.
Huson DH, Rupp R, Scornavacca C. Phylogenetic networks. Cambridge: Cambridge University Press; 2009.

23.
Huson DH, Bryant D. Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006;23:254–67.
CAS  PubMed  Article  Google Scholar 

24.
Solís-Lemus C, Yang M, Ané C. Inconsistency of species tree methods under gene flow. Syst Biol. 2016;65:843–51.
PubMed  Article  Google Scholar 

25.
Deiner K, Bik HM, Mächler E, Seymour M, Lacoursière-Roussel A, Altermatt F, et al. Environmental DNA metabarcoding: Transforming how we survey animal and plant communities. Mol Ecol. 2017;26:5872–95.
PubMed  Article  Google Scholar 

26.
Pawlowski J, Audic S, Adl S, Bass D, Belbahri L, Berney C, et al. CBOL protist working group: barcoding eukaryotic richness beyond the animal, plant, and fungal kingdoms. PLoS Biol. 2012;10:e1001419.
CAS  PubMed  PubMed Central  Article  Google Scholar 

27.
Trobajo R, Mann DG, Clavero E, Evans KM, Vanormelingen P, McGregor RC. The use of partial cox1, rbcL and LSU rDNA sequences for phylogenetics and species identification within the Nitzschia palea species complex (Bacillariophyceae). Eur J Phycol. 2010;45:413–25.
CAS  Article  Google Scholar 

28.
Decelle J, Suzuki N, Mahé F, De Vargas C, Not F. Molecular phylogeny and morphological evolution of the acantharia (Radiolaria). Protist. 2012;163:435–50.
PubMed  Article  Google Scholar 

29.
Stoeck T, Przybos E, Dunthorn M. The D1-D2 region of the large subunit ribosomal DNA as barcode for ciliates. Mol Ecol Resour. 2014;14:458–68.
CAS  PubMed  Article  Google Scholar 

30.
Moniz MBJ, Kaczmarska I. Barcoding of diatoms: nuclear encoded ITS revisited. Protist. 2010;161:7–34.
CAS  PubMed  Article  Google Scholar 

31.
Gile GH, Stern RF, James ER, Keeling PJ. DNA barcoding of chlorarachniophytes using nucleomorph ITS sequences. J Phycol. 2010;46:743–50.
CAS  Article  Google Scholar 

32.
Stern RF, Andersen RA, Jameson I, Küpper FC, Coffroth M-A, Vaulot D, et al. Evaluating the ribosomal internal transcribed spacer (ITS) as a candidate dinoflagellate barcode marker. PLoS ONE. 2012;7:e42780.
CAS  PubMed  PubMed Central  Article  Google Scholar 

33.
Saunders GW. Applying DNA barcoding to red macroalgae: a preliminary appraisal holds promise for future applications. Philos Trans R Soc B Biol Sci. 2005;360:1879–88.

34.
MacGillivary ML, Kaczmarska I. Survey of the efficacy of a short fragment of the rbcL gene as a supplemental DNA barcode for diatoms. J Eukaryot Microbiol. 2011;58:529–36.
CAS  PubMed  Article  Google Scholar 

35.
Zimmermann J, Jahn R, Gemeinholzer B. Barcoding diatoms: evaluation of the V4 subregion on the 18S rRNA gene, including new primers and protocols. Org Divers Evol. 2011;11:173–92.
Article  Google Scholar 

36.
Piredda R, Tomasino MP, D’Erchia AM, Manzari C, Pesole G, Montresor M, et al. Diversity and temporal patterns of planktonic protist assemblages at a Mediterranean Long Term Ecological Research site. FEMS Microbiol Ecol. 2016;93:fiw200.
PubMed  Article  CAS  Google Scholar 

37.
Pawlowski J, Lecroq B. Short rDNA barcodes for species identification in foraminifera. J Eukaryot Microbiol. 2010;57:197–205.
CAS  PubMed  Article  Google Scholar 

38.
Mordret S, Piredda R, Vaulot D, Montresor M. Kooistra WHCF, Sarno D. dinoref: a curated dinoflagellate (Dinophyceae) reference database for the 18S rRNA gene. Mol Ecol Resour. 2018;18:974–87.
CAS  Article  Google Scholar 

39.
Gaonkar CC, Piredda R, Minucci C, Mann DG, Montresor M, Sarno D, et al. Annotated 18S and 28S rDNA reference sequences of taxa in the planktonic diatom family Chaetocerotaceae. PLoS ONE. 2018;13:e0208929.
PubMed  PubMed Central  Article  Google Scholar 

40.
Balzano S, Percopo I, Siano R, Gourvil P, Chanoine M, Marie D, et al. Morphological and genetic diversity of Beaufort Sea diatoms with high contributions from the Chaetoceros neogracilis species complex. J Phycol. 2017;53:161–87.
CAS  PubMed  Article  Google Scholar 

41.
Kopf A, Bicak M, Kottmann R, Schnetzer J, Kostadinov I, Lehmann K, et al. The ocean sampling day consortium. Gigascience. 2015;4. https://doi.org/10.1186/s13742-015-0066-5.

42.
Pesant S, Not F, Picheral M, Kandels-Lewis S, Le Bescot N, Gorsky G, et al. Open science resources for the discovery and analysis of Tara Oceans data. Sci Data. 2015;2:150023.
CAS  PubMed  PubMed Central  Article  Google Scholar 

43.
Yau S, Lopes dos Santos A, Eikrem W, Gérikas Ribeiro C, Gourvil P, Balzano S, et al. Mantoniella beaufortii and Mantoniella baffinensis sp. nov. (Mamiellales, Mamiellophyceae), two new green algal species from the high arctic. J Phycol. 2020;56:37–51.
PubMed  Article  Google Scholar 

44.
Lopes Dos Santos A, Gourvil P, Tragin M, Noël M-H, Decelle J, Romac S, et al. Diversity and oceanic distribution of prasinophytes clade VII, the dominant group of green algae in oceanic waters. ISME J. 2017;11:512–28.
PubMed  Article  Google Scholar 

45.
Kuwata A, Yamada K, Ichinomiya M, Yoshikawa S, Tragin M, Vaulot D, et al. Bolidophyceae, a sister picoplanktonic group of diatoms—a review. Front Mar Sci. 2018;5:370.
Article  Google Scholar 

46.
Segawa T, Matsuzaki R, Takeuchi N, Akiyoshi A, Navarro F, Sugiyama S, et al. Bipolar dispersal of red-snow algae. Nat Commun. 2018;9:1–8.
CAS  Article  Google Scholar 

47.
Ichinomiya M, Dos Santos AL, Gourvil P, Yoshikawa S, Kamiya M, Ohki K, et al. Diversity and oceanic distribution of the Parmales (Bolidophyceae), a picoplanktonic group closely related to diatoms. ISME J. 2016;10:2419–34.
PubMed  PubMed Central  Article  Google Scholar 

48.
Tragin M, Vaulot D. Novel diversity within marine Mamiellophyceae (Chlorophyta) unveiled by metabarcoding. Sci Rep. 2019;9:1–14.
CAS  Article  Google Scholar 

49.
Morard R, Vollmar NM, Greco M, Kucera M. Unassigned diversity of planktonic foraminifera from environmental sequencing revealed as known but neglected species. PLoS ONE. 2019;14:e0213936.
CAS  PubMed  PubMed Central  Article  Google Scholar 

50.
Pinseel E, Janssens SB, Verleyen E, Vanormelingen P, Kohler TJ, Biersma EM, et al. Global radiation in a rare biosphere soil diatom. Nat Commun. 2020;11:1–12.
Article  CAS  Google Scholar 

51.
Hasle GR, Syvertsen EE. Marine diatoms. In: Tomas CR, editor. Identifying marine phytoplankton. San Diego: Academic Press; 1997. pp 5–385.

52.
Kooistra WHCF, Sarno D, Hernández-Becerril DU, Assmy P, Di Prisco C, Montresor M. Comparative molecular and morphological phylogenetic analyses of taxa in the Chaetocerotaceae (Bacillariophyta). Phycologia. 2010;49:471–500.
Article  Google Scholar 

53.
De Luca D, Sarno D, Piredda R, Kooistra WHCF. A multigene phylogeny to infer the evolutionary history of Chaetocerotaceae (Bacillariophyta). Mol Phylogenet Evol. 2019;140:106575.
PubMed  Article  Google Scholar 

54.
Longhurst AR. Toward and ecological geography of the sea. In: Longhurst AR, editor. Ecological geography of the sea. 2nd ed. Cambridge: Academic Press; 2007. pp 1–17.

55.
Stoeck T, Bass D, Nebel M, Christen R, Jones MDM, Breiner H-W, et al. Multiple marker parallel tag environmental DNA sequencing reveals a highly complex eukaryotic community in marine anoxic water. Mol Ecol. 2010;19:21–31.
CAS  PubMed  Article  Google Scholar 

56.
Amaral-Zettler LA, McCliment EA, Ducklow HW, Huse SM. A method for studying protistan diversity using massively parallel sequencing of V9 hypervariable regions of small-subunit ribosomal RNA genes. PLoS ONE. 2009;4:e6372.
PubMed  PubMed Central  Article  CAS  Google Scholar 

57.
Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009;75:7537–41.
CAS  PubMed  PubMed Central  Article  Google Scholar 

58.
De Vargas C, Audic S, Tara Oceans Consortium C, Tara Oceans Expedition P. Total V9 rDNA information organized at the metabarcode level for the Tara Oceans Expedition (2009–12). 2017. PANGAEA. https://doi.org/10.1594/PANGAEA.873277.

59.
Ibarbalz FM, Henry N, Brandão MC, Martini S, Busseni G, Byrne H, et al. Global trends in marine plankton diversity across kingdoms of life. Cell. 2019;179:1084–97.
CAS  PubMed  PubMed Central  Article  Google Scholar 

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

61.
Katoh K, Rozewicki J, Yamada KD. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief Bioinform. 2019;20:1160–6.
CAS  PubMed  PubMed Central  Article  Google Scholar 

62.
Price MN, Dehal PS, Arkin AP. FastTree 2—approximately maximum-likelihood trees for large alignments. PLoS ONE. 2010;5:e9490.
PubMed  PubMed Central  Article  CAS  Google Scholar 

63.
Han MV, Zmasek CM. phyloXML: XML for evolutionary biology and comparative genomics. BMC Bioinform. 2009;10:356.
Article  CAS  Google Scholar 

64.
Templeton AR, Crandall KA, Sing CF. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics. 1992;132:619–33.
CAS  PubMed  PubMed Central  Article  Google Scholar 

65.
Clement M, Posada D, Crandall KA. TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000;9:1657–9.
CAS  PubMed  Article  Google Scholar 

66.
Leigh JW, Bryant D. popart: full‐feature software for haplotype network construction. Methods Ecol Evol. 2015;6:1110–6.
Article  Google Scholar 

67.
Nguyen L-T, Schmidt HA, von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32:268–74.
CAS  Article  Google Scholar 

68.
Kalyaanamoorthy S, Minh BQ, Wong TKF, von Haeseler A, Jermiin LS. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat Methods. 2017;14:587–9.
CAS  PubMed  PubMed Central  Article  Google Scholar 

69.
Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30:2725–9.
CAS  PubMed  PubMed Central  Article  Google Scholar 

70.
Jukes TH, Cantor CR. Evolution of protein molecules. Mamm Protein Metab. 1969;3:21–132.
CAS  Article  Google Scholar 

71.
Meyer CP, Paulay G. DNA barcoding: error rates based on comprehensive sampling. PLoS Biol. 2005;3:e422.
PubMed  PubMed Central  Article  CAS  Google Scholar 

72.
R Core Team. R: a language and environment for statistical computing. 2019. Vienna, Austria: R Foundation for Statistical Computing; 2019.

73.
McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE. 2013;8:e61217.
CAS  PubMed  PubMed Central  Article  Google Scholar 

74.
Wickham H. ggplot2: elegant graphics for data analysis. New York: Springer-Verlag; 2016. https://ggplot2.tidyverse.org.

75.
Becker A, Wilks AR. Maps: draw geographical maps. 2018. https://CRAN.R-project.org/package=maps.

76.
Markmann M, Tautz D. Reverse taxonomy: an approach towards determining the diversity of meiobenthic organisms based on ribosomal RNA signature sequences. Philos Trans R Soc Lond B Biol Sci. 2005;360:1917–24.
CAS  PubMed  PubMed Central  Article  Google Scholar 

77.
López-Escardó D, Paps J, de Vargas C, Massana R, Ruiz-Trillo I, Del Campo J. Metabarcoding analysis on European coastal samples reveals new molecular metazoan diversity. Sci Rep. 2018;8:9106.
PubMed  PubMed Central  Article  CAS  Google Scholar 

78.
Álvarez I, Wendel JF. Ribosomal ITS sequences and plant phylogenetic inference. Mol Phylogenet Evol. 2003;29:417–34.

79.
Alverson AJ, Kolnick L. Intragenomic nucleotide polymorphism among small subunit (18S) rDNA paralogs in the diatom genus Skeletonema (Bacillariophyta). J Phycol. 2005;41:1248–57.
CAS  Article  Google Scholar 

80.
Gaonkar CC, Piredda R, Sarno D, Zingone A, Montresor M, Kooistra WHCF. Species detection and delineation in the marine planktonic diatoms Chaetoceros and Bacteriastrum through metabarcoding: making biological sense of haplotype diversity. Environ Microbiol. 2020;22:1917–29.
CAS  PubMed  Article  Google Scholar 

81.
Cleve PT. Pelagisk Diatomeer från Kattegat. In: Petersen CGJ, editor. Det Videnskabelige Udbytte af Kanonbaaden ‘Hauchs’ Togter i de Danske Have Indefor Skagen, I. Aarene 1883–86. Kjøbenhavn: Andr. Fred. Høst & Sons Forlag; 1889. pp 53–56.

82.
Gran HH. Den Norske Nordhaus-Expedition 1876-1878. Botanik, Protophyta: Diatomaceae, Silicoflagellata og Cilioflagellata. Christiania: Grøndal & Søns; 1897.

83.
De Luca D, Kooistra WHCF, Sarno D, Gaonkar CC, Piredda R. Global distribution and diversity of Chaetoceros (Bacillariophyta, Mediophyceae): integration of classical and novel strategies. PeerJ. 2019;7:e7410.
PubMed  PubMed Central  Article  Google Scholar 

84.
Wang J, Wu J. Occurrence and potential risks of harmful algal blooms in the East China Sea. Sci Total Environ. 2009;407:4012–21.
CAS  PubMed  Article  Google Scholar 

85.
Zhen Y, Mi T, Yu Z. Detection of several harmful algal species by sandwich hybridization integrated with a nuclease protection assay. Harmful Algae. 2009;8:651–7.
CAS  Article  Google Scholar 

86.
Richter DJ, Watteaux R, Vannier T, Leconte J, Frémont P, Reygondeau G, et al. Genomic evidence for global ocean plankton biogeography shaped by large-scale current systems. bioRxiv. 2019. https://doi.org/10.1101/867739.

87.
Sarno D, Kooistra WHCF, Balzano S, Hargraves PE, Zingone A. Diversity in the genus Skeletonema (Bacillariophyceae). III. Phylogenetic position and morphological variability of Skeletonema costatum and Skeletonema grevillei, with the description of Skeletonema ardens sp. nov. J Phycol. 2007;43:156–70.
CAS  Article  Google Scholar 

88.
Hasle GR. The biogeography of some marine planktonic diatoms. Deep Sea Res Oceanogr Abstr. 1976;23:319–338, IN1-IN6.

89.
Pargana A. Functional and molecular diversity of the diatom family Leptocylindraceae. 2017. PhD Thesis, The Open University, Milton Keynes, UK.

90.
Novis PM. Taxonomy of Klebsormidium (Klebsormidiales, Charophyceae) in New Zealand streams and the significance of low-pH habitats. Phycologia. 2006;45:293–301.
Article  Google Scholar 

91.
Rindi F, Guiry MD, López-Bautista JM. Distribution, morphology, and phylogeny of Klebsormidium (Klebsormidiales, Charophyceae) in urban environments in Europe. J Phycol. 2008;44:1529–40.
PubMed  Article  Google Scholar 

92.
Rindi F, Mikhailyuk TI, Sluiman HJ, Friedl T, López-Bautista JM. Phylogenetic relationships in Interfilum and Klebsormidium (Klebsormidiophyceae, Streptophyta). Mol Phylogenet Evol. 2011;58:218–31.
PubMed  Article  Google Scholar 

93.
Škaloud P, Rindi F. Ecological differentiation of cryptic species within an asexual protist morphospecies: a case study of filamentous green alga Klebsormidium (Streptophyta). J Eukaryot Microbiol. 2013;60:350–62.
PubMed  Article  CAS  Google Scholar 

94.
Baas Becking LGM. Geobiologie of Inleiding tot de Milieukunde. The Hague: Van Stockum & Zoon; 1934).

95.
Shapiro BJ, Leducq J-B, Mallet J. What is speciation? PLoS Genet. 2016;12:e1005860.
PubMed  PubMed Central  Article  CAS  Google Scholar 

96.
Godhe A, Rynearson T. The role of intraspecific variation in the ecological and evolutionary success of diatoms in changing environments. Philos Trans R Soc Lond B Biol Sci. 2017;372:20160399.
PubMed  PubMed Central  Article  Google Scholar 

97.
de Vargas C, Norris R, Zaninetti L, Gibb SW, Pawlowski J. Molecular evidence of cryptic speciation in planktonic foraminifers and their relation to oceanic provinces. Proc Natl Acad Sci USA. 1999;96:2864–8.
PubMed  Article  Google Scholar 

98.
Amato A, Kooistra WHCF, Levialdi Ghiron JH, Mann DG, Pröschold T, Montresor M. Reproductive isolation among sympatric cryptic species in marine diatoms. Protist. 2007;158:193–207.
CAS  PubMed  Article  Google Scholar 

99.
Weisse T. Distribution and diversity of aquatic protists: an evolutionary and ecological perspective. Biodivers Conserv. 2007;17:243–59.
Article  Google Scholar 

100.
Vanelslander B, Créach V, Vanormelingen P, Ernst A, Chepurnov VA, Sahan E, et al. Ecological differentiation between sympatric pseudocryptic species in the estuarine benthic diatom Navicula phyllepta (Bacillariophyceae). J Phycol. 2009;45:1278–89.
CAS  PubMed  Article  Google Scholar  More

1.
Haken, J. Transnational Crime in the Developing World (Global Financial Integrity, 2011).
2.
Patel, N. G. et al. Quantitative methods of identifying the key nodes in the illegal wildlife trade network. Proc. Natl Acad. Sci. USA 112, 7948–7953 (2015).
CAS  PubMed  Article  PubMed Central  Google Scholar 

3.
Bager Olsen, M. T. et al. Thirty-six years of legal and illegal wildlife trade entering the USA. Oryx https://doi.org/10.1017/S0030605319000541 (2019).

4.
Tittensor, D. P. et al. Evaluating the relationships between the legal and illegal international wildlife trades. Conserv. Lett. https://doi.org/10.1111/conl.12724 (2020).

5.
Harfoot, M. et al. Unveiling the patterns and trends in 40 years of global trade in CITES-listed wildlife. Biol. Conserv. 223, 47–57 (2018).
Article  Google Scholar 

6.
Scheffers, B. R., Oliveira, B. F., Lamb, I. & Edwards, D. P. Global wildlife trade across the tree of life. Science 76, 71–76 (2019).
Article  CAS  Google Scholar 

7.
Maxwell, S. L., Fuller, R. A., Brooks, T. M. & Watson, J. E. M. Biodiversity: the ravages of guns, nets and bulldozers. Nature 536, 143–145 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

8.
Nielsen, M. R., Meilby, H., Smith-Hall, C., Pouliot, M. & Treue, T. The importance of wild meat in the global south. Ecol. Econ. 146, 696–705 (2018).
Article  Google Scholar 

9.
’t Sas-Rolfes, M., Challender, D. W. S., Hinsley, A., Veríssimo, D. & Milner-Gulland, E. J. Illegal wildlife trade: scale, processes, and governance. Annu. Rev. Environ. Resour. 44, 201–228 (2019).
Article  Google Scholar 

10.
Cooney, R. et al. From poachers to protectors: engaging local communities in solutions to illegal wildlife trade. Conserv. Lett. 10, 367–374 (2017).
Article  Google Scholar 

11.
Bodmer, R. E. & Lozano, E. P. Rural development and sustainable wildlife use in Peru. Conserv. Biol. 15, 1163–1170 (2001).
Article  Google Scholar 

12.
McClenachan, L., Cooper, A. B. & Dulvy, N. K. Rethinking trade-driven extinction risk in marine and terrestrial megafauna. Curr. Biol. 26, 1640–1646 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

13.
Wittemyer, G. et al. Illegal killing for ivory drives global decline in African elephants. Proc. Natl Acad. Sci. USA 111, 13117–13121 (2014).
CAS  PubMed  Article  PubMed Central  Google Scholar 

14.
Brook, S., Van Coeverden De Groot, P. J., Mahood, S. & Long, B. Extinction of the Javan Rhinoceros (Rhinoceros sondaicus) from Vietnam (WWF, 2011).

15.
Heinrich, S. et al. Where did all the pangolins go? International CITES trade in pangolin species. Glob. Ecol. Conserv. 8, 241–253 (2016).
Article  Google Scholar 

16.
Cowlishaw, G., Mendelson, S. & Rowcliffe, J. M. Evidence for post-depletion sustainability in a mature bushmeat market. J. Appl. Ecol. 42, 460–468 (2005).
Article  Google Scholar 

17.
Hutton, J. M. & Webb, G. in The Trade in Wildlife: Regulation for Conservation (ed. Oldfield, S.) Ch. 11 (Earthscan, 2003).

18.
Harris, J. B. C. et al. Using market data and expert opinion to identify overexploited species in the wild bird trade. Biol. Conserv. 187, 51–60 (2015).
Article  Google Scholar 

19.
Milner-Gulland, E. J. & Clayton, L. The trade in babirusas and wild pigs in North Sulawesi, Indonesia. Ecol. Econ. 42, 165–183 (2002).
Article  Google Scholar 

20.
Thuiller, W. et al. Vulnerability of African mammals to anthropogenic climate change under conservative land transformation assumptions. Glob. Change Biol. 12, 424–440 (2006).
Article  Google Scholar 

21.
Benítez-López, A. et al. The impact of hunting on tropical mammal and bird populations. Science 356, 180–183 (2017).
PubMed  Article  CAS  PubMed Central  Google Scholar 

22.
Gibson, L. et al. Primary forests are irreplaceable for sustaining tropical biodiversity. Nature 478, 378–381 (2011).
CAS  PubMed  Article  PubMed Central  Google Scholar 

23.
CITES Trade Statistics Derived from the CITES Trade Database (CITES, 2020).

24.
Phelps, J. & Webb, E. L. ‘Invisible’ wildlife trades: Southeast Asia’s undocumented illegal trade in wild ornamental plants. Biol. Conserv. 186, 296–305 (2015).
Article  Google Scholar 

25.
Davies, G., Schulte-Herbrüggen, B., Kümpel, N. F. & Mendelson, S. Hunting and trapping in Gola Forests, south-eastern Sierra Leone: bushmeat from farm, fallow and forest. Bushmeat Livelihoods Wildl. Manage. Poverty Reduct. https://doi.org/10.1002/9780470692592.ch1 (2008).

26.
Linder, J. M. & Oates, J. F. Differential impact of bushmeat hunting on monkey species and implications for primate conservation in Korup National Park, Cameroon. Biol. Conserv. 144, 738–745 (2011).
Article  Google Scholar 

27.
Gilroy, J. J. & Edwards, D. P. Source–sink dynamics: a neglected problem for landscape-scale biodiversity conservation in the tropics. Curr. Landsc. Ecol. Rep. 2, 51–60 (2017).
Article  Google Scholar 

28.
Watson, J. E. M. et al. Catastrophic declines in wilderness areas undermine global environment targets. Curr. Biol. 26, 2929–2934 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

29.
Marshall, H. et al. Spatio-temporal dynamics of consumer demand driving the Asian songbird crisis. Biol. Conserv. 241, 108237 (2020).
Article  Google Scholar 

30.
Harris, J. B. C. et al. Measuring the impact of the pet trade on Indonesian birds. Conserv. Biol. 31, 394–405 (2017).
PubMed  Article  PubMed Central  Google Scholar 

31.
Carrasco, L. R., Chan, J., Mcgrath, F. L. & Nghiem, L. T. P. Biodiversity conservation in a telecoupled world. Ecol. Soc. 22, 24 (2017).
Article  Google Scholar 

32.
Blundell, A. G. & Mascia, M. B. Discrepancies in reported levels of international wildlife trade. Conserv. Biol. 19, 2020–2025 (2005).
Article  Google Scholar 

33.
Nelson, A. et al. A suite of global accessibility indicators. Sci. Data 6, 266 (2019).
PubMed  PubMed Central  Article  Google Scholar 

34.
Rao, M., Zaw, T., Htun, S. & Myint, T. Hunting for a living: wildlife trade, rural livelihoods and declining wildlife in the Hkakaborazi National Park, North Myanmar. Environ. Manage. 48, 158–167 (2011).
PubMed  Article  PubMed Central  Google Scholar 

35.
Symes, W. S., Edwards, D. P., Miettinen, J., Rheindt, F. E. & Carrasco, L. R. Combined impacts of deforestation and wildlife trade on tropical biodiversity are severely underestimated. Nat. Commun. 9, 4052 (2018).
PubMed  PubMed Central  Article  CAS  Google Scholar 

36.
Hinsley, A. et al. Building sustainability into the Belt and Road Initiative’s traditional Chinese medicine trade. Nat. Sustain. 3, 96–100 (2020).
Article  Google Scholar 

37.
Lechner, A. M., Chan, F. K. S. & Campos-Arceiz, A. Biodiversity conservation should be a core value of China’s Belt and Road Initiative. Nat. Ecol. Evol. 2, 408–409 (2018).
PubMed  Article  PubMed Central  Google Scholar 

38.
Farhadinia, M. S. et al. Belt and Road Initiative may create new supplies for illegal wildlife trade in large carnivores. Nat. Ecol. Evol. 3, 1267–1268 (2019).
PubMed  Article  Google Scholar 

39.
Courchamp, F. et al. Rarity value and species extinction: the anthropogenic Allee effect. PLoS Biol. 4, 2405–2410 (2006).
CAS  Article  Google Scholar 

40.
Jetz, W. & Freckleton, R. P. Towards a general framework for predicting threat status of data-deficient species from phylogenetic, spatial and environmental information. Phil. Trans. R. Soc. B 370, 20140016 (2015).
PubMed  Article  PubMed Central  Google Scholar 

41.
Dulac, J. Global Land Transport Infrastructure Requirements: Estimating Road and Railway Infrastructure Capacity and Costs to 2050 (IEA, 2013).

42.
Vilela, T. et al. A better Amazon road network for people and the environment. Proc. Natl Acad. Sci. USA 117, 7095–7102 (2020).
CAS  PubMed  Article  PubMed Central  Google Scholar 

43.
Challender, D. W. S., Harrop, S. R. & MacMillan, D. C. Towards informed and multi-faceted wildlife trade interventions. Glob. Ecol. Conserv. 3, 129–148 (2015).
Article  Google Scholar 

44.
Papworth, S., Milner-Gulland, E. J. & Slocombe, K. Hunted woolly monkeys (Lagothrix poeppigii) show threat-sensitive responses to human presence. PLoS ONE 8, e62000 (2013).
CAS  PubMed  PubMed Central  Article  Google Scholar 

45.
Toledo-Aceves, T., Garcia-Franco, J. G. & Lopez-Barrera, F. Bromeliad rain: an opportunity for cloud forest management. Ecol. Manage. 329, 129–136 (2014).
Article  Google Scholar 

46.
Challender, D. W. S. et al. Mischaracterization of wildlife trade threat. Science (30 October 2019).

47.
Leung, B. et al. Clustered versus catastrophic global vertebrate declines. Nature https://doi.org/10.1038/s41586-020-2920-6 (2020).

48.
Tierney, M. et al. Use it or lose it: measuring trends in wild species subject to substantial use. Oryx 48, 420–429 (2014).
Article  Google Scholar 

49.
Jachmann, H. Monitoring law-enforcement performance in nine protected areas in Ghana. Biol. Conserv. 141, 89–99 (2008).
Article  Google Scholar 

50.
Cardador, L., Tella, J. L., Anadón, J. D., Abellán, P. & Carrete, M. The European trade ban on wild birds reduced invasion risks. Conserv. Lett. https://doi.org/10.1111/conl.12631 (2019).

51.
Grames, E. M., Stillman, A. N., Tingley, M. W. & Elphick, C. S. An automated approach to identifying search terms for systematic reviews using keyword co-occurrence networks. Methods Ecol. Evol. 10, 1645–1654 (2019).
Article  Google Scholar 

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

53.
Espinosa-Andrade, S. R. Road development, bushmeat extraction and jaguar conservation in Yasuni Biosphere Reserve Ecuador PhD Thesis, Univ. Florida (2012); https://ufdc.ufl.edu/UFE0044108/00001/citation

54.
Schoppe, S., Matillano, J., Cervancia, M. & Acosta, D. Conservation needs of the critically endangered Philippine forest turtle, Siebenrockiella leytensis, in Palawan, Philippines. Chelonian Conserv. Biol. 9, 145–153 (2010).
Article  Google Scholar 

55.
Wilman, H. et al. EltonTraits 1.0: species-level foraging attributes of the world’s birds and mammals. Ecology 95, 2027 (2014).
Article  Google Scholar 

56.
Slavenko, A., Tallowin, O. J. S., Itescu, Y., Raia, P. & Meiri, S. Late Quaternary reptile extinctions: size matters, insularity dominates. Glob. Ecol. Biogeogr. 25, 1308–1320 (2016).
Article  Google Scholar 

57.
Aquino, R. & Calle, A. Evaluation of the conservation status of the game mammals: a comparative model in communities of the Pacaya Samira National Reserve (Loreto, Peru). Rev. Peru. Biol. 10, 163–174 (2003).
Google Scholar 

58.
Aquino, R., Lopez, L., Garcia, G., Charpentier, E. & Arevalo, I. Conservation status and threats to atelids in the northeastern Peruvian Amazon. Primate Conserv. 30, 21–29 (2016).
Google Scholar 

59.
Carrillo, E., Wong, G. & Cuarón, A. D. Monitoring mammal populations in Costa Rican protected areas under different hunting restrictions. Conserv. Biol. 14, 1580–1591 (2000).
Article  Google Scholar 

60.
Cronin, D. T. The Impact of Bushmeat Hunting on the Primates of Bioko Island, Equatorial Guinea (Drexel Univ., 2013).

61.
Dasgupta, S. & Hilaluddin Differential effects of hunting on populations of hornbills and imperial pigeons in the rainforests of the eastern Indian Himalaya. Indian Forester 138, 902–909 (2012).
Google Scholar 

62.
De Thoisy, B., Renoux, F. & Julliot, C. Hunting in northern French Guiana and its impact on primate communities. Oryx 39, 149–157 (2005).
Article  Google Scholar 

63.
Fay, J. M. An elephant (Loxodonta africana) survey using dung counts in the forests of the Central African Republic. J. Trop. Ecol. 7, 25–36 (1991).
Article  Google Scholar 

64.
Gamble, T. & Simons, A. M. Comparison of harvested and nonharvested painted turtle populations. Wildl. Soc. Bull. 32, 1269–1277 (2004).
Article  Google Scholar 

65.
Gonzalez, J. A. Harvesting, local trade, and conservation of parrots in the northeastern Peruvian Amazon. Biol. Conserv. 114, 437–446 (2003).
Article  Google Scholar 

66.
Gray, T. N. E. & Phan, C. Habitat preferences and activity patterns of the larger mammal community in Phnom Prich Wildlife Sanctuary, Cambodia. Raffles Bull. Zool. 59, 311–318 (2011).
Google Scholar 

67.
Hall, J. S. et al. A survey of elephants (Loxodonta africana) in the Kahuzi-Biega National Park lowland sector and adjacent forest in eastern Zaire. Afr. J. Ecol. 35, 213–223 (1997).
Article  Google Scholar 

68.
Klemens, M. W. & Moll, D. An assessment of the effects of commercial exploitation on the pancake tortoise, Malacochersus tornieri, in Tanzania. Chelonian Conserv. Biol. 1, 197–206 (1995).
Google Scholar 

69.
Kümpel, N. F., Milner-Gulland, E. J., Rowcliffe, J. M. & Cowlishaw, G. Impact of gun-hunting on diurnal primates in continental Equatorial Guinea. Int. J. Primatol. 29, 1065–1082 (2008).
Article  Google Scholar 

70.
Magige, F. J., Holmern, T., Stokke, S., Mlingwa, C. & Røskaft, E. Does illegal hunting affect density and behaviour of African grassland birds? A case study on ostrich (Struthio camelus). Biodivers. Conserv. 18, 1361–1373 (2009).
Article  Google Scholar 

71.
Maldonado, A. M. & Peck, M. R. Research and in situ conservation of owl monkeys enhances environmental law enforcement at the Colombian–Peruvian border. Am. J. Primatol. 76, 658–669 (2014).
PubMed  Article  PubMed Central  Google Scholar 

72.
Maldonado, A. M., Nijman, V. & Bearder, S. K. Trade in night monkeys Aotus spp. in the Brazil–Colombia–Peru tri-border area: international wildlife trade regulations are ineffectively enforced. Endanger. Species Res. 9, 143–149 (2009).
Article  Google Scholar 

73.
Muchaal, P. K. & Ngandjui, G. Impact of village hunting on wildlife populations in the Western Dja Reserve, Cameroon. Conserv. Biol. 13, 385–396 (1999).
Article  Google Scholar 

74.
Nuñez-iturril, G. & Howe, H. F. Bushmeat and the fate of trees with seeds dispersed by large primates in a lowland rain forest in western Amazonia. Biotropica 39, 348–354 (2007).
Article  Google Scholar 

75.
O’Brien, S. et al. Decline of the Madagascar radiated tortoise Geochelone radiata due to overexploitation. Oryx 37, 338–343 (2003).
Google Scholar 

76.
Patrick, D. A., Shirk, P., Vonesh, J. R., Harper, E. B. & Howell, K. M. Abundance and roosting ecology of chameleons in the Eastern Arc Mountains of Tanzania and potential effects of harvesting. Herpetol. Conserv. Biol. 6, 422–431 (2011).
Google Scholar 

77.
Poulsen, J. R., Clark, C. J. & Bolker, B. M. Decoupling the effect of logging and hunting on an Afrotropical animal community. Ecol. Appl. 21, 1819–1836 (2011).
CAS  PubMed  Article  Google Scholar 

78.
Remis, M. J. & Kpanou, J. B. Primate and ungulate abundance in response to multi-use zoning and human extractive activities in a central African reserve. Afr. J. Ecol. 49, 70–80 (2010).
Article  Google Scholar 

79.
Rist, J., Milner-Gulland, E. J., Cowlishaw, G. & Rowcliffe, J. M. The importance of hunting and habitat in determining the abundance of tropical forest species in Equatorial Guinea. Biotropica 41, 700–710 (2009).
Article  Google Scholar 

80.
Rovero, F., Mtui, A. S., Kitegile, A. S. & Nielsen, M. R. Hunting or habitat degradation? Decline of primate populations in Udzungwa Mountains, Tanzania: an analysis of threats. Biol. Conserv. 146, 89–96 (2012).
Article  Google Scholar 

81.
Segura, A. & Acevedo, P. The importance of protected and unprotected areas for the Mediterranean spur-thighed tortoise demography in northwest Morocco. Amphib. Reptilia https://doi.org/10.6084/m9.figshare.7751783.v1 (2019).

82.
Sung, Y.-H., Karraker, N. E. & Hau, B. C. H. Demographic evidence of illegal harvesting of an endangered Asian turtle. Conserv. Biol. 27, 1421–1428 (2013).
PubMed  Article  PubMed Central  Google Scholar 

83.
Topp-Jorgensen, E., Nielsen, M. R., Marshall, A. R. & Pedersen, U. Relative densities of mammals in response to different levels of bushmeat hunting in the Udzungwa Mountains, Tanzania. Trop. Conserv. Sci. 2, 70–87 (2009).
Article  Google Scholar 

84.
Yasuoka, H. The sustainability of duiker (Cephalophus spp.) hunting for the Baka hunter–gatherers in southeastern Cameroon. Afr. Study Monogr. 33, 95–120 (2006).
Google Scholar 

85.
Protected Planet: The World Database on Protected Areas (WDPA) (UNEP-WCMC and IUCN, 2020).

86.
QGIS Development Team QGIS Geographic Information System v.3.12.0 (Open Source Geospatial Foundation Project, 2020).

87.
User Manual for the World Database on Protected Areas and World Database on Other Effective Area-Based Conservation Measures: 1.6 (UNEP-WCMC, 2019).

88.
Hedges, L. V., Gurevitch, J. & Curtis, P. S. The meta-analysis of response ratios in experimental ecology. Ecology 80, 1150–1156 (1999).
Article  Google Scholar 

89.
Smithson, M. & Verkuilen, J. A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables. Psychol. Methods 11, 54–71 (2006).
PubMed  Article  PubMed Central  Google Scholar 

90.
Lajeunesse, M. J. in Handbook of Meta-analysis in Ecology and Evolution (eds Koricheva, J. et al.) Ch. 13 (Princeton Univ. Press, 2013); https://doi.org/10.23943/princeton/9780691137285.003.0013

91.
Rubin, D. B. & Schenker, N. Multiple imputation in health-care databases: an overview and some applications. Stat. Med. 10, 585–598 (1991).
CAS  PubMed  Article  PubMed Central  Google Scholar 

92.
Sweeting, M. J., Sutton, A. J. & Lambert, P. C. What to add to nothing? Use and avoidance of continuity corrections in meta-analysis of sparse data. Stat. Med. 23, 1351–1375 (2004).
PubMed  Article  PubMed Central  Google Scholar 

93.
Viechtbauer, W. Conducting meta-analyses in R with the metafor. J. Stat. Softw. 36, 1–48 (2010).
Article  Google Scholar 

94.
Hothorn, T. et al. multcomp: Simultaneous inference in general parametric models. R package version 1.4-15 (2016).

95.
Cochran, W. G. The combination of estimates from different experiments. Biometrics 10, 101–129 (1954).
Article  Google Scholar 

96.
Hoaglin, D. C. Misunderstandings about Q and ‘Cochran’s Q test’ in meta-analysis. Stat. Med. 35, 485–495 (2016).
PubMed  Article  PubMed Central  Google Scholar 

97.
McFadden, D. Quantitative Methods for Analyzing Travel Behaviour of Individuals: Some Recent Developments Discussion Paper No. 474 (Cowles Foundation, 1977).

98.
Geary, R. C. The frequency distribution of the quotient of two normal variates. J. R. Stat. Soc. 93, 442–446 (1930).
Article  Google Scholar 

99.
Lajeunesse, M. J. Bias and correction for the log response ratio in ecological meta-analysis. Ecol. Soc. Am. 96, 2056–2063 (2015).
Google Scholar 

100.
Rosenthal, R. The file drawer problem and tolerance for null results. Psychol. Bull. 86, 638–641 (1979).
Article  Google Scholar 

101.
Rosenberg, M. S. The file-drawer problem revisited: a general weighted method for calculating fail-safe number in meta-analysis. Evolution 59, 464–468 (2005).
PubMed  Article  PubMed Central  Google Scholar  More

Nonnutritive sweeteners promote conjugative transfer
To test the effect of nonnutritive sweeteners on the conjugative transfer of ARGs, both intra- and intergenus-transfer experiments (model I) were first conducted, in which the bacteria were exposed to various concentrations of four commonly used nonnutritive sweeteners (SAC, SUC, ASP, and ACE-K) for 8 h at room temperature. Notably, in both mating systems, the whole concentration range (from 0.03 to 300 mg/L) of three sweeteners (SUC, ASP, and ACE-K) caused a significant concentration-dependent increase (p = 0.00017 ~ 0.047, Fig. S1a, b); Pearson correlation analysis was shown in Table S3 in conjugative transfer compared to the control (Fig. 1a, b). The intragenus (donor Escherichia coli K-12 LE392 and recipient E. coli K-12 MG1655) spontaneous conjugative transfer frequency was (1.9 ± 0.2) × 10−4 transconjugants per recipient cell (Fig. S2). However, the conjugative transfer frequencies were significantly enhanced by the sweeteners SUC, ASP, and ACE-K at 0.3 mg/L or above. For example, SUC, ASP, and ACE-K at 300 mg/L enhanced the conjugative frequencies by 1.5- (p = 0.00027), 4.1- (p = 0.000000089), and 3.4-fold (p = 0.0000020), respectively (Fig. 1a). In contrast, SAC did not significantly change the conjugative transfer frequency in the conjugation system (p = 0.200 ~ 0.670, Fig. 1a). For intergenus conjugation (donor E. coli K-12 LE392 and recipient Pseudomonas alloputida), all sweeteners at concentrations of 3 mg/L or higher (except for SAC) were seen to promote the conjugative transfer of the donor RP4 plasmid to recipients of different genera (p = 0.000047 ~0.042, Fig. 1b). SUC, ASP, and ACE-K at 300 mg/L caused a great increase in conjugative transfer, by 2.6- (p = 0.0000020), 4.1- (p = 0.000036), and 4.2-fold (p = 0.000019), respectively (Fig. 1b). It should be noted that the enhanced transfer frequencies were associated with the increased number of colonies on selective transconjugant plates, rather than decreased recipient numbers (Fig. S3).
Fig. 1: Nonnutritive sweeteners (SAC, SUC, ASP, and ACE-K) promoted RP4 plasmid-mediated conjugative transfer.

a Fold changes in conjugative ARG transfer within genera. At high concentrations ( >0.3 mg/L), all tested sweeteners (except for SAC) promoted conjugation (N = 6; ANOVA, p  More