Ecology

1.FOREST EUROPE. State of Europe’s Forests (Ministerial Conference on the Protection of Forests in Europe, Bratislava, 2020).2.Gamfeldt, L. et al. Higher levels of multiple ecosystem services are found in forests with more tree species. Nat. Commun. 4, 1–8 (2013).
Google Scholar 
3.Brockerhoff, E. G. et al. Forest biodiversity, ecosystem functioning and the provision of ecosystem services. Biodiv. Conserv. 26, 3005–3035 (2017).
Google Scholar 
4.Mori, A. S., Lertzman, K. P. & Gustafsson, L. Biodiversity and ecosystem services in forest ecosystems: a research agenda for applied forest ecology. J. Appl. Ecol. 54, 12–27 (2017).
Google Scholar 
5.Forzieri, G. et al. Emergent vulnerability to climate-driven disturbances in European forests. Nat. Commun. 12, 1–12 (2021).
Google Scholar 
6.Senf, C. & Seidl, R. Mapping the forest disturbance regimes of Europe. Nat. Sustain. 4, 63–70 (2021).
Google Scholar 
7.Talluto, M. V., Boulangeat, I., Vissault, S., Thuiller, W. & Gravel, D. Extinction debt and colonization credit delay range shifts of eastern North American trees. Nat. Ecol. Evol. 1, 1–6 (2017).
Google Scholar 
8.Zhu, K., Woodall, C. W. & Clark, J. S. Failure to migrate: lack of tree range expansion in response to climate change. Glob. Change Biol. 18, 1042–1052 (2012).ADS 

Google Scholar 
9.Williams, J. W., Ordonez, A. & Svenning, J.-C. A unifying framework for studying and managing climate-driven rates of ecological change. Nat. Ecol. Evol. 5, 17–26 (2021).PubMed 

Google Scholar 
10.Jump, A. S. & Penuelas, J. Running to stand still: adaptation and the response of plants to rapid climate change. Ecol. Lett. 8, 1010–1020 (2005).PubMed 

Google Scholar 
11.Saltré, F. et al. Climate or migration: what limited European beech post-glacial colonization? Glob. Ecol. Biogeogr. 22, 1217–1227 (2013).
Google Scholar 
12.Svenning, J.-C. & Skov, F. Limited filling of the potential range in European tree species. Ecol. Lett. 7, 565–573 (2004).
Google Scholar 
13.Pedlar, J. H. et al. Placing forestry in the assisted migration debate. BioScience 62, 835–842 (2012).
Google Scholar 
14.Overpeck, J. T. & Breshears, D. D. The growing challenge of vegetation change. Science 372, 786–787 (2021).ADS 
CAS 
PubMed 

Google Scholar 
15.Strona, G. et al. Far from naturalness: How much does spatial ecological structure of European tree assemblages depart from potential natural vegetation? Plos One 11, e0165178 (2016).PubMed 
PubMed Central 

Google Scholar 
16.Giesecke, T. et al. Postglacial change of the floristic diversity gradient in Europe. Nat. Commun. 10, 1–7 (2019).CAS 

Google Scholar 
17.Kaplan, J. O., Krumhardt, K. M. & Zimmermann, N. The prehistoric and preindustrial deforestation of Europe. Quat. Sci. Rev. 28, 3016–3034 (2009).ADS 

Google Scholar 
18.Sabatini, F. M. et al. Where are Europe’s last primary forests? Divers. Distrib. 24, 1426–1439 (2018).
Google Scholar 
19.Nabuurs, G.-J. et al. Next-generation information to support a sustainable course for European forests. Nat. Sustain. 2, 815–818 (2019).
Google Scholar 
20.Williams, J. W., Jackson, S. T. & Kutzbach, J. E. Projected distributions of novel and disappearing climates by 2100 AD. Proc. Natl. Acad. Sci. 104, 5738–5742 (2007).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
21.Hoegh-Guldberg, O. et al. Assisted colonization and rapid climate change. Science 321, 345–346 (2008).CAS 
PubMed 

Google Scholar 
22.Jandl, R., Spathelf, P., Bolte, A. & Prescott, C. E. Forest adaptation to climate change – is non-management an option? Ann. For. Sci. 76, 1–13 (2019).
Google Scholar 
23.Dyderski, M. K., Paź, S., Frelich, L. E. & Jagodziński, A. M. How much does climate change threaten European forest tree species distributions? Glob. Change Biol. 24, 1150–1163 (2018).ADS 

Google Scholar 
24.Hanewinkel, M., Cullmann, D. A., Schelhaas, M.-J., Nabuurs, G.-J. & Zimmermann, N. E. Climate change may cause severe loss in the economic value of European forest land. Nat. Clim. Change 3, 203–207 (2013).ADS 

Google Scholar 
25.Thurm, E. A. et al. Alternative tree species under climate warming in managed European forests. For. Ecol. Manag. 430, 485–497 (2018).
Google Scholar 
26.Thuiller, W., Lavorel, S., Araújo, M. B., Sykes, M. T. & Prentice, I. C. Climate change threats to plant diversity in Europe. Proc. Natl. Acad. Sci. 102, 8245–8250 (2005).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
27.Isbell, F. et al. Biodiversity increases the resistance of ecosystem productivity to climate extremes. Nature 526, 574–577 (2015).ADS 
CAS 
PubMed 

Google Scholar 
28.Morin, X. et al. Long-term response of forest productivity to climate change is mostly driven by change in tree species composition. Sci. Rep. 8, 1–12 (2018).ADS 

Google Scholar 
29.Hisano, M., Searle, E. B. & Chen, H. Y. Biodiversity as a solution to mitigate climate change impacts on the functioning of forest ecosystems. Biol. Rev. 93, 439–456 (2018).PubMed 

Google Scholar 
30.Messier, C. et al. The functional complex network approach to foster forest resilience to global changes. For. Ecosyst. 6, 1–16 (2019).
Google Scholar 
31.Di Sacco, A. et al. Ten golden rules for reforestation to optimize carbon sequestration, biodiversity recovery and livelihood benefits. Glob. Change Biol. 27, 1328–1348 (2021).ADS 

Google Scholar 
32.Jacob, D. et al. EURO-CORDEX: new high-resolution climate change projections for European impact research. Reg. Environ. Change 14, 563–578 (2014).
Google Scholar 
33.Buras, A. & Menzel, A. Projecting tree species composition changes of European forests for 2061–2090 under RCP 4.5 and RCP 8.5 scenarios. Front. Plant Sci. 9, 1986 (2019).PubMed 
PubMed Central 

Google Scholar 
34.Chakraborty, D., Móricz, N., Rasztovits, E., Dobor, L. & Schueler, S. Provisioning forest and conservation science with high-resolution maps of potential distribution of major European tree species under climate change. Ann. For. Sci. 78, 1–18 (2021).
Google Scholar 
35.Noce, S., Collalti, A. & Santini, M. Likelihood of changes in forest species suitability, distribution, and diversity under future climate: The case of Southern Europe. Ecol. Evol. 7, 9358–9375 (2017).PubMed 
PubMed Central 

Google Scholar 
36.Hickler, T. et al. Projecting the future distribution of European potential natural vegetation zones with a generalized, tree species-based dynamic vegetation model. Glob. Ecol. Biogeogr. 21, 50–63 (2012).
Google Scholar 
37.Takolander, A., Hickler, T., Meller, L. & Cabeza, M. Comparing future shifts in tree species distributions across Europe projected by statistical and dynamic process-based models. Reg. Environ. Change 19, 251–266 (2019).
Google Scholar 
38.Chen, M. et al. Global land use for 2015–2100 at 0.05 resolution under diverse socioeconomic and climate scenarios. Sci. Data 7, 1–11 (2020).ADS 

Google Scholar 
39.Mauri, A., Strona, G. & San-Miguel-Ayanz, J. EU-Forest, a high-resolution tree occurrence dataset for Europe. Sci. Data 4, 1–8 (2017).
Google Scholar 
40.Strona, G., Mauri, A. & San-Miguel-Ayanz, J. A high-resolution pan-European tree occurrence dataset. Figshare https://doi.org/10.6084/m9.figshare.c.3288407.v1 (2016).41.Benito-Garzón, M. & Fernández-Manjarrés, J. F. Testing scenarios for assisted migration of forest trees in Europe. New For. 46, 979–994 (2015).
Google Scholar 
42.Thuiller, W., Lavorel, S., Sykes, M. T. & Araújo, M. B. Using niche-based modelling to assess the impact of climate change on tree functional diversity in Europe. Divers. Distrib. 12, 49–60 (2006).
Google Scholar 
43.Robinet, C. et al. A suite of models to support the quantitative assessment of spread in pest risk analysis. PLoS ONE 7, 10 (2012).
Google Scholar 
44.European Commission. The European Green Deal. (Publications office of the European Union, 2019).45.European Commission. EU Biodiversity Strategy for 2030, Bringing nature back into our lives. (Publications office of the European Union, 2020).46.European Commission. A sustainable bioeconomy for Europe: strengthening the connection between economy, society and the environment. (Publications office of the European Union, 2018).47.European Commission. New EU Forest Strategy for 2030. (Publications office of the European Union, 2021).48.Thuiller, W., Lafourcade, B., Engler, R. & Araújo, M. B. BIOMOD–a platform for ensemble forecasting of species distributions. Ecography 32, 369–373 (2009).
Google Scholar 
49.ICP Forests. International Co-operative Programme on Assessment and Monitoring of Air Pollution Effects on Forests. http://icp-forests.net/ (2019).50.Zając, A., Zając, M., Tertil, R. & Harman, I. Atlas rozmieszczenia roślin naczyniowych w Polsce–Distribution Atlas of Vascular Plants in Poland. (Nakladem Pracowni Chorologii Komputerowej Instytutu Botaniki Uniwersytetu – Laboratory of Computer Corology – Institute of Botany – Jagiellonian University, 2001).51.Gschwantner, T. et al. Common tree definitions for national forest inventories in Europe. Silva Fennica 43, 303–321 (2009).
Google Scholar 
52.Rivers, M. et al. European Red List of Trees. (International Union for Conservation of Nature and Natural Resources, 2019).53.Rocchini, D. et al. Anticipating species distributions: Handling sampling effort bias under a Bayesian framework. Sci. Total Environ. 584, 282–290 (2017).ADS 
PubMed 

Google Scholar 
54.Bartlein, P. J., Prentice, I. C. & Webb III, T. Climatic response surfaces from pollen data for some eastern North American taxa. J. Biogeogr. 35–57 (1986).55.Woodward, F. I. & Woodward, F. Climate and plant distribution. (Cambridge University Press, 1987).56.Harrison, S. et al. Towards a global scheme of plant functional types for ecosystem modelling, palaeoecology and climate impact research. J Veg Sci 21, 300–317 (2009).
Google Scholar 
57.Thuiller, W. BIOMOD–optimizing predictions of species distributions and projecting potential future shifts under global change. Glob. Change Biol. 9, 1353–1362 (2003).ADS 

Google Scholar 
58.Prentice, I. C. et al. Special paper: a global biome model based on plant physiology and dominance, soil properties and climate. J. Biogeogr. 117–134 (1992).59.Pouteau, R. et al. Potential alien ranges of European plants will shrink in the future, but less so for already naturalized than for not yet naturalized species. Divers. Distrib. 27, 2063–2076 (2021).
Google Scholar 
60.Naimi, B. USDM: Uncertainty analysis for species distribution models. https://www.rdocumentation.org/packages/usdm/versions/ (2015).61.Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. J. R. Meteorol. Soc. 25, 1965–1978 (2005).
Google Scholar 
62.Title, P. O. & Bemmels, J. B. ENVIREM: an expanded set of bioclimatic and topographic variables increases flexibility and improves performance of ecological niche modeling. Ecography 41, 291–307 (2018).
Google Scholar 
63.Teutschbein, C. & Seibert, J. Bias correction of regional climate model simulations for hydrological climate-change impact studies: Review and evaluation of different methods. J. Hydrol. 456, 12–29 (2012).ADS 

Google Scholar 
64.Ekström, M., Grose, M. R. & Whetton, P. H. An appraisal of downscaling methods used in climate change research. Wiley Interdiscip. Rev. Clim. Change 6, 301–319 (2015).
Google Scholar 
65.Beck, H. E. et al. Present and future Köppen-Geiger climate classification maps at 1-km resolution. Sci. Data 5, 1–12 (2018).ADS 

Google Scholar 
66.Baker, B., Diaz, H., Hargrove, W. & Hoffman, F. Use of the Köppen–Trewartha climate classification to evaluate climatic refugia in statistically derived ecoregions for the People’s Republic of China. Clim. Change 98, 113–131 (2010).ADS 

Google Scholar 
67.Barredo, J. I., Caudullo, G. & Dosio, A. Mediterranean habitat loss under future climate conditions: Assessing impacts on the Natura 2000 protected area network. Appl. Geogr. 75, 83–92 (2016).
Google Scholar 
68.Klausmeyer, K. R. & Shaw, M. R. Climate change, habitat loss, protected areas and the climate adaptation potential of species in Mediterranean ecosystems worldwide. PloS One 4, e6392 (2009).ADS 
PubMed 
PubMed Central 

Google Scholar 
69.Tabor, K. & Williams, J. W. Globally downscaled climate projections for assessing the conservation impacts of climate change. Ecol. Appl. 20, 554–565 (2010).PubMed 

Google Scholar 
70.Collins, M. et al. Long-term climate change: projections, commitments and irreversibility. in Climate Change 2013-The Physical Science Basis: Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change 1029–1136 (Cambridge University Press, 2013).71.Hengl, T. et al. SoilGrids250m: Global gridded soil information based on machine learning. PLoS One 12, e0169748 (2017).PubMed 
PubMed Central 

Google Scholar 
72.Zhang, L. et al. Consensus forecasting of species distributions: The effects of niche model performance and niche properties. PloS One 10, e0120056 (2015).PubMed 
PubMed Central 

Google Scholar 
73.Merow, C., Smith, M. J. & Silander, J. A. Jr. A practical guide to MaxEnt for modeling species’ distributions: what it does, and why inputs and settings matter. Ecography 36, 1058–1069 (2013).
Google Scholar 
74.Barbet-Massin, M., Jiguet, F., Albert, C. H. & Thuiller, W. Selecting pseudo-absences for species distribution models: how, where and how many? Methods Ecol. Evol. 3, 327–338 (2012).
Google Scholar 
75.De Jong, R., Verbesselt, J., Zeileis, A. & Schaepman, M. E. Shifts in global vegetation activity trends. Remote Sens. 5, 1117–1133 (2013).ADS 

Google Scholar 
76.Engler, R. & Guisan, A. MigClim: predicting plant distribution and dispersal in a changing climate. Divers. Distrib. 15, 590–601 (2009).
Google Scholar 
77.Engler, R., Hordijk, W. & Guisan, A. The MIGCLIM R package–seamless integration of dispersal constraints into projections of species distribution models. Ecography 35, 872–878 (2012).
Google Scholar 
78.Merow, C., Wilson, A. M. & Jetz, W. Integrating occurrence data and expert maps for improved species range predictions. Glob. Ecol. Biogeogr. 26, 243–258 (2017).
Google Scholar 
79.Caudullo, G., Welk, E. & San-Miguel-Ayanz, J. Chorological maps for the main European woody species. Data Brief 12, 662–666 (2017).PubMed 
PubMed Central 

Google Scholar 
80.Euro+Med. Euro+Med PlantBase – the information resource for Euro-Mediterranean plant diversity. http://ww2.bgbm.org/EuroPlusMed/ (2019).81.Summers, D. M., Bryan, B. A., Crossman, N. D. & Meyer, W. S. Species vulnerability to climate change: impacts on spatial conservation priorities and species representation. Glob. Change Biol. 18, 2335–2348 (2012).ADS 

Google Scholar 
82.García-Valdés, R., Zavala, M. A., Araujo, M. B. & Purves, D. W. Chasing a moving target: Projecting climate change-induced shifts in non-equilibrial tree species distributions. J. Ecol. 101, 441–453 (2013).
Google Scholar 
83.Lischke, H., Zimmermann, N. E., Bolliger, J., Rickebusch, S. & Löffler, T. J. TreeMig: a forest-landscape model for simulating spatio-temporal patterns from stand to landscape scale. Ecol. Model. 199, 409–420 (2006).
Google Scholar 
84.Tamme, R. et al. Predicting species’ maximum dispersal distances from simple plant traits. Ecology 95, 505–513 (2014).PubMed 

Google Scholar 
85.Thomson, F. J., Letten, A. D., Tamme, R., Edwards, W. & Moles, A. T. Can dispersal investment explain why tall plant species achieve longer dispersal distances than short plant species? New Phytol. 217, 407–415 (2018).PubMed 

Google Scholar 
86.Kattge, J. et al. TRY plant trait database–enhanced coverage and open access. Glob. Change Biol. 26, 119–188 (2020).ADS 

Google Scholar 
87.Mauri, A., Girardello, M. & Strona, G. EU-Trees4F. A dataset on the future distribution of European tree species, figshare, https://doi.org/10.6084/m9.figshare.c.5525688 (2021).88.Vítková, M., Müllerová, J., Sádlo, J., Pergl, J. & Pyšek, P. Black locust (Robinia pseudoacacia) beloved and despised: A story of an invasive tree in Central Europe. For. Ecol. Manag. 384, 287–302 (2017).
Google Scholar 
89.Muscarella, R. et al. ENM eval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods Ecol. Evol. 5, 1198–1205 (2014).
Google Scholar 
90.Allouche, O., Tsoar, A. & Kadmon, R. Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). J. Appl. Ecol. 43, 1223–1232 (2006).
Google Scholar 
91.Fitzpatrick, M. C. & Hargrove, W. W. The projection of species distribution models and the problem of non-analog climate. Biodivers. Conserv. 18, 2255–2261 (2009).
Google Scholar 
92.Elith, J. et al. A statistical explanation of MaxEnt for ecologists. Divers. Distrib. 17, 43–57 (2011).
Google Scholar 
93.R Core Team. R: A language and environment for statistical computing. (2020). More

Study designThe primary data source for this study is the cross-sectional 2013–2014 Demographic and Health Surveys for the DRC which is joined with remote-sensed environmental measures and land use data for mining and logging concessions extracted to DHS survey cluster locations. The DHS was administered using a multi-stage cluster survey design to represent the population of the DRC26. Briefly, survey clusters were selected to be representative of all 26 DRC provinces. Within clusters, households were randomly selected proportional to the population size, and within each household, adults ages 15–59 years were consented, interviewed, and asked to provide a dried blood spot (DBS) sample. Only adults who provided a DBS and consented for biospecimen use in future studies were included in this analysis. The outcome of prevalent malaria infections in the DRC was measured through PCR detection of the P. falciparum lactate dehydrogenase gene from DBS samples collected during DHS administration as described previously12.The main exposures were residence within 15 km of a mining concession and residence within 15 km of a logging concession. Additional covariates included individual-level variables for participant age, sex, use of a long-lasting insecticidal net (LLIN), education, and occupation; household variables for wealth, house roofing material, and the ratio of the number of household members using a bed-net to the total number of household members; and cluster variables for elevation, temperature, precipitation, vegetation, percentage of land cover identified as cropland, grassland, forest, and flooded/swamp land. All individual and household variables were obtained through the DHS. Occupation was recoded such that the manual labor and army category included laborers in mining and logging industries. Cluster variables were extracted from various satellite imagery platforms and other spatial datasets; the methods are described in more detail in the “Appendix”. The main exposures were extracted from geographic data sources as described below.Mining and logging concession data were obtained from the Global Forest Watch online repository27. Mining concessions were subset to only include operations that were active or in remediation spanning the DHS study years (2013–2014); logging concessions only included active operations during 2013. Distance to a mining or logging concession was measured from each cluster location to the boundary of a concession. Clusters were considered exposed to mining or logging if they were located within 15 km of a concession. This distance was chosen to account for the estimated 10 km maximum flight distance of a blood-fed mosquito5, with an additional 5 km to compensate for boundaries and non-residential land near the concessions. This range also accounts for the 5–10 km random spatial offset implemented by the DHS. Locations of mining and logging concessions along with cluster locations were mapped across the DRC. All maps were created in ArcGIS version 10.7.1, shapefiles for administrative boundaries were obtained from GADM.org.Data analysisCharacteristics of the study population were evaluated across quantiles of P. falciparum cluster prevalence and grouped by individual, household, and cluster level variables. To further examine distributions of malaria interventions and risk factors such as age, sex, LLIN use, occupation, household wealth, and household roof materials by mining and logging exposure, we compared mining exposed and logging exposed clusters with mining and logging doubly unexposed clusters stratified by urban and rural residence.We then modelled the prevalence odds of malaria across the DRC using hierarchical logistic regression models to account for the nested structure of the DHS data and to allow for inclusion of spatially varying effects. Models were implemented in a Bayesian framework using Integrated Nested Laplace Approximation (INLA) and stochastic partial differential equations for spatial effects28. In all models, we included two separate indicator terms for proximity to a mining concession and to a logging concession; since these areas are non-overlapping, the referent condition for each of these exposures is therefore locations exposed neither to mining nor to logging.The model fitting process followed two approaches. The first approach evaluated population-level effects of mining and logging on malaria prevalence adjusting for covariates and accounting for cluster-level random effects, which were assumed to vary independently across clusters. The second approach retained covariates and the cluster-level random intercept from the first model and additionally incorporated a spatial field to account for confounding due to space. For the spatial approach, two models were constructed. The first included a spatially varying intercept which borrowed information from neighboring cluster locations assuming a Gaussian random field. The second spatial model explored possible residual confounding due to environmental covariates by allowing spatially varying slopes for temperature, precipitation, vegetation, elevation, and land cover classes while including both independently and spatially varying intercepts across clusters. We introduced spatially varying slopes to account for the unobserved vector population across the DRC. Temperature, precipitation, vegetation, elevation, and various land cover classes have been shown to influence vector composition, survival, and competence for P. falciparum5,23,25, and associations with these covariates may vary due to their effects on the unobserved vector population. Using the spatial modelling approach, we also constructed a smoothed predicted prevalence map of malaria across the DRC, additional details are in the “Appendix”.For all models, confounding variables were selected based on a directed acyclic graph analysis and retained for adjustment if the 95% uncertainty interval (UI) of the variable excluded the null. Variables were coded as they were presented in the DHS with the exception of collapsing wealth into moderate or higher versus low wealth and recategorization of occupation as: professional, sales, or services; not working; manual labor or army; and agricultural work. All environmental variables were coded as continuous and scaled. Land cover variables were coded in intervals of 10 percentage points. Model comparison was done using Deviance Information Criterion (DIC), with the best fitting model having the smallest DIC29. All models were run using the ‘INLA’ package in R version 4.0.428, additional details are described in the “Appendix”.Differences in urban and rural residence were considered an important potential source of bias. Urban residence has been associated with lower prevalence of malaria due to many factors including different vector habitats, better access to healthcare, improved housing construction, and overall higher wealth4,12. To address possible bias introduced by urban residence, we stratified all models by urban and rural residence based on the DHS classification of clusters as urban or rural.A discrete set of confounding variables was identified from fixed effect models for mining and logging in rural and urban areas. The final adjustment set included age, sex, LLIN use, household wealth, temperature, precipitation, vegetation, and elevation. These variables had statistical or substantive significance and were adjusted for in all consecutive analyses.Ethical approval for this study was obtained from the University of North Carolina Institutional Review Board (UNC IRB# 20-3175) and the Kinshasa School of Public Health. Informed consent was obtained from all participants and all methods were conducted in accordance with guidelines and regulations set forth by the UNC IRB and the Kinshasa School of Public Health. More

As stated above this study aims to shed light on the deep uncertainties that are associated with the climate change phenomenon. The seasonally-averaged surface air temperature, hereafter simply referred to as temperature, was selected as the non-conditional climatic variable to be monitored within the Karkheh River Basin, Iran, during the baseline period (1975–2005). The CORDEX datasets (RCP 8.5) were employed to make climate-change projections.The first step in the proposed framework is to identify the most suitable theoretical distribution function to represent the stochastic behavior patterns of both historical and climate change data sets. Such identification considered the following theoretical distributions: normal, lognormal, exponential, Weibull, 3-parameter Weibull, extreme value, gamma, logistic, and loglogistic. It is important to note here that the primary strategy in this study is to analyze the data from a numeric standpoint without any presumption about the stochastic structure of the data44. As such, the study would opt for any distribution that is deemed fittest to describe the data. A summary of the fitted distributions to represent the prior distributions and likelihood functions is found in Tables S1 through S4 (see the Appendix). Furthermore, the climate-change period was divided into three mutually exclusive time frames which are short-term (2010–2039), mid-term (2040–2069), and long-term (2070–2099) future to gain a better understanding of the evoluton of future temperature changes.With Bayes’ theorem in mind, a Markov Chain Monte Carlo (MCMC) method was then applied to merge the prior distributions and the likelihood functions and to generate a sample set from the posterior distribution set. After a series of trials-and-errors, the sample size for the MCMC algorithm was set to be 1000 (n = 1000). These generated sample sets were then used to specify the most suitable theoretical pdfs to represent the posterior distribution functions. Figure 3, for instance, illustrates the most appropriate theoretical distribution that could represent the posterior distribution for the Seimareh sub-basin during spring under the short-term period.Figure 3The step-by-step process of computing the posterior pdf: (a) the prior distribution of Seimareh sub-basin during spring and the likelihood function of this sub-basin in the short-term future; (b) the histogram of the generated samples; and (c) the posterior distribution.Full size imageFigure 4 demonstrates the frequency with which each specific theoretical distribution functions was deemed the most suitable to characterize the prior, likelihood, and posterior distributions. Analyzing the fitted pdfs in Fig. 4 reveals an important point about the nature of RCMs’ raw projections. Specifically, the most frequently chosen distribution function for prior and posterior distributions is the 3-parameter Weibull. As for the likelihood function, however, it was the normal distribution that outperformed other available alternatives. Furthermore, the type of selected theoretical distribution for prior and posterior pdfs seems far more diverse compared to those from the likelihood functions. In fact, the likelihood functions were only limited to three types of distributions, most of which are normal distributions. Keep in mind that these functions are the most suitable pdfs that were fitted to the RCMs projected results. The cause behind this notion might be traced back to the nature of RCMs’ projections. RCMs operate at a finer horizontal resolution than GCMs, and thus they provide localized and high-resolution detailed climatic information that can be of importance for many management purposes, especially in regions with complex topography. However, the analyzed data revealed that among the distributions fitted for the likelihood function the normal distribution was found to be the best distribution to describe the data 70% of the time. This could be interpreted as signaling that employing RCMs’ raw projections, especially for regions that have considerable volatility in their climatic variables, should be used with caution, and further adjustment to the raw projected data may be required in some cases. Note that from a statistical standpoint, the normal distribution is not heavy-tailed, and as such, may not be the best way to portray this data. The fact that, in most cases, it has been selected as the best way to portray the stochastic nature of the likelihood function (i.e., RCM’s projections) means that innate characteristics of these data might prevent them to truly represent these types of variables on their own.Figure 4The frequency of using each individual theoretical pdfs as prior, likelihood function, and posterior distributions.Full size imageFigure 5 provides additional information regarding the frequency in which each individual theoretical distributions were deemed suitable to represent the posterior pdfs. While posterior distribution sets are, indeed, the most diverse in terms of the number of different types of distributions, a significant proportion of fitted pdfs (approximately 52%), however, are fitted by the 3-parameter Weibull distribution. Further information regarding the fitted distributions to represent the posterior pdfs is found in Tables S5 to S7 (Appendix).Figure 5The frequency of using different theoretical distributions as posterior pdfs.Full size imageThe computed posterior distribution functions can be interpreted as modified representations of the stochastic behavior of temperature variable concerning the short-term, mid-term, and long-term climate change projections. In that spirit, employing the confidence interval of 95%, the average temperature of the entire basin is depicted in Fig. 6 associated with historical and climate change conditions. Two sets of behavior patterns are observed. The first one is a broad trend in summer. The second pattern describes the rest of the seasons. In summer (Fig. 6b) the presence of a mild, yet, steady positive trend (upward) is detected. Here, one can expect the average temperature of the basin to increase steadily with the passage of time. As for the rest of the seasons, while it seems that the average temperature of the basin would experience a mild drop in the short-term, the temperature would begin to rise with a steady trend with time. In spring (Fig. 6a) and autumn (Fig. 6c) time series, it is projected that the expected average temperature in the basin would eventually surpass those that had been experienced in the baseline condition in the mid-term and long-term future. Concerning winter temperature it is seen in Fig. 6d that it is projected to increase over time. Yet, it has been estimated that it might not reach the observed average temperature of the basin in neither of the expressed time frames. Of course, given the upward trend in the data, this temperature would indeed be reached in a longer timeframe. It is worth noting that these patterns are in line with the idea that the earlier impacts of climate change are to amplify the historical patterns in climatic variables. That is why the data show a slight drop in colder seasons and an uptick in the warmer ones. That is, of course, until eventually, a new climatic equilibrium is reached on a global scale. At this point, the temperature as shown here would start to increase gradually.Figure 6The historical and simulated average temperatures of the entire basin with the 95%confidence interval in (a) spring; (b) summer; (c) autumn; and (d) winter.Full size imageThe other notable implication that can be understood by analyzing Fig. 6 is the variation in the width of the confidence intervals under baseline and climate change conditions. In comparison to the baseline condition, the length of the 95% confidence intervals would dramatically decrease under climate change conditions. This shrinking indicates that the generated results are more densely surrounding the central tendency measure herein chosen as the mean (μ) of the data. This notion is still in line with the idea that RCM’s projections mostly resemble the stochastic characteristic of normal distributions. A normal distribution is by nature not a heavily-tailed distribution, meaning that it rarely generates tail values. Even though the MCMC framework has mitigated this effect to some extent, they inevitably inherit this stochastic property from the likelihood functions.Again, to truly understand the obtained results here, one must first acknowledge how Bayesian models work. The main idea behind a Bayesian-based framework is to adjust the prior assumptions about a stochastic phenomenon through observed samples. In this case, the prior information represents the historical data, and the likelihood function (i.e., the samples) is obtained from RCM projects. As can be seen here, while RCMs’ projections might be perfectly capable of portraying the normal behavior of a variable under climate change conditions, which is usually sufficient for most lumped evaluation of climate change impact assessments, they might not be suitable to study extreme hydro-climatic events. The main problem with the raw RCM projections is that they follow a normal distribution, which is a symmetric distribution. Figure 4 suggests that while the MCMC framework here is mitigating this impact the final projections inherit this property from the likelihood functions. This simply means that while any RCM-based projection is perfectly suitable to understand the general outline of the climate change impacts, they are not the best option to study extreme events because even by modifying their pdfs, they rarely generate truly extreme values. The average temperatures in all sub-basins under baseline and climate change conditions are summarized in Table 1.Table 1 The average surface air temperature in all sub-basins under baseline and climate change conditions (°C).Full size tableAs for the impact of climate change, it is clear that these data are associated with deep uncertainty; that is, the parameters used to describe the stochastic behavior of a variable may be subjected to some degree of uncertainty. These parameters, μ for one, may also be represented by a pdf of their own. This study focuses on highlighting this type of deep uncertainty that might interfere with the central tendency measure μ.The deep uncertainty in this instance dictates that the recorded parameters for each posterior distribution are not deterministic values. While for a given prior distribution and likelihood function the MCMC would lead to a specific type of posterior pdf, the parameters that are used to define this pdf (e.g., μ), could vary each time the algorithm is used. If this variation is mild, there is more certainty about the nature of the variable’s stochastic behavior pattern (i.e., the posterior distribution function). If it is determined that the parameters are experiencing severe variations then the deep uncertain environment would leave the decision-makers unsure about the variables’ stochastic behavior pattern.With that idea in mind the combination of prior distributions and likelihood functions was executed for 100 times, and in each iteration the mean of 1000 samples was recorded. A theoretical distribution function was then fitted to the recorded values. Naturally, if the recorded values are generally close to one another numerically, the parameters of the computed posterior pdfs are less subjected to deep uncertainty. If, however, these values show significant fluctuation then the deep uncertainty of climate change would impede predictions of the stochastic behavior pattern of temperature. Figure 7, for example, portrays the uncertainty of the computed μ parameter for Seimareh sub-basin in spring under short-term future condition.Figure 7The uncertainty of the computed μ parameter for Seimareh sub-basin in spring under short-term future condition demonstrated by (a) a histogram and (b) a probability distribution function.Full size imageFigure 8 demonstrates the number of times each theoretical distribution was chosen to portray the stochastic behavior of the μ parameter. As can be seen here, the normal and lognormal distributions are the most common pdfs used to describe the variation in the μ parameter. One should also note the fact that about 65% of the distributions used to describe the future condition are normal distributions. The list of fitted pdfs is summarized in Tables S8 to S10.Figure 8The frequency with which each theoretical distribution was found suitable to describe the stochastic distribution of the μ parameter.Full size imageTable 2 summarizes the variation in the computed μ parameter in each given sub-basin. It is seen in Table 2 the 95% confidence interval of the μ parameter in all cases ranges between ± 0.1 and ± 0.3 °C. In 55% of the cases, this interval was found not to be more than ± 0.1 °C, and, furthermore, in 97% of them the interval was less than ± 0.2 °C. Needless to say, a widened confidence interval for the μ parameter can only signal that the deep uncertainty has a more pronounced impact on the temperature’s stochastic behavior. As for the case of the spring data set of the Seimareh sub-basin under the short-term condition, or the case of the Gharesou sub-basin’s winter data series under short-term period, the confidence interval for the μ parameter is estimated to be ± 0.3 °C wide. This indicates that compared to other projected posterior pdfs there is less certainty about the predicted stochastic behavior pattern of temperature variable for these particular cases. As shown in Table 2 in some cases, the variation in the projected μ temperature’ posterior pdfs is decreasing over time (for a given season over different timeframes). As discussed earlier, this was interpreted as the deep uncertainties of the climate change projections, meaning that lower volatility in this measure indicates that the said variable is less affected by the deep uncertainty of the climate-change phenomenon. This observation is in line with the general belief that, in the near future, the climate change phenomenon is most likely to intensify the historical patterns in climatic variable, but gradually we expect to see an upward trend in temperature in the longer run45. In this case, there is more volatility in the earlier time frames, but as time progresses, this volatility seems to decrease in some cases. This means that the obtained projections are showing less uncertainty about the outline behavior of the parameter for the long-term future as the models that are simulating the climatic behavior under climate change conditions have already reached a new equilibrium by that point.Table 2 The variation in the computed μ of the temperature’ posterior pdfs.Full size table More

ModelTo investigate the evolution of cyclically competing species with intraspecific interaction which sensitively plays to the territory awareness, we employ the spatial RPS model11,19,20,23. At the microscopic level, the model can be demonstrated on a lattice system, and for convenience, we consider a square lattice of size N with periodic boundary conditions where all sites have von Neumann neighbors. Each site can be occupied by an individual from one of the three species (referred to as A, B, and C, respectively) or left empty(E), and thus the system describes a limited carrying capacity. In addition, to explore the effect of territory awareness on intraspecific interaction, we assume that the given lattice is divided into two areas of the same size which may possibly realize different territorial ranges. Here we simply divide the two regions into the top and bottom halves of the given square lattice. To reflect the territorial awareness on intraspecific interaction, we distribute population in each group into two sub-networks randomly, and denote species (X_1) for the top and (X_2) for the bottom ((X in {A, B, C})) to distinguish the emergence of intraspecific interaction between individuals who lie on different domains. The distribution of all species with respect to the separation of the domain is illustrated in Fig. 1a.Figure 1Schematic diagrams of network structure and the invasion rules among species. (a) Each circle represents a node, and individuals of species A, B, and C are evenly and randomly distributed on each node. To realize territorial awareness, the lattice is divided into two regions of equal size: the top and the bottom where the dashed line indicates the regional boundary. Two genera of the same species are distributed in different regions, and different color markers represent different species types. Nodes without color markers are empty nodes. (b) Interspecific interaction among three species A, B, and C (indicated by three boxes) occurs cyclically with a rate (p_1). A box of each group describes the intraspecific interaction between individuals who belong to different territories where the interaction is regulated by territorial consciousness. Here intraspecific interaction in each group occurs with a rate (k_i cdot p_2) ((i in {A, B, C})).Full size imageUnder the given assumption for the lattice, all the interactions between individuals occur within nearest neighboring sites by the following set of rules (see Fig. 1b):$$begin{aligned}&A_{i}+B_{j}overset{p_{1}}{longrightarrow }A_{i}+E,quad B_{i}+C_{j}overset{p_{1}}{longrightarrow }B_{i}+E,quad C_{i}+A_{j}overset{p_{1}}{longrightarrow }C_{i}+E, end{aligned}$$
(1)
$$begin{aligned}&A_{i}+A_{j}overset{k_{A} cdot p_{2}}{longrightarrow }A_{i}+E,mathrm{or},E+A_{j}, quad B_{i}+B_{j}overset{k_{B} cdot p_{2}}{longrightarrow }B_{i}+E,mathrm{or},E+B_{j}, quad C_{i}+C_{j}overset{k_{C} cdot p_{2}}{longrightarrow }C_{i}+E,mathrm{or},E+C_{j},quad ine j, end{aligned}$$
(2)
$$begin{aligned}&A_{i}+Eoverset{r}{longrightarrow }A_{i}+A_{i},quad B_{i}+Eoverset{r}{longrightarrow }B_{i}+B_{i},quad C_{i}+Eoverset{r}{longrightarrow }C_{i}+C_{i}, end{aligned}$$
(3)
$$begin{aligned}&A_{i}+otimes overset{m}{longrightarrow }otimes +A_{i},quad B_{i}+otimes overset{m}{longrightarrow }otimes +B_{i},quad C_{i}+otimes overset{m}{longrightarrow }otimes +C_{i}, end{aligned}$$
(4)
where (i, j=1, 2). The mark (otimes) stands for any species or empty sites. Relation (1) describes interspecific interaction among three species which occurs cyclically with a rate (p_1): (A_{i}) dominates (B_{i}), (B_{i}) dominates (C_{i}), and (C_{i}) dominates (A_{i}) ((i=1, 2)). The defeated individual dies and the site becomes an empty site. Relation (2) demonstrates the intraspecific interaction which will sensitively depend on territorial awareness. Since we assume the intraspecific interaction is related to the territorial consciousness, the rate in each species may be defined by (k_{A} cdot p_{2}), (k_{B} cdot p_{2}), (k_{C} cdot p_{2}) for species A, B, C, respectively, where (p_{2}) is the given rate of interaction, k is the the sensitive parameter to territorial awareness. Similar to previous works, the result of intraspecific interaction eventually results in a death of one individual at random with a 1/2 chance. Relation (3) stands for the reproduction with a rate r which is allowed when an empty site in neighbors is selected, and migration defined by an exchange between two neighboring sites is denoted in Relation (4). Based on the theory of random walks38, it occurs with a rate (m=2MN) where M and N indicate individuals’ mobility and a system size, respectively, as usual to previous works. Thus, an actual time step is defined when each individuals has interacted with others once on average, i.e., N pairwise interactions will occur in one actual time step unit. In order to make an unbiased comparison with previous works15,19,20,21 and for the convenience of interpretations, we assume parameters as (p_{1}=p_{2}=r=1) and (k_{A}=k_{B}=k_{C}=k) (see the Methods for the meanings of specific parameters) in our simulations. Three species are divided into two types to distinguish distributions on different regions: (A_{1,2}), (B_{1,2}), and (C_{1,2}), and randomly distributed initially on a square lattice of size (N= 300 times 300). In addition, in all our simulations, species coexistence refers to the coexistence of (A_i), (B_j), and (C_k) for any combination of (i,,j,,k in left{ 1,,2 right}).Biodiversity under territorial awarenessWe first consider the effect of territorial awareness on species biodiversity. In general, it is well-known that, the spatial RPS game exhibits a transition of survival states from coexistence to extinction (which is presented by the uniform state) as individuals’ mobility increases. The phase transition occurs when M exceed a certain value, referred to as a critical mobility (M_c = (4.5 pm 0.5) times 10^{-4}), which is identified in Ref.11. To address the effect of territorial awareness, we consider two different mobility values (M=1 times 10^{-5}) and (M=1 times 10^{-3}) which eventually yield different survival states: coexistence and extinction, respectively, for different sensitivity parameter k.In general, the total simulation time T in classic spatial RPS games is considered as (T = N) which can yield the extinction for the critical mobility (M_c)11. In this regard, using the time (T=N) may yield different results for species evolution and corresponding survival states due to stochastic events, and such behaviors may be induced by the choice of mobility. In our simulations, since we consider two different mobility values where the one (Fig. 2a–c) is quite lower and the other (Fig. 2d–f) is higher than (M_c), we thus consider different simulation times at (M=1 times 10^{-5}) and (M=1 times 10^{-3}): more than 490, 000 and 180, 000 steps, respectively, to obtain robust features on species survival states. The time dependent evolution of densities are illustrated in Fig. 2 where the top and bottom panels are obtained from simulations with the first 250, 000 and 140, 000 steps, respectively.Figure 2Time dependent evolution of densities in the system for different M and k. Top and bottom panels are obtained with (M=1times 10^{-5}) and (M=1times 10^{-3}), respectively, and the sensitivity parameter k is given by (k=5), 10, and 20 from the left to right in each row. (a–c) Regardless of the choice k, the low mobility still leads species coexistence as usual. (d–f) At high mobility regimes, the system also always exhibit the extinction state.Full size imageEven if different k are considered, the panels in Fig. 2 show features similar to previous works11,14,15,16,18,20: coexistence and extinction for tops and bottoms, respectively. At the low mobility (M=1times 10^{-5}) as shown in Fig. 2a–c, even if the individuals located in different domains in each group disappear, the spatial RPS game eventually exhibits coexistence as k increases since some of individuals in A, B, C are survived. For instance, in our simulations, coexistence can be presented by survival of species (A_1), (B_1), (C_1) (Fig. 2a–b) or (A_2), (B_2), (C_1) (Fig. 2c). Since the typical waiting time for extinction is exponentially increasing to the size N at low mobility11, there will be extinction and eventually only one species will dominate the system after extremely long times. Thus, within the finite time steps, one type of each species will disappear slowly with the increase of k and the system exhibits coexistence.On the other hand, the high mobility (M=1times 10^{-3}) leads the extinction and only one species dominate the whole domain. As shown in Fig. 2d–f, the extinction that is defined by the two types of the species disappear occurs and the only one species finally dominate the system [e.g., (C_2), (B_2), and (C_1) for (k=5), 10, and 20, respectively]. In this case, the increase in k has little effect on the disappearance of one of the species, but has a tendency to accelerate the complete extinction of the second species. Take Fig. 2d for example, when species (A_2), (B_1) and (B_2) became extinct, (A_1), (C_1) and (C_2) is left in the system, the density of (A_1) in the system had an absolute advantage, while (C_1) and (C_2) had intraspecific interaction. Since the intensity of intraspecific interaction sensitive to territorial awareness was greater than that of interspecific interaction, the interaction was mainly intraspecific interaction between (C_1) and (C_2), then (C_1) was defeated into extinction, (C_2) preyed on the only specie (A_2) remaining in the system and eventually occupied the whole system. As k value affects the intraspecific interaction intensity, it determines the waiting time for the extinction of two species in the system. For example, we found that the larger the k value is, the shorter the waiting time for the extinction of two species is, as shown in Fig. 2e–f. But this is only the observation result of a single simulation. Due to the randomness of the simulation, this phenomenon needs further verification, so we give specific results about the effect of k value on the average extinction time in the next section. Due to stochastic events during Monte-Carlo simulations, the combination of survival species for coexistence and extinction at the final step can be different, but the such states at two mobility regimes will be still maintained. Fig. 2 may impose the follows: territorial awareness on intraspecific interaction can eventually yield similar feature to previous works in a broad aspect, but the composition of the surviving species type for each state may vary.Average extinction time versus territorial awarenessWhile survival states in both cases are consistent with previous works on the effects of species migration in Fig. 2, we found an interesting feature that the evolutionary time when some type of species disappear is changed depending on k. To be concrete, at (M=1 times 10^{-5}), we found that one type of each species (A_1), (B_1), (C_1) (Fig. 2a) will eventually coexist while their companion species (A_2), (B_2), (C_2) are extinct as t exceeds (t approx 50{,}000). As k is increasing, the time point when one genus of each species disappears shows an increasing pattern as presented in Fig. 2b,c. The opposite trend can be captured at high mobility (M=1 times 10^{-3}), that is, the increase of k seems to shorten the evolution time of two species extinction in the system. Based on these observations, we may assume that the critical time for such disappearance phenomena has a certain relationship with k and the relation may differ to the choice of M.To answer the issue, we measure the average extinction time T. In classic RPS games, traditionally, the extinction state on spatially extended systems has been identified by the uniform state that only one species dominates whole domain11,14,15,16,18,20. As shown in Fig. 2a–c which ultimately describe a coexistence state in a finite time, however, any one of type in each species disappeared and the time associated with the phenomena is changed by the strength of k. In a slightly different aspect to the classic meaning of extinction, we here define the average extinction time T with respect to the regime of mobility: (a) the evolutionary time when one genus of each species disappears for low mobility and (b) the time when two of the three species disappear completely for high mobility. In this consideration, for both given cases of M in Fig. 2, the average extinction time T in each k is measured from 30 independent realizations and presented in Fig. 3.Figure 3The average extinction time T as a function of the territorial sensitive parameter k. (a) Two cases of fixed mobility in territorial sensitive intraspecific interaction. For low mobility (M=1times 10^{-5}), the time T which is measured by detecting the time when one genus of each species disappears tends to increase with the increase of k, i.e., the high sensitivity of territorial awareness has the effect of delaying the waiting time for extinction. Similarly, it can be seen that at high mobility (M=1times 10^{-3}), an increase in k value will also delay the waiting time for extinction, but the effect is much more gentle. (b) At low mobility value (M=1times 10^{-5}), traditional intraspecific interaction (i.e., intraspecific interaction among all individuals of the same species, regardless of territorial residence) was compared with territorial sensitive intraspecific interaction. Here, for the traditional case, k represents intraspecific interaction intensity, and the running time of the simulation is 810, 000 steps. In the case that the final steady state has not occurred before the end of the simulation, we take the maximum time step ((t=810{,}000)) as the extinction time T value, which causes the blue line to become gentle when (k >14). Compared with the traditional situation, the intraspecific interaction affected by territorial awareness significantly reduced the average extinction time, that is, accelerated the damage of species diversity in the system. The results were averaged from 30 independent simulations, and error bars (using standard errors, which defined as the sample standard deviation divided by the square root of the number of samples) are shown in the figure.Full size imageAs shown in Fig. 3a, we find clearly that the average extinction time is obviously affected by the strength of sensitivity coefficient k, especially, when the mobility is low. When species has no consciousness on territories ((k=0)), the system becomes exactly the classic RPS model11 since intraspecific interaction is undefined, and the waiting time T generally tends to increase exponentially to the choice of M. However, our simulation shows the T is approximately measured at (T=110{,}000) at (k=0). Traditionally, it is well known that the average waiting time for extinction in the classic RPS game is taken (T=N) near the critical mobility regime ((M approx M_c)), and the coexistence duration is exponentially increasing as M decreases from (M_c). Within this knowledge, our simulation results may seem inconsistent with the general concept of extinction time. In our model, however, the definition of extinction is different at the low mobility regime, and the change into a single RPS system as one genus of each individual disappears may have a similar meaning to the previous definition of extinction in some sense, the above result can be said to be reasonable.The important point is actually addressed for (k >0). In this case, species can allow intraspecific interaction and the strength of intraspecific interaction is also increasing since the territorial awareness is intensified. As a result, it is found that the average extinction time T shows a tendency to gradually increase with the increase of k at (M=1 times 10^{-5}). In addition, this trend can also be observed at (M = 1 times 10^{-3}), but it is more gradual. To investigate whether the tendency to prolong the waiting time for extinction time at low migration rates is caused by territorial awareness or the existence of intraspecific interaction, we compared traditional intraspecific interaction (i.e., intraspecific interaction among all individuals of the same species, regardless of territorial residence, which equivalent to removing the condition (ine j) from Relation (2)) with territorial-sensitive intraspecific interaction in our model, the results are shown in Fig. 3b. We found that in the presence of intraspecific interaction, the average extinction time increased with the intensity of intraspecific interaction. Specifically, the stronger the intraspecific interaction, the slower the loss of species diversity. However, compared with the traditional situation, intraspecific interaction influenced by territorial consciousness controlled the delay of extinction to a certain extent. Even if our simulations have been carried on for two specific M, it is obvious that the territorial awareness can affect the average extinction time, and we suggest that a strong sense of territoriality can also delay species extinction and lead to long-term coexistence of systems at low mobility regimes, although the introduction of territoriality leads to faster damage to species diversity than is traditionally the case, while it does not affect significantly on the extinction time and the biodiversity (which eventually appears as extinction) at high mobility regimes.Evolution of the interface between territoriesFrom the investigation on the average extinction time in Fig. 3, we know that the territorial awareness can affect not only species survival but also the maintenance period of survival states. Here, we may wonder why the territorial awareness can affect the waiting time to extinction. In order to investigate such an issue, we observe evolution of the spatial system, in particular invasion between species near the border on two territories, i.e., the evolution of the interface. To capture the phenomena in detail, we consider pattern formations associated with the given two mobility values at the initial state of the evolution (e.g. (t=1000)) which are represented in Fig. 4.Figure 4The typical snapshots of evolution on patterns at (t=1000) for different k: 0.5 for (a) and (e), 2 for (b) and (f), 10 for (c) and (g), and 20 for (d) and (h), where the mobility is considered as (M = 1times 10^{-5}) for tops and (M=1 times 10^{-3}) for bottoms. Different colors correspond to different species types, as shown in Figs. 1 and 2, with white indicating vacancy. As k increases, the invasion among species between two territories occurs more gradually, and such phenomena are clearly observed for the high mobility as shown in the panel (h).Full size imageThe top and bottom panels in Fig. 4 exhibit spatial patterns for the low and high mobility values, respectively. For (M=1 times 10^{-5}), when the value k is small such as (k=0.5) (see Fig. 4a), interspecific interaction can occur more frequently than intraspecific interaction among all pairwise reactions (1)–(4). The system can exhibit similar pattern formations to the classic RPS game11. Three species, even if they are distinguished into six subgroups, are spirally entangled with clearly exhibiting spiral waves which are appeared in both two territories. Since the given lattice has periodic boundaries, species in both territories can migrate to the other region each other, but such migration is weak because the normalized probability for migration (Relation (4)) is small at the low mobility. Thus, when the system exhibits coexistence, it may be possible to predict that the top and bottom territories present dominance of species (X_1) and (X_2) ((X in {A, B, C})), respectively, while our simulations only present spatial patterns at the first 1,000 steps which may be too short to lead phase transitions.We also find that the spiral-wave patterns are getting to fuzzy as k increases. In particular, such fuzzy patterns are conspicuous near the boundary between the two territories at the large k (see Fig. 4c,d). The increase of k directly means the intensification of intraspecific interaction, and according to the setting on the initial distribution of population, intraspecific interaction will have many chances to occur in the vicinity of the boundary than near the top and bottom periodic boundaries. Frequent intraspecific interaction can provide as many chances to allow reproduction as possible, and high intraspecific interaction rate can dominate on pairwise invasions than interspecific interaction.In the vicinity of the border between two territories, the occurrence of intraspecific interaction is observed more prominently at (M=1 times 10^{-3}), and such features are clear as k increases. To be concrete, compared with figures among Fig. 4e–h, we found that empty sites are produced near the border and their presence is clear for high strength k such as 10 and 20 (Fig. 4g–h). In this case, the two domains appear to be more clearly divided and each domain is dominated by a single RPS system. Each single RPS system shows extinction state (only one genus survives) at high mobility, and eventually shows extinction state through interspecific or intraspecific interaction depending on the type of surviving genus. This is in good agreement with the results we obtained in Fig. 2. However, it can be seen from Fig. 3 that the time for each domain system to reach extinction at high mobility is very short compared to that at low mobility, and this has no relation with the degree of territorial awareness in interspecific interaction.From our simulations, we find that: the relationship between territorial awareness and the average extinction time is particularly prominent at the low mobility, and the likelihood of intraspecific interaction is relatively high near territorial boundaries. Under these considerations, we may expect a new relationship between the delay of the extinction time and boundary of two territories. To uncover this veil, we try to quantify a width for occurrence of intraspecific interaction near the border between two area with respect to the territorial awareness. Specifically, we give each node on a two-dimensional grid a coordinate, defined by its row and column position. For each column (j=1,ldots ,L), calculate the interface width, defined as I:$$begin{aligned} I_{j}= {left{ begin{array}{ll} P(1,j)-P(2,j),&{}quad text {if } P(1,j) >P(2,j)\ 0,&{}quad text {if } P(1,j) More

1.Gibson, L. et al. Primary forests are irreplaceable for sustaining tropical biodiversity. Nature 478(7369), 378–381. https://doi.org/10.1038/nature10425 (2011).CAS 
Article 
ADS 

Google Scholar 
2.Newbold, T. et al. A global model of the response of tropical and sub-tropical forest biodiversity to anthropogenic pressures. Proc. R. Soc. B. 281, 20141371. https://doi.org/10.1098/rspb.2014.1371 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
3.Hansen, M. C. et al. The fate of tropical forest fragments. Sci. Adv. 6(11), eaax8574. https://doi.org/10.1126/sciadv.aax8574 (2020).Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
4.Peres, C. A. et al. Biodiversity conservation in human-modified Amazonian Forest landscapes. Biol. Conserv. 143, 2314–2327. https://doi.org/10.1016/j.biocon.2010.01.021 (2010).Article 

Google Scholar 
5.PRODES INPE. Monitoring Deforestation of the Brazilian Amazon Forest by Satellite. TerraBrasilis (inpe.br) (accessed in october 2020, 2020).6.Barlow, J. et al. Quantifying the biodiversity value of tropical primary, secondary, and plantation forests. Proc. Natl. Acad. Sci. 104, 18555–18560. https://doi.org/10.1073/pnas.0703333104 (2007).Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
7.Peres, C. A., Barlow, J. & Laurance, W. F. Detecting anthropogenic disturbance in tropical forests. Trends Ecol. Evol. 21, 227–229. https://doi.org/10.1016/j.tree.2006.03.007 (2006).Article 
PubMed 

Google Scholar 
8.Arroyo-Rodríguez, V. et al. Designing optimal human-modified landscapes for forest biodiversity conservation. Ecol. Lett. 23, 1404–1420. https://doi.org/10.1111/ele.13535 (2020).Article 
PubMed 

Google Scholar 
9.Gardner, T. A. et al. Prospects for tropical forest biodiversity in a human-modified world. Ecol. Lett. 12, 1–21. https://doi.org/10.1111/j.1461-0248.2009.01294.x (2009).Article 

Google Scholar 
10.Hardwick, S. R. et al. The relationship between leaf area index and microclimate in tropical forest and oil palm plantation: Forest disturbance drives changes in microclimate. Agric. For. Meteorol. 201, 187–195. https://doi.org/10.1016/j.agrformet.2014.11.010 (2015).Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
11.Sambuichi, R. H. et al. Cabruca agroforests in southern Bahia, Brazil: Tree component, management practices and tree species conservation. Biodivers. Conserv. 21, 1055–1077. https://doi.org/10.1007/s10531-012-0240-3 (2012).Article 

Google Scholar 
12.Devictor, V., Julliard, R. & Jiguet, F. Distribution of specialist and generalist species along spatial gradients of habitat disturbance and fragmentation. Oikos 117, 507–514. https://doi.org/10.1111/j.0030-1299.2008.16215.x (2008).Article 

Google Scholar 
13.Banks-Leite, C. Using ecological thresholds to evaluate the costs and benefits of set-asides in a biodiversity hotspot. Science 345, 1041–1045. https://doi.org/10.1126/science.1255768 (2014).CAS 
Article 
PubMed 
ADS 

Google Scholar 
14.Newbold, T. et al. Global patterns of terrestrial assemblage turnover within and among land uses. Ecography 39, 1151–1163. https://doi.org/10.1111/ecog.01932 (2016).Article 

Google Scholar 
15.Paglia, A. P. et al. Annotated checklist of Brazilian mammals. Occas. Pap. Conserv. Int. 6, 1–82 (2012).
Google Scholar 
16.Dirzo, R. et al. Defaunation in the anthropocene. Science 345, 401–406. https://doi.org/10.1126/science.1251817 (2014).CAS 
Article 
PubMed 
ADS 

Google Scholar 
17.Estrada, A. et al. Impending extinction crisis of the world’s primates: Why primates matter. Sci. Adv. 3, e1600946. https://doi.org/10.1126/sciadv.1600946 (2017).Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
18.Newbold, T. et al. Global effects of land use on local terrestrial biodiversity. Nature 520, 45–50. https://doi.org/10.1038/nature14324 (2015).CAS 
Article 
PubMed 
ADS 

Google Scholar 
19.Phillips, H. R., Newbold, T. & Purvis, A. Land-use effects on local biodiversity in tropical forests vary between continents. Biodivers. Conserv. 26, 2251–2270. https://doi.org/10.1007/s10531-017-1356-2 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
20.Teixeira, D. F., Guillera-Arroita, G., Hilário, R. R., Fonseca, C. & Rosalino, L. M. Influence of life-history traits on the occurrence of carnivores within exotic Eucalyptus plantations. Divers. Distrib. 26, 1071–1082. https://doi.org/10.1111/ddi.13114 (2020).Article 

Google Scholar 
21.Asner, G. P. et al. Selective logging in the Brazilian Amazon. Science 310, 480–482. https://doi.org/10.1126/science.1118051 (2005).CAS 
Article 
PubMed 
ADS 

Google Scholar 
22.Robinson, J. G. & Redford, K. H. Body size, diet, and population density of neotropical forest mammals. Am. Nat. 128, 665–680. https://doi.org/10.1086/284596 (1986).Article 

Google Scholar 
23.Cardillo, M. et al. Multiple causes of high extinction risk in large mammal species. Science 309, 1239–1241. https://doi.org/10.1890/05-0112 (2005).CAS 
Article 
PubMed 
ADS 

Google Scholar 
24.Almeida-Maués, P.C.R. Efeitos antropogênicos sobre a diversidade de mamíferos de médio e grande porte na Amazônia Oriental. PhD. Thesis, Graduate Program in Ecology, Federal University of Pará, Belém, Pará, Brazil (2019).25.Parry, L., Barlow, J. & Peres, C. A. Large-vertebrate assemblages of primary and secondary forests in the Brazilian Amazon. J. Trop. Ecol. 23, 653–662. https://doi.org/10.1017/S0266467407004506 (2007).Article 

Google Scholar 
26.Mendes-Oliveira, A. C. et al. Oil palm monoculture induces drastic erosion of an Amazonian forest mammal fauna. PLoS ONE 12, e0187650. https://doi.org/10.1371/journal.pone.0187650 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
27.Coelho, M., Juen, L. & Mendes-Oliveira, A. C. The role of remnants of Amazon savanna for the conservation of Neotropical mammal communities in eucalyptus plantations. Biodivers. Conserv. 23, 3171–3184. https://doi.org/10.1007/s10531-014-0772-9 (2014).Article 

Google Scholar 
28.Bicknell, J. E., Struebig, M. J. & Davies, Z. G. Reconciling timber extraction with biodiversity conservation in tropical forests using reduced-impact logging. J. Appl. Ecol. 52, 379–388. https://doi.org/10.1111/1365-2664.12391 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
29.Chazdon, R. L. et al. The potential for species conservation in tropical secondary forests. Conserv. Biol. 23, 1406–1417. https://doi.org/10.1111/j.1523-1739.2009.01338.x (2009).Article 
PubMed 

Google Scholar 
30.Koh, L. P. & Wilcove, D. S. Is oil palm agriculture really destroying tropical biodiversity?. Conserv. Lett. 1, 60–64. https://doi.org/10.1111/j.1755-263X.2008.00011.x (2008).Article 

Google Scholar 
31.Putz, F. E. & Pinard, M. A. Reduced-impact logging as a carbon-offset method. Conserv. Biol. 7, 755–757. https://doi.org/10.1046/j.1523-1739.1993.7407551.x (1993).Article 

Google Scholar 
32.Pinard, M. A. & Putz, F. E. Retaining forest biomass by reducing logging damage. Biotropica 28, 278–295. https://doi.org/10.2307/2389193 (1996).Article 

Google Scholar 
33.Prudente, B. S., Pompeu, P. S., Juen, L. & Montag, L. F. A. Effects of reduced-impact logging on physical habitat and fish assemblages in streams of Eastern Amazonia. Freshw. Biol. 62, 303–316. https://doi.org/10.1111/fwb.12868 (2017).Article 

Google Scholar 
34.Kanowski, J., Catterall, C. P. & Wardell-Johnson, G. W. Consequences of broadscale timber plantations for biodiversity in cleared rainforest landscapes of tropical and subtropical Australia. For. Ecol. Manage. 208, 359–372. https://doi.org/10.1016/j.foreco.2005.01.018 (2005).Article 

Google Scholar 
35.Correa, F. S., Juen, L., Rodrigues, L. C., Silva-Filho, H. F. & Santos-Costa, M. C. Effects of oil palm plantations on anuran diversity in the eastern Amazon. Anim. Biol. 65, 321–335. https://doi.org/10.1163/15707563-00002481 (2015).Article 

Google Scholar 
36.Peres, C. A. & Cunha, A. A. Line-Transect Censuses of Large-Bodied Tropical Forest Vertebrates: A Handbook (Wildlife Conservation Society, 2011).
Google Scholar 
37.Chao, A. & Jost, L. Coverage-based rarefaction and extrapolation: Standardizing samples by completeness rather than size. Ecology 93, 2533–2547. https://doi.org/10.1890/11-1952.1 (2012).Article 
PubMed 

Google Scholar 
38.Oksanen, J. F. et al. vegan: Community Ecology Package. R package version 2.5–6. https://CRAN.R-project.org/package=vegan (2019).39.Ceballos, G., Ehrlich, P. R. & Dirzo, R. Biological annihilation via the ongoing sixth mass extinction signaled by vertebrate population losses and declines. Proc. Natl. Acad. Sci. 114, 6089–6096. https://doi.org/10.1073/pnas.1704949114 (2017).CAS 
Article 

Google Scholar 
40.Kricher, J. Tropical Ecology 632 (Princeton University Press, 2011).
Google Scholar 
41.Edwards, D. P. et al. Reduced-impact logging and biodiversity conservation: A case study from Borneo. Ecol. Appl. 22, 561–571. https://doi.org/10.1890/11-1362.1 (2012).Article 
PubMed 

Google Scholar 
42.Melo, F. P. L., Arroyo-Rodríguez, V., Fahrig, L., Martínez-Ramos, M. & Tabarelli, M. On the hope for biodiversity friendly tropical landscapes. Trends Ecol. Evol. 28, 462–468. https://doi.org/10.1016/j.tree.2013.01.001 (2013).Article 
PubMed 

Google Scholar 
43.Benton, T. G., Vickery, J. A. & Wilson, J. D. Farmland biodiversity: Is habitat heterogeneity the key?. Trends Ecol. Evol. 18, 182–188. https://doi.org/10.1016/S0169-5347(03)00011-9 (2003).Article 

Google Scholar 
44.Almeida-Rocha, J. M., Peres, C. A. & Oliveira, L. C. Primate responses to anthropogenic habitat disturbance: A pantropical meta-analysis. Biol. Conserv. 215, 30–38. https://doi.org/10.1016/j.biocon.2017.08.018 (2017).Article 

Google Scholar 
45.Palmeirim, A. F., Vieira, M. V. & Peres, C. A. Herpetofaunal responses to anthropogenic forest habitat modification across the neotropics: Insights from partitioning β-diversity. Biodivers. Conserv. 26, 2877–2891. https://doi.org/10.1007/s10531-017-1394-9 (2017).Article 

Google Scholar 
46.Christie, A. P. et al. Quantifying and addressing the prevalence and bias of study designs in the environmental and social sciences. Nat. Commun. 11, 6377. https://doi.org/10.1038/s41467-020-20142-y (2020).CAS 
Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
47.Whitworth, A. et al. Human disturbance impacts on rainforest mammals are most notable in the canopy, especially for larger-bodied species. Divers. Distrib. 25, 1166–1178. https://doi.org/10.1111/ddi.12930 (2019).Article 

Google Scholar 
48.Johns, A. D. & Skorupa, J. P. Responses of rain-forest primates to habitat disturbance: A review. Int. J. Primatol. 8, 157–191. https://doi.org/10.1007/BF02735162 (1987).Article 

Google Scholar 
49.Wearn, O. R. et al. Mammalian species abundance across a gradient of tropical land-use intensity: A hierarchical multi-species modelling approach. Biol. Conserv. 212, 162–171. https://doi.org/10.1016/j.biocon.2017.05.007 (2017).Article 

Google Scholar 
50.Benchimol, M. & Peres, C. A. Determinants of population persistence and abundance of terrestrial and arboreal vertebrates stranded in tropical forest land-bridge islands. Conserv. Biol. 35(3), 870–883. https://doi.org/10.1111/cobi.13619 (2020).Article 
PubMed 

Google Scholar 
51.Gittleman, J. L. & Harvey, P. H. Carnivore home-range size, metabolic needs and Ecology. Behav. Ecol. Sociobiol. 10(1), 57–63. https://doi.org/10.1007/BF00296396 (1982).Article 

Google Scholar 
52.Edwards, D. P., Tobias, J. A., Sheil, D., Meijaard, E. & Laurance, W. F. Maintaining ecosystem function and services in logged tropical forests. Trends Ecol. Evol. 29, 511–520. https://doi.org/10.1016/j.tree.2014.07.003 (2014).Article 
PubMed 

Google Scholar 
53.Mollinari, M. M., Peres, C. A. & Edwards, D. P. Rapid recovery of thermal environment after selective logging in the Amazon. Agric. For. Meteorol. 278, 107637. https://doi.org/10.1016/j.agrformet.2019.107637 (2019).Article 
ADS 

Google Scholar 
54.Azevedo-Ramos, C., de Carvalho, O. & de Amaral, B. D. Short-term effects of reduced-impact logging on eastern Amazon fauna. For. Ecol. Manag. 232, 26–35. https://doi.org/10.1016/j.foreco.2006.05.025 (2006).Article 

Google Scholar 
55.Bicknell, J. E. & Peres, C. A. Vertebrate population responses to reduced-impact logging in a neotropical forest. For. Ecol. Manage. 259, 2267–2275. https://doi.org/10.1016/j.foreco.2010.02.027 (2010).Article 

Google Scholar 
56.Laufer, J., Michalski, F. & Peres, C. A. Effects of reduced-impact logging on medium and large-bodied forest vertebrates in eastern Amazonia. Biota Neotrop. 15, e20140131. https://doi.org/10.1590/1676-06032015013114 (2015).Article 

Google Scholar 
57.Carvalho Jr, E. A. R., Mendonça, E. N., Martins, A. & Haugaasen, T. Effects of illegal logging on Amazonian medium and large-sized terrestrial vertebrates. For. Ecol. Manage. 466, 118105. https://doi.org/10.1016/j.foreco.2020.118105 (2020).Article 

Google Scholar 
58.Kuussaari, M. et al. Extinction debt: A challenge for biodiversity conservation. Trends Ecol. Evol. 24, 564–571. https://doi.org/10.1016/j.tree.2009.04.011 (2009).Article 
PubMed 

Google Scholar 
59.Richardson, V. A. & Peres, C. A. Temporal decay in timber species composition and value in Amazonian logging concessions. PLoS ONE 11, e0159035. https://doi.org/10.1371/journal.pone.0159035 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
60.Chazdon, R. L. Second Growth: The Promise of Tropical Forest Regeneration in an Age of Deforestation (University of Chicago Press, 2014).Book 

Google Scholar 
61.Acevedo-Charry, O. & Aide, T. M. Recovery of amphibian, reptile, bird and mammal diversity during secondary forest succession in the tropics. Oikos 128, 1065–1078. https://doi.org/10.1111/oik.06252 (2019).Article 

Google Scholar 
62.Sodhi, N. S. et al. Conserving Southeast Asian forest biodiversity in human-modified landscapes. Biol. Conserv. 143, 2375–2384. https://doi.org/10.1016/j.biocon.2009.12.029 (2010).Article 

Google Scholar 
63.Dunn, R. R. Recovery of faunal communities during tropical forest regeneration. Conserv. Biol. 18, 302–309. https://doi.org/10.1111/J.1523-1739.2004.00151.X (2004).Article 

Google Scholar 
64.Luskin, M. S. & Potts, M. D. Microclimate and habitat heterogeneity through the oil palm lifecycle. Basic Appl. Ecol. 12, 540–551. https://doi.org/10.1016/j.baae.2011.06.004 (2011).Article 

Google Scholar 
65.Fitzherbert, E. B. et al. How will oil palm expansion affect biodiversity?. Trends Ecol. Evol. 23(10), 538–545. https://doi.org/10.1016/j.tree.2008.06.012 (2008).Article 
PubMed 

Google Scholar 
66.Martello, F. et al. Homogenization and impoverishment of taxonomic and functional diversity of ants in Eucalyptus plantations. Sci. Rep. 8, 3266. https://doi.org/10.1038/s41598-018-20823-1 (2018).CAS 
Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
67.da Rocha, P. L. B. What is the value of eucalyptus monocultures for the biodiversity of the Atlantic Forest? A multitaxa study in southern Bahia, Brazil. J. For. Res. 24, 263–272. https://doi.org/10.1007/s11676-012-0311-z (2013).Article 

Google Scholar 
68.Martin, P. S., Gheler-Costa, C., Lopes, P. C., Rosalino, L. M. & Verdade, L. M. Terrestrial non-volant small mammals in agro-silvicultural landscapes of Southeastern Brazil. For. Ecol. Manag. 282, 185–195. https://doi.org/10.1016/j.foreco.2012.07.002 (2012).Article 

Google Scholar 
69.Fayle, T. M. et al. Oil palm expansion into rain forest greatly reduces ant biodiversity in canopy, epiphytes and leaf-litter. Basic Appl. Ecol. 11, 337–345. https://doi.org/10.1016/j.baae.2009.12.009 (2010).Article 

Google Scholar 
70.Koh, L. P. Can oil palm plantations be made more hospitable for forest butterflies and birds?. J. Appl. Ecol. 45, 1002–1009. https://doi.org/10.1007/s10531-009-9760-x (2008).Article 

Google Scholar 
71.Martins, C. A. & Júnior, A. P. P. Production of biodiesel: Source strategies and efficiency in the Brazilian energy matrix. Energy Sour. Part A Recov. Util. Environ. Eff. 38, 277–285. https://doi.org/10.1080/15567036.2012.716139 (2016).CAS 
Article 

Google Scholar 
72.Peres, C. A. Why we need megareserves in Amazonia. Cons. Biol. 19, 728–733. https://doi.org/10.1111/j.1523-1739.2005.00691.x (2005).Article 

Google Scholar  More

Sampling site and sediment collectionPaleolimnological analysis by using sediment core samples were applied to reconstruct historical variations of Daphnia eDNA concentrations in Lake Biwa (Fig. 1b). Lake Biwa is the largest monomictic and mesotrophic lake in Japan. In this lake, during the last several decades the industrial revolution, multiple stressors of human origins impacted this ecosystem and the resident biological communities34,35,58. In our study, four 30-cm-long gravity core samples (namely LB1, LB2, LB4, and LB7; Fig. 1a–e) were collected on 17 August 2017 at the anchoring site of Hasu, a Center of Ecological Research boat from Kyoto University (Supplementary Fig. S2). A gravity corer with an inner diameter of 10.9 cm and a length of 30 cm was used to obtain the core samples. LB7 core was analyzed for chronological and reconstruction of temporal variation in Daphnia remain abundance (Fig. 1a,b). LB2 and LB1, LB4 cores were analyzed for reconstruction of temporal variation in sedimentary Daphnia DNA concentrations and resting egg production, respectively (Fig. 1b). Additionally, two 30-cm-long gravity core samples (namely IM1 and IM8; Fig. 1c,f) were collected at a pelagic site in the northern basin of Lake Biwa in August 2019 (Supplementary Fig. S2). The collected cores were sectioned at intervals of 1-cm thickness using a vertical extruder with a cutting apparatus, except for core number IM8, which was sectioned at 5-cm intervals (Fig. 1f). During the sectioning process, several millimeters of outer edge in each layer disturbed during the splitting process were carefully removed from the entire samples using a knife. After sectioning, each sliced sample were homogenized by shaking and then, all subsamples were taken from each homogenized sample. The pipes, knives, and cutting apparatus were cleaned with 0.6% sodium hypochlorite, tap water, and Milli-Q water to avoid DNA cross-contamination. Each sliced sample was transferred to lightproof bags and frozen at − 80 °C until further analysis.To examine the contamination due to core splitting, DNA extraction, and qPCR analysis, control water samples were inserted at a depth of 14.5–29.5 cm in the sediment cores, and the water samples for core IM1 were used as the negative control (Fig. 1c).Chronology of sediment coresSediment chronology was performed for the LB7 core based on the constant rate of supply (CRS) method of 210Pb dating59 and verified using the 137Cs peak traced in the period 1962 to 196360. Details of the chronological method have been reported elsewhere61. Briefly, dried samples were sealed in holders for a month to allow 222Rn and its short-lived decay product (214Pb) to equilibrate, which were determined by gamma counting using a germanium detector (GXM25P; EG & G ORTEC, Tokyo, Japan) equipped with a multi-channel analyzer (MCA7700; SEIKO EG & G, Tokyo, Japan) at the Center for Marine Environmental Studies, Ehime University. The activity of supported 210Pb was estimated by measuring the activity of 214Pb, whereas that of 210Pbexcess was determined according to the difference between the total and the supported 210Pb (210Pbexcess = 210Pbtotal − 214Pb). The age and age error of the remaining cores (LB1, LB2, and LB4) were indirectly estimated using stratigraphic correlations between the cores based on chronological controls in chlorophyll pigments and magnetic susceptibilities of the chronological LB7 core61. To compare these proxies, the marked peak or trough layers were used as reference layers (Supplementary Fig. S3).DNA extraction and purificationDNA extraction in the sediment samples was performed according to methods described in previous studies45,62. In brief, 9 g of each sediment sample was incubated at 94 °C for 50 min in a 9 mL alkaline solution comprising 6 mL of 0.33 M sodium hydroxide and a 3 mL Tris–EDTA buffer (pH 6.7). After centrifugation at 10,000×g for 60 min, 7.5 mL of the supernatant of the alkalized mixture was neutralized with 7.5 mL of 1 M Tris–HCl (pH 6.7). After adding 1.5 mL of 3 M sodium acetate (pH 5.2) and 30 mL absolute ethanol, the solution was preserved at − 20 °C for more than 1 h and then centrifuged at 10,000×g for 60 min. The pellet was transferred into a power bead tube that was installed in a fecal-soil DNA extraction kit (Power Soil DNA Isolation Kit, Qiagen, Germany). The ‘Experienced User Protocol 3 to 22’ of the Power Soil DNA Isolation Kit was followed. Finally, 200 μL of the DNA solution was obtained and stored at − 20 °C until qPCR analysis.12S rRNA gene primer-prove development for Daphnia geleata and Daphnia pulicaria
As the primer–probe for Daphnia galeata and D. pulicaria in qPCR analysis were not purchased by a company, thus we developed them for the two Daphnia species (see Supplementary Table S1). We preliminary obtained the mitochondrial 12S, 16S and COI gene of Daphnia genus from the National Center for Biotechnology Information (NCBI, http://www.ncbi.nlm.nih.gov/) and compared among them. From the preliminary results, we decided to use 12S because of the variability of sequences among Daphnia genus. Then we obtained the 12S sequences of Daphnia genus and other inhabiting plankton species in Lake Biwa, including Copepoda. We designed the primer–probe using Primer3plus (https://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi). The reference sequences for the targeted gene regions are queried for potential amplicons between 50–150 bp using NCBI primer blast. The specificities of the primers and probes were then assessed in silico with homologous sequences from other Daphnia species in Japan using NCBI targeting 154 bp of the mitochondrial 12S rRNA gene. Once suitable amplicons are found the respective primers and probes are tested against template DNA originating from the species of D. galeata and D. pulicaria to verify amplification. During the in silico screening for specificity, we performed Primer-BLAST (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). We checked all species from Japan of the order Daphnia. Using the D. galeata primer-set, we did not detect any Daphnia species. However, the D. pulicaria primer can amplify D. pulex DNA, as these species are known to have very similar sequences47. In Lake Biwa, another subgenus Daphnia (D. pulex group) different from D. pulicaria was temporally found around during the 1920s, although thereafter it was never reported47. Thus, the D. pulicaria primer may temporally detect another subgenus Daphnia (D. pulex group). However, their appearance time do not overlap, therefore we used the primer for our measurement to detect D. pulicaria during the last several decades.Quantitative PCRThe DNA samples were quantified by real-time TaqMan quantitative PCR using the PikoReal Real-Time PCR system (Thermo Fisher Scientific, Waltham, MA, USA). The primer–probe sets for the two Daphnia species were used for qPCR (Supplementary Table S1). The TaqMan reaction contained 900 nM of each forward and reverse primer, 125 nM TaqMan-Probe, 5 μL qPCR master mix (TaqPath; Thermo Fisher Scientific), and 2.0 μL sedimentary DNA solution. The final volume of the PCR was 10 μL after adding distilled water (DW). The qPCR conditions were as follows: 50 °C for 2 min and 95 °C for 10 min, followed by 50 cycles of 95 °C for 15 s and 60 °C for 60 s. We used a dilution series of 10,000, 1000, 100, and 10 copies per PCR reaction (n = 4) for the standard curve using the target DNA cloned into a plasmid. The R2 values of the standard curves ranged from 0.988 to 0.996 (PCR efficiencies = 93.1–102.0%). The quantitative data of the DNA copies (copies g−1 dry sed.) were reported by mean values ± standard deviation, which were calculated from DNA copies µL−1 PCR reaction with four replicates including zero (i.e., no detection). We also performed four replicates for each sample and an NTC (n = 4). No positives were detected from the NTC and the negative control of DNA extraction, confirming that there was no cross-contamination in any of the DNA measurements.To confirm primer specificity, an in vivo test for the primer/probe set was performed using the extracted DNA (10 pg per PCR reaction, n = 4) of D. galeata and D. pulicaria. In addition, qPCR amplicons were sequenced directly from a positive PCR from each site (n = 21) after treatment with ExoSAP-IT (USB Corporation, Cleveland, OH, USA). Sequences were determined using a commercial sequencing service (Eurofins Genomics, Tokyo, Japan).Inhibitor testSpike tests were performed for the LB2 core sample to evaluate the PCR inhibition effect of several substances and minerals in the sediment samples (Fig. 1c). For the spike test, 1 µL plasmid, including the internal positive control (IPC, 207-bp, Nippon Gene Co. Ltd., Tokyo, Japan;100 copies per PCR reaction), was added to the PCR template with 1.6 µL DNA-free DW. We used the primer and probe sets for IPC as follows:

IPC1-5′: CCGAGCTTACAAGGCAGGTT

IPC1-3′: TGGCTCGTACACCAGCATACTAG

IPC1-Taq: [FAM] TAGCTTCAAGCATCTGGCTGTCGGC [TAMRA].

To measure the relative degree of PCR inhibition in the samples, the Ct shift was compared between the samples and controls with the same number of known target DNA copies. The presence of PCR inhibitors was evaluated as ΔCt = Ct sample − Ct positive control. ΔCt ≥ 3 cycles was considered evidence of inhibition63 because the presence of PCR inhibitors will delay the Ct with a given quantity of template DNA.
Daphnia abundance and resting egg production as potential sources of Daphnia DNA archived in sedimentsTo unveil the potential source of sedimentary DNA of Daphnia, we reconstructed the historical variation in Daphnia abundance by counting remains of the post abdominal claw for LB7 core. There are two dominant Daphnia species: D. galeata Sars (Hyalodaphnia) and D. pulicaria Forbes (Daphnia)47,61, which have different post-abdominal claw characteristics64 and are known to be preserved in centuries-old sediments65. The post-abdominal claw remains were counted for core LB7 from the surface to a depth layer of 21.5 cm and additionally 23.5 cm, 25.5 cm, and 29.5 cm, totaling 25 samples, though each layer was expressed as mid-depth; e.g., 0.5 cm for the 0–1 cm depth layer. The enumeration method was based on a simplified standard method65 as previously reported29.Daphnia resting eggs enveloped by thickened carapaces, referred to as ephippial cases, and these ephippia can be preserved in sediments for decades to centuries29,30,33. In Lake Biwa, Daphnia species in Lake Biwa are distinguished on the basis of the size of the ephippium, with a boundary length of approximately 860 μm between them61. We collected ephippia from the surface to a depth layer of 29.5 cm for cores LB1 and LB4 (except for several layers of the LB1 core), totaling 56 samples (Supplementary Table S4). A detailed method for collecting ephippia is described in a previous study61. The total number of collected ephippia with an almost perfect shape, namely complete formation, or with a partial body constituting more than half of the original shape, namely incomplete formation, are shown in Supplementary Table S4. In our study, at least 16 ephippia in each sample were measured by photographs taken by a digital camera, excluding those from the samples in which fewer than 16 complete ephippia were detected (Supplementary Fig. S4). Species identification was then performed based on length.To determine whether the Daphnia sedimentary DNA concentrations were regulated by DNA derived directly from Daphnia remains or ephippia included in the analytical sediment, we divided the sediment sample into two fractions to exclude the remains and ephippia (Fig. 1d). The minimum size of Daphnia remains in this lake was approximately 55 μm (Tsugeki et al., in preparation). The analytical sediments for DNA extraction were divided into particles  38 μm using 38-μm mesh sieves on three-layer samples (specifically, LB2-5; 4.5 cm, LB2-7; 6.5 cm, and LB2-17; 16.5 cm expressed in middle depth of each sample) for core sample LB2, whose layers were known to include abundant ephippia and Daphnia remains. Furthermore, to test the possibility of the vertical movement of Daphnia sedimentary DNA through pore waters, we examined the sedimentary DNA concentration in pore water and its residual sediment by qPCR analysis (Fig. 1e). All DNA extractions were evaluated for sediment with and without sieves, and pore waters and the associated residual sediment samples were evaluated according to previous studies45.Measurement of DNA concentration in sediment ephippiaTo determine the potential source of sedimentary Daphnia DNA, we quantified the DNA concentration extracted from several ephippia obtained from the 0–5 cm and 5–10 cm layers of core IM8 using qPCR analysis (Fig. 1f). We selected 34 and 23 ephippia for D. galeata and D. pulicaria, respectively. We then measured the ephippial lengths and determined whether they contained resting eggs using a microscope. Among the selected ephippia, the well-preserved 17 ephippia with almost complete formation were set aside and grouped into 6 samples together in two or three ephippia for DNA analysis (Supplementary Table S5). Grouping was performed because of the low DNA concentrations typically associated with individual ephippium61.Possible factors regulating sedimentary Daphnia DNATo explore potential factors regulating temporal variation in sedimentary DNA concentrations, we analyzed chlorophyll pigments and algal remains. Sedimentary pigments of chlorophyll a were investigated for the LB 2 core, and algal remains were investigated for the LB7 core (Fig. 1a). Details of the method used for chlorophyll-a and algal remains are described in previous study61. In short, the concentrations of chlorophyll-a and phaeopigments were calculated according to the method66 and the diatom remains were analyzed according to the simplified method67. Green algae, Micrasterias hardyi, Staurastrum dorsidentiferum, S. arctiscon, S. limneticum, S. pingue, and Pediastrum biwae, were enumerated in a Sedgewick–Rafter chamber, following the method of zooplankton enumeration.Data analysisRegression models along with the standardized major axis method were used to determine the relationship between the sedimentary DNA concentration obtained from qPCR analysis and abundance or resting egg production in the sediment layers. Since qPCR (LB2), remains (LB7), and ephippia (LB1, LB4) analyses were performed on different cores, the chronological age of each analytical sample differed slightly. Therefore, prior to performing the statistical analysis, the sedimentary DNA (LB2) and ephippia data (LB1, LB4) in each chronological age were converted to annual data by linear interpolation and averaged for the year corresponding to the period in each sample of the chronology core (LB7). This conversion was possible because the time resolution at 1-cm intervals represented several years, depending on the sediment depth29,61. We employed the Gaussian type II model because our preliminary evaluation showed higher R2 values for type II regression models with a Gaussian distribution than for those with a logarithmic distribution, in all cases. All statistical analyses were performed using R ver. 4.0.3 (R Core Team 2020) with the package “smatr” ver. 3.4-8 for type II regressions. The significant criteria of all analyses were set as α = 0.05. In addition, to explore the potential environmental factors driving temporal variation in sedimentary DNA concentrations, we performed Pearson’s correlation analysis among the sedimentary DNA concentrations, chlorophyll a concentration, and algal remains using the SPSS version 20.0 statistical package. More

Changepoint detection methodAlthough most of the reflectance time series used in the BinSeg–Normal–Mean and BinSeg–Normal–MeanVar algorithms had a normal distribution, several lagoons had distributions that were skewed or did not follow a normal distribution (Fig. S1). However, results suggested that the accuracy of the detected changepoints were not sensitive to the normality assumption or distributional characteristics.The BinSeg-Normal-Mean algorithm had the highest performance (81% of the 340 validation sites) in detecting the correct year of swine waste lagoon construction, followed by BinSeg-Normal-MeanVar (77%). The two algorithms did not detect the same year of construction for 19 waste lagoons; of these 19, the BinSeg-Normal-Mean detected the correct year for 84% of them, while the BinSeg-Normal-MeanVar detected the correct year for only 16%. Therefore, the BinSeg-Normal-MeanVar algorithm was abandoned given it did not provide additional useful information relative to the BinSeg-Normal-Mean algorithm.Despite good performance, the BinSeg-Normal-Mean algorithm consistently detected a changepoint during the period of record for all sites included in the 10% validation set (n = 340 swine waste lagoons). However, 58 of the 340 swine waste lagoons were constructed prior to 1986, before the period of record suitable for detecting an accurate changepoint. Changepoints before 1986 either (1) detected the correct construction year, or (2) incorrectly detected a changepoint due to artifact signals identified on the images taken in 1984, probably associated with the initial satellite commissioning. In the latter circumstance, if the algorithm detected a changepoint due to this signal, it meant that no land-use change was detected after 1986. Therefore, these waste lagoons were estimated as having been constructed before 1986. In some conditions, when a large number of images was available for the year 1985 and 1986, the algorithm was able to detect the changepoint occurring for the years 1985 or 1986. Further, the BinSeg-Normal-Mean algorithm detected a false year of construction for swine waste lagoons for which the mean of the segment after the changepoint (S2) had a greater average than the segment before the changepoint (S1).To increase algorithm performance, we developed a workflow to address some of the aforementioned caveats (Fig. 4). In this workflow, the BinSeg-Normal-Mean algorithm is applied to a B4 reflectance time series at location j. If the BinSeg-Normal-Mean changepoint is identified for a time in or prior to 1986 (Fig. 4a,i,b,i) we assume that the lagoon was constructed in or prior to 1986. Similarly, a lagoon is assumed to be constructed in or prior to 1986 if a BinSeg-Normal-Mean changepoint is identified after 1986 and the mean of S2 is greater than the mean of S1 (Fig. 4a,ii,b,ii). If a changepoint occurred after 1986 and the mean of S1 was greater than S2, then the changepoint was estimated as having occurred between 1987 and 2010 (Fig. 4a,iii,b,iii).Figure 4Changepoint detection algorithm for determining the year of construction of swine waste lagoons. Panel (a) summarizes the algorithm workflow, while panel (b) illustrates specific examples corresponding to each step (i–iii) in the workflow.Full size imageThe performance of the workflow was evaluated using the validation set composed of 10% of the total number of swine waste lagoons (n = 340). With the new approach, 94% of the swine waste lagoon construction years (+ /- one year) were accurately retrieved. A tolerance of + /− 1 year was chosen to account for a lack of images in some years due to issues with image quality (e.g. high cloud cover) (e.g., Fig. 5a), or because construction spanned at least a year (e.g., Fig. 5b). The changepoint detection workflow incorrectly estimated the construction years for 19 of the 340 swine waste lagoons in the validation set; the differences between the observed and predicted years of construction of these lagoons ranged from 2 to 26 years with a median of 8 years.Figure 5Examples of limitations to the changepoint detection algorithm. In some cases, an insufficient number of high-quality Landsat 5 images were available to capture the year of construction of an individual swine waste lagoon (a), resulting in errors of + /− 1 year. In other cases, the changepoint algorithms detected the start of the construction of the swine waste lagoon but the swine waste lagoon was not fully operational until later years due to prolonged construction timelines (b).Full size imageBy visually inspecting historical Google Earth images for each of the lagoon sites for which the model incorrectly estimated construction year, we identified that model errors were associated with swine waste lagoon expansion, pixel transitions to land-use classes other than swine waste lagoons, or issues with pixels being partly covered by clouds or incompletely covered by the lagoon (i.e., narrow and small waste lagoons that do not entirely cover a pixel).Estimating swine waste lagoon construction yearsUsing the newly developed algorithm (Fig. 4), construction years were estimated for each swine waste lagoon in the NC Coastal Plain (Fig. 6); the years of construction for each swine waste lagoon are included in the supplementary material. Most swine waste lagoons were built in the early 90s and prior to the moratorium of 1997. More specifically, 80% of the swine waste lagoons (n = 2,736) were built between 1987 and 1997. Sixteen percent of the swine waste lagoons were constructed in or prior to 1986. A large decrease in the construction of swine waste lagoons occurred after the moratorium of 1997, with only 3.7% of swine waste lagoons being constructed after the moratorium. These results suggest that the 1997 moratorium did not completely halt the construction of lagoons, but dramatically slowed the rate of expansion.Figure 6Spatiotemporal distribution of swine waste lagoon construction (+/- 1 year) across the HUC6 watersheds. This figure was produced using QGIS version QGIS 3.18.3 (https://www.qgis.org/).Full size imageWith regards to hydrological boundaries (Fig. 7a–h), the Cape Fear River watershed had the highest number of swine waste lagoons (i.e., 56%; Fig. 7b), followed by the Neuse River (i.e., 23%; Fig. 7d), the Lower Pee Dee River (i.e., 9%; Fig. 7c) watersheds. The Albemarle-Chowan (Fig. 7a), Onslow Bay (Fig. 7e), Pamlico (Fig. 7f), Roanoke (Fig. 7g), and Upper Pee Dee (Fig. 7h) watersheds all had less than 9% of the total lagoons within the study area.Figure 7Year of construction of the swine waste lagoons (+ /− 1 year) for the HUC6 watersheds. The y-axis scales are unequal between the plots to improve readability. The dashed red lines correspond to the establishment of the moratorium in 1997.Full size imageResults suggested that the Cape Fear River watershed was the center of the historical growth of the swine industry, where over 300 swine waste lagoons were built prior to 1987. The Cape Fear River watershed experienced a steady increase in the number of swine waste lagoons from 1987 to 1990, with an average of 46 swine waste lagoons being built annually. However, after 1991, the pace of swine waste lagoon construction increased dramatically with an average of 192 swine waste lagoons built annually between 1991 and 1997. The highest construction rate occurred in 1994, with 242 swine waste lagoons built. However, after the 1997 moratorium, the construction rate decreased dramatically; in 1997, 153 swine waste lagoons were constructed, and this number dropped to 23 in 1998. After 1998, the annual average number of swine waste lagoons constructed plunged to 5. Although the swine waste lagoon construction rate fell considerably after the 1997 moratorium, the decrease had already started in 1995. The same pattern was observed for the Neuse, Pamlico, Albemarle-Pamlico, and Onslow Bay watersheds.The spatiotemporal distribution of swine waste lagoons at the HUC12 watershed scale emphasized the historical clustering of the swine industry in the NC Coastal Plain. After the moratorium, swine waste lagoons were present within 436 HUC12 watersheds. However, before 1986, they were spread across only 197 HUC12 watersheds (Fig. 8). Before 1986, the density of waste lagoons was relatively low with an average of 3.38 swine waste lagoons per 100 km2 and a maximum of 15.13 swine waste lagoons per 100 km2 (i.e., Clayroot Swamp-Swift Creek watershed) (Fig. 8). In the 90s, swine waste lagoon construction expanded and continued to intensify in the region. After the moratorium of 1997, the average density of waste lagoons per HUC12 watersheds was 10 per 100 km2 with a maximum of 78 waste lagoons per 100 km2 identified in the Maxwell Creek-Stocking Head Creek basin. After 1997, 16 of 436 HUC12 watersheds had a swine waste lagoon density greater than 40 per 100 km2 (Fig. 8).Figure 8Cumulative swine waste lagoon density per 100 km2 reported at the HUC12 watershed scale; HUC6 watersheds shown in gray for reference. This figure was produced using QGIS version QGIS 3.18.3 (https://www.qgis.org/).Full size imageSpatiotemporal distribution of swine waste lagoons in relation to water resourcesDistance of swine waste lagoon sites to the nearest water feature (i.e., reservoir, canal/ditch, lake/pond, stream/river, estuary) were assessed using the NHD. The analysis revealed that over 150 swine waste lagoons were misclassified by the NHD and were documented in the NHD as lake/pond (n = 102) or swamp/marsh (n = 46). Further, we observed that some NHD water features were misclassified as other non-water features (e.g., forest, pasture), and most of these misclassifications were for polygons with an area less than 0.05 km2. Therefore, NHD water features with areas less than 0.05 km2 were removed from subsequent analyses. Distances between swine waste lagoons and waterways were computed from the NHD without features with areas less than 0.05 km2. The new analysis revealed that 3 swine waste lagoons remained misclassified as lake/pond (n = 1) and swamp/marsh (n = 2). Canal/Ditch, lake/pond, stream/river, and swamp/marsh were identified as the NHD features that were most commonly near swine waste lagoons (Fig. 9). Two swine waste lagoons were near a reservoir in which one was identified as a treatment-sewage pond by the NHD.Figure 9Nearest water features distance to swine waste lagoons.Full size imageThe average and median distance of all swine waste lagoons (including those built early and late in the period of record) to the nearest water features were 234 and 177 m, respectively. Further, 92% of the swine waste lagoons were less than 500 m from the nearest waterways. The Mann–Kendall results revealed a significant upward trend over time of swine waste lagoon distances to the nearest water features (alpha = 0.05, p-value = 0.01). A slight increase over time of swine waste lagoon distances to the nearest water feature is also documented in Table 1.Table 1 Temporal average and median of nearest distance (m) of swine waste lagoons to water features. NA indicated that the water feature was not the closest waterway to any of the studied swine waste lagoons for the time period.Full size table More

1.MacLachlan, N. J. & Guthrie, A. J. Re-emergence of bluetongue, African horse sickness, and other Orbivirus diseases. Vet. Res. 41, 35 (2010).Article 

Google Scholar 
2.Zientara, S., Weyer, C. T. & Lecollinet, S. African horse sickness. OIE Revue Sci. Tech. 34, 315–327 (2015).CAS 
Article 

Google Scholar 
3.Ayelet, G. et al. Outbreak investigation and molecular characterization of African horse sickness virus circulating in selected areas of Ethiopia. Acta Trop. 127, 91–96 (2013).Article 

Google Scholar 
4.Diarra, M. et al. Spatial distribution modelling of Culicoides (Diptera: Ceratopogonidae) biting midges, potential vectors of African horse sickness and bluetongue viruses in Senegal. Parasit. Vectors 11, 1–15 (2018).Article 

Google Scholar 
5.Karamalla, S. T. et al. Sero-epidemioloical survey on African horse sickness virus among horses in Khartoum State, Central Sudan. BMC Vet. Res. 14, 1–6 (2018).Article 

Google Scholar 
6.Escobar, L. E. Ecological Niche modeling: An introduction for veterinarians and epidemiologists. Front. Vet. Sci. 7, 519059. https://doi.org/10.3389/fvets.2020.519059 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
7.Okely, M., Anan, R., Gad-Allah, S. & Samy, A. M. Mapping the environmental suitability of etiological agent and tick vectors of Crimean-Congo hemorrhagic fever. Acta Trop. 203, 105319 (2020).CAS 
Article 

Google Scholar 
8.Chavy, A. et al. Ecological niche modelling for predicting the risk of cutaneous leishmaniasis in the Neotropical moist forest biome. PLoS Negl. Trop. Diseases 13, e0007629 (2019).Article 

Google Scholar 
9.Sloyer, K. E. et al. Ecological niche modeling the potential geographic distribution of four Culicoides species of veterinary significance in Florida, USA. PLoS ONE 14, e0206648 (2019).CAS 
Article 

Google Scholar 
10.Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
11.Cao, Z., Jin, Y., Shen, T., Xu, F. & Li, Y. Risk factors and distribution for peste des petits ruminants (PPR) in Mainland China. Small Rumin. Res. 162, 12–16 (2018).Article 

Google Scholar 
12.Naimi, B. & Araújo, M. B. sdm: a reproducible and extensible R platform for species distribution modelling. Ecography 39, 368–375 (2016).Article 

Google Scholar 
13.Naimi, B., Hamm, N. A. S., Groen, T. A., Skidmore, A. K. & Toxopeus, A. G. Where is positional uncertainty a problem for species distribution modelling. undefined 37, 191–203 (2014).
Google Scholar 
14.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, (2020).15.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 
16.Araújo, M. B. & New, M. Ensemble forecasting of species distributions. Trends Ecol. Evol. 22, 42–47 (2007).Article 

Google Scholar 
17.Uusitalo, R. et al. Predicting spatial patterns of sindbis virus (Sinv) infection risk in finland using vector, host and environmental data. Int. J. Environ. Res. Public Health 18, 7064 (2021).Article 

Google Scholar 
18.Raffini, F. et al. From nucleotides to satellite imagery: Approaches to identify and manage the invasive pathogen Xylella fastidiosa and its insect vectors in Europe. Sustainability (Switzerland) 12, 4508 (2020).CAS 
Article 

Google Scholar 
19.Phillips, S. B., Aneja, V. P., Kang, D. & Arya, S. P. Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231–259 (2006).Article 

Google Scholar 
20.Barbet-Massin, M., Jiguet, F., Albert, C. H. & Thuiller, W. Selecting pseudo-absences for species distribution models: how, where and how many?. Methods Ecol. Evol. 3, 327–338 (2012).Article 

Google Scholar 
21.Hernández-Urcera, J., Murillo, F. J., Regueira, M., Cabanellas-Reboredo, M. & Planas, M. Preferential habitats prediction in syngnathids using species distribution models. Marine Environ. Res. 172, 105488 (2021).Article 

Google Scholar 
22.Smeraldo, S. et al. Generalists yet different: distributional responses to climate change may vary in opportunistic bat species sharing similar ecological traits. Mammal Rev. 51, 571–584 (2021).Article 

Google Scholar 
23.Thomson, A. M. et al. RCP4.5: A pathway for stabilization of radiative forcing by 2100. Clim. Change 109, 77–94 (2011).ADS 
CAS 
Article 

Google Scholar 
24.QGIS Development Team. QGIS Geographic Information System. Open-Source Geospatial Foundation Project. (2020).25.Ramirez-Reyes, C. et al. Embracing ensemble species distribution models to inform at-risk species status assessments. J. Fish Wildl. Manag. 12, 98–111 (2021).Article 

Google Scholar 
26.Stephenson, F. et al. Presence-only habitat suitability models for vulnerable marine ecosystem indicator taxa in the South Pacific have reached their predictive limit. ICES J. Mar. Sci. 78, 2830–2843 (2021).Article 

Google Scholar 
27.Zhu, G., Fan, J. & Peterson, A. T. Cautions in weighting individual ecological niche models in ensemble forecasting. Ecol. Modelling 448, 109502 (2021).Article 

Google Scholar 
28.Leta, S. et al. Modeling the global distribution of Culicoides imicola: an Ensemble approach. Sci. Rep. 9, 1–9 (2019).ADS 
CAS 
Article 

Google Scholar 
29.Onyango, M. G. et al. Delineation of the population genetic structure of Culicoides imicola in East and South Africa. Parasit. Vectors 8, 660 (2015).Article 

Google Scholar 
30.Carpenter, S., Mellor, P. S., Fall, A. G., Garros, C. & Venter, G. J. African horse sickness virus: history. Transm. Curr. Status. 62, 343–358. https://doi.org/10.1146/annurev-ento-031616-035010 (2017).CAS 
Article 

Google Scholar 
31.Carpenter, S., Mellor, P. S., Fall, A. G., Garros, C. & Venter, G. J. African Horse Sickness Virus: History, Transmission, and Current Status. Annu. Rev. Entomol. 62, 343–358 (2017).CAS 
Article 

Google Scholar 
32.Fall, M. et al. Culicoides (Diptera: Ceratopogonidae) midges, the vectors of African horse sickness virus—a host/vector contact study in the Niayes area of Senegal. Parasit. Vectors 8, 1–13 (2015).Article 

Google Scholar 
33.Mellor, P. S. Epizootiology and vectors of African horse sickness virus. Comp. Immunol. Microbiol. Infect. Dis. 17, 287–296 (1994).CAS 
Article 

Google Scholar 
34.Wu, X., Lu, Y., Zhou, S., Chen, L. & Xu, B. Impact of climate change on human infectious diseases: Empirical evidence and human adaptation. Environ. Int. 86, 14–23 (2016).Article 

Google Scholar 
35.Nosrat, C. et al. Impact of recent climate extremes on mosquito-borne disease transmission in Kenya. PLOS Negl. Trop. Diseases 15, e0009182 (2021).CAS 
Article 

Google Scholar 
36.Abiodun, G. J., Maharaj, R., Witbooi, P. & Okosun, K. O. Modelling the influence of temperature and rainfall on the population dynamics of Anopheles arabiensis. Malar. J. 15, 1–15 (2016).Article 

Google Scholar  More