Philippa Kaur

1.Arora, V. K. & Boer, G. J. A parameterization of leaf phenology for the terrestrial ecosystem component of climate models. Glob. Change Biol. 11, 39–59 (2005).
Google Scholar 
2.Richardson, A. D. et al. Terrestrial biosphere models need better representation of vegetation phenology: results from the North American Carbon Program Site Synthesis. Glob. Change Biol. 18, 566–584 (2012).
Google Scholar 
3.Peñuelas, J., Rutishauser, T. & Filella, I. Phenology feedbacks on climate change. Science 324, 887–888 (2009).
Google Scholar 
4.Richardson, A. D. et al. Influence of spring and autumn phenological transitions on forest ecosystem productivity. Philos. Trans. R. Soc. Lond. B Biol. Sci. 365, 3227–3246 (2010).
Google Scholar 
5.Diez, J. M. et al. Forecasting phenology: from species variability to community patterns. Ecol. Lett. 15, 545–553 (2012).
Google Scholar 
6.Hegland, S. J., Nielsen, A., Lazaro, A., Bjerknes, A. L. & Totland, O. How does climate warming affect plant-pollinator interactions? Ecol. Lett. 12, 184–195 (2009).
Google Scholar 
7.Fu, Y. H. et al. Declining global warming effects on the phenology of spring leaf unfolding. Nature 526, 104–107 (2015).CAS 

Google Scholar 
8.Zhang, H., Yuan, W., Liu, S. & Dong, W. Divergent responses of leaf phenology to changing temperature among plant species and geographical regions. Ecosphere 6, art250 (2015).
Google Scholar 
9.Zhang, G., Zhang, Y., Dong, J. & Xiao, X. Green-up dates in the Tibetan Plateau have continuously advanced from 1982 to 2011. Proc. Natl Acad. Sci. USA 110, 4309–4314 (2013).CAS 

Google Scholar 
10.Menzel, A. et al. European phenological response to climate change matches the warming pattern. Glob. Change Biol. 12, 1969–1976 (2006).
Google Scholar 
11.Cleland, E. E., Chuine, I., Menzel, A., Mooney, H. A. & Schwartz, M. D. Shifting plant phenology in response to global change. Trends Ecol. Evol. 22, 357–365 (2007).
Google Scholar 
12.Menzel, A., Sparks, T. H., Estrella, N. & Roy, D. B. Altered geographic and temporal variability in phenology in response to climate change. Glob. Ecol. Biogeogr. 15, 498–504 (2006).
Google Scholar 
13.Zhang, X., Tarpley, D. & Sullivan, J. T. Diverse responses of vegetation phenology to a warming climate. Geophys. Res. Lett. https://doi.org/10.1029/2007gl031447 (2007).14.Fitter, A. H. & Fitter, R. S. Rapid changes in flowering time in British plants. Science 296, 1689–1691 (2002).CAS 

Google Scholar 
15.Primack, R. B. et al. Spatial and interspecific variability in phenological responses to warming temperatures. Biol. Conserv. 142, 2569–2577 (2009).
Google Scholar 
16.Cleland, E. E., Chiariello, N. R., Loarie, S. R., Mooney, H. A. & Field, C. B. Diverse responses of phenology to global changes in a grassland ecosystem. Proc. Natl Acad. Sci. USA 103, 13740–13744 (2006).CAS 

Google Scholar 
17.Wang, H., Dai, J., Zheng, J. & Ge, Q. Temperature sensitivity of plant phenology in temperate and subtropical regions of China from 1850 to 2009. Int. J. Climatol. 35, 913–922 (2015).
Google Scholar 
18.Chuine, I. M., Morin, X. & Bugmann, H. Warming, photoperiods, and tree phenology. Science 329, 277–278 (2010).
Google Scholar 
19.Chuine, I. A unified model for budburst of trees. J. Theor. Biol. 207, 337–347 (2000).CAS 

Google Scholar 
20.Murray, M., Cannell, M. G. R. & Smith, R. I. Date of budburst of fifteen tree species in Britain following climatic warming. J. Appl. Ecol. 26, 693–700 (1989).
Google Scholar 
21.Man, R., Lu, P. & Dang, Q. L. Insufficient chilling effects vary among boreal tree species and chilling duration. Front. Plant Sci. 8, 1354 (2017).
Google Scholar 
22.Cannell, M. G. R. & Smith, R. I. L. Thermal time, chill days and prediction of budburst in Picea sitchensis. J. Appl. Ecol. 20, 951–963 (1983).
Google Scholar 
23.Fu, Y. H. et al. Increased heat requirement for leaf flushing in temperate woody species over 1980-2012: effects of chilling, precipitation and insolation. Glob. Change Biol. 21, 2687–2697 (2015).
Google Scholar 
24.Zhang, H., Liu, S., Regnier, P. & Yuan, W. New insights on plant phenological response to temperature revealed from long-term widespread observations in China. Glob. Change Biol. 24, 2066–2078 (2018).
Google Scholar 
25.Yu, H., Luedeling, E. & Xu, J. Winter and spring warming result in delayed spring phenology on the Tibetan Plateau. Proc. Natl Acad. Sci. USA 107, 22151–22156 (2010).CAS 

Google Scholar 
26.Asse, D. et al. Warmer winters reduce the advance of tree spring phenology induced by warmer springs in the Alps. Agric. For. Meteorol. 252, 220–230 (2018).
Google Scholar 
27.Ettinger, A. K. et al. Winter temperatures predominate in spring phenological responses to warming. Nat. Clim. Change 10, 1137–1142 (2020).
Google Scholar 
28.Chuine, I. & Régnière, J. Process-based models of phenology for plants and animals. Annu. Rev. Ecol. Evol. Syst. 48, 159–182 (2017).
Google Scholar 
29.Caffarra, A., Donnelly, A., Chuine, I. & Jones, M. B. Modelling the timing of Betula pubescens budburst. I. Temperature and photoperiod: a conceptual model. Clim. Res. 46, 147–157 (2011).
Google Scholar 
30.Luterbacher, J., Dietrich, D., Xoplaki, E., Grosjean, M. & Wanner, H. European seasonal and annual temperature variability, trends, and extremes since 1500. Science 303, 1499–1503 (2004).CAS 

Google Scholar 
31.Ciais, P. et al. 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 (eds Stocker, T. F. et al.) (Cambridge Univ. Press, 2013).32.Fu, Y. H. et al. Daylength helps temperate deciduous trees to leaf-out at the optimal time. Glob. Change Biol. 25, 2410–2418 (2019).
Google Scholar 
33.Wolkovich, E. M. et al. A simple explanation for declining temperature sensitivity with warming. Glob. Change Biol. 27, 4947–4949 (2021).CAS 

Google Scholar 
34.Templ, B. et al. Pan European Phenological database (PEP725): a single point of access for European data. Int. J. Biometeorol. 62, 1109–1113 (2018).
Google Scholar 
35.Kramer, K. Selecting a model to predict the onset of growth of Fagus sylvatica. J. Appl. Ecol. 31, 172–181 (1994).
Google Scholar 
36.Chuine, I., Cour, P. & Rousseau, D.-D. Selecting models to predict the timing of flowering of temperate trees: implications for tree phenology modelling. Plant Cell Environ. 22, 1–13 (1999).37.Savas, R. Investigations on the annual cycle of development of forest trees. II. Autumn dormancy and winter dormancy https://eurekamag.com/research/000/414/000414639.php (1974).38.Hänninen, H. Modelling bud dormancy release in trees from cool and temperate regions. Acta. Fenn. 14, 499–454 (1990).
Google Scholar 
39.Harrington, C. A., Gould, P. J. & St. Clair, J. B. Modeling the effects of winter environment on dormancy release of Douglas-fir. Ecol. Manag. 259, 798–808 (2010).
Google Scholar 
40.Zhang, H., Yuan, W., Liu, S., Dong, W. & Fu, Y. Sensitivity of flowering phenology to changing temperature in China. J. Geophys. Res. Biogeosci. 120, 1658–1665 (2015).
Google Scholar 
41.Richardson, A. D. et al. Influence of spring phenology on seasonal and annual carbon balance in two contrasting New England forests. Tree Physiol. 29, 321–331 (2009).CAS 

Google Scholar 
42.Piao, S. et al. Plant phenology and global climate change: current progresses and challenges. Glob. Change Biol. 25, 1922–1940 (2019).
Google Scholar 
43.Körner, C. & Basler, D. Phenology under global warming. Science 327, 1461–1462 (2010).
Google Scholar 
44.Zohner, C. M. & Renner, S. S. Common garden comparison of the leaf-out phenology of woody species from different native climates, combined with herbarium records, forecasts long-term change. Ecol. Lett. 17, 1016–1025 (2014).
Google Scholar 
45.Vitasse, Y. & Basler, D. What role for photoperiod in the bud burst phenology of European beech. Eur. J. For. Res. 132, 1–8 (2012).
Google Scholar 
46.Lenz, A., Hoch, G., Körner, C. & Vitasse, Y. Convergence of leaf-out towards minimum risk of freezing damage in temperate trees. Funct. Ecol. 30, 1480–1490 (2016).
Google Scholar 
47.Wang, Y. et al. Forest controls spring phenology of juvenile Smith fir along elevational gradients on the southeastern Tibetan Plateau. Int. J. Biometeorol. 63, 963–972 (2019).
Google Scholar 
48.Marquis, B., Bergeron, Y., Simard, M. & Tremblay, F. Probability of sping frosts, not growing degree-days, drives onset of spruce bud burst in plantations at the boreal-temperate forest ecotone. Front. Plant Sci. 11, 1031 (2020).
Google Scholar 
49.Shen, M., Piao, S., Cong, N., Zhang, G. & Jassens, I. A. Precipitation impacts on vegetation spring phenology on the Tiberan Plateau. Glob. Change Biol. 21, 3647–3656 (2015).
Google Scholar 
50.Liu et al. Temperature, precipitation, and insolation effects on autumn vegetation phenology in temperate China. Glob. Change Biol. 22, 644–655 (2016).CAS 

Google Scholar 
51.Minder, J. R., Mote, P. W. & Lundquist, J. D. Surface temperature lapse rates over complex terrain: lessons from the Cascade Mountains. J. Geophys. Res. 115, D14122 (2010).
Google Scholar 
52.Navarro-Serrano et al. Elevation effects on air temperature in a topographically complex mountain valley in the Spanish Pyrenees. Atmosphere 11, 656 (2020).
Google Scholar 
53.Chen, L. et al. Leaf senescence exhibits stronger climatic responses during warm than during cold autumns. Nat. Clim. Change 10, 777–780 (2020).CAS 

Google Scholar 
54.Leys, C., Ley, C., Klein, O., Bernard, P. & Licata, L. Detecting outliers: do not use standard deviation around the mean, use absolute deviation around the median. J. Exp. Soc. Psychol. 49, 764–766 (2013).
Google Scholar 
55.Beer, C. et al. Harmonized European long-term climate data for assessing the effect of changing temporal variability on land–atmosphere CO2 fluxes. J. Clim. 27, 4815–4834 (2014).
Google Scholar 
56.Olsson, C. & Jönsson, A. M. Process-based models not always better than empirical models for simulating budburst of Norway spruce and birch in Europe. Glob. Change Biol. 20, 3492–3507 (2014).
Google Scholar 
57.Duan, Q., Sorooshian, S. & Gupta, V. K. Optimal use of the SCE-UA global optimization method for calibrating watershed models. J. Hydrol. 158, 265–284 (1994).
Google Scholar 
58.Bluemel, K. & Chmielewski, F. Shortcomings of classical phenological forcing models and a way to overcome them. Agric. For. Meteorol. 164, 10–19 (2012).
Google Scholar  More

1.Brooks, B. W. et al. Are harmful algal blooms becoming the greatest inland water quality threat to public health and aquatic ecosystems? Environ. Toxicol. Chem. 35, 6–13 (2016).
Google Scholar 
2.Lopez, C., Jewett, E., Dortch, Q., Walton, B. & Hudnell, H. Scientific Assessment of Freshwater Harmful Algal Blooms (United States National Ocean Service, 2008)3.Huisman, J. et al. Cyanobacterial blooms. Nat. Rev. Microbiol. 16, 471–483 (2018).
Google Scholar 
4.Paerl, H. W. & Paul, V. J. Climate change: links to global expansion of harmful cyanobacteria. Water Res. 46, 1349–1363 (2012).
Google Scholar 
5.Carmichael, W. W. The toxins of cyanobacteria. Sci. Am. 270, 78–86 (1994).
Google Scholar 
6.Carmichael, W. W. et al. Human fatalities from cyanobacteria: chemical and biological evidence for cyanotoxins. Environ. Health Persp. 109, 663–668 (2001).
Google Scholar 
7.Botswana: mystery elephant deaths caused by cyanobacteria. BBC News https://www.bbc.com/news/world-africa-54234396 (2020).8.Paerl, H. W. & Huisman, J. Blooms like it hot. Science 320, 57–58 (2008).
Google Scholar 
9.O’Neil, J. M., Davis, T. W., Burford, M. A. & Gobler, C. J. The rise of harmful cyanobacteria blooms: the potential roles of eutrophication and climate change. Harmful Algae 14, 313–334 (2012).
Google Scholar 
10.Kutser, T. Quantitative detection of chlorophyll in cyanobacterial blooms by satellite remote sensing. Limnol. Oceanogr. 49, 2179–2189 (2004).
Google Scholar 
11.Kutser, T., Metsamaa, L., Strömbeck, N. & Vahtmäe, E. Monitoring cyanobacterial blooms by satellite remote sensing. Estuar. Coast. Shelf Sci. 67, 303–312 (2006).
Google Scholar 
12.Binding, C. E., Pizzolato, L. & Zeng, C. EOLakeWatch; delivering a comprehensive suite of remote sensing algal bloom indices for enhanced monitoring of Canadian eutrophic lakes. Ecol. Indic. 121, 106999 (2021).
Google Scholar 
13.Stumpf, R. P. et al. Challenges for mapping cyanotoxin patterns from remote sensing of cyanobacteria. Harmful Algae 54, 160–173 (2016).
Google Scholar 
14.Matthews, M. W. Eutrophication and cyanobacterial blooms in South African inland waters: 10 years of MERIS observations. Remote Sens. Environ. 155, 161–177 (2014).
Google Scholar 
15.Mishra, S. et al. Measurement of cyanobacterial bloom magnitude using satellite remote sensing. Sci. Rep. 9, 18310 (2019).
Google Scholar 
16.Hu, C. et al. Moderate Resolution Imaging Spectroradiometer (MODIS) observations of cyanobacteria blooms in Taihu Lake, China. J. Geophys. Res. 115, C04002 (2010).
Google Scholar 
17.Song, K. et al. Climatic versus anthropogenic controls of decadal trends (1983–2017) in algal blooms in lakes and reservoirs across China. Environ. Sci. Technol. 55, 2929–2938 (2021).
Google Scholar 
18.Coffer, M. M., Schaeffer, B. A., Darling, J. A., Urquhart, E. A. & Salls, W. B. Quantifying national and regional cyanobacterial occurrence in US lakes using satellite remote sensing. Ecol. Indic. 111, 105976 (2020).
Google Scholar 
19.Ho, J., Michalak, A. & Pahlevan, N. Widespread global increase in intense lake phytoplankton blooms since the 1980s. Nature 574, 667–670 (2019).
Google Scholar 
20.Dierssen, H. M., Kudela, R. M., Ryan, J. P. & Zimmerman, R. C. Red and black tides: quantitative analysis of water-leaving radiance and perceived color for phytoplankton, colored dissolved organic matter, and suspended sediments. Limnol. Oceanogr. 51, 2646–2659 (2006).
Google Scholar 
21.Michalak, A. M. et al. Record-setting algal bloom in Lake Erie caused by agricultural and meteorological trends consistent with expected future conditions. Proc. Natl Acad. Sci. USA 110, 6448–6452 (2013).
Google Scholar 
22.Binding, C., Greenberg, T., McCullough, G., Watson, S. & Page, E. An analysis of satellite-derived chlorophyll and algal bloom indices on Lake Winnipeg. J. Great Lakes Res. 44, 436–446 (2018).
Google Scholar 
23.Guo, L. Doing battle with the green monster of Taihu Lake. Science 317, 1166–1166 (2007).
Google Scholar 
24.Moradi, M. Comparison of the efficacy of MODIS and MERIS data for detecting cyanobacterial blooms in the southern Caspian Sea. Mar. Pollut. Bull. 87, 311–322 (2014).
Google Scholar 
25.Schindler, D. W. Eutrophication and recovery in experimental lakes: implications for lake management. Science 184, 897–899 (1974).
Google Scholar 
26.Qin, B. et al. Water depth underpins the relative roles and fates of nitrogen and phosphorus in lakes. Environ. Sci. Technol. 54, 3191–3198 (2020).
Google Scholar 
27.Beman, J. M., Arrigo, K. R. & Matson, P. A. Agricultural runoff fuels large phytoplankton blooms in vulnerable areas of the ocean. Nature 434, 211–214 (2005).
Google Scholar 
28.Yu, C. et al. Managing nitrogen to restore water quality in China. Nature 567, 516–520 (2019).
Google Scholar 
29.Zhang, X. et al. Managing nitrogen for sustainable development. Nature 528, 51–59 (2015).
Google Scholar 
30.Hobbie, S. E. et al. Contrasting nitrogen and phosphorus budgets in urban watersheds and implications for managing urban water pollution. Proc. Natl Acad. Sci. USA 114, 4177–4182 (2017).
Google Scholar 
31.Wang, Z. China’s wastewater treatment goals. Science 338, 604–604 (2012).
Google Scholar 
32.Sutton, M. A. et al. The European Nitrogen Assessment: Sources, Effects and Policy Perspectives (Cambridge Univ. Press, 2011).33.Litke, D. W. Review of Phosphorus Control Measures in the United States and Their Effects on Water Quality (US Geological Survey, 1999).34.Kosten, S. et al. Warmer climates boost cyanobacterial dominance in shallow lakes. Glob. Change Biol. 18, 118–126 (2012).
Google Scholar 
35.Carey, C. C., Ibelings, B. W., Hoffmann, E. P., Hamilton, D. P. & Brookes, J. D. Eco-physiological adaptations that favour freshwater cyanobacteria in a changing climate. Water Res. 46, 1394–1407 (2012).
Google Scholar 
36.Wells, M. L. et al. Harmful algal blooms and climate change: learning from the past and present to forecast the future. Harmful Algae 49, 68–93 (2015).
Google Scholar 
37.Elliott, J. A. The seasonal sensitivity of cyanobacteria and other phytoplankton to changes in flushing rate and water temperature. Glob. Change Biol 16, 864–876 (2010).
Google Scholar 
38.Jeppesen, E. et al. in Shallow Lakes ’95 (eds Kufel, L. et al.) 151–164 (Springer, 1997).39.O’Reilly, C. M. et al. Rapid and highly variable warming of lake surface waters around the globe. Geophys. Res. Lett. 42, 773–710,781 (2015).
Google Scholar 
40.Janssen, A. B. G. et al. How to model algal blooms in any lake on earth. Curr. Opin. Environ. Sustain 36, 1–10 (2019).
Google Scholar 
41.Woodcock, C. E. et al. Free access to Landsat imagery. Science 320, 1011 (2008).
Google Scholar 
42.Zhu, Z. & Woodcock, C. E. Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sens. Environ. 118, 83–94 (2012).
Google Scholar 
43.Masek, J. G. et al. A Landsat surface reflectance dataset for North America, 1990–2000. IEEE Geosci. Remote Sens. Lett. 3, 68–72 (2006).
Google Scholar 
44.Vermote, E., Justice, C., Claverie, M. & Franch, B. Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product. Remote Sens. Environ. 185, 46–56 (2016).
Google Scholar 
45.Irish, R. R. Landsat 7 Science Data Users Handbook 415–430 (US Geological Survey, 2000).46.Messager, M. L., Lehner, B., Grill, G., Nedeva, I. & Schmitt, O. Estimating the volume and age of water stored in global lakes using a geo-statistical approach. Nat. Commun. 7, 13603 (2016).
Google Scholar 
47.Wang, J. et al. Recent global decline in endorheic basin water storages. Nat. Geosci. 11, 926–932 (2018).
Google Scholar 
48.Pekel, J.-F., Cottam, A., Gorelick, N. & Belward, A. S. High-resolution mapping of global surface water and its long-term changes. Nature 540, 418–422 (2016).
Google Scholar 
49.McNally, A. et al. A land data assimilation system for sub-Saharan Africa food and water security applications. Sci. Data 4, 170012 (2017).
Google Scholar 
50.CIESIN Gridded Population of the World v.4 (NASA SEDAC, 2018).51.Bouwman, L. et al. Exploring global changes in nitrogen and phosphorus cycles in agriculture induced by livestock production over the 1900–2050 period. Proc. Natl Acad. Sci. USA 110, 20882–20887 (2013).
Google Scholar 
52.Pickens, A. H. et al. Mapping and sampling to characterize global inland water dynamics from 1999 to 2018 with full Landsat time-series. Remote Sens. Environ. 243, 111792 (2020).
Google Scholar 
53.Feng, L. & Hu, C. Land adjacency effects on MODIS Aqua top-of-atmosphere radiance in the shortwave infrared: statistical assessment and correction. J. Geophys. Res. Oceans 122, 4802–4818 (2017).
Google Scholar 
54.Walsh, S. E. et al. Global patterns of lake ice phenology and climate: model simulations and observations. J. Geophys. Res. Atmos. 103, 28825–28837 (1998).
Google Scholar 
55.Yang, X., Pavelsky, T. M. & Allen, G. H. The past and future of global river ice. Nature 577, 69–73 (2020).
Google Scholar 
56.Hu, C. et al. Dynamic range and sensitivity requirements of satellite ocean color sensors: learning from the past. Appl. Opt. 51, 6045–6062 (2012).
Google Scholar 
57.Kuhn, C. & Butman, D. Declining greenness in Arctic-boreal lakes. Proc. Natl Acad. Sci. USA 118, e2021219118 (2021).
Google Scholar 
58.Kirillin, G. et al. Physics of seasonally ice-covered lakes: a review. Aquat. Sci. 74, 659–682 (2012).
Google Scholar 
59.Kotovirta, V., Toivanen, T., Järvinen, M., Lindholm, M. & Kallio, K. Participatory surface algal bloom monitoring in Finland in 2011–2013. Environ. Syst. Res. 3, 24 (2014).
Google Scholar 
60.Cronberg, G., Annadotter, H. & Lawton, L. A. The occurrence of toxic blue-green algae in Lake Ringsjön, southern Sweden, despite nutrient reduction and fish biomanipulation. Hydrobiologia 404, 123–129 (1999).
Google Scholar 
61.Romarheim, A. T. & Riise, G. Development of Cyanobacteria in Årungen (Norsk vannforening, 2009)62.Robertson, A. R. The CIE 1976 color‐difference formulae. Color Res. Appl. 2, 7–11 (1977).
Google Scholar 
63.Mouw, C. B. et al. Aquatic color radiometry remote sensing of coastal and inland waters: challenges and recommendations for future satellite missions. Remote Sens. Environ. 160, 15–30 (2015).
Google Scholar 
64.Wasmund, N., Nausch, G. & Matthäus, W. Phytoplankton spring blooms in the southern Baltic Sea—spatio-temporal development and long-term trends. J. Plankton Res. 20, 1099–1117 (1998).
Google Scholar 
65.Hu, C. A novel ocean color index to detect floating algae in the global oceans. Remote Sens. Environ. 113, 2118–2129 (2009).
Google Scholar 
66.Fairman, H. S., Brill, M. H. & Hemmendinger, H. How the CIE 1931 color-matching functions were derived from Wright-Guild data. Color Res. Appl. 22, 11–23 (1997).
Google Scholar 
67.Chander, G., Markham, B. L. & Helder, D. L. Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote Sens. Environ. 113, 893–903 (2009).
Google Scholar 
68.Feng, L. et al. Radiometric cross-calibration of Gaofen-1 WFV cameras using Landsat-8 OLI images: a solution for large view angle associated problems. Remote Sens. Environ. 174, 56–68 (2016).
Google Scholar 
69.Yu, X. et al. An empirical algorithm to seamlessly retrieve the concentration of suspended particulate matter from water color across ocean to turbid river mouths. Remote Sens. Environ. 235, 111491 (2019).
Google Scholar 
70.Hou, X., Feng, L., Chen, X. & Zhang, Y. Dynamics of the wetland vegetation in large lakes of the Yangtze Plain in response to both fertilizer consumption and climatic changes. ISPRS J. Photogramm. Remote Sens. 141, 148–160 (2018).
Google Scholar 
71.Lee, Z., Pahlevan, N., Ahn, Y.-H., Greb, S. & O’Donnell, D. Robust approach to directly measuring water-leaving radiance in the field. Appl. Opt. 52, 1693–1701 (2013).
Google Scholar 
72.Liu, L., Peng, W., Wu, L. & Liu, L. Water quality assessment of Danjiangkou Reservoir and its tributaries in China. IOP Conf. Ser. Earth Environ. Sci. 112, 012008 (2018).
Google Scholar 
73.Li, X. et al. The color formation mechanism of the blue karst lakes in Jiuzhaigou Nature Reserve, Sichuan, China. Water 12, 771 (2020).
Google Scholar 
74.Wurtsbaugh, W. & Marcarelli, A. Eutrophication in Farmington Bay, Great Salt Lake, Utah 2005 Annual Report (Utah State Univ., 2006).75.Hammer, U. T. Saline Lake Ecosystems of the World Vol. 59 (Springer, 1986). More

1.Macdonald IR, Guinasso NL, Sassen R, Brooks JM, Lee L, Scott KT. Gas hydrate that breaches the sea-floor on the continental-slope of the Gulf-of-Mexico. Geology. 1994;22:699–702.CAS 

Google Scholar 
2.Kvenvolden KA. A review of the geochemistry of methane in natural gas hydrate. Org Geochem. 1995;23:997–1008.CAS 

Google Scholar 
3.Milkov AV. Molecular and stable isotope compositions of natural gas hydrates: a revised global dataset and basic interpretations in the context of geological settings. Org Geochem. 2005;36:681–702.CAS 

Google Scholar 
4.Cragg BA, Parkes RJ, Fry JC, Weightman AJ, Rochelle PA, Maxwell JR. Bacterial populations and processes in sediments containing gas hydrates (ODP Leg 146: Cascadia Margin). Earth Planet Sc Lett. 1996;139:497–507.CAS 

Google Scholar 
5.Yoshioka H, Maruyama A, Nakamura T, Higashi Y, Fuse H, Sakata S, et al. Activities and distribution of methanogenic and methane-oxidizing microbes in marine sediments from the Cascadia Margin. Geobiology. 2010;8:223–33.CAS 
PubMed 

Google Scholar 
6.Yoshioka H, Sakata S, Cragg BA, Parkes RJ, Fujii T. Microbial methane production rates in gas hydrate-bearing sediments from the eastern Nankai Trough, off central Japan. Geochem J. 2009;43:315–21.CAS 

Google Scholar 
7.Heuer VB, Inagaki F, Morono Y, Kubo Y, Spivack AJ, Viehweger B, et al. Temperature limits to deep subseafloor life in the Nankai Trough subduction zone. Science. 2020;370:1230–4.CAS 
PubMed 

Google Scholar 
8.Wellsbury P, Goodman K, Cragg BA, Parkes RJ. The geomicrobiology of deep marine sediments from Blake Ridge containing methane hydrate (sites 994, 995 and 997). Proc Ocean Drill Program Sci results. 2000;164:379–91.
Google Scholar 
9.Bidle KA, Kastner M, Bartlett DH. A phylogenetic analysis of microbial communities associated with methane hydrate containing marine fluids and sediments in the Cascadia margin (ODP site 892B). Fems Microbiol Lett. 1999;177:101–8.CAS 
PubMed 

Google Scholar 
10.Reed DW, Fujita Y, Delwiche ME, Blackwelder DB, Sheridan PP, Uchida T, et al. Microbial communities from methane hydrate-bearing deep marine sediments in a forearc basin. Appl Environ Microbiol. 2002;68:3759–70.CAS 
PubMed 
PubMed Central 

Google Scholar 
11.Briggs BR, Inagaki F, Morono Y, Futagami T, Huguet C, Rosell-Mele A, et al. Bacterial dominance in subseafloor sediments characterized by methane hydrates. FEMS Microbiol Ecol. 2012;81:88–98.CAS 
PubMed 

Google Scholar 
12.Kendall MM, Boone DR. Cultivation of methanogens from shallow marine sediments at Hydrate Ridge, Oregon. Archaea. 2006;2:31–38.CAS 
PubMed 
PubMed Central 

Google Scholar 
13.Fry JC, Parkes RJ, Cragg BA, Weightman AJ, Webster G. Prokaryotic biodiversity and activity in the deep subseafloor biosphere. FEMS Microbiol Ecol. 2008;66:181–96.CAS 
PubMed 

Google Scholar 
14.Nunoura T, Takaki Y, Shimamura S, Kakuta J, Kazama H, Hirai M, et al. Variance and potential niche separation of microbial communities in subseafloor sediments off Shimokita Peninsula, Japan. Environ Microbiol. 2016;18:1889–906.CAS 
PubMed 

Google Scholar 
15.Mikucki JA, Liu Y, Delwiche M, Colwell FS, Boone DR. Isolation of a methanogen from deep marine sediments that contain methane hydrates, and description of Methanoculleus submarinus sp. nov. Appl Environ Microbiol. 2003;69:3311–6.CAS 
PubMed 
PubMed Central 

Google Scholar 
16.Weng C-Y, Chen S-C, Lai M-C, Wu S-Y, Lin S, Yang TF, et al. Methanoculleus taiwanensis sp. nov., a methanogen isolated from deep marine sediment at the deformation front area near Taiwan. Int J Syst Evol Micr. 2015;65:1044–9.CAS 

Google Scholar 
17.Kendall MM, Liu Y, Sieprawska-Lupa M, Stetter KO, Whitman WB, Boone DR. Methanococcus aeolicus sp. nov., a mesophilic, methanogenic archaeon from shallow and deep marine sediments. Int J Syst Evol Microbiol. 2006;56:1525–9.CAS 
PubMed 

Google Scholar 
18.Strąpoć D, Ashby M, Wood L, Levinson R, Huizinga B. Significant contribution of methyl/methanol-utilising methanogenic pathway in a subsurface biogas environment. In: Skovhus T, Whitby C, editors. Applied microbiology and molecular biology in oilfield systems. Dordrecht: Springer; 2010. p. 211–6.19.Guo H, Liu R, Yu Z, Zhang H, Yun J, Li Y, et al. Pyrosequencing reveals the dominance of methylotrophic methanogenesis in a coal bed methane reservoir associated with Eastern Ordos Basin in China. Int J Coal Geol. 2012;93:56–61.CAS 

Google Scholar 
20.Katayama T, Yoshioka H, Muramoto Y, Usami J, Fujiwara K, Yoshida S, et al. Physicochemical impacts associated with natural gas development on methanogenesis in deep sand aquifers. ISME J. 2015;9:436–46.CAS 
PubMed 

Google Scholar 
21.Yanagawa K, Tani A, Yamamoto N, Hachikubo A, Kano A, Matsumoto R, et al. Biogeochemical cycle of methanol in anoxic deep-sea sediments. Microbes Environ. 2016;31:190–3.PubMed 
PubMed Central 

Google Scholar 
22.Colwell F, Matsumoto R, Reed D. A review of the gas hydrates, geology, and biology of the Nankai Trough. Chem Geol. 2004;205:391–404.CAS 

Google Scholar 
23.Uchida T, Waseda A, Namikawa T. Methane accumulation and high concentration of gas hydrate in marine and terrestrial sandy sediments. In: Collett T, Johnson A, Knapp C, Boswell R, editors. Natural gas hydrates: energy resource potential and associated geologic hazards. Tulsa: American Association of Petroleum Geologists Memoir 89; 2009. p. 401–13.24.Katayama T, Yoshioka H, Takahashi HA, Amo M, Fujii T, Sakata S. Changes in microbial communities associated with gas hydrates in subseafloor sediments from the Nankai Trough. FEMS Microbiol Ecol. 2016;92:fiw093.PubMed 

Google Scholar 
25.Oba M, Sakata S, Fujii T. Archaeal polar lipids in subseafloor sediments from the Nankai Trough: Implications for the distribution of methanogens in the deep marine subsurface. Org Geochem. 2015;78:153–60.CAS 

Google Scholar 
26.Noguchi S, Shimoda N, Takano O, Oikawa N, Inamori T, Saeki T, et al. 3-D internal architecture of methane hydrate-bearing turbidite channels in the eastern Nankai Trough, Japan. Mar Pet Geol. 2011;28:1817–28.
Google Scholar 
27.Fujii T, Suzuki K, Takayama T, Tamaki M, Komatsu Y, Konno Y, et al. Geological setting and characterization of a methane hydrate reservoir distributed at the first offshore production test site on the Daini-Atsumi Knoll in the eastern Nankai Trough, Japan. Mar Pet Geol. 2015;66:310–22.CAS 

Google Scholar 
28.Kanno T, Fukuhara M, Osawa O, Chee S, Takekoshi M, Wang X, et al. Estimation of geothermal gradient in marine gas-hydrate-bearing formation in the Eastern Nankai Trough. Beijing, China: Proceedings of the 8th International Conference on Gas Hydrates (ICGH8–2014); 2014.29.Kaneko M, Takano Y, Ogawa NO, Sato Y, Yoshida N, Ohkouchi N. Estimation of methanogenesis by quantification of coenzyme F430 in marine sediments. Geochem J. 2016;50:453–60.CAS 

Google Scholar 
30.Kaneko M, Takano Y, Chikaraishi Y, Ogawa NO, Asakawa S, Watanabe T, et al. Quantitative analysis of coenzyme F430 in environmental samples: a new diagnostic tool for methanogenesis and anaerobic methane oxidation. Anal Chem. 2014;86:3633–8.CAS 
PubMed 

Google Scholar 
31.Katayama T, Kamagata Y Cultivation of Methanogens. Hydrocarbon and lipid microbiology protocols. In: McGenity T, Timmis K, Nogales B, editors. Springer protocols handbooks. Berlin, Heidelberg: Springer; 2016. p. 177–95.32.Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009;75:7537–41.CAS 
PubMed 
PubMed Central 

Google Scholar 
33.Pruesse E, Quast C, Knittel K, Fuchs BM, Ludwig W, Peplies J, et al. SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB. Nucleic Acids Res. 2007;35:7188–96.CAS 
PubMed 
PubMed Central 

Google Scholar 
34.Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011;28:2731–9.CAS 
PubMed 
PubMed Central 

Google Scholar 
35.Jobb G, von Haeseler A, Strimmer K. TREEFINDER: a powerful graphical analysis environment for molecular phylogenetics. BMC Evol Biol. 2004;4:18.PubMed 
PubMed Central 

Google Scholar 
36.Whiticar MJ. Carbon and hydrogen isotope systematics of bacterial formation and oxidation of methane. Chem Geol. 1999;161:291–314.CAS 

Google Scholar 
37.Scheller S, Goenrich M, Thauer RK, Jaun B. Methyl-coenzyme M reductase from methanogenic archaea: Isotope effects on the formation and anaerobic oxidation of methane. J Am Chem Soc. 2013;135:14975–84.CAS 
PubMed 

Google Scholar 
38.Diekert G, Konheiser U, Piechulla K, Thauer RK. Nickel requirement and factor F430 content of methanogenic bacteria. J Bacteriol. 1981;148:459–64.CAS 
PubMed 
PubMed Central 

Google Scholar 
39.Mayr S, Latkoczy C, Krüger M, Günther D, Shima S, Thauer RK, et al. Structure of an F430 variant from archaea associated with anaerobic oxidation of methane. J Am Chem Soc. 2008;130:10758–67.CAS 
PubMed 

Google Scholar 
40.House CH, Orphan VJ, Turk KA, Thomas B, Pernthaler A, Vrentas JM, et al. Extensive carbon isotopic heterogeneity among methane seep microbiota. Environ Microbiol. 2009;11:2207–15.CAS 
PubMed 

Google Scholar 
41.Lloyd KG, Alperin MJ, Teske A. Environmental evidence for net methane production and oxidation in putative ANaerobic MEthanotrophic (ANME) archaea. Environ Microbiol. 2011;13:2548–64.CAS 
PubMed 

Google Scholar 
42.Laso-Pérez R, Wegener G, Knittel K, Widdel F, Harding KJ, Krukenberg V, et al. Thermophilic archaea activate butane via alkyl-coenzyme M formation. Nature. 2016;539:396–401.PubMed 

Google Scholar 
43.Inagaki F, Nunoura T, Nakagawa S, Teske A, Lever M, Lauer A, et al. Biogeographical distribution and diversity of microbes in methane hydrate-bearing deep marine sediments on the Pacific Ocean Margin. Proc Natl Acad Sci USA. 2006;103:2815–20.CAS 
PubMed 
PubMed Central 

Google Scholar 
44.Marchesi JR, Weightman AJ, Cragg BA, Parkes RJ, Fry JC. Methanogen and bacterial diversity and distribution in deep gas hydrate sediments from the Cascadia Margin as revealed by 16S rRNA molecular analysis. FEMS Microbiol Ecol. 2001;34:221–8.CAS 
PubMed 

Google Scholar 
45.Nunoura T, Inagaki F, Delwiche ME, Colwell FS, Takai K. Subseafloor microbial communities in methane hydrate-bearing sediment at two distinct locations (ODP Leg 204) in the Cascadia Margin. Microbes Environ. 2008;23:317–25.PubMed 

Google Scholar 
46.Cord-Ruwisch R, Ollivier B. Interspecific hydrogen transfer during methanol degradation by Sporomusa acidovorans and hydrogenophilic anaerobes. Arch Microbiol. 1986;144:163–5.CAS 

Google Scholar 
47.Heijthuijsen JHFG, Hansen TA. Interspecies hydrogen transfer in co-cultures of methanol-utilizing acidogens and sulfate-reducing or methanogenic bacteria. FEMS Microbiol Ecol. 1986;2:57–64.
Google Scholar 
48.Eichler B, Schink B. Oxidation of primary aliphatic alcohols by Acetobacterium carbinolicum sp. nov., a homoacetogenic anaerobe. Arch Microbiol. 1984;140:147–52.CAS 

Google Scholar 
49.Parkes RJ, Cragg B, Roussel E, Webster G, Weightman A, Sass H. A review of prokaryotic populations and processes in sub-seafloor sediments, including biosphere: geosphere interactions. Mar Geol. 2014;352:409–25.CAS 

Google Scholar 
50.Imachi H, Aoi K, Tasumi E, Saito Y, Yamanaka Y, Saito Y, et al. Cultivation of methanogenic community from subseafloor sediments using a continuous-flow bioreactor. ISME J. 2011;5:1913–25.CAS 
PubMed 
PubMed Central 

Google Scholar 
51.Newberry CJ, Webster G, Cragg BA, Parkes RJ, Weightman AJ, Fry JC. Diversity of prokaryotes and methanogenesis in deep subsurface sediments from the Nankai Trough, Ocean Drilling Program Leg 190. Environ Microbiol. 2004;6:274–87.PubMed 

Google Scholar 
52.Orsi WD, Edgcomb VP, Christman GD, Biddle JF. Gene expression in the deep biosphere. Nature 2013;499:205–8.CAS 
PubMed 

Google Scholar 
53.Vigneron A, L’Haridon S, Godfroy A, Roussel EG, Cragg BA, Parkes RJ, et al. Evidence of active methanogen communities in shallow sediments of the sonora margin cold seeps. Appl Environ Microbiol. 2015;81:3451–9.CAS 
PubMed 
PubMed Central 

Google Scholar  More

The Levant region, the major land bridge connecting Africa with Eurasia, was a significant dispersal route for Hominins and fauna during the early Pleistocene1,2,3. But while there are numerous Eurasian early Pleistocene sites, fossil hominin remains are rare and present only at four localities dated between 1.1 and 1.9 Mya4,5,6,7,8,9,10,11: Dmanisi (Georgia), Venta Micena (Orce, Granada), Modjokerto and Sangiran (Java, Indonesia), and Sima De Elefante (Atapuerca, Spain) (Supplementary 2: Table 1; Fig. 1a). In contrast, early Pleistocene east African sites containing Homo cranial remains are more abundant, but postcranial remains are scarcer, and the best-preserved skeleton is Nariokotome KNM-WT 1500012,13.Figure 1‘Ubeidya site locality. (a) Map of Africa and Eurasia with major Pleistocene paleoanthropological sites. Black circles denote sites with no osteological remains; red circles denote sites with human osteological remains. (b) The location of the site of ‘Ubeidiya, south of lake Kineret (Sea of Galilee), on the western banks of the Jordan Valley (red circle) (c) aerial photograph of the excavation plan of ‘Ubeidiya with the location of layer II-23 where UB 10749 was recovered.Full size imageIn the Levant, the only site from this time-period with hominin remains is ‘Ubeidiya at the western escarpment of the Jordan Valley which is a part of the broader Rift Valley (Supplementary 1: Fig. 1b,c). The fossil remains include cranial fragments (UB 1703, 1704, 1705, and 1706), two incisor (UB 1700, UB 335) and a molar (UB 1701), identified as Homo cf. erectus/ergaster14,15,16,17,18. It is important to note that some of these fragments were bulldozed out of the ground preceding the first season, while others are considered intrusive and younger than the surroundings deposits17.In 2018 during a reanalysis of the faunal assemblages done by two of the authors (A. B, and M. B.) a complete vertebral body (UB 10749) with hominin characteristics was found. This is the first hominin postcranial remain found at ‘Ubeidiya securely assigned to early Pleistocene deposits (See “Materials and methods”).Here we assess the taxonomic affinity of UB 10749, its serial location along the spinal column, its chronological and physiological age at death, estimate the specimen’s height and weight, and detect any pathological or taphonomic changes. Based on our findings, we explore the unique developmental characteristics of the UB 10749 within the context of early Homo paleobiology and its implications for hominin dispersal out of Africa.Description of the findingUB 10749 is a complete vertebral body (Fig. 2). The superior plate of the vertebra is oval, with an uneven surface, indicating non-ossification of the vertebral endplate. Similarly, the inferior plate is also oval with marked postero-lateral edges. A small pit is found in the center of both superior and inferior plates. The inferior plate is bilaterally wider than the superior plate. The anterior and lateral walls are smooth and slightly concave i.e., their superior and inferior edges are more prominent than the center. There is no evidence of rib attachment to the body on the lateral walls. The posterior wall can be divided into three parts, the center and right and left lateral thirds. The central part is smooth with two nutritional foramina. The two lateral thirds are located at the junction between the vertebral body and the pedicles. Their surface is uneven, indicating that the pedicle had not yet ossified to the vertebral body. In a lateral view, the vertebra shows a lordotic wedging as the height of the anterior wall is greater than that of the posterior wall (Supplementary 2: Table 2). The oval shape of the vertebral body, the concavity of the inferior plate, the lordotic wedging, and the lack of rib bearing facets all indicate a lower lumbar vertebra, i.e., presacral (PS) 1, PS2, or PS3 (corresponding to L5, L4, and L3 in modern humans).Figure 2UB 10749 vertebral body. (a) Superior view; (b) posterior view; (c) inferior view; (d) anterior view.Full size imageA micro-CT (µCT) scan of UB 10749 (Fig. 3) reveals a well-developed cortical bone on the anterior and lateral walls and the central part of the posterior wall. The cancellous bone at the superior and inferior plates is very thin, as is the bone at the lateral thirds of the posterior wall, indicating that these were not yet ossified. The µCT scan also reveals well-developed canals within the vertebral body –Bastons’ venous plexus19 (Fig. 3c). A small pit at the superior and inferior plates is seen in the mid-sagittal and coronal planes of the CT scan (Fig. 3a, b). A thin vertical region that appears black on the µCT connects the two pits, indicating that this area was not yet ossified.Figure 3µCT scan of UB 10749. (a) Midsagittal section. (b) Coronal section. (c) Horizontal section.Full size imageTaxonomic identificationWe compared UB 10749 to a range of mammalian species from, but not limited to, those present in ‘Ubeidiya such as carnivores (e.g., Ursus, Hyeana, Panthera), Artiodactyla (e.g., Hippopotamus, Praemegaceros), Perissodactyla (Rhinocertidae, Equidae), Proboscidea (Mamuthus, Elephas), and Primates (Homo, Pongo, Gorilla, Theropithecus and Papio).UB 10749 lacks the inward indentation on the posterior wall distinctive of Ursus and is short cranio-caudally, as opposed to the longer vertebral bodies of ungulates.The size, the large vertebral plate, and the relatively short vertebral body of UB 10749 indicates that it belongs to hominoidea. The lordotic wedging and the concavity of the inferior plate further suggests that this is a hominin vertebra20,21.To narrow the taxonomic identification, we compared UB 10749 to a range of extant and extinct hominin species, and to Pan as an outgroup (Supplementary 2: Table 3). The analysis revealed that the best index to which best differentiates between lumbar vertebral bodies of Homo and Pan is ‘superior length to posterior height’ (Fig. 4; Supplementary 2: Table 4). This index also differentiates between Homo and Australopithecus22. Compared to the three presacral vertebrae (PS1–PS3) of hominins and Pan, UB 10749 falls within the range of Homo and outside the range of Pan or Australopithecus. It falls near the position of the vertebrae of KNM-WT-15000, an early Pleistocene sub adult specimen from east Africa. Therefore, we conclude that the vertebra at hand most probably belongs to an early Pleistocene Homo.Figure 4Comparison of UB 10749 to other hominoids. Vertebral body ratio (superior length to posterior height) of each of the lower 3 presacral vertebra in modern humans, neandertals, australopith, chimpanzees, KNM-WT 15000 and UB 10749. Note that UB 10749 is consistently falls within the range of Homo, and beyond the range of chimpanzees and australopith.Full size imageSerial allocation of the vertebral bodyIt is well known, especially in Hominoidea, that there is a vast overlap in the shape of adjacent lumbar vertebral bodies23. We conducted three separate analyses to address this problem: (1) Vertebral wedging i.e., the ratio of posterior height/anterior height which significantly separates the vertebral segments PS1, PS2, and PS3 of modern humans (Supplementary 2: Fig. 1; Supplementary 2: Table 4), (2) A principal component analysis (PCA) of vertebral linear indices (Fig. 5a; Supplementary 2: Table 4), and (3) Geometric morphometrics (GM) shape analysis (Fig. 5b). Vertebral wedging sets UB 10749 as PS2. The vertebral linear indices PCA sets the UB 10749 as either PS2 or PS3, and the GM shape analysis sets the vertebra as either PS1 or PS2. Based on these results, we estimate that the serial allocation of UB 10749 is most likely PS2.Figure 5Serial allocation of UB 10749. (a) PCA of vertebral body ratios of modern humans, KNM-WT 15000, STS 14, and UB 10749 (see Supplementary Table 4). Note the overlap between adults and juvenile in each presacral vertebra. UB 10749 falls within the range of PS2–PS3. Note that KNM-WT 15000 and STS 14 follow the same pattern as modern humans. (b) PCA results for GM shape analysis. UB 10749 falls within the range of PS1, but with proximity to PS2. Note that KNM-WT 15000 and STS 14 follow the same pattern as modern humans. In both analysis: Circles denotes juvenile; Squares denotes adults. Red denotes PS1; Blue denotes PS2; Green denotes PS3.Full size imageAge at deathAge at death is estimated based on level of ossification, relative vertebral size, or vertebral shape. The lack of neural canal ossification in UB 10749 indicates an approximate age of 3–6-years-old compared to modern humans24,25 (Supplementary 2: Fig. 2), although it is important to note that several authors report high variability in pedicle ossification, up to 16-years-old26,27. The absence of vertebral endplate ossification also supports the young age of UB 10749, indicating that the vertebra belongs to an individual that had not reached puberty28.In contrast, based on its size alone, UB 10749 would be assigned an older age, probably between 11 and 15-year-old modern human (Fig. 6a: Supplementary 2: Table 5). However, vertebral size is highly variable with age, and we cannot rule out either a younger or older age. Finally, geometric morphometric principal component shape analysis suggests that UB 10749 falls within the range of 6–10-years-old modern humans (Fig. 6b). This is confirmed by a linear discriminant analysis which also places UB 10749 well within the 6–10 years old group (Supplementary 2: Fig. 4). Considering all the above information, we estimate that the age at death for UB 10749 is between 6 and 12-years-old.Figure 6Age at death of UB 10749. (a) Vertebral body size (combined sample of PS1–PS3) of modern humans, KNM-WT 15000 and UB 10749 (see Supplementary 2: Table 5). UB 10749 falls within the range of 11–15 years or the lower end of adults. (b) PCA results for GM shape analysis of modern human, KNM-WT 15000, STS 14, and UB 10749 vertebral bodies. UB 10749 falls within the range of the 6–10 age group. In both analyses: Red, 0–5 years old; Green, 6–10 years old; Blue, 11–15 years old; Brown, 16-adults.Full size imageHeight and weight estimationHeight (stature) and weight at death is estimated based on a range of equations and growth charts for modern humans (Supplementary 2: Tables 6–8). The estimated average height at death of UB 10749, points to a height of 155 cm. This height is comparable to a 13 years-old boy or a 12.5 years-old girl, based on CDC growth charts. A height of 155 cm is above the 95 percentile of 10 years old and above the 75 percentile for 12 years old modern humans29. As the age estimation for UB 10749 is 6–12 years, it seems that this individual was tall for its age.Weight is estimated based on height or based on chronological age. Based on height, UB 10749 was 45–55 kg, while based on age, the weight of UB 10749 was 20–43 kg (Supplementary 2: Table 7). Since height is a stronger predictor for weight than age30, we estimate the weight at death at about 45–50 kg.A single juvenile vertebral body is not a definitive predictor for adult height and weight. Even more so, the growth pattern of early Pleistocene hominins is unknown. Thus, to cautiously estimate the adult height and weight of UB 10749, calculations were based on several methods: modern American (CDC growth charts), modern Sudanese population31, and chimpanzees32.Assuming UB 10749 was 6–12 years old, based on chimpanzees’ growth charts, it would have reached adult height of 155–192 cm and weighted 50–101 kg. Based on modern American and Sudanese growth charts, UB 10749 displays a range of a height between 168 and 247 cm and a weight between 62 and 173 kg (Supplementary 2: Table 8). The average height and weight predication for adult size of UB 10749 is 198 cm and 100 kg. Although we cannot rule out any of the estimations, based on its size at death and predicated adult size, UB 10749 was most likely a large-bodied hominin33,34,35.TaphonomyVery thin fluviatile deposits are evident on the surface of the vertebra, despite being cleaned during the excavation. Aside from that, there is no apparent taphonomic alteration or post-depositional breakage.PathologyThe completeness of the vertebral body and its bilateral symmetry do not suggest pathological processes or developmental deformities that may have affected the vertebra, such as osteoarthritis, disc herniation, spondylosis, tuberculosis, brucellosis, or scoliosis36. However, in the absence of the vertebral arch, we cannot rule out anterior slippage of the vertebral body, i.e., spondylolisthesis or facet joint deformities. The discrepancy between the size of the vertebral body and the level of ossification is puzzling. The size of UB 10749 is equivalent to an 11–15-year-old modern human, and the level of ossification is equivalent to a 3–6-year-old modern human child. This discrepancy might result from several factors, including developmental or pathological conditions, such as: persistent notochondral canal; hypopituitarism; androgen deficiency; genetic mutation24,37,38 (see Supplementary 2 for discussion regarding possible pathology). While these conditions are rare in modern humans, they cannot be ruled out. Another possibility is that UB 10749 displays a different ossification pattern than observed in modern humans or great apes25,39. 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