Ecology

Banks, N. C., Paini, D. R., Bayliss, K. L. & Hodda, M. The role of global trade and transport network topology in the human-mediated dispersal of alien species. Ecol. Lett. 18, 188–199 (2015).
Google Scholar 
Epanchin-Niell, R. et al. Controlling invasive species in complex social landscapes. Front. Ecol. Environ. 8, 210–216 (2009).
Google Scholar 
Charles, H. & Dukes, J. S. in Biological Invasions (ed. Nentwig, W.) 217–237 (Springer, 2007). https://doi.org/10.1007/978-3-540-36920-2_13Gallardo, B., Clavero, M., Sánchez, M. & Vilà, M. Global ecological impacts of invasive species in aquatic ecosystems. Glob. Change Biol. 22, 151–163 (2016).
Google Scholar 
Diagne, C. et al. High and rising economic costs of biological invasions worldwide. Nature 592, 571–576 (2021).CAS 

Google Scholar 
Sardain, A., Sardain, E. & Leung, B. Global forecasts of shipping traffic and biological invasions to 2050. Nat. Sustain. 2, 274–282 (2019).
Google Scholar 
Epanchin-Niell, R. S. & Hastings, A. Controlling established invaders: integrating economics and spread dynamics to determine optimal management. Ecol. Lett. 13, 528–541 (2010).
Google Scholar 
Chades, I. et al. General rules for managing and surveying networks of pests, diseases, and endangered species. Proc. Natl. Acad. Sci. USA 108, 8323–8328 (2011).CAS 

Google Scholar 
Epanchin-Niell, R. S. & Wilen, J. E. Optimal spatial control of biological invasions. J. Environ. Econ. Manag. 63, 260–270 (2012).
Google Scholar 
Epanchin-Niell, R. S. & Wilen, J. E. Individual and cooperative management of invasive species in human-mediated landscapes. Am. J. Agric. Econ. 97, 180–198 (2015).
Google Scholar 
Aadland, D., Sims, C. & Finnoff, D. Spatial dynamics of optimal management in bioeconomic systems. Comput. Econ. 45, 545–577 (2015).
Google Scholar 
Baker, C. M. Target the source: optimal spatiotemporal resource allocation for invasive species control. Conserv. Lett. 10, 41–48 (2017).
Google Scholar 
Bushaj, S., Büyüktahtakın, İ. E., Yemshanov, D. & Haight, R. G. Optimizing surveillance and management of emerald ash borer in urban environments. Nat. Res. Model. 34, e12267 (2021).
Google Scholar 
Fischer, S. M., Beck, M., Herborg, L.-M. & Lewis, M. A. Managing aquatic invasions: optimal locations and operating times for watercraft inspection stations. J. Environ. Manag. 283, 111923 (2021).
Google Scholar 
Büyüktahtakın, İ. E. & Haight, R. G. A review of operations research models in invasive species management: state of the art, challenges, and future directions. Ann. Oper. Res. 271, 357–403 (2018).
Google Scholar 
Epanchin-Niell, R. S. Economics of invasive species policy and management. Biol. Invasions 19, 3333–3354 (2017).
Google Scholar 
Bodin, Ö. et al. Improving network approaches to the study of complex social–ecological interdependencies. Nat. Sustain. 2, 551–559 (2019).CAS 

Google Scholar 
Nowzari, C., Precaido, V. M. & Pappas, G. J. Analysis and control of epidemics: a survey of spreading processes on complex networks. IEEE Control Syst. 36, 26–46 (2016).
Google Scholar 
Newman, M. E. J. Spread of epidemic disease on networks. Phys. Rev. E 66, 016128 (2002).CAS 

Google Scholar 
Kempe, D., Kleinberg, J. & Tardos, E. Maximizing the spread of influence through a social network. In Proc. 9th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining 137–146 (ACM Press, 2003).Pastor-Satorras, R. & Vespignani, A. Immunization of complex networks. Phys. Rev. E 65, 036104 (2002).
Google Scholar 
Pastor-Satorras, R., Castellano, C., Van Mieghem, P. & Vespignani, A. Epidemic processes in complex networks. Rev. Mod. Phys. 87, 925–979 (2015).
Google Scholar 
Holme, P., Kim, B. J., Yoon, C. N. & Han, S. K. Attack vulnerability of complex networks. Phys. Rev. E 65, 056109 (2002).
Google Scholar 
Muirhead, J. R. & Macisaac, H. J. Development of inland lakes as hubs in an invasion network. J. Appl. Ecol. 42, 80–90 (2005).
Google Scholar 
de la Fuente, B., Saura, S. & Beck, P. S. Predicting the spread of an invasive tree pest: the pine wood nematode in southern europe. J. Appl. Ecol. 55, 2374–2385 (2018).
Google Scholar 
Minor, E. S. & Urban, D. L. A graph-theory framework for evaluating landscape connectivity and conservation planning. Conserv. Biol. 22, 297–307 (2008).
Google Scholar 
Morel-Journel, T., Assa, C. R., Mailleret, L. & Vercken, E. Its all about connections: hubs and invasion in habitat networks. Ecol. Lett. 22, 313–321 (2019).
Google Scholar 
Perry, G. L. W., Moloney, K. A. & Etherington, T. R. Using network connectivity to prioritise sites for the control of invasive species. J. Appl. Ecol. 54, 1238–1250 (2017).
Google Scholar 
Kvistad, J. T., Chadderton, W. L. & Bossenbroek, J. M. Network centrality as a potential method for prioritizing ports for aquatic invasive species surveillance and response in the Laurentian Great Lakes. Manag. Biol. Invasions 10, 403 (2019).
Google Scholar 
Haight, R. G., Kinsley, A. C., Kao, S.-Y., Yemshanov, D. & Phelps, N. B. Optimizing the location of watercraft inspection stations to slow the spread of aquatic invasive species. Biol. Invasions 23, 3907–3919 (2021).
Google Scholar 
McEachran, M. C. et al. Stable isotopes indicate that zebra mussels (Dreissena polymorpha) increase dependence of lake food webs on littoral energy sources. Freshw, Biol. 64, 183–196 (2019).CAS 

Google Scholar 
Karatayev, A. Y., Burlakova, L. E. & Padilla, D. K. in Invasive Aquatic Species of Europe. Distribution, Impacts and Management (eds Leppäkoski, E. et al.) 433–446 (Springer, 2002).Prescott, T. H., Claudi, R. & Prescott, K. L. Impact of Dreissenid mussels on the infrastructure of dams and hydroelectric power plants. In Quagga and Zebra Mussels (eds Nalepa, T. F. & Schloesser, D. W.) 243–258 (CRC Press, 2013).Invasive Species of Aquatic Plants and Wild Animals in Minnesota: Annual Report for 2020 (Minnesota Department of Natural Resources, 2020).Kanankege, K. S., Alkhamis, M. A., Phelps, N. B. & Perez, A. M. A probability co-kriging model to account for reporting bias and recognize areas at high risk for zebra mussels and eurasian watermilfoil invasions in Minnesota. Front. Vet. Sci. 4, 231 (2018).
Google Scholar 
Mallez, S. & McCartney, M. Dispersal mechanisms for zebra mussels: population genetics supports clustered invasions over spread from hub lakes in Minnesota. Biol. Invasions 20, 2461–2484 (2018).
Google Scholar 
Kao, S.-Y. Z. et al. Network connectivity of Minnesota waterbodies and implications for aquatic invasive species prevention. Biol. Invasions 23, 3231–3242 (2021).
Google Scholar 
Kleinberg, J. M. Authoritative sources in a hyperlinked environment. In Proc. 9th Annual ACM-SIAM Symposium on Discrete Algorithms 668–677 (1998).McDonald-Madden, E. et al. Using food-web theory to conserve ecosystems. Nat. Commun. 7, 10245 (2016).CAS 

Google Scholar 
Bossenbroek, J. M., Kraft, C. E. & Nekola, J. C. Prediction of long-distance dispersal using gravity models: zebra mussel invasion of inland lakes. Ecol. Appl. 11, 1778–1788 (2001).
Google Scholar 
Leung, B., Bossenbroek, J. M. & Lodge, D. M. Boats, pathways, and aquatic biological invasions: estimating dispersal potential with gravity models. Biol. Invasions 8, 241–254 (2006).
Google Scholar 
Beger, M. et al. Integrating regional conservation priorities for multiple objectives into national policy. Nat. Commun. 6, 8208 (2015).Runting, R. K. et al. Larger gains from improved management over sparing–sharing for tropical forests. Nat. Sustain. 2, 53–61 (2019).
Google Scholar 
Kinsley, A. C. et al. AIS Explorer: prioritization for watercraft inspections. A decision-support tool for aquatic invasive species management. J. Environ. Manage. 314, 115037 (2022).
Google Scholar 
Vander Zanden, M. J. & Olden, J. D. A management framework for preventing the secondary spread of aquatic invasive species. Can. J. Fish. Aquat. Sci. 65, 1512–1522 (2008).
Google Scholar 
Kanankege, K. S. et al. Lessons learned from the stakeholder engagement in research: application of spatial analytical tools in one health problems. Front. Vet. Sci. 7, 254 (2020).
Google Scholar 
Kroetz, K. & Sanchirico, J. The bioeconomics of spatial-dynamic systems in natural resource management. Annu. Rev. Resour. Econ. 7, 189–207 (2015).
Google Scholar 
Cade, B. S. & Noon, B. R. A gentle introduction to quantile regression for ecologists. Front. Ecol. Environ. 1, 412–420 (2003).
Google Scholar 
Koenker, R. in Asymptotic Statistics (eds Mandl, P. & Hušková, M.) 349–359 (Springer, 1994).Ashander, J. Analysis code and data for ‘Guiding large-scale management of invasive species using network metrics’. figshare https://doi.org/10.6084/m9.figshare.14402447 (2021). More

Climate and land cover dataOur study of tropical dry season dynamics required climatic variables with high temporal resolution (i.e., daily) and full coverage of tropic regions. To reduce uncertainties associated with the choice of precipitation (P) and evapotranspiration (Ep or E) datasets, we used an ensemble of eight precipitation products, three reanalysis-based products for Ep, and one satellite-based land E product. These precipitation datasets were derived four gauge-based or satellite observation (CHIRPS58, GPCC59, CPC-U60 and PERSIANN-CDR61), three reanalyses (ERA-562, MERRA-263, and PGF64) and a multi-source weighted ensemble product (MSWEP v2.865). The potential evapotranspiration (Ep) was calculated using the FAO Penman–Monteith equation66 (Eqs. (1, 2)), which requires meteorological inputs of wind speed, net radiation, air temperature, specific humidity, and surface pressure. We derived these meteorological variables from the three reanalysis products (ERA-5, MERRA-2, and GLDAS-2.067). Since PGF reanalysis lacked upward short- and long-wave radiation output and thus net radiation, we used available meteorological outputs from GLDAS-2.0 instead, which was forced entirely with the PGF input data.$${Ep}=frac{0.408cdot triangle cdot left({R}_{n}-Gright)+gamma cdot frac{900}{T+273}cdot {u}_{2}cdot left({e}_{s}-{e}_{a}right)}{triangle +{{{{{rm{gamma }}}}}}cdot left(1+0.34cdot {u}_{2}right)}$$
(1)
$${VPD}={e}_{s}-{e}_{a}=0.6108cdot {e}^{frac{17.27cdot T}{T+237.3}}cdot left(1-frac{{RH}}{100}right)$$
(2)
Where Ep is the potential evapotranspiration (mm day−1). Rn is net radiation at the surface (MJ m−2 day−1), T is mean daily air temperature at 2 m height (°C), ({u}_{2}) is wind speed at 2 m height (m s−1), ((,{e}_{s}-{e}_{a})) is the vapor pressure deficit of the air (kPa), ({RH}) is the relative air humidity near surface (%), ∆ is the slope of the saturation vapor pressure-temperature relationship (kPa °C−1), γ is the psychrometric constant (kPa °C−1), G is the soil heat flux (MJ m−2 day−1, is often ignored for daily time steps G ≈ 0).We derived the daily evapotranspiration data from the Global Land Evaporation Amsterdam Model (GLEAM v3.3a68), which is a set of algorithms dedicated to developing terrestrial evaporation and root-zone soil moisture data. GLEAM fully assimilated the satellite-based soil moisture estimates from ESA CCI, microwave L-band vegetation optical depth (VOD), reanalysis-based temperature and radiation, and multi-source precipitation forcings. The direct assimilation of observed soil moisture allowed us to detect true soil moisture dynamic and its impacts on evapotranspiration. Besides, the incorporation of VOD, which is closely linked to vegetation water content69,70, allowed us to detect the effect of water stress, heat stress, and vegetation phenological constraints on evaporation. Other observation-driven ET products from remote-sensing physical estimation and flux-tower are not included due to their low temporal resolution (i.e., monthly)71 or short duration72,73. ET outputs of reanalysis products are not considered in our analysis, because the assimilation systems lack explicit representation of inter-annual variability of vegetation activities and thus may not fully capture hydrological response to vegetation changes62,63,67.We used land cover maps for the year 2001 from the Moderate-Resolution Imaging Spectroradiometer (MODIS, MCD12C1 C574) based on the IGBP classification scheme to exclude water-dominated and sparely-vegetated pixels (like Sahara, Arabian Peninsula). All climate and land cover datasets mentioned above were remapped to a common 0.25° × 0.25° grid and unified to daily resolution. The main characteristics of the datasets mentioned above are summarized in Supplementary Table 1.Outputs of CMIP6 simulationsTo understand how modeled dry season changes compare with observed changes, we analyzed outputs from the “historical” (1983-2014) runs of 34 coupled models participating in the 6th Coupled Model Inter-comparison Project75 (CMIP6, Supplementary Table 3). We used these models because they offered daily outputs of all climatic variables needed for our analysis, including precipitation, latent heat (convert to E), and multiple meteorological variables for Ep (air temperature, surface specific humidity, wind speed, and net radiation). All outputs were remapped to a common 1.0° × 1.0° grid and unified to daily resolution.Defining dry season length and timingFor each grid cell and each dry season definition (P  More

The management of introduced species, whether kudzu or zebra mussels, is costly and complex. Now, a paper reports a workable, effective solution that harnesses network analyses of ecological phenomena.Invasive species can pose severe economic and environmental problems, costing more than US$1 trillion worldwide since 1970 (ref. 1). Yet managing this human-driven issue is difficult in itself. The regions involved can be vast — entire continents or countries, for instance — while budgets are typically limited. As well, the sites potentially affected and management options can be numerous. Real systems (for example, all the lakes in the United States) can have thousands of locations that could potentially be infested. By contrast, considering just 40 locations means dealing theoretically with over 1 trillion unique combinations (240) where management could be applied (for instance, to reduce the number of invasive species leaving infested areas or entering uninfested ones). Given these constraints, a key problem is how and where to deploy control measures such as invasive-species removal. While sophisticated optimization approaches exist2, which use mathematical rules to exclude most suboptimal combinations and quickly zoom in to which locations should be managed to minimize new invasions, these algorithms are generally unfeasible for very large systems. Now, writing in Nature Sustainability, Ashander et al.3 demonstrate that simpler network metrics revealing linkages between patches can provide solutions that are often comparable to the more complex optimization algorithms. More

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain
the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in
Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles
and JavaScript. More

Gause, G. F. Experimental analysis of Vito Volterra’s mathematical theory of the struggle for existence. Science 79, 16–17 (1934).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Amarasekare, P. Competitive coexistence in spatially structured environments: A synthesis. Ecol. Lett. 6, 1109–1122 (2003).Article 

Google Scholar 
HilleRisLambers, J., Adler, P. B., Harpole, W. S., Levine, J. M. & Mayfield, M. M. Rethinking community assembly through the lens of coexistence theory. Annu. Rev. Ecol. Evol. Syst. 43, 227–248 (2012).Article 

Google Scholar 
Wisz, M. S. et al. The role of biotic interactions in shaping distributions and realised assemblages of species: Implications for species distribution modelling. Biol. Rev. 88, 15–30 (2013).PubMed 
Article 

Google Scholar 
Frey, S., Fisher, J. T., Burton, A. C. & Volpe, J. P. Investigating animal activity patterns and temporal niche partitioning using camera-trap data: Challenges and opportunities. Remote Sens. Ecol. Conserv. 3, 123–132 (2017).Article 

Google Scholar 
Davis, C. L. et al. Ecological correlates of the spatial co-occurrence of sympatric mammalian carnivores worldwide. Ecol. Lett. 21, 1401–1412 (2018).PubMed 
Article 

Google Scholar 
Durant, S. M. Competition refuges and coexistence: An example from Serengeti carnivores. J. Anim. Ecol. 67, 370–386 (1998).Article 

Google Scholar 
Fedriani, J. M., Fuller, T. K., Sauvajot, R. M. & York, E. C. Competition and intraguild predation among three sympatric carnivores. Oecologia 125, 258–270 (2000).ADS 
PubMed 
Article 

Google Scholar 
Kamler, J. F., Ballard, W. B., Gilliland, R. L. & Mote, K. Spatial relationships between swift foxes and coyotes in northwestern Texas. Can. J. Zool. 81, 168–172 (2003).Article 

Google Scholar 
Vanak, A. T. et al. Moving to stay in place: Behavioral mechanisms for coexistence of African large carnivores. Ecology 94, 2619–2631 (2013).PubMed 
Article 

Google Scholar 
Donadio, E. & Buskirk, S. W. Diet, morphology, and interspecific killing in carnivora. Am. Nat. 167, 524–536 (2006).PubMed 
Article 

Google Scholar 
Tsunoda, H. et al. Food niche segregation between sympatric golden jackals and red foxes in central Bulgaria. J. Zool. 303, 64–71 (2017).Article 

Google Scholar 
Palomares, F. & Caro, T. M. Interspecific killing among mammalian carnivores. Am. Nat. 153, 492–508 (1999).CAS 
PubMed 
Article 

Google Scholar 
Linnell, J. D. C. & Strand, O. Interference interactions, co-existence and conservation of mammalian carnivores. Divers. Distrib. 6, 169–176 (2000).Article 

Google Scholar 
Kamler, J. F., Stenkewitz, U., Klare, U., Jacobsen, N. F. & MacDonald, D. W. Resource partitioning among cape foxes, bat-eared foxes, and black-backed jackals in South Africa. J. Wildl. Manag. 76, 1241–1253 (2012).Article 

Google Scholar 
Di Bitetti, M. S., Di Blanco, Y. E., Pereira, J. A., Paviolo, A. & Pírez, I. J. Time Partitioning favors the coexistence of sympatric crab-eating foxes (Cerdocyon thous) and Pampas Foxes (Lycalopex gymnocercus). J. Mammal. 90, 479–490 (2009).Article 

Google Scholar 
Lesmeister, D. B., Nielsen, C. K., Schauber, E. M. & Hellgren, E. C. Spatial and temporal structure of a mesocarnivore guild in Midwestern North America. Wildl. Monogr. 191, 1–61 (2015).Article 

Google Scholar 
Di Bitetti, M. S., De Angelo, C. D., Di Blanco, Y. E. & Paviolo, A. Niche partitioning and species coexistence in a Neotropical felid assemblage. Acta Oecologica 36, 403–412 (2010).ADS 
Article 

Google Scholar 
Monterroso, P., Alves, P. C. & Ferreras, P. Plasticity in circadian activity patterns of mesocarnivores in southwestern Europe: Implications for species coexistence. Behav. Ecol. Sociobiol. 68, 1403–1417 (2014).Article 

Google Scholar 
Tsunoda, H., Ito, K., Peeva, S., Raichev, E. & Kaneko, Y. Spatial and temporal separation between the golden jackal and three sympatric carnivores in a human-modified landscape in central Bulgaria. Zool. Ecol. 28, 172–179 (2018).Article 

Google Scholar 
Tsunoda, H. et al. Spatio-temporal partitioning facilitates mesocarnivore sympatry in the Stara Planina Mountains, Bulgaria. Zoology 141, 125801 (2020).PubMed 
Article 

Google Scholar 
Ramesh, T., Kalle, R., Sankar, K. & Qureshi, Q. Spatio-temporal partitioning among large carnivores in relation to major prey species in Western Ghats. J. Zool. 287, 269–275 (2012).Article 

Google Scholar 
Gómez-Ortiz, Y., Monroy-Vilchis, O. & Castro-Arellano, I. Temporal coexistence in a carnivore assemblage from central Mexico: Temporal-domain dependence. Mammal Res. 64, 333–342 (2019).Article 

Google Scholar 
Ridout, M. S. & Linkie, M. Estimating overlap of daily activity patterns from camera trap data. J. Agric. Biol. Environ. Stat. 14, 322–337 (2009).MathSciNet 
MATH 
Article 

Google Scholar 
Meredith, M. & Ridout, M. Overlap: Estimates of coefficient of overlapping for animal activity patterns. https://cran.r-project.org/web/packages/overlaphttps://cran.r-project.org/web/packages/overlap/index.html (2018).Marinho, P. H., Fonseca, C. R., Sarmento, P., Fonseca, C. & Venticinque, E. M. Temporal niche overlap among mesocarnivores in a Caatinga dry forest. Eur. J. Wildl. Res. 66, 1–13 (2020).Article 

Google Scholar 
Vilella, M., Ferrandiz-Rovira, M. & Sayol, F. Coexistence of predators in time: Effects of season and prey availability on species activity within a Mediterranean carnivore guild. Ecol. Evol. 10, 11408–11422 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Zhao, G. et al. Spatio-temporal coexistence of sympatric mesocarnivores with a single apex carnivore in a fine-scale landscape. Glob. Ecol. Conserv. 21, e00897 (2020).Article 

Google Scholar 
Farmer, M. J., Allen, M. L., Olson, E. R., Van Stappen, J. & Van Deelen, T. R. Agonistic interactions and island biogeography as drivers of carnivore spatial and temporal activity at multiple scales. Can. J. Zool. 99, 309–317 (2021).Article 

Google Scholar 
Watabe, R. & Saito, M. U. Diel activity patterns of three sympatric medium-sized carnivores during winter and spring in a heavy snowfall area in northeastern Japan. Mammal Study 46, 69–75 (2021).Article 

Google Scholar 
Lashley, M. A. et al. Estimating wildlife activity curves: comparison of methods and sample size. Sci. Rep. 8, 4173 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Niedballa, J., Wilting, A., Sollmann, R., Hofer, H. & Courtiol, A. Assessing analytical methods for detecting spatiotemporal interactions between species from camera trapping data. Remote Sens. Ecol. Conserv. 5, 272–285 (2019).Article 

Google Scholar 
Karanth, K. U. et al. Spatio-temporal interactions facilitate large carnivore sympatry across a resource gradient. Proc. R. Soc. B Biol. Sci. 284, 20161860 (2017).Article 

Google Scholar 
Cusack, J. J. et al. Revealing kleptoparasitic and predatory tendencies in an African mammal community using camera traps: A comparison of spatiotemporal approaches. Oikos 126, 812–822 (2017).Article 

Google Scholar 
Balme, G. et al. Big cats at large: density, structure, and spatio-temporal patterns of a leopard population free of anthropogenic mortality. Popul. Ecol. 61, 256–267 (2019).Article 

Google Scholar 
Li, Z. et al. Coexistence of two sympatric flagship carnivores in the human-dominated forest landscapes of Northeast Asia. Landsc. Ecol. 34, 291–305 (2019).Article 

Google Scholar 
Lahkar, D., Ahmed, M. F., Begum, R. H., Das, S. K. & Harihar, A. Inferring patterns of sympatry among large carnivores in Manas National Park: A prey-rich habitat influenced by anthropogenic disturbances. Anim. Conserv. 24, 589–601 (2021).Article 

Google Scholar 
Paúl, M. J., Layna, J. F., Monterroso, P. & Álvares, F. Resource partitioning of sympatric African Wolves (Canis lupaster) and side-striped jackals (Canis adustus) in an arid environment from West Africa. Diversity 12, 477 (2020).Article 

Google Scholar 
Prat-Guitart, M., Onorato, D. P., Hines, J. E. & Oli, M. K. Spatiotemporal pattern of interactions between an apex predator and sympatric species. J. Mammal. 101, 1279–1288 (2020).Article 

Google Scholar 
Stone, L. & Roberts, A. The checkerboard score and species distributions. Oecologia 85, 74–79 (1990).ADS 
PubMed 
Article 

Google Scholar 
Griffith, D. M., Veech, J. A. & Marsh, C. J. Cooccur: Probabilistic species co-occurrence analysis in r. J. Stat. Softw. 69, 1–17 (2016).Article 

Google Scholar 
Noor, A., Mir, Z. R., Veeraswami, G. G. & Habib, B. Activity patterns and spatial co-occurrence of sympatric mammals in the moist temperate forest of the Kashmir Himalaya, India. Folia Zool. 66, 231–241 (2017).Article 

Google Scholar 
de Satgé, J., Teichman, K. & Cristescu, B. Competition and coexistence in a small carnivore guild. Oecologia 184, 873–884 (2017).ADS 
PubMed 
Article 

Google Scholar 
Kass, J. M., Tingley, M. W., Tetsuya, T. & Koike, F. Co-occurrence of invasive and native carnivorans affects occupancy patterns across environmental gradients. Biol. Invasions 22, 2251–2266 (2020).Article 

Google Scholar 
Louppe, V., Herrel, A., Pisanu, B., Grouard, S. & Veron, G. Assessing occupancy and activity of two invasive carnivores in two Caribbean islands: implications for insular ecosystems. J. Zool. 313, 182–194 (2020).Article 

Google Scholar 
Proulx, G. et al. World distribution and status of the genus Martes in 20. In Martens and Fishers (Martes) in Human-Altered Environments (eds Harrison, D. J. et al.) 21–76 (Springer, Berlin, 2005). https://doi.org/10.1007/b99487.Chapter 

Google Scholar 
Ohdachi, S. D., Ishibashi, Y., Iwasa, M., Fukuki, D. & Saitoh, T. The Wild Mammals of Japan 2nd edn. (Shokadoh Book Seller, Kyoto, 2015).
Google Scholar 
Kauhala, K. & Saeki, M. Nyctereutes procyonoides. The IUCN Red List of Threatened Species. https://www.iucnredlist.org/species/14925/85658776 (2016).Yamamoto, Y. Comparative analyses on food habits of Japanese marten, red fox, badger and raccoon dog in the Mt. Nyugasa, Nagano Prefecture, Japan. Nat. Environ. Sci. Res. 7, 45–52 (1994) (in Japanese with English summary).
Google Scholar 
Hisano, M. et al. A comparison of visual and genetic techniques for identifying Japanese marten scats enabling diet examination in relation to seasonal food availability in a sub-alpine area of Japan. Zool. Sci. 34, 137–146 (2017).Article 

Google Scholar 
Lindstrom, E. R., Brainerd, S. M., Helldin, J. O. & Overskaug, K. Pine marten-red fox interactions: A case of intraguild predation?. Ann. Zool. Fenn. 32, 123–130 (1995).
Google Scholar 
Waggershauser, C. N., Ruffino, L., Kortland, K. & Lambin, X. Lethal interactions among forest-grouse predators are numerous, motivated by hunger and carcasses, and their impacts determined by the demographic value of the victims. Ecol. Evol. 11, 7164–7186 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Watabe, R., Saito, M. U., Enari, H. S. & Enari, H. Mammalian fauna of the Kaminagawa Experimental Forest of Yamagata University detected by camera traps. Tohoku J. For. Sci. 25, 37–40 (2020) (in Japanese).
Google Scholar 
Hofmeester, T. R., Rowcliffe, J. M. & Jansen, P. A. A simple method for estimating the effective detection distance of camera traps. Remote Sens. Ecol. Conserv. 3, 81–89 (2017).Article 

Google Scholar 
Di Bitetti, M. S., Paviolo, A. & De Angelo, C. Camera trap photographic rates on roads vs. off roads: Location does matter. Mastozoología Neotrop. 21, 37–46 (2014).
Google Scholar 
Borcard, D. & Legendre, P. Is the Mantel correlogram powerful enough to be useful in ecological analysis? A simulation study. Ecology 93, 1473–1481 (2012).PubMed 
Article 

Google Scholar 
Oksanen, J. et al. Vegan: community ecology package. https://cran.r-project.org/web/packages/veganhttps://cran.r-project.org/web/packages/vegan/index.html (2019).R Core Team. R: a language environment for statistical computing. r foundation for statistical computing, Vienna, Austria. https://www.r-project.org/https://www.r-project.org/ (2021).Linkie, M. & Ridout, M. S. Assessing tiger-prey interactions in Sumatran rainforests. J. Zool. 284, 224–229 (2011).Article 

Google Scholar 
Watabe, R. & Saito, M. U. Effects of vehicle-passing frequency on forest roads on the activity patterns of carnivores. Landsc. Ecol. Eng. 17, 225–231 (2021).Article 

Google Scholar 
Furukawa, G. genkiFurukawa/rSetDayNightAttr documentation. https://rdrr.io/github/genkiFurukawa/rSetDayNightAhttps://rdrr.io/github/genkiFurukawa/rSetDayNightAttr/ (2019).Mielke, P. W., Berry, K. J. & Johnson, E. S. Multi-response permutation procedures for a priori classifications. Commun. Stat. Theory Methods 5, 1409–1424 (1976).MATH 
Article 

Google Scholar 
Kronfeld-Schor, N. & Dayan, T. Partitioning of time as an ecological resource. Annu. Rev. Ecol. Evol. Syst. 34, 153–181 (2003).Article 

Google Scholar 
Monterroso, P., Alves, P. C. & Ferreras, P. Catch me if you can: Diel activity patterns of mammalian prey and predators. Ethology 119, 1044–1056 (2013).Article 

Google Scholar 
Hendrichsen, D. K. & Tyler, N. J. C. How the timing of weather events influences early development in a large mammal. Ecology 95, 1737–1745 (2014).CAS 
PubMed 
Article 

Google Scholar 
Herfindal, I. et al. Weather affects temporal niche partitioning between moose and livestock. Wildlife Biol. https://doi.org/10.2981/wlb.00275 (2017).Article 

Google Scholar 
Haswell, P. M., Jones, K. A., Kusak, J. & Hayward, M. W. Fear, foraging and olfaction: How mesopredators avoid costly interactions with apex predators. Oecologia 187, 573–583 (2018).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Barrull, J. et al. Factors and mechanisms that explain coexistence in a Mediterranean carnivore assemblage: An integrated study based on camera trapping and diet. Mamm. Biol. 79, 123–131 (2014).Article 

Google Scholar 
Tattersall, E. R., Burgar, J. M., Fisher, J. T. & Burton, A. C. Boreal predator co-occurrences reveal shared use of seismic lines in a working landscape. Ecol. Evol. 10, 1678–1691 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Moll, R. J. et al. Humans and urban development mediate the sympatry of competing carnivores. Urban Ecosyst. 21, 765–778 (2018).Article 

Google Scholar 
McCreadie, J. W. & Bedwell, C. R. Patterns of co-occurrence of stream insects and an examination of a causal mechanism: Ecological checkerboard or habitat checkerboard?. Insect Conserv. Divers. 6, 105–113 (2013).Article 

Google Scholar  More

Palaeogeographical modelWe use palaeogeographical reconstructions describing Earth’s palaeotopography and palaeobathymetry for a series of time slices from 541 Ma to the present day. The reconstructions merge existing models from two published global reconstruction datasets—those of ref. 32 and ref. 33 (https://doi.org/10.5281/zenodo.5348492), which themselves are syntheses of a wealth of previous work.For continental regions, estimates of palaeoelevation and continental flooding rely on a diverse range of geological evidence, such as sedimentary depositional environments and the spatiotemporal distribution of volcanic activity. For a full description, see a recent review34. Together, these data can be used to define the past locations of mountain ranges and palaeoshorelines34. For this part of our reconstruction, we used the compilation of ref. 33 with updated palaeoshorelines based on depositional environment information in current fossil databases35. This compilation comprises 82 palaeotopography maps covering the entire Phanerozoic. Note that each palaeogeographical map is a time slice representing the concatenation of geological data over several million years36.We quantified the impact of using the original compilation of ref. 33 on our model results and found only small changes with respect to using the reconstructions with updated palaeoshorelines (Extended Data Fig. 3a–c). Similarly, eustatic sea level is thought to have varied by around 100 m at timescales much shorter than the duration of the time-slices throughout the Phanerozoic37, such that the extent of continental flooding could have varied within each time slice by an amount significant for our analysis. For this reason, and to assess the uncertainty of our results to continental palaeogeography in general, we computed additional maps of continental flooding in which the sea level is raised or lowered by 100 m compared with the original palaeo–digital elevation model grids of ref. 33 (Extended Data Fig. 3d–f).For deep-ocean regions, the primary control on seafloor depth is the age of the seafloor, so reconstructing palaeobathymetry relies on constructing maps of seafloor age back in time38. As a consequence, we rely on reconstruction models that incorporate a continuous network of plate boundaries. For this part, we used the reconstruction of ref. 32 and derived maps of seafloor age from the plate tectonic model using the method of ref. 39 for which source code is available at GitHub (https://github.com/siwill22/agegrid-0.1). Palaeobathymetry is derived from the seafloor age maps following the steps outlined in ref. 38. It is important to note that seafloor age maps for most of the Phanerozoic (that is, pre-Pangaea times) are not directly constrained by data due to recycling of oceanic crust at subduction zones. Rather, they are model predictions generated by constructing plate motions and plate boundary configurations from the geological and palaeomagnetic record of the continents. Nonetheless, the first-order trends in ocean-basin volume and mean seafloor age are consistent with independent estimates for at least the last 410 million years (Myr)39.The reconstructions of refs. 32,33 differ in the precise locations of the continents through time. To resolve this discrepancy, we reverse reconstructed the continental palaeoelevation model of ref. 33 to present-day coordinates using their rotation parameters, then reconstructed them back in time using the rotations of ref. 32. Owing to the differences in how the continents are divided into different tectonic units, this process leads to some gaps and overlaps in the results40, which we resolved primarily through a combination of data interpolation and averaging. Manual adjustments were made to ensure that the flooding history remained consistent with the original palaeotopography in areas in which interpolation gives a noticeably different history of seafloor ages. The resulting palaeotopography maps are therefore defined in palaeomagnetic reference frame32 appropriate for use in Earth system models.For the biodiversity modelling, we generate estimates of the age of the seafloor for discrete points within the oceans and flooded continents, and track these ages through the lifetime of each point (Supplementary Video 5). For the oceans, this is achieved using the method described in ref. 39 in which the seafloor is represented by points that are incrementally generated at the mid-ocean ridges for a series of time steps 1 Myr apart, with each point tracked through subsequent time steps based on Euler poles of rotation until either the present-day is reached, or they arrive at a subduction zone and are considered to be destroyed.For the continents, tracking the location of discrete points is generally simpler as most crust is conserved throughout the timespan of the reconstruction. In contrast to the deep oceans (where we assume that crust is at all times submerged), we model the ‘age’ of the seafloor from the history of continental flooding and emergence within the palaeogeographical interpretation33. The continents are seeded with uniformly distributed points at the oldest timeslice (541 Ma) at which they are assigned an age of zero. These points are tracked to subsequent time slices of which the palaeogeography is used to determine whether the point lies within a flooded or emergent region. Points within flooded regions of continents are considered to be seafloor, and the age of this seafloor is accumulated across consecutive time slices where a given point lies within a flooded region. When a point is within an emergent region, the seafloor age is reset to zero. Following this approach, individual points within stable continents may undergo several cycles of seafloor age increasing from zero before being reset. At the continental margins formed during the Pangaea breakup, the age of the seafloor continuously grows from the onset of rifting. Intraoceanic island arcs represent an additional case, which can appear as new tectonic units with the reconstructions at various times. In these cases, we assume that the seafloor has a zero-age at the time at which the intraoceanic arc first develops, then remains predominantly underwater for the rest of its lifetime.Thus, for each of the 82 palaeogeographical reconstructions, we annotate 0.5° by 0.5° grids as continental, flooded continental shelf or oceanic for later use in model coupling and production of regional diversity maps.Palaeoenvironmental conditions under the cGENIE Earth system modelWe use cGENIE41, an Earth system model of intermediate complexity, to simulate palaeoenvironmental conditions of seawater temperature and organic carbon export production (as a surrogate for food supply) throughout the Phanerozoic (from 541 Ma to the present day).cGENIE is based on a three-dimensional (3D) ocean circulation model coupled to a 2D energy–moisture balance atmospheric component and a sea-ice module. We configured the model on a 36 × 36 (latitude, longitude) equal area grid with 17 unevenly spaced vertical levels in depth, down to a maximum ocean depth of 5,900 m. The cycling of carbon and associated tracers in the ocean is based on a size-structured plankton ecosystem model with a single (phosphate) nutrient42,43, and adopts an Arrhenius-type temperature-dependent scheme for the remineralization of organic matter exported to the ocean interior44.cGENIE provides a spatially resolved representation of ocean physics and biogeochemistry, which is a prerequisite for the present study to be able to reconstruct the spatial patterns of biodiversity in deep time. However, owing to the computational impracticality of generating a single transient simulation of physics (that is, temperature) and biogeochemistry (that is, export production) over the entire Phanerozoic, we therefore generate 30 model equilibria at regular time intervals throughout the Phanerozoic that are subsequently used as inputs for the regional diversification model (see the ‘Model coupling’ section).We used 30 Phanerozoic palaeogeographical reconstructions through time (~20 Myr evenly spaced time intervals) to represent key time periods. For each continental configuration corresponding to a given age in Earth history, we generate idealized 2D (but zonally averaged) wind speed and wind stress, and 1D zonally averaged albedo forcing fields45 required by the cGENIE model using the ‘muffingen’ open-source software (see the ‘Code availability’ section). For each palaeogeographical reconstruction, the climatic forcing (that is, solar irradiance and carbon dioxide concentration) is adapted to match the corresponding geological time interval. The partial pressure of CO2 is taken from the recent update of the GEOCARB model46. Solar luminosity is calculated using the model of stellar physics of ref. 47. We impose modern-day orbital parameters (obliquity, eccentricity and precession). The simulations are initialized with a sea-ice-free ocean, homogeneous oceanic temperature (5 °C) and salinity (34.9‰). As variations in the oceanic concentration of bio-available phosphate remain challenging to reconstruct in the geological past48,49, we impose a present-day mean ocean phosphate concentration (2.159 μmol kg−1) in our baseline simulations. We quantify the impact of this uncertainty on our model results by conducting additional simulations using half and twice the present-day ocean phosphate concentration (Extended Data Fig. 3g–i). For each ocean phosphate scenario (that is, 0.5×, 1× and 2× the present-day value), each of the 30 model simulations is then run for 20,000 years, a duration ensuring that deep-ocean temperature and geochemistry reach equilibrium. For each model simulation, the results of the mean annual values of the last simulated year are used for the analysis. Note that, although cGENIE makes projections of the distribution of dissolved oxygen ([O2]) in the ocean, our diversification model does not currently consider oxygenation to be a limit on diversity. Thus, we assumed a modern atmospheric partial pressure of O2 in all 30 palaeo simulations and did not use the resulting projected [O2] fields.Regional diversification modelWe tested two models of diversification—the logistic model and the exponential model—describing the dynamics of regional diversity over time. In both models, the net diversification rate (ρ), with units of inverse time (Myr−1), varies within a pre-fixed range of values as a function of seawater temperature and food availability. The net diversification rate is then calculated for a given location and time according to the following equation:$$rho ={rho }_{max }-({rho }_{max }-{rho }_{min })(1-{Q}_{{rm{temp}}}{Q}_{{rm{food}}})$$
(1)
where ρmin and ρmax set the lower and upper net diversification rate limits within which ρ is allowed to vary, and Qtemp and Qfood are non-dimensional limitation terms with values between 0 and 1 that define the dependence of ρ on temperature and food, respectively (Extended Data Table 1).The model considers a direct relationship between seawater temperature, food supply and the rate of net diversification on the basis of the theoretical control that temperature and food supply exert on the rates of origination and extinction (Supplementary Fig. 1). Temperature rise is expected to accelerate the biochemical kinetics of metabolism50 and shorten the development times of individuals51, leading to higher rates of mutation and origination. Greater food availability increases population sizes, which increases the rates of mutation and reduces the probability of extinction52. Furthermore, a large body of observations shows the existence of a positive relationship between resource availability (that is, food supply) and the standing stock of species in marine and terrestrial communities53,54. A larger food supply would support a greater number of individuals. A greater diversity of food resources could also lead to a finer partitioning of available resources55.The temperature dependence of ρ is calculated using the following equation:$${Q}_{{rm{temp}}}=frac{{Q}_{10}^{frac{T-Tmin }{10}}}{{Q}_{10}^{frac{Tmax -Tmin }{10}}}$$
(2)
where the Q10 coefficient measures the temperature sensitivity of the origination rate. In equation (2) above, T is the seawater temperature (in °C) at a given location and time, and Tmin and Tmax are the 0.01 percentile and the 0.99 percentile, respectively, of the temperature frequency distribution in each time interval. In the model, the values of Tmin and Tmax used to calculate Qtemp are therefore recomputed for every time interval (~5 Myr) according to the temperature frequency distribution of the corresponding time interval. This enables us to use updated Tmin and Tmax values in each Phanerozoic time interval and to account for the thermal adaptation of organisms to ever changing climate conditions.The food limitation term is parameterized using a Michaelis–Menten formulation as follows:$${Q}_{{rm{food}}}=frac{text{POC flux}}{left({K}_{{rm{food}}},+,text{POC flux}right)}$$
(3)
where POC flux (mol m−2 yr−1) is the particulate organic carbon export flux, which is used as a surrogate for food availability, at a given location and time of the simulated seafloor. The parameter Kfood (mol m−2 yr−1) in equation (3) is the half-saturation constant, that is, the POC flux at which the diversification rate is half its maximum value, provided that other factors were not limiting. These temperature and food supply limitation terms vary in space and time as a result of changes in seawater temperature and particulate organic carbon export rate, respectively, thereby controlling the spatial and temporal variability of ρ (Supplementary Video 6).The net diversification rate becomes negative (1) in the event of mass extinctions or (2) in response to regional-scale processes, such as sea-level fall and/or seafloor deformation along convergent plate boundaries. Mass extinction events are imposed as external perturbations to the diversification model by imputing negative net diversification rates to all active seafloor points (ocean points and flooded continental points) and assuming non-selective extinction. The percentage of diversity loss as well as the starting time and duration of mass extinctions are extracted from three fossil diversity curves of reference20,21,22 (Source Data for Fig. 1). Each of these fossil diversity curves provides different insights into the Phanerozoic history of marine animal diversity based on uncorrected range-through genus richness estimates20,22 and sampling standardized estimates21. Regional-scale processes—such as sea level fall during marine regressions and/or seafloor destruction at plate boundaries, either by subduction or uplift—are simulated by the combined plate tectonic–palaeoelevation model, and constrain the time that seafloor habitats have to accumulate diversity.The model assumes non-selective extinction during mass extinction events (that is, the field of bullets model of extinction; everything is equally likely to die, no matter the age of the clade and regardless of adaptation)56. However, there is much fossil evidence supporting extinction selectivity57,58. It could be argued that higher extinction rates at diversity hotspots would have delayed their subsequent recovery, flattening global diversity trends. This argument is difficult to reconcile with Sepkoski’s genus-level global diversity curve but could be consistent with the standardized diversity curve of ref. 21. Similarly, the model is also not suitable for reproducing the explosive radiations of certain taxonomic groups after mass extinctions, which could explain the offset between the model and fossil observations in the early Mesozoic (Fig. 1).Letting D represent regional diversity (number of genera within a given seafloor point) and t represent time, the logistic model is formalized by the following differential equation:$$frac{partial Dleft(tright)}{partial t}=rho Dleft[1-frac{D}{{K}_{{rm{eff}}}}right]$$
(4)
where D(t) is the number of genera at time t and Keff is the effective carrying capacity or maximum number of genera that a given seafloor point (that is, grid cell area after gridding) can carry at that time, t. In our logistic model, Keff is allowed to vary within a fixed range of values (from Kmin to Kmax) as a positive linear function of the POC flux at a given location and time as follows:$${K}_{{rm{eff}}}={K}_{max }-left({K}_{max }-{K}_{min }right)frac{{text{POC flux}}_{max }-text{POC flux}}{{text{POC flux}}_{max }-{text{POC flux}}_{min }}$$
(5)
where POC fluxmin and POC fluxmax correspond to the 0.01 and 0.99 quantiles of the POC flux range in the whole Phanerozoic dataset.In the logistic model, the net diversification rate decreases as regional diversity approaches its Keff. The exponential model is a particular case of the logistic model when Keff approaches infinity and, therefore, neither the origination rate nor the extinction rate depend on the standing diversities. In this scenario, diversity grows in an unlimited manner over time only truncated by the effect of mass extinctions and/or by the dynamics of the seafloor (creation versus destruction). Thus, the exponential model is as follows:$$frac{partial Dleft(tright)}{partial t}=rho D$$
(6)
where the rate of change of diversity (the time derivative) is proportional to the standing diversity D such that the regional diversity will follow an exponential increase in time at a speed controlled by the temperature- and food-dependent net diversification rate. Even if analytical solutions exist for the steady-state equilibrium of the logistic and exponential functions, we solved the ordinary differential equations (4) and (6) using numerical methods with a time lag of 1 Myr to account for the spatially and temporally varying environmental constraints, seafloor dynamics and mass extinction events.As the analysis of global fossil diversity curves is unable to discern the causes of diversity loss during mass extinctions, our imputation of negative diversification rates could have overestimated the loss of diversity in those cases in which sea level fall, a factor already accounted for by our model, contributed to mass extinction. This effect was particularly recognizable across the Permian–Triassic mass extinction (Extended Data Fig. 6d–f), and supports previous claims that the decline in the global area of the shallow water shelf exacerbated the severity of the end-Permian mass extinction34.Model couplingAs stated above, the coupled plate tectonic–palaeoelevation (palaeogeographical) model corresponds to a tracer-based model (a Lagrangian-based approach) that simulates and tracks the spatiotemporal dynamics of ocean and flooded continental points. The diversification models start at time 541 Ma with all active points having a D0 = 1 (one single genus everywhere) and we let points accumulate diversity heterogeneously with time according to seafloor age distributions (for ocean points) and the time that continents have been underwater (for flooded continental points). The ocean points are created at mid-ocean ridges and disappear primarily at subduction zones. Between their origin and demise, the points move following plate tectonic motions and we trace their positions while accumulating diversity. The flooded continental points begin to accumulate diversity from the moment that they are submerged, starting with a D value equal to the nearest neighbour flooded continental point with D  > 1, thereby simulating a process of coastal recolonization (or immigration). The diversification process remains active while the seafloor points remain underwater, but it is interrupted, and D set to 0, in those continental points that emerge above sea level. Similarly, seafloor points corresponding to ocean domains disappear in subduction zones, and their diversity is lost. We track the geographical position of the ocean and flooded continental points approximately every 5 Myr, from 541 Ma to the present. Each and every one of the tracked points accumulates diversity over time at a different rate, which is modulated by the environmental history (seawater temperature and food availability) of each point, as described in equations (1)–(3). When a point arrives in an environment with a carrying capacity lower than the diversity it has accumulated through time, we reset the diversity of the point to the value of the carrying capacity, thereby simulating local extinction.Seawater temperature (T) and food availability (POC flux) are provided by the cGENIE model, which has a spatial and temporal resolution coarser than the palaeogeographical model. The cGENIE model provides average seawater T and POC flux values in a 36 × 36 equal area grid (grid cell area equivalent to 2° latitude by 10° longitude at the equator) and 30 time slices or snapshots (from 541 Ma to the present: each ~20 Myr time intervals). To have environmental inputs for the 82 time slices of the plate tectonic–palaeoelevation model, we first interpolate the cGENIE original model output data on a 0.5° by 0.5° grid to match the annotated grids provided by the plate tectonic–palaeoelevation model. As the relatively coarse spatial resolution of the cGENIE model prevents rendering the coast–ocean gradients, we assign surface T and POC flux at the base of the euphotic zone to the flooded continental shelf grid cells, and deep ocean T and POC flux at the bottom of the ocean to the ocean grid cells. As there are time slices without input data of seawater T and POC flux, we interpolate/extrapolate seawater T and POC flux values into the 0.5° by 0.5° flooded continental shelf and ocean grids independently. Finally, we interpolate values from these 0.5° by 0.5° flooded continental shelf and ocean grids into the exact point locations in each time frame. Thus, each active point is tracked with its associated time-varying T and POC flux values throughout its lifetime. On average, 6,000 flooded continental points and 44,000 oceanic points were actively accumulating diversity in each time frame. The model cannot simulate the singularities of relatively small enclosed seas for which the spatial resolution of the palaeogeographical and Earth system models is insufficient to capture relevant features (such as palaeobathymetry, seawater temperature) in detail. The method is also likely to underestimate the diversity of epeiric (inland) seas due to the difficulty of simulating immigration, a process that is strongly influenced by the effect of surface ocean currents and is not considered here. However, as stated above, the model considers recolonization of recently submerged areas by marine biota from nearby coastal environments, which partially explains coastal immigration.Estimation of global diversity from regional diversityOur regional diversity maps are generated by separately interpolating ocean point diversity and flooded continental point diversity into the 0.5° by 0.5° annotated grids provided by the palaeogeographical model. We calculate global diversity at each time step from each of the regional diversity maps following a series of steps to integrate diversity along line transects from diversity peaks (maxima) to diversity troughs (minima) (Extended Data Fig. 1). To select the transects, first, we identify on each of the regional diversity maps the geographical position of the diversity peaks. We identify local maxima (that is, grid cells with diversity greater than their neighbour cells), and define the peaks as those local maxima with diversity greater than the 0.75 quantile of diversity values in all local maxima in the map. In the case of grid cells with equal neighbour diversity, the peak is assigned to the grid cell in the middle. We subsequently identify the geographical position of the diversity troughs, which are defined as newly formed ocean grid cells (age = 0 Myr) and, therefore, with diversities equal to one. The troughs are mostly located at mid-ocean ridges.On each of the 82 spatial diversity maps, we trace a line transect from each diversity peak to its closest trough, provided that the transect does not cross land in more than 20% of the grid cells along the linear path (Supplementary Video 7). On average, for each spatial diversity map, we trace 400 (σ = ±75) linear transects. This sampling design gives rise to transects of different lengths, which may bias the estimates of global diversity. To minimize this bias, we cut the tail of the transects to have a length of 555 km equivalent to 5° at the equator. We tested an alternative cut-off threshold, 1,110 km, and the results do not alter the study’s conclusions.We apply Bresenham’s line algorithm59 to detect the grid cells crossed by the transects and annotate their diversity. To integrate regional diversity along the transects, we developed a method to simplify the scenario of peaks and troughs heterogeneously distributed on the 2D diversity maps. The method requires (1) a vector (the transect) of genus richness (αn) at n different locations (grids) arranged in a line (1D) of L grids, and (2) a coefficient of similarity (Vn,n + 1) between each two neighbouring locations, n and n + 1. Vn,n + 1, the coefficient of similarity, follows a decreasing exponential function with distance between locations. The number of shared genera between n and n + 1 is Vn,n + 1 × min(αn; αn + 1). We integrate diversity from peaks to troughs and assume that, along the transect, αn + 1 is lower than αn. We further assume that the genera present in n and n + 2 cannot be absent from n + 1. Using this method, we integrate the transect’s diversity (γi) using the following equation:$${gamma }_{i}={ {mathbf{upalpha}}}_{1}+{sum }_{n=1}^{L-1}left(1-{V}_{n,n+1}right){ {mathbf{upalpha}}}_{n+1}$$
(7)
To integrate the diversity of all transects (γi) on each 2D diversity map (or time slice), we apply the same procedure as described above (Extended Data Fig. 1). We first sort the transects in descending order from the highest to the lowest diversity. We then assume that the number of shared genera between transect i and the rest of the transects with greater diversity {1, 2, …, i − 1} is given by the distance of its peak to the nearest neighbour peak (NN(i)) of those already integrated {1, 2, …, i − 1}. Thus, we perform a zigzag integration of transects’ diversities down gradient, from the greatest to the poorest, weighted by the nearest neighbour distance among the peaks already integrated. As a result, the contribution of each transect to global diversity will depend on its diversity and its distance to the closest transect out of all those transects already integrated. Using this method, we linearize the problem to simplify the cumbersome procedure of passing from a 2D regional diversity map to a global diversity estimate without knowing the identity (taxonomic affiliation) of the genera. If γtotal is the global diversity at time t:$${gamma }_{{rm{total}}}={gamma }_{1}+{sum }_{i=2}^{j}left(1-{V}_{{rm{NN}}left(iright),i}right){gamma }_{i}$$
(8)
Finally, the resulting global estimates are plotted against the midpoint value of the corresponding time interval to generate a synthetic global diversity curve. To compare the global diversity curves produced by the diversification models with those composed from the fossil record, Lin’s CCC60 is applied to the data normalized to the min–max values of each time series (that is, rescaled within the range 0–1). Lin’s CCC combines measures of both precision and accuracy to determine how far the observed data deviate from the line of perfect concordance (that is, the 1:1 line). Lin’s CCC increases in value as a function of the nearness of the data’s reduced major axis to the line of perfect concordance (the accuracy of the data) and of the tightness of the data around its reduced major axis (the precision of the data).The time series of global diversity generated from the fossil record and from the diversification model exhibit serial correlation and the resulting CCCs are therefore inflated. The use of methods for analysing non-zero autocorrelation time series data, such as first differencing or generalised least squares regression, enables high-frequency variations along the time series to be taken into account. However, the relative simplicity of our model, which was designed to reproduce the main Phanerozoic trends in global diversity, coupled with the fact that biases in the fossil data would introduce uncertainty into the analysis, leads us to focus our analysis on the long-term trends, obviating the effect of autocorrelation.Model parameterization and calibrationThe diversification models are parameterized assuming a range of values that constrain the lower and upper limits of the genus-level net diversification rate (ρmin and ρmax, respectively) (Extended Data Table 1) according to previously reported estimates from fossil records (figures 8 and 11 of ref. 5). A range of realistic values is assigned for the parameters Q10 and Kfood, determining, respectively, the thermal sensitivity and food dependence of the net diversification rate. We test a total of 40 different parameter combinations (Extended Data Table 2). The resulting estimates of diversity are then compared against the fossil diversity curves of ref. 20, ref. 21 or ref. 22, and the 15 parameter combinations providing the highest CCCs are selected.The results of the logistic diversification model rely on the values of the minimum and maximum carrying capacities (Kmin and Kmax, respectively) within which the spatially resolved effective carrying capacities (Keff) are allowed to vary. The values of Kmin and Kmax are therefore calibrated by running 28 simulations of pair-wise Kmin and Kmax combinations increasing in a geometric sequence of base 2, from 2 to 256 genera (Extended Data Fig. 4). We perform these simulations independently for each of the 15 parameter settings selected previously (Extended Data Fig. 4 and Extended Data Table 2). Each combination of Kmin and Kmax produces a global diversity curve, which is evaluated as described above using Lin’s CCC.Calculating estimates of global diversity from regional diversity maps in the absence of information on genus-level taxonomic identities requires that we assume a spatial turnover of taxa with geographical distance (Extended Data Fig. 1). Distance-decay curves are routinely fitted by calculating the ecological similarity (for example, the Jaccard similarity index) between each pair of sampling sites, and fitting an exponential decay function to the points on a scatter plot of similarity (y axis) versus distance (x axis). Following this method, we fit an exponential decay function to the distance–decay curves reported in ref. 61, depicting the decrease in the Jaccard similarity index (J) of fossil genera with geographical distance (great circle distance) at different Phanerozoic time intervals:$$J={J}_{{rm{o}}{rm{f}}{rm{f}}}+(,{J}_{max}text{-}{J}_{{rm{o}}{rm{f}}{rm{f}}}){{rm{e}}}^{-lambda times {rm{d}}{rm{i}}{rm{s}}{rm{t}}{rm{a}}{rm{n}}{rm{c}}{rm{e}}}$$
(9)
where Joff = 0.06 (n.d.) is a small offset, Jmax = 1.0 (n.d.) is the maximum value of the genus-based Jaccard similarity index and λ = 0.0024 (km−1) is the distance-decay rate.The Jaccard similarity index (J) between consecutive points n and n + 1 is bounded between 0 and min(αn; αn + 1)/max(αn; αn + 1). A larger value for J would mean that there are more shared genera between the two communities than there are genera within the least diverse community, which is ecologically absurd. However, using a single similarity decay function can lead the computed value of J to be locally larger than min(αn; αn + 1)/max(αn; αn + 1). To prevent this artefact, we use the Simpson similarity index or ‘overlap coefficient’ (V) instead of J. V corresponds to the percentage of shared genera with respect to the least diverse community (min(αn ; αn + 1)). V is bounded between 0 and 1, whatever the ratio of diversities. As the pre-existing estimates of similarity are expressed using J (ref. 61), we perform the conversion from J to V using the algebraic expression V = (1 + R) × J/(1 + J) where R = max(αn; αn + 1)/min(αn; αn + 1) (Supplementary Note 1). In the cases in which J exceeds the min(αn; αn + 1)/max(αn; αn + 1), V becomes >1 and, in those cases, we force V to be More

Brandes, N. & Linial, M. Giant viruses—big surprises. Viruses 11, 404 (2019).CAS 
Article 

Google Scholar 
Jover, L. F., Effler, T. C., Buchan, A., Wilhelm, S. W. & Weitz, J. S. The elemental composition of virus particles: implications for marine biogeochemical cycles. Nat. Rev. Microbiol. 12, 519–528 (2014).CAS 
Article 

Google Scholar 
Madsen, E. L. Microorganisms and their roles in fundamental biogeochemical cycles. Curr. opinion biotechnology 22, 456–464 (2011).CAS 
Article 

Google Scholar 
Gummow, B. Challenges posed by new and re-emerging infectious diseases in livestock production, wildlife and humans. Livest. Sci. 130, 41–46 (2010).CAS 
Article 

Google Scholar 
Becker, D. J., Streicker, D. G. & Altizer, S. Linking anthropogenic resources to wildlife–pathogen dynamics: a review and meta-analysis. Ecol. letters 18, 483–495 (2015).Article 

Google Scholar 
Woolhouse, M. E. & Gowtage-Sequeria, S. Host range and emerging and reemerging pathogens. Emerg. infectious diseases 11, 1842 (2005).Article 

Google Scholar 
Foster, R. et al. Pathogens co-transported with invasive non-native aquatic species: implications for risk analysis and legislation. NeoBiota 69, 79–102 (2021).Article 

Google Scholar 
Brasier, C. The biosecurity threat to the uk and global environment from international trade in plants. Plant Pathol. 57, 792–808 (2008).Article 

Google Scholar 
Ruiz, G. M. et al. Global spread of microorganisms by ships. Nature 408, 49–50 (2000).ADS 
CAS 
Article 

Google Scholar 
Essl, F. et al. Which taxa are alien? criteria, applications, and uncertainties. BioScience 68, 496–509 (2018).Article 

Google Scholar 
Blackburn, T. M., Bellard, C. & Ricciardi, A. Alien versus native species as drivers of recent extinctions. Front. Ecol. Environ. 17, 203–207 (2019).Article 

Google Scholar 
Hawkins, C. L. et al. Framework and guidelines for implementing the proposed iucn environmental impact classification for alien taxa (eicat). Divers. Distributions 21, 1360–1363 (2015).Article 

Google Scholar 
Corrales, X. et al. Advances and challenges in modelling the impacts of invasive alien species on aquatic ecosystems. Biol. Invasions 22, 907–934 (2020).Article 

Google Scholar 
Regulation, E. Regulation (eu) no 1143/2014 of the European Parliament and of the Council of 22 October 2014 on the prevention and management of the introduction and spread of invasive alien species. Off. J. Eur. Union 57, 35–55 (2014).
Google Scholar 
EU. Regulation (eu) 2016/2031 of the European Parliament of the Council of 26 October 2016 on protective measures against pests of plants, amending regulations (eu) 228/2013,(eu) 652/2014 and (eu) 1143/2014 and repealing council directives 69/464/eec, 74/647/eec, 93/85/eec, 98/57/ec, 2000/29/ec, 2006/91/ec and 2007/33/ec. Off. J. 317, 4–104 (2016).
Google Scholar 
Murtaugh, M. P. et al. The science behind one health: at the interface of humans, animals, and the environment. Tech. Rep. (2017).Ogden, N. H. et al. Emerging infectious diseases and biological invasions: a call for a one health collaboration in science and management. Royal Soc. open science 6, 181577 (2019).ADS 
Article 

Google Scholar 
Roy, H. E. et al. Alien pathogens on the horizon: Opportunities for predicting their threat to wildlife. Conserv. Lett. 10, 477–484 (2017).Article 

Google Scholar 
Ikner, L. A., Gerba, C. P. & Bright, K. R. Concentration and recovery of viruses from water: a comprehensive review. Food Environ. Virol. 4, 41–67 (2012).
Google Scholar 
Taylor, M. W. Introduction: A short history of virology. In Viruses and Man: A History of Interactions, 1–22 (Springer, 2014).Thakur, M. P., Van der Putten, W. H., Cobben, M. M., van Kleunen, M. & Geisen, S. Microbial invasions in terrestrial ecosystems. Nat. Rev. Microbiol. 17, 621–631 (2019).CAS 
Article 

Google Scholar 
Desprez-Loustau, M.-L. et al. The fungal dimension of biological invasions. Trends ecology & evolution 22, 472–480 (2007).Article 

Google Scholar 
Rivett, D. W. et al. Elevated success of multispecies bacterial invasions impacts community composition during ecological succession. Ecol. Lett. 21, 516–524 (2018).Article 

Google Scholar 
Dunn, A. M. & Hatcher, M. J. Parasites and biological invasions: parallels, interactions, and control. TRENDS Parasitol. 31, 189–199 (2015).Article 

Google Scholar 
Pyšek, P. et al. Macroecological framework for invasive aliens (mafia): disentangling large-scale context dependence in biological invasions. (2020).Hulme, P. E. et al. Blurring alien introduction pathways risks losing the focus on invasive species policy. Conserv. Lett. 10, 265–266 (2017).Article 

Google Scholar 
Gilroy, J. J., Avery, J. D. & Lockwood, J. L. Seeking international agreement on what it means to be “native”. Conserv. Lett. 10, 238–247 (2017).Article 

Google Scholar 
Webber, B. L. & Scott, J. K. Rapid global change: implications for defining natives and aliens. Glob. Ecol. Biogeogr. 21, 305–311 (2012).Article 

Google Scholar 
CBD Secretariat. Decision VI/23: Alien species that threaten ecosystems, habitats and species. Document UNEP/CBD/COP/6/23 (2002).World Health Organization. A brief guide to emerging infectious diseases and zoonoses. Tech. Rep. https://apps.who.int/iris/handle/10665/204722 (2014).Firrao, G. et al. Candidatus phytoplasma’, a taxon for the wall-less, non-helical prokaryotes that colonize plant phloem and insects. Int. J. Syst. Evol. Microbiol. 54, 1243–1255 (2004).CAS 
Article 

Google Scholar 
CBD. Pathways of introduction of invasive species, their prioritization and management (Secretariat of the Convention on Biological Diversity Montreal, 2014).OIE. Terrestrial Animal Health Code 2021 (OIE, 2021).Magliozzi, C. et al. bacteria and viruses traits and species-related factors. figshare https://doi.org/10.6084/m9.figshare.18550907.v2 (2022).Katsanevakis, S. et al. Implementing the European policies for alien species: networking, science, and partnership in a complex environment. Manag. Biol. Invasions 4, 3–6 (2013).Article 

Google Scholar 
Tsiamis, K. et al. The EASIN Editorial Board: quality assurance, exchange and sharing of alien species information in europe. Manag. Biol. invasions 7, 321–328 (2016).Article 

Google Scholar 
Wieczorek, J. et al. Darwin core: an evolving community-developed biodiversity data standard. PloS one 7, e29715 (2012).ADS 
CAS 
Article 

Google Scholar 
Darwin Core. Darwin Core quick reference guide. https://dwc.tdwg.org/terms/ (2018).R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, https://www.R-project.org/ (2020).Wickham, H. ggplot2: Elegant Graphics for Data Analysis, https://ggplot2.tidyverse.org (Springer-Verlag New York, 2016).Schwarzl, T. ggBubbles: Mini Bubble Plots for Comparison of Discrete Data with ‘ggplot2’ R package version 0.1.4 (2019).Moon, K. R statistics and graphs for medical papers (Hannarae Seoul, 2015).Current, C. Invasive species compendium. Wallingford, UK: CAB Int. Available online: www.cabi.org/isc (accessed on 19 August 2020) (2011).Adams, M. J. & Antoniw, J. F. Dpvweb: An open access internet resource on plant viruses and virus diseases. Outlooks on Pest Manag. 16, 268 (2005).Article 

Google Scholar 
Adams, M. J. & Antoniw, J. F. Dpvweb: a comprehensive database of plant and fungal virus genes and genomes. Nucleic acids research 34, D382–D385 (2006).CAS 
Article 

Google Scholar 
Benson, D. A. et al. Genbank. Nucleic acids research 41, D36–D42 (2012).Article 

Google Scholar  More

Sakamoto, M. Primary production by phytoplankton community in some Japanese lakes and its dependence on lake depth. Archiv für Hydrobilogie. 62, 1–28 (1966).
Google Scholar 
Vollenweider, R. A. Scientific fundamentals of the eutrophication of lakes and flowing waters, with particular reference to nitrogen and phosphorus as factors in eutrophication (Organisation for Economic Co-operation and Development, 1968).
Google Scholar 
Edmondson, W. T. Phosphorus, nitrogen, and algae in Lake Washington after diversion of sewage. Science 169, 690–691 (1970).ADS 
CAS 
Article 

Google Scholar 
Dillon, P. J. & Rigler, F. H. The phosphorus-chlorophyll relationship in lakes. Limnol. Oceanogr. 19, 767–773 (1974).ADS 
CAS 
Article 

Google Scholar 
Jones, J. R. & Bachmann, R. W. Prediction of phosphorus and chlorophyll levels in lakes. J. Water Pollut. Control Feder. 48, 2176–2182 (1976).CAS 

Google Scholar 
Schindler, D. W. Evolution of phosphorus limitation in lakes. Science 195, 260–262 (1977).ADS 
CAS 
Article 

Google Scholar 
Filstrup, C. T. & Downing, J. A. Relationship of chlorophyll to phosphorus and nitrogen in nutrient-rich lakes. Inland Waters. 7, 385–400 (2017).CAS 
Article 

Google Scholar 
Schindler, D. W. Recent advances in the understanding and management of eutrophication. Limnol. Oceanogr. 51, 356–363 (2006).ADS 
Article 

Google Scholar 
Quinlan, R. et al. Relationships of total phosphorus and chlorophyll in lakes worldwide. Limnol. Oceanogr. 66, 392–404 (2020).ADS 
Article 

Google Scholar 
Yuan, L. L. & Jones, J. R. Rethinking phosphorus–chlorophyll relationships in lakes. Limnol. Oceanogr. 65, 1847–1857 (2020).ADS 
CAS 
Article 

Google Scholar 
Carlson, R. E. A trophic state index for lakes. Limnol. Oceanogr. 11, 361–369 (1977).ADS 
Article 

Google Scholar 
Neveux, J. et al. Comparison of chlorophyll and phaeopigment determinations by spectrophotometric, fluorometric, spectrofluorometric and HPLC methods. Mar. Microb. Food Webs 4, 217–238 (1990).
Google Scholar 
Lampert, W. & Sommer, U. Limnoecology (Oxford University, 2007).
Google Scholar 
Kovalevskaya, R. Z., Zhukava, H. A. & Adamovich, B. V. Modification of the method of spectrophotometric determination of chlorophyll a in the suspended matter of water bodies. J. Appl. Spectrosc. 87, 72–78 (2020).ADS 
CAS 
Article 

Google Scholar 
Søndergaard, M., Lauridsen, T. L., Johansson, L. S. & Jeppesen, E. Nitrogen or phosphorus limitation in lakes and its impact on phytoplankton biomass and submerged macrophyte cover. Hydrobiologia 795, 35–48 (2017).Article 

Google Scholar 
Søndergaard, M., Jensen, J. P., Jeppesen, E. & Møller. P. H. Seasonal dynamics in the concentrations and retention of phosphorus in shallow Danish lakes after reduced loading. Aquat. Ecosyst. Health Manag. 5(1), 19–29 (2002).Magumba, D., Atsushi, M., Michiko, T., Akira, K. & Masao, K. Relationships between Chlorophyll-a, phosphorus and nitrogen as fundamentals for controlling phytoplankton biomass in lakes. Environ. Control. Biol. 51(4), 179–185 (2013).CAS 
Article 

Google Scholar 
Smith, V. H. & Shapiro, J. Chlorophyll-phosphorus relations in individual lakes. Their importance to lake restoration strategies. Environ. Sci. Technol. 15(4), 444–451 (1981).Pothoven, S. A. & Vanderploeg, H. A. Seasonal patterns for Secchi depth, chlorophyll a, total phosphorus, and nutrient limitation differ between nearshore and offshore in Lake Michigan. J. Great Lakes Res. 46, 519–527 (2020).CAS 
Article 

Google Scholar 
Søndergaard, M. & Jeppesen, E. Lake Søbygaard, Denmark: phosphorus dynamics during the first 35 years after an external loading reduction. In: Internal Phosphorus Loading: Causes, Case Studies, and Management (ed. Steinman, A.D. & Spears, B. M.) 285–299 (J. Ross, Plantation, 2020).Guildford, S. J. & Hecky, R. E. Total nitrogen, total phosphorus, and nutrient limitation in lakes and oceans: Is there a common relationship?. Limnol. Oceanogr. 45, 1213–1223 (2000).ADS 
CAS 
Article 

Google Scholar 
Jones, J.R. et al. Nutrients, seston, and transparency of Missouri reservoirs and oxbow lakes. An analysis of regional limnology. Lake Reser. Manag. 24, 155–180 (2008).Pikovsky, A., Rosenblum, M. & Kurths, J. Synchronization. A universal concept in nonlinear sciences (Cambridge University, 2001).Book 

Google Scholar 
Kuramoto, Y. Chemical Oscillations, Waves, and Turbulence (Springer, 1984).Book 

Google Scholar 
Sazonov, A. V. et al. An investigation of the phase locking index for measuring of interdependency of cortical source signals recorded in the EEG. Biol. Cybern. 100, 129–146 (2009).Article 

Google Scholar 
Medvinsky, A. B. et al. Temperature as a factor affecting fluctuations and predictability of the abundance of lake bacterioplankton. Ecol. Complex. 32, 90–98 (2017).Article 

Google Scholar 
Zhukova, T. V. & Ostapenya, A. P. Estimation of efficiency of nature protection measures in water catchment area of the Naroch lakes. Natural Resources. 3, 68–73 (2000) ((in Russian)).
Google Scholar 
Burlakova, L. E., Karatayev, A. Y. & Padilla, D. K. Changes in the distribution and abundance of Dreissena polymorpha within lakes through time. Hydrobiologia 571, 133–146 (2006).Article 

Google Scholar 
Ostapenia, A. P. et al. Bentification of lake ecosystem: causes, mechanisms, possible consequences, prospect for future research. Trudy BGU. 7, 135–148 (2012) ((in Russian)).
Google Scholar 
Karatayev, A.Y., Burlakova, L.E. & Padilla, D.K. Impacts of Zebra Mussels on aquatic communities and their role as ecosystem engineers. In: Leppäkoski, E., Gollasch, S., Olenin, S. (eds) Invasive Aquatic Species of Europe. Distribution, Impacts and Management (Springer, Dordrecht, 2002).Adamovich, B. V. et al. The divergence of chlorophyll dynamics in the Naroch Lakes. Biophysics 60, 632–638 (2015).CAS 
Article 

Google Scholar 
Zhukova, T. V. et al. Long-term dynamics of suspended matter in Naroch Lakes: Trend or intervation. Inland Water Biol. 10, 250–257 (2017).Article 

Google Scholar 
Adamovich, B. V. et al. Eutrophication, oligotrophication, and benthiphication in Naroch Lakes: 40 years of monitoring. J. Siber. Federal Univ. Biol. 10, 379–394 (2017).Article 

Google Scholar 
Ostapenya A.P. et al. Ecological passport of Lake Myastro (EcoMir, Minsk, 1994) (in Russian).Kantz, H. & Schreiber, T. Nonlinear time series analysis (Cambridge University, 1997).MATH 

Google Scholar 
Kot, M. Elements of mathematical ecology (Cambridge University, 2001).Book 

Google Scholar 
Turchin, P. Complex population dynamics. A Theoretical/Empirical Synthesis (Princeton University, Princeton, 2003).MATH 

Google Scholar 
Cazelles, B. & Stone, L. Detection of imperfect population synchrony in an uncertain world. J. Anim. Ecol. 72, 953–968 (2003).Article 

Google Scholar 
Karatayev, A. Y., Burlakova, L. & Padilla, D. K. The effects of Dreissena polymorpha (Pallas) invasion on aquatic communities in Eastern Europe. J. Shellfish Res. 16, 187–203 (1997).
Google Scholar 
Lia, J. et al. Benthic invaders control the phosphorus cycle in the world’s largest freshwater ecosystem. PNAS 118(6), e2008223118. https://doi.org/10.1073/pnas.2008223118 (2021).CAS 
Article 

Google Scholar 
Mikheyeva, T. M. et al. The dynamics of freshwater phytoplankton stability in the Naroch Lakes (Belarus). Ecol. Ind. 81, 481–490 (2017).Article 

Google Scholar 
Harris, P. H. Phytoplankton ecology. Structure, functioning and flucttuation (Chapman & Hall, London, New York, 1986).Jeppesen, E., Jensen, J. P., Søndergaard, M. & Lauridsen, T. L. Response of fish and plankton to nutrient loading reduction in eight shallow Danish lakes with special emphasis on seasonal dynamics. Freshw. Biol. 50, 1616–1627 (2005).CAS 
Article 

Google Scholar 
Nezlin, N.P. & Li, B-L. Time-series analysis of remote-sensed chlorophyll and environmental factors in the Santa Monica–San Pedro Basin off Southern California. J. Mar. Syst. 39, 185–202 (2003).French, T. D. & Petticrew, E. L. Chlorophyll a seasonality in four shallow eutrophic lakes (northern British Columbia, Canada) and the critical roles of internal phosphorus loading and temperature. Hydrobiologia 575, 285–299 (2007).CAS 
Article 

Google Scholar 
SCOR-UNESCO Working Group no. 17. Determination of photosynthetic pigments in sea-water. Monographs on Oceanologic Methodology 9–18 (UNESSCO, Paris, 1966).Semenov, A. D. Guide on the chemical analysis of continental surface waters (Gidrometeoizdat, 1977) ((in Russian)).
Google Scholar 
Wetzel, R. G. & Likens, G. E. Limnological analysis (Springer, 2000).Book 

Google Scholar 
Steffen, M. & Bartz-Beielstein, T. imputeTS: time series missing value imputation in R. R J. 9(1), 207–218 (2017).Article 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, 2020). More