Ecology

1.
Charbonneau, M. R. et al. A microbial perspective of human developmental biology. Nature 535, 48–55 (2016).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 
2.
Bäckhed, F. et al. Dynamics and stabilization of the human gut microbiome during the first year of life. Cell Host Microbe 17, 690–703 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

3.
Stewart, C. J. et al. Temporal development of the gut microbiome in early childhood from the TEDDY study. Nature 562, 583–588 (2018).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

4.
Derrien, M., Alvarez, A. S. & de Vos, W. M. The gut microbiota in the first decade of life. Trends Microbiol. 27, 997–1010 (2019).
CAS  PubMed  Article  Google Scholar 

5.
Yatsunenko, T. et al. Human gut microbiome viewed across age and geography. Nature 486, 222–227 (2012).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

6.
Lim, E. S. et al. Early life dynamics of the human gut virome and bacterial microbiome in infants. Nat. Med. 21, 1228–1234 (2015).
CAS  PubMed  PubMed Central  Article  Google Scholar 

7.
Palmer, C., Bik, E. M., DiGiulio, D. B., Relman, D. A. & Brown, P. O. Development of the human infant intestinal microbiota. PLoS Biol. 5, e177 (2007).
PubMed  PubMed Central  Article  CAS  Google Scholar 

8.
Lynch, S. V. & Pedersen, O. The human intestinal microbiome in health and disease. N. Engl. J. Med. 375, 2369–2379 (2016).
CAS  PubMed  Article  Google Scholar 

9.
Honda, K. & Littman, D. R. The microbiota in adaptive immune homeostasis and disease. Nature 535, 75–84 (2016).
ADS  CAS  PubMed  Article  Google Scholar 

10.
Fischbach, M. A. Microbiome: focus on causation and mechanism. Cell 174, 785–790 (2018).
CAS  PubMed  PubMed Central  Article  Google Scholar 

11.
Widder, S. et al. Challenges in microbial ecology: building predictive understanding of community function and dynamics. ISME J. 10, 2557–2568 (2016).
MathSciNet  PubMed  PubMed Central  Article  Google Scholar 

12.
Vrancken, G., Gregory, A. C., Huys, G. R. B., Faust, K. & Raes, J. Synthetic ecology of the human gut microbiota. Nat. Rev. Microbiol. 17, 754–763 (2019).
CAS  PubMed  Article  Google Scholar 

13.
Walter, J., Armet, A. M., Finlay, B. B. & Shanahan, F. Establishing or exaggerating causality for the gut microbiome: lessons from human microbiota-associated rodents. Cell 180, 221–232 (2020).
CAS  PubMed  Article  Google Scholar 

14.
Wolfe, B. E., Button, J. E., Santarelli, M. & Dutton, R. J. Cheese rind communities provide tractable systems for in situ and in vitro studies of microbial diversity. Cell 158, 422–433 (2014).
CAS  PubMed  PubMed Central  Article  Google Scholar 

15.
Connell, J. H. & Slatyer, R. O. Mechanisms of succession in natural communities and their role in community stability and organization. Am. Nat. 111, 1119–1144 (1977).
Article  Google Scholar 

16.
Bertness, M. D. & Callaway, R. Positive interactions in communities. Trends Ecol. Evol. 9, 191–193 (1994).
CAS  PubMed  Article  Google Scholar 

17.
Shade, A. et al. Macroecology to unite all life, large and small. Trends Ecol. Evol. 33, 731–744 (2018).
PubMed  Article  Google Scholar 

18.
Gregory, K. E. et al. Influence of maternal breast milk ingestion on acquisition of the intestinal microbiome in preterm infants. Microbiome 4, 68 (2016).
PubMed  PubMed Central  Article  Google Scholar 

19.
Gibson, M. K. et al. Developmental dynamics of the preterm infant gut microbiota and antibiotic resistome. Nat. Microbiol. 1, 16024 (2016).
CAS  PubMed  PubMed Central  Article  Google Scholar 

20.
DiBartolomeo, M. E. & Claud, E. C. The developing microbiome of the preterm infant. Clin. Ther. 38, 733–739 (2016).
PubMed  PubMed Central  Article  Google Scholar 

21.
La Rosa, P. S. et al. Patterned progression of bacterial populations in the premature infant gut. Proc. Natl Acad. Sci. USA 111, 12522–12527 (2014).
ADS  PubMed  Article  CAS  Google Scholar 

22.
Costello, E. K., Carlisle, E. M., Bik, E. M., Morowitz, M. J. & Relman, D. A. Microbiome assembly across multiple body sites in low-birthweight infants. MBio 4, e00782-13 (2013).
PubMed  PubMed Central  Article  CAS  Google Scholar 

23.
Stewart, C. J. et al. Temporal bacterial and metabolic development of the preterm gut reveals specific signatures in health and disease. Microbiome 4, 67 (2016).
PubMed  PubMed Central  Article  Google Scholar 

24.
Pammi, M. et al. Intestinal dysbiosis in preterm infants preceding necrotizing enterocolitis: a systematic review and meta-analysis. Microbiome 5, 31 (2017).
PubMed  PubMed Central  Article  Google Scholar 

25.
Gasparrini, A. J. et al. Persistent metagenomic signatures of early-life hospitalization and antibiotic treatment in the infant gut microbiota and resistome. Nat. Microbiol. 4, 2285–2297 (2019).
PubMed  PubMed Central  Article  CAS  Google Scholar 

26.
Reynolds, L. A. & Finlay, B. B. Early life factors that affect allergy development. Nat. Rev. Immunol. 17, 518–528 (2017).
CAS  PubMed  Article  Google Scholar 

27.
Gensollen, T., Iyer, S. S., Kasper, D. L. & Blumberg, R. S. How colonization by microbiota in early life shapes the immune system. Science 352, 539–544 (2016).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

28.
Bokulich, N. A. et al. Antibiotics, birth mode, and diet shape microbiome maturation during early life. Sci. Transl. Med. 8, 343ra82 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

29.
Shao, Y. et al. Stunted microbiota and opportunistic pathogen colonization in caesarean-section birth. Nature 574, 117–121 (2019).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

30.
Pasolli, E. et al. Extensive unexplored human microbiome diversity revealed by over 150,000 genomes from metagenomes spanning age, geography, and lifestyle. Cell 176, 649–662 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

31.
The Integrative HMP (iHMP) Research Network Consortium. The Integrative Human Microbiome Project. Nature 569, 641–648 (2019).
ADS  Article  CAS  Google Scholar 

32.
Nash, A. K. et al. The gut mycobiome of the Human Microbiome Project healthy cohort. Microbiome 5, 153 (2017).
PubMed  PubMed Central  Article  Google Scholar 

33.
Limon, J. J., Skalski, J. H. & Underhill, D. M. Commensal fungi in health and disease. Cell Host Microbe 22, 156–165 (2017).
CAS  PubMed  PubMed Central  Article  Google Scholar 

34.
Koskinen, K. et al. First insights into the diverse human archaeome: specific detection of Archaea in the gastrointestinal tract, lung, and nose and on skin. MBio 8, 00824-17 (2017).
Article  Google Scholar 

35.
Durán, P. et al. Microbial interkingdom interactions in roots promote Arabidopsis survival. Cell 175, 973–983 (2018).
PubMed  PubMed Central  Article  CAS  Google Scholar 

36.
Carr, A., Diener, C., Baliga, N. S. & Gibbons, S. M. Use and abuse of correlation analyses in microbial ecology. ISME J. 13, 2647–2655 (2019).
PubMed  PubMed Central  Article  Google Scholar 

37.
Contijoch, E. J. et al. Gut microbiota density influences host physiology and is shaped by host and microbial factors. eLife 8, e40553 (2019).
PubMed  PubMed Central  Article  Google Scholar 

38.
Vandeputte, D. et al. Quantitative microbiome profiling links gut community variation to microbial load. Nature 551, 507–511 (2017).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

39.
Stämmler, F. et al. Adjusting microbiome profiles for differences in microbial load by spike-in bacteria. Microbiome 4, 28 (2016).
PubMed  PubMed Central  Article  Google Scholar 

40.
Ishwaran, H. & Rao, J. S. Spike and slab variable selection: frequentist and Bayesian strategies. Ann. Stat. 33, 730–773 (2005).
MathSciNet  MATH  Article  Google Scholar 

41.
Gonze, D., Coyte, K. Z., Lahti, L. & Faust, K. Microbial communities as dynamical systems. Curr. Opin. Microbiol. 44, 41–49 (2018).
PubMed  Article  Google Scholar 

42.
Bucci, V. et al. MDSINE: Microbial Dynamical Systems INference Engine for microbiome time-series analyses. Genome Biol. 17, 121 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

43.
Freilich, M. A., Wieters, E., Broitman, B. R., Marquet, P. A. & Navarrete, S. A. Species co-occurrence networks: can they reveal trophic and non-trophic interactions in ecological communities? Ecology 99, 690–699 (2018).
PubMed  Article  Google Scholar 

44.
Friedman, J. & Alm, E. J. Inferring correlation networks from genomic survey data. PLOS Comput. Biol. 8, e1002687 (2012).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

45.
Faust, K. et al. Microbial co-occurrence relationships in the human microbiome. PLOS Comput. Biol. 8, e1002606 (2012).
CAS  PubMed  PubMed Central  Article  Google Scholar 

46.
Watts, S. C., Ritchie, S. C., Inouye, M. & Holt, K. E. FastSpar: rapid and scalable correlation estimation for compositional data. Bioinformatics 35, 1064–1066 (2019).
CAS  PubMed  Article  Google Scholar 

47.
Stein, R. R. et al. Ecological modeling from time-series inference: insight into dynamics and stability of intestinal microbiota. PLOS Comput. Biol. 9, e1003388 (2013).
PubMed  PubMed Central  Article  CAS  Google Scholar 

48.
Fisher, C. K. & Mehta, P. Identifying keystone species in the human gut microbiome from metagenomic timeseries using sparse linear regression. PLoS ONE 9, e102451 (2014).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

49.
Pammi, M., Liang, R., Hicks, J., Mistretta, T. A. & Versalovic, J. Biofilm extracellular DNA enhances mixed species biofilms of Staphylococcus epidermidis and Candida albicans. BMC Microbiol. 13, 257 (2013).
PubMed  PubMed Central  Article  CAS  Google Scholar  More

1.
BPS. Statistik Indonesia 2018 (Badan Pusat Statistik, Jakarta, 2018).
Google Scholar 
2.
KLHK. Pedoman penilaian efektivitas pengelolaan kawasan konservasi di Indonesia. (Direktorat Kawasan Konservasi dan Direktorat Jenderal Konservasi Sumber Daya Alam dan Ekosistem, Kementerian Lingkungan Hidup dan Kehutanan, Jakarta, 2015).

3.
Ditjen KSDAE. Statistik Direktorat Jenderal KSDAE 2017 (Kementerian Lingkungan Hidup dan Kehutanan Direktorat Jenderal Konservasi Sumber Daya Alam dan Ekosistem, Jakarta, 2018).

4.
Mulyana, A. et al. Kebijakan pengelolaan zona khusus. Dapatkah meretas kebuntuan dalam menata ruang Taman Nasional di Indonesia? Brief, 1 (2010).

5.
CBD. Convention on Biological Diversity https://www.cbd.int/sp/ (2019).

6.
Murninngtyas, E. et al. (eds) Indonesian Biodiversity Strategy and Action Plan (IBSAP) 2015–2020 (Ministry of the National Development Planning/ BAPPENAS, Jakarta, 2016).
Google Scholar 

7.
Wiratno. Sepuluh cara baru kelola kawasan konservasi di Indonesia: membangun “organisasi pembelajar” (Direktorat Jenderal KSDAE Kementrian Lingkungan Hidup dan Kehutanan, Jakarta, 2018).

8.
Coad, L. et al. Measuring impact of protected area management interventions: current and future use of the Global Database of Protected Area Management Effectiveness. Philos. Trans. R. Soc. B 370, 20140281. https://doi.org/10.1098/rstb.2014.0281 (2015).
Article  Google Scholar 

9.
Geldmann, J. Evaluating the effectiveness of protected areas for maintaining biodiversity, securing habitats, and reducing threats. PhD thesis. (PhD School of the Faculty of Science, University of Copenhagen, 2013).

10.
Supriatna, J., Dwiyahreni, A. A., Winarni, N., Mariati, S. & Margules, C. Deforestation of primate habitat on Sumatra and adjacent islands Indonesia. Primate Conserv. 31, 71–82 (2017).
Google Scholar 

11.
Gray, C. L. et al. Local biodiversity is higher inside than outside terrestrial protected areas worldwide. Nat. Commun. 7(12306), 1–7. https://doi.org/10.1038/ncomms12306 (2016).
CAS  Article  Google Scholar 

12.
Eklund, J. & Cabeza, M. Quality of governance and effectiveness of protected areas: crucial concepts for conservation planning. Ann. N. Y. Acad. Sci. https://doi.org/10.1111/nyas.13284 (2016).
Article  PubMed  Google Scholar 

13.
Jones, K. R. et al. One-third of global protected land is under intense human pressure. Science 360, 788–789. https://doi.org/10.1126/science.aap9565 (2018).
CAS  Article  PubMed  Google Scholar 

14.
Wiratno. Tersesat di jalan yang benar. Seribu hari mengelola Leuser (Direktorat PKPS, Jakarta, 2012).

15.
Gaveau, D. L. A. et al. Examining protected area effectiveness in Sumatra: importance of regulations governing unprotected lands. Conserv. Lett. 5, 142–148. https://doi.org/10.1111/j.1755-263X.2011.00220.x (2012).
Article  Google Scholar 

16.
Azmi, W. & Gunaryadi, D. Current status of Asian elephants in Indonesia. Gajah 35, 55–61 (2011).
Google Scholar 

17.
O’Brien, T. G. & Kinnaird, M. F. Changing populations of birds and mammals in North Sulawesi. Oryx 30, 150–156 (1996).
Article  Google Scholar 

18.
Wheeler, P. & Dwiyahreni, A. Large mammal monitoring in Lambusango. Interim Progress Report January 2007 (University of Hull, UK, 2007).

19.
Gaveau, D. L. A., Wich, S. A. & Marshall, A. J. Are protected areas conserving primate habitat in Indonesia? In An introduction to primate conservation (eds Wich, S. A. & Marshal, A. J.) 193–200 (Oxford University Press, Oxford, 2016).
Google Scholar 

20.
Wibisono, H. T. & Pusparini, W. Sumatran tiger (Panthera tigris sumatrae): a review of conservation status. Integr. Zool. 5, 313–323. https://doi.org/10.1111/j.1749-4877.2010.00219.x (2010).
Article  PubMed  Google Scholar 

21.
Dwiyahreni, A.A. et al. Forest cover changes in Indonesia’s terrestrial national parks between 2012 and 2017. Biodiversitas, 22(3), 1235–1242. https://doi.org/10.13057/biodiv/d220320 (2021).  

22.
Sanderson, E. W. et al. The human footprint and the last of the wild. Bioscience 52(10), 891–904. https://doi.org/10.1641/0006-3568(2002)052[0891:THFATL]2.0.CO;2 (2002).
Article  Google Scholar 

23.
CIESIN. United Nation’s World Population Prospects (UN WPP)-Adjusted Population Density https://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-density-adjusted-to-2015-unwpp-country-totals-rev11 (2018).

24.
Wan, J. Z., Wang, C. J. & Yuc, F. H. Human footprint and climate disappearance in vulnerable ecoregions of protected areas. Global Planet. Change 170, 260–268. https://doi.org/10.1016/j.gloplacha.2018.09.002 (2018).
ADS  Article  Google Scholar 

25.
Leroux, S. J. et al. Global protected areas and IUCN designations: do the categories match the conditions?. Biol. Conserv. 143, 609–616. https://doi.org/10.1016/j.biocon.2009.11.018 (2010).
Article  Google Scholar 

26.
Ayram, C. A. C., Mendoza, M. E., Etter, A. & Salicrup, D. R. P. Anthropogenic impact on habitat connectivity: a multidimensional human footprint index evaluated in a highly biodiverse landscape of Mexico. Ecol. Ind. 72, 895–909. https://doi.org/10.1016/j.ecolind.2016.09.007 (2017).
Article  Google Scholar 

27.
Tapia-Armijos, M. F., Homeier, J. & Muntac, D. D. Spatio-temporal analysis of the human footprint in South Ecuador: influence of human pressure on ecosystems and effectiveness of protected areas. Appl. Geogr. 78, 22–32. https://doi.org/10.1016/j.apgeog.2016.10.007 (2017).
Article  Google Scholar 

28.
Li, S., Wu, J., Gong, J. & Li, S. Human footprint in Tibet: assessing the spatial layout and effectiveness of nature reserves. Sci. Total Environ. 621, 18–29. https://doi.org/10.1016/j.scitotenv.2017.11.216 (2018).
ADS  CAS  Article  PubMed  Google Scholar 

29.
Geldmann, J. et al. Effectiveness of terrestrial protected areas in reducing habitat loss and population declines. Biol. Conserv. 161, 230–238. https://doi.org/10.1016/j.biocon.2013.02.018 (2013).
Article  Google Scholar 

30.
Anderson, E. & Mammides, C. The role of protected areas in mitigating human impact in the world’s last wilderness areas. Ambio https://doi.org/10.1007/s13280-019-01213-x (2019).
Article  PubMed  PubMed Central  Google Scholar 

31.
Mulyana, A., Kosmaryandi, N., Hakim, N., Suryadi, S. & Suwito. Ruang adaptif: refleksi penataan zona/blok di kawasan konservasi (Direktorat Pemolaan dan Informasi Konservasi Alam dan Direktorat Jenderal Konservasi Sumberdaya Alam dan Ekosistem, Kementerian Lingkungan Hidup dan Kehutanan, Bogor, 2019).

32.
Landuse thematic map 2012 and 2017 (Indonesian Ministry of Environment and Forestry)

33.
Binamarga Map 2014 (Indonesian Ministry of Public Works)

34.
Rusmin, N., Edy, S., Robby, A. & Dwi, K. R. Study of the potential expansion of new rice fields in Central Maluku District to support food security in Maluku Province. IOP Conf. Ser. Earth Environ. Sci. 334, 012067. https://doi.org/10.1088/1755-1315/334/1/012067 (2019).
Article  Google Scholar 

35.
Harrison, M. E., Capilla, B. R., Thornton, S. A., Cattau, M. E., & Page, S.E. Impacts of the 2015 fire season on peat-swamp forest biodiversity in Indonesian Borneo. In 15th International Peat Congress 2016, 713–717 (2016).

36.
Laurance, W. F. et al. Averting biodiversity collapse in tropical forest protected areas. Nature 00, 1–5. https://doi.org/10.1038/nature11318 (2012).
CAS  Article  Google Scholar 

37.
Opdam, P. Exploring the role of science in sustainable landscape management: an introduction to the special issue. Sustainability 10, 1–6. https://doi.org/10.3390/su10020331 (2018).
Article  Google Scholar 

38.
Field, D. R. Symbiotic relationships between national parks and neighboring social-biological regions in National Parks and rural development. In Practice, Policy in the United States (eds Machlis, G. E. & Field, D. R.) 211–218 (Island Press, Washington, DC, 2000).
Google Scholar 

39.
IUCN. Protected area categories. Category II: National Park. https://www.iucn.org/theme/protected-areas/about/protected-areas-categories/category-ii-national-park (2017).

40.
Alamgir, M. et al. High-risk infrastructure projects pose imminent threats to forests in Indonesian Borneo. Sci. Rep. https://doi.org/10.1038/s41598-018-36594-8 (2019).
Article  PubMed  PubMed Central  Google Scholar 

41.
Sloan, S. et al. Transnational conservation and infrastructure development in the heart of Borneo. PLoS ONE 14(9), e0221947. https://doi.org/10.1371/journal.pone.0221947 (2019).
CAS  Article  PubMed  PubMed Central  Google Scholar 

42.
Sloan, S., Alamgir, M., Campbell, M. J., Setyawati, T. & Laurance, W. F. Development corridors and remnant-forest conservation in Sumatra Indonesia. Trop. Conserv. Sci. 12, 1–9. https://doi.org/10.1177/1940082919889509 (2019).
Article  Google Scholar 

43.
Healey, R. M. et al. Road mortality threatens endemic species in a national park in Sulawesi Indonesia. Glob. Ecol. Conserv. 24, e01281. https://doi.org/10.1016/j.gecco.2020.e01281 (2020).
Article  Google Scholar 

44.
Du, W., Penabaz-Wiley, S. M., Njeru, A. M. & Kinoshita, I. Models and approaches for integrating protected areas with their surroundings: a review of the literature. Sustainability 7, 8151–8177. https://doi.org/10.3390/su7078151 (2015).
Article  Google Scholar 

45.
Verma, M. et al. Severe human pressures in the Sundaland biodiversity hotspot. Conserv. Sci. Pract. 2(e169), 2020. https://doi.org/10.1111/csp2.169 (2020).
Article  Google Scholar 

46.
DeFries, R., Hansen, A., Newton, A. C. & Hansen, M. C. Increasing isolation of protected areas in tropical forests over the past twenty years. Ecol. Appl. 15(1), 19–26 (2005).
Article  Google Scholar 

47.
Brun, C. et al. Analysis of deforestation and protected area effectiveness in Indonesia: a comparison of Bayesian spatial models. Glob. Environ. Change 31, 285–295. https://doi.org/10.1016/j.gloenvcha.2015.02.004 (2015).
Article  Google Scholar 

48.
Mariati, S., Kusnoputranto, H., Supriatna, J. & Koestoer, R. H. Habitat Loss of Sumatran elephants (Elephas maximus sumatranus) in Tesso Nilo Forest, Riau, Indonesia. Aust. J. Basic Appl. Sci. 8(2), 248–255 (2014).
Google Scholar 

49.
Busch, J. & Ferretti-Gallon, K. What drives deforestation and what stops it? A meta-analysis. Rev. Environ. Econ. Policy 11(1), 3–23. https://doi.org/10.1093/reep/rew013 (2017).
Article  Google Scholar 

50.
Tacconi, L., Rodriguesa, R. J. & Maryudi, A. Law enforcement and deforestation: lessons for Indonesia from Brazil. For. Policy Econ. 108, 101943. https://doi.org/10.1016/j.forpol.2019.05.029 (2019).
Article  Google Scholar 

51.
Bruner, A., Gullison, G., Rice, R. E., da Fonseca, R. E. & Gustavo, A. B. Effectiveness of parks in protecting tropical biodiversity. Science 291, 125–128. https://doi.org/10.1126/science.291.5501.125 (2001).
ADS  CAS  Article  PubMed  Google Scholar 

52.
Adams, V. M., Iacona, G. D. & Possingham, H. P. Weighing the benefits of expanding protected areas versus managing existing ones. Nat. Sustain. https://doi.org/10.1038/s41893-019-0275-5 (2019).
Article  Google Scholar 

53.
Bickford, D. et al. In Indonesia’s protected areas need more protection: suggestions from island examples in biodiversity and human livelihoods in protected areas: case studies from the Malay Archipelago (eds Sodhi, N. S. et al.) 53–77 (Cambridge University Press, Cambridge, 2008).
Google Scholar 

54.
Baier, M. The Kayan Mentarang National Park: Indonesia’s new national park in North Central Borneo bordering Northern Sarawak and Sabah. Borneo Res. Bull. 40, 297 (2009).
Google Scholar 

55.
Anau, N., Hakim, A., Lekson, A. S. & Setyowati, E. Local wisdom practices of Dayak indigenous people in the management of Tana’ Ulen in the Kayan Mentarang National Park of Malinau Regency, North Kalimantan Province Indonesia. RJOAS 7(91), 156–167. https://doi.org/10.18551/rjoas.2019-07.16 (2019).
Article  Google Scholar 

56.
Susanti, R. & Zuhud, E. A. M. Traditional ecological knowledge and biodiversity conservation: the medicinal plants of the Dayak Krayan people in Kayan Mentarang National Park Indonesia. Biodiversitas 20(9), 2764–2779. https://doi.org/10.13057/biodiv/d200943 (2019).
Article  Google Scholar 

57.
Blankespoor, B., Dasgupta, S. & Wheeler, D. Protected areas and deforestation: new results from high-resolution panel data. Nat. Resour. Forum 41, 55–68. https://doi.org/10.1111/1477-8947.12118 (2017).
Article  Google Scholar 

58.
Sanderson, E. W., Walston, J. & Robinson, J. G. From bottleneck to breakthrough: urbanization and the future of biodiversity conservation. Bioscience 20, 1–15 (2018).
Google Scholar 

59.
Immanuel, G. et al. Indonesia. Pathways to sustainable land-use and food systems. FABLE Report (2019).

60.
Allan, J. R. et al. Recent increases in human pressure and forest loss threaten many natural world heritage sites. Biol. Conserv. 206, 47–55. https://doi.org/10.1016/j.biocon.2016.12.011 (2017).
Article  Google Scholar 

61.
Damayanti, E.K. Legality of National Parks and Involvement of Local People: Case Studies in Java, Indonesia and Kerala, India. Thesis. Ph.D. in Agricultural Science. (University of Tsukuba, Japan, 2008).

62.
Smiet, A. C. Forest ecology on Java: human impact and vegetation of montane forest. J. Trop. Ecol. 8, 129–152 (1992).
Article  Google Scholar 

63.
Wijaya, I. K. M. The semiotics of banyan trees spaces in Denpasar. Bali. LivaS Int. J. Livable 4(2), 48–59. https://doi.org/10.25105/livas.v4i2.5564 (2019).
Article  Google Scholar 

64.
Wijaya, M. H. & Sutrisni, K. How can the existence of customary laws protect the water preservation in the Cau Belayu (Tabanan) traditional village?. Yustisia 7(3), 600–613. https://doi.org/10.20961/yustisia.v7i3.21560 (2018).
Article  Google Scholar 

65.
Swandi, I. W. Kearifan lokal Bali untuk pelestarian alam: kajian wacana kartun-kartun majalah “Bog-Bog”. Jurnal Kajian Bali 7(2), 229–248. https://doi.org/10.24843/JKB.2017.v07.i02.p12 (2017).
Article  Google Scholar 

66.
Sobirin, S. Pranata Mangsa dan budaya kearifan lingkungan. Jurnal Budaya Nusantara 2(1), 250–264. https://doi.org/10.36456/b.nusantara.vol2.no1.a1719 (2018).
Article  Google Scholar 

67.
Apriando, T. Brubuh, kearifan masyarakat Jawa menjaga hutan. Mongabay. https://www.mongabay.co.id/2015/02/12/brubuh-kearifan-masyarakat-jawa-menjaga-hutan/ (2015).

68.
Sulfiantono, A., Hermawan, M. T. T. & Maluyi, A. Comparison of Effectiveness of the management of conservation areas of China and Indonesia. Int. J. Sci. 11, 73–82 (2013).
Google Scholar 

69.
Supriatna, J. Berwisata di Taman Nasional (Yayasan Obor, Jakarta, 2014).
Google Scholar 

70.
Venter, O. et al. Global terrestrial human footprint maps for 1993 and 2009. Sci. Data 3, 160067. https://doi.org/10.1038/sdata.2016.67 (2016).
Article  PubMed  PubMed Central  Google Scholar 

71.
Willmott, C. J. & Matsuura, K. Advantages of the mean absolute error (MAE) over the root mean square error (RMSE) in assessing average model performance. Clim. Res. 30, 79. https://doi.org/10.3354/CR030079 (2005).
Article  Google Scholar 

72.
Viera, A. J. & Garrett, J. M. Understanding interobserver agreement: the kappa statistic. Fam. Med. 37, 360–363 (2005).
PubMed  Google Scholar 

73.
R Development Core Team. R. A Language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, Austria, 2014). More

Tropical cyclone-induced bloom in NIO
We have analysed the tropical cyclones that occurred over the bay from 1997 to 2019 in accordance with the availability of satellite Chl-a measurements. Out of 51 storm events, 30 are identified as the phytoplankton bloom events (i.e. the Chl-a values greater than of 0.2 mg/m3) in BoB and 18 in AS across all seasons35,36,37. In the case of BoB, we have divided our analyses for pre-monsoon and post-monsoon, as the cyclone occurrences are rare in other seasons (e.g. winter and monsoon). Over AS, some cyclones also occur in the beginning of monsoon season. The spatio-temporal variability of cyclones is closely connected to the seasonal changes in the monsoon trough35,38. In pre-monsoon season, the trough passes over the northern BoB, but it passes through the central bay with an east-west orientation in the post-monsoon, and facilitates the formation of more number of TCs during the period39. The big seasonal change in wind shear and relative vorticity are the reasons for the lower number of cyclones in the pre-monsoon season. In general, the upwelling driven nutrient influx to the surface together with sunlight leads to the enhancement of Chl-a or phytoplankton bloom after the passage of cyclones in the open ocean40.
The Bay of Bengal cyclones
During pre-monsoon, the bay is least productive, but the western boundary current helps more production in the coastal regions41. The higher wind speed associated with TCs deepens (about 30 m) Mixed Layer Depth (MLD), and rupture the pycnocline and pumps nutrients to the surface24. In general, TCs occurring during pre-monsoon move northwards and pass north-eastern coast of India or Bangladesh42 (Fig. 1). To estimate the cyclone-induced phytoplankton bloom, we performed spatial analyses for each cyclone during the period 1997–2019 and selected cases are shown in Supplementary Fig. 1 for BOB01 in 2003, BOB01 in 2004 and Mala in 2006. The maps of Chl-a concentrations superimposed with Sea Level Anomaly (SLA) for the same period are shown in Supplementary Fig. 1. The cyclone BOB01 was a category 1 storm with a maximum sustained wind (MSW) of 39 m/s, which occurred during 10–19 May 2003. The Chl-a remained well below 0.2 mg/m3 prior to occurrence of the cyclone, which enhanced to 0.5 mg/m3 with a small area of about 1 mg/m3 at 9°–10° N, 86°–87° E, just after passage of the cyclone. This is consistent with its high Ekman Pumping Velocity (EPV) and small Translational Speed (TS) in that period. In addition, the Chl-a enhancement was higher for the Category 1 cyclone BOB01 that occurred during 14–19 May 2004 and was about 0.5 mg/m3 after passage of the cyclone. The bloom sustained for the next 5 days as also shown in Supplementary Fig. 1. On the other hand, Mala was the strongest cyclone occurred during the pre-monsoon period (24–28 April 2006) in the last two decades over BoB with a MSW speed of 61 m/s. The Chl-a increased from 0.1 to 1.0 mg/m3 after the cyclone passage in five days, with a small area of Chl-a about 1.0 mg/m3 on the immediate left of its track. The Chl-a concentration remained close to 0.5 mg/m3 in the next 5 days. Both EPV and TS were favourable for sustained upwelling and Chl-a bloom in the case of cyclone Mala. In all three cases, the closed contours of Sea surface Height Anomaly (SSHA) is negative, which indicate the presence of cold-core (cyclonic) eddies that triggered turbulent mixing and sustained Chl-a bloom. We have not used any specific eddy detection method but used the composite of SSHA to identify the presence of eddies, as done by Girishkumar et al.43.
Fig. 1: The cyclone tracks and pre-cyclone background Chl-a value.

a The study region North Indian Ocean (NIO) and the tracks for the pre- (green) and post-monsoon (yellow) cyclones for BoB and AS (red), as analysed for all cyclones occurred during the study period (1997–2019).

Full size image

We applied the same method to identify the phytoplankton bloom that occurred for each cyclone event in BoB after 1996 and estimated the corresponding EPV and TS for diagnosing the physical mechanisms that made different scales of blooms. The results are presented in Table 1 and Fig. 2. The analyses show that the bloom was comparatively larger for the cyclone BOB01 in 2004, about 3.28 mg/m3. The increase in Chl-a is negatively correlated with TS and is in agreement with the intensity of cyclone with a statistically significant correlation value of −0.30 (at the 95% confidence interval as per the P-test44). The TS is lower and bloom is larger for BOB01 in 2003, as the faster moving storms tend to intensify rapidly when compared to slower moving storms (i.e. wind speed 14 m/s). In contrast, the slow-moving storms expend more time over the ocean and thereby, increases the magnitude of upwelling to enhance the Chl-a over the region45. The estimated EPV is about 1.8 × 10−4 m/s for BOB01 in 2003 and is consistent with the observed Chl-a concentrations, whereas the EPV is about 1 × 10−4 m/s and Chl-a concentration is about 1.87 mg/m3 for the cyclone Mala. It suggests that TS has a prominent role in cyclone-induced upwelling and associated phytoplankton bloom.
Table 1 The cyclone-induced Chl-a blooms in the pre-monsoon seasons since 1997 in BoB.
Full size table

Fig. 2: The cyclone-induced Chl-a in BoB.

The observed enhancement in Chl-a (mg/m3) following the cyclone passage (5-day average) in the post and pre-monsoon seasons in BoB. The translational speed of the closest track points where the bloom occurred in (m/s), the ratio of wind speed to translational speed (WS/TS), and EPV (Ekman Pumping Velocity) as the cyclones reach their maximum intensity (m/s) are also shown in the lower panels.

Full size image

Post-monsoon is the active storm season over BoB, and about 25 cyclones with significant enhancement in Chl-a concentration are identified during the 1997–2019 period. As the haline stratification is stronger in BoB due to the monsoon rain and river water influx, the presence of BL increases SST, which fuel the storms over the bay46. Presence of BL weakens the impact of cooling in the mixed layer driven by cyclones and favours the intensification of post-monsoon cyclones2. We have also analysed the variability in BL, MLD, isothermal layer depth (ILD) and Chl-a for selected storms passed over the Argo Floats (see next section). The cyclones either form or develop further over the southeast BoB, but some move west northwest and cross the peninsular coast. Some cyclones recurve towards the west central bay and pass the central and northeast coast of India, but some hit Bangladesh and upper Burma coast47, as illustrated in Fig. 1.
Figure 3 presents a closer look at the bloom and its spatial distribution for selected cyclones during the post-monsoon season; e.g. the cyclones Sidr, Madi and Vardah overlaid with SSHA contours. Sidr, a category 4 cyclone occurred during 11–16 November 2007 with a MSW of about 44 m/s. The Chl-a is about 0.5 mg/m3 during the cyclone period at the right side of the track, but the bloom has spread to a wider area with values close to 0.5 mg/m3 just after the passage of cyclone. The analyses of SSHA further provide evidence for the eddy-mediated phytoplankton bloom. The phytoplankton bloom also sustained for another 5 days. This is also in agreement with that reported in other analyses, although bloom values were estimated for 19 November in the other studies48,49. The cyclone Madi, occurred during 6–13 December 2013, showed an enhancement of about 0.5 mg/m3 during the cyclone period with a region of 1 mg/m3 in the left side of the track. Some regions with 2–3 mg/m3 are also observed at the right and left sides of cyclone track, and the bloom sustained for the next 5 days with values of about 1 mg/m3 in the adjacent areas. The closed contours of negative SSHA suggest the presence of cyclonic eddies there. The phytoplankton bloom during this particular period is also contributed by the cyclone Lehar that occurred a week before, in 23–28 November; demonstrating the impact of occurrences of consecutive storms over the same oceanic region. Nevertheless, the cyclone Vardah showed an enhancement of about 1.92 mg/m3, which is higher than that of Sidr due to the higher EPV of the former. As for Lehar and Madi, there was another cyclone Nada that appeared during the period 29 November–2 December 2016, just before the appearance of Vardah, and that storm might have also contributed to the Chl-a bloom during the period of Vardah.
Fig. 3: The spatial extension of Chl-a bloom for selected post-monsoon cyclones in BoB.

The Chl-a averaged, overlaid with SSHA contours (solid—positive and dashed—negative), for 5 days before cyclone, during the entire cyclone period, five days just after the passage of cyclone and next 5 days for a Sidr (2007), b Madi (2013), and c Vardah (2016). The tracks of the respective cyclones are also shown.

Full size image

Figure 4 shows the time evolution of physical and biological observations during the period of TCs Phailin, Hudhud and Vardah. The biogeochemical Argo float WMO ID 2902086 was closer to the track of TC Phailin, and the float WMO ID 2902114 was near the tracks of Hudhud and Vardah. Supplementary Table 1 shows the name of cyclones, Argo IDs, and distance between the float and nearest track point of respective cyclones. Figure 4 (right) represents the subsurface temperature up to 200 m with MLD, ILD, BLT and D23 (23° isotherm) for selected cyclones occurred over BoB. Figure 4 (left) represents the subsurface Chl-a concentration up to 200 m driven by the same cyclones. Prior to the passage of cyclones, the profiles represent typical hydrographic state of the oceans with warm waters near the surface and cold waters in the subsurface. The MLD was shallower about 20 m and the ocean was warm from September to December in 2013 during the passage of cyclone Phailin. The other cyclones occurred in 2013 were Helen, Lehar and Madi. Similar situation was observed in September–December of 2014 for Hudhud, but a colder and deeper MLD is observed in September–December of 2016 for Vardah. These are consistent with the climatological oceanic characteristics observed in BoB during the post-monsoon seasons. The time-depth cross-section of Chl-a reveals that the Chl-a concentration remains small in the surface, but about 0.8–1.0 mg/m3 at 40–60 m for the cyclone Phailin. The Chl-a concentration is about 3 mg/m3 for the cyclone Hudhud and about 1.5 mg/m3 for Vardah.
Fig. 4: Bio-Argo measurements.

Temporal evolution (right panel) of depth-time section up to 200 m of temperature of some selected cyclones in the BoB. The MLD (cyan), ILD (blue), BLT (green), and D23 (black) are also indicated in the figure. The vertical black lines indicate the cyclone period. The subsurface Chl-a concentration (left panel) up to 200 m for some selected cyclones in BoB. The vertical black lines indicate the cyclone period.

Full size image

To identify the differential oceanic response of the cyclones at the Argo float locations, we further examined the presence of eddies that play a major role in regulating the physical and biogeochemical processes. The analysis of 7-day SSHA composite before and during the cyclone period at the location of Argo float shows the presence of cold-core (negative SSHA) eddies before the passage of cyclone Phailin, Madi and Hudhud, whereas a warm-core (positive SSHA) eddy prior to the passage of Vardah (Fig. 5). The lower TS and a cold-core eddy during the cyclone Hudhud, and higher TS and a warm-core eddy during the cyclone Vardah produce contrasting oceanic response43. The temperature measurements during the periods of Hudhud and Vardah exhibit comparable response to cold and warm-core eddies, respectively. Another feature of cold-core eddies is trapping the near inertial oscillations in the mixed layer50, which accelerates the entrainment at the bottom of mixed layer and vertical shear as observed during the period of Hudhud. Conversely, the proximity of warm-core eddies triggers rapid vertical dispersion of near inertial energy, which suppresses the mixing and shear as for Vardah50.
Fig. 5: Eddies and primary productivity.

The 7-day composite of sea level anomaly (m) before and during the cyclones in BoB. The black solid lines represent the cyclone tracks, the stars represent the genesis location of the cyclone and the box represents the Argo float location.

Full size image

Vardah was a category 1 cyclone that occurred during 6–13 December 2016. The Chl-a amount before the passage of cyclone was about 0.23 mg/m3, but it escalated to 1.92 mg/m3 in 5 days after the passage of cyclone. The bloom continued to exist for the next 5 days as shown in Fig. 3. Unlike the pre-monsoon cases, for which the Chl-a is restored back to open ocean values in 10 days after the passage of cyclones, the bloom continued to persist even longer periods for the post-monsoon cases. The changes in Chl-a concentrations before and after the passage of cyclones in all three cases are greater than 0.2 mg/m3 and are higher for the lower category tropical storms. These analyses are consistent with the frequent occurrence of cyclones over the south east BoB during this season, as shown in Supplementary Table 2. It is also compelling to note that higher intensity cyclones occur over the north as compared to south BoB, which may be due to the presence of BL in the northern BoB as BL does not exist or insignificantly shallow in the south BoB in any season. The barrier layer in turn favours intensification of tropical cyclones whereas the absence of BL favours the storm-induced upwelling that eventually makes the Chl-a blooms45. Note that the stratification is very strong in northern BoB due to the river water input there51.
Supplementary Table 2 and Fig. 2 also show the results of Chl–bloom events in the post-monsoon seasons in 1997–2019. For instance, the Chl-a increased from 0.53 to 1.13 mg/m3 in the northern and southwestern BoB after the super cyclone of 1999 (25 October–3 November), as also shown by Madhu et al.52. Similar enhancements in Chl-a are estimated for BOB08 (1997), BOB05 (2000) and Madi (2013), about 0.6–3 mg/m3, depending on the cyclones. Although the cyclones BOB06 (1999), Sidr (2007), Giri (2010) and Phailin (2013) were category 4 or 5 cyclones, the high TS and lower EPV did not magnify the Chl-a concentrations to the level of bloom initiated by other cyclones. The enhancement of Chl-a estimated in our study during the cyclone Phailin is in agreement with the reported value of 0.9 mg/m3 for the period 16–24 October 2013 by Vidya et al.21. They also computed the bloom associated with Thane, about 0.7 mg/m3 in the post-cyclone period (1–8 January 2012) at 10°–13° N and 82°–86° E. Nevertheless, we have estimated about 1.8 mg/m3 for the post-cyclone period for Thane. The cyclone Hudhud produced a Chl-a bloom of up to 2.8 mg/m3 in 8–15 October 2014 along the track, as analysed by Chacko27 using the MODIS data, which is very close to our estimate of 2.8 mg/m3. We find similar enhancements in Chl-a that reported by Rao et al.26 for BOB05 in 16–23 November 2000, about 1.2 mg/m3 as deduced from the MODIS data. The other cyclones show moderate bloom values, below 1 mg/m3 as listed in Supplementary Table 2. Nevertheless, the low intensity cyclones such as BOB08 (1997) and Thane (2011) exhibit notable increment in Chl-a following the passage of cyclone, about 1.25–1.8 mg/m3, which is in agreement with their comparatively lower TS and higher EPV during the cyclone period. It also attests the impact and significance of TS in deciding the amplitude of phytoplankton bloom; suggesting sustained low intensity winds trigger strong upwelling to cause intense bloom events.
The Arabian Sea cyclones
In Arabian Sea, about 18 out of 33 cyclones are identified as phytoplankton bloom events (55%) during the study period 1997–2019, in which one occurred in pre-monsoon, three in monsoon and nine in post-monsoon seasons. Since the frequency of occurrences is very small, we have not separated the analyses into seasons or regions of landfall, but a selected case is presented in Supplementary Fig 2. For a better understanding of the behaviour of cyclones, we have selected three cyclones, ARB01 (2001), Mukda (2006) and Megh (2015), one in each season for this discussion. Supplementary Fig 3 illustrates the spatial distribution of Chl-a superimposed with SSHA for the selected cyclones. The ARB01 (2001) was a category 3 cyclone that occurred during 21–28 May 2001 with a MSW of about 60 m/s. The surface Chl-a concentration before the cyclone appearance was below 0.2 mg/m3 due to the profound heating in May53,54. It was one of the strongest cyclones appeared over AS, but measurements were sparse during the cyclone period, and thus, only a small area of about 1 mg/m3 is observed at 16°–17° N, 69°–72° E after passage of the cyclone. Analysis of IRS–P4 measurements by Subramanyam et al.29 found a very large bloom of about 5–8 mg/m3 at 17° N, from 67° E to 71° E. However, our analyses show the bloom of about 2.07 mg/m3 in the region 67°–68° E, 16°–17° N. The difference in bloom values could be due to the difference in datasets, region and period of analyses. As found in the case of BoB, there are closed contours of negative SSHA, which strengthens the observed cyclone-induced and eddy-mediated phytoplankton bloom.
Mukda was a tropical storm that occurred during 21–24 September 2006 with a MSW of 28 m/s. The Chl-a along the right side of the track was above 0.5 mg/m3 even before the cyclone period. The cyclone Megh was considered as the worst to hit Yemen and it occurred just after the passage of another cyclone Chapala over the same region. Megh was a category 3 cyclone with a MSW of 57 m/s. The Chl-a was about 0.7 mg/m3 in the post-cyclone stage, but a small region of about 1.0 mg/m3 was also observed at the right end of the track. Table 2 and Fig. 6 show the analysis of phytoplankton bloom occurrences in AS during the study period (1997–2019). It shows higher Chl-a concentrations in connection with the cyclones ARB01 and ARB02 in 2001, Mukda in 2006 and Megh in 2015, and are consistent with their lower TS. The situation in 2015 was also similar, in which the Category 4 cyclone Chapala showed higher bloom than that of the category 3 cyclone Megh. Similarly, the Category 4 cyclone Kyaar triggered higher bloom than that of the Category 3 cyclone Maha in 2019.
Table 2 The cyclone-induced Chl-a blooms in the post-monsoon seasons since 1997 in AS.
Full size table

Fig. 6: The Chl-a blooms associated with cyclones in Arabian Sea.

The observed enhancement in Chl-a (mg/m3) following the cyclone passage (5-day average) over AS. The translational speed of the closest track points where the bloom occurred in (m/s), the ratio of wind speed to translational speed (WS/TS) and EPV when cyclone reached its maximum intensity (m/s) are also shown.

Full size image

Supplementary Fig. 3 shows the time evolution of biological and physical observations from the Argo float (WMO ID 2902120) in the period of TCs Nilofar, Chapala and Megh. The distance between float location and nearest track point is provided in Supplementary Table 1. The temperature profiles show warm waters near the surface and cold waters in the subsurface before the cyclone passage, as for a typical oceanic state (Supplementary Fig. 4). The MLD is shallower about 20 m and the ocean is warm throughout September–December 2014 and October–December 2015 during the passage of cyclones Nilofar, Megh and Chapala. The time-depth cross-section of chlorophyll shows about 0.8–1 mg/m3 at 40–60 m. However, the Chl-a values remain about 1 mg/m3 for cyclones Megh and Chapala close to the surface; supporting the satellite measurements. The analysis of 7-day SSHA composite shows a cold-core eddy before the passage of cyclone Nilofar whereas warm-core eddies before the passage of Chapala and Megh at the buoy location (Supplementary Fig. 5). The warm-core eddies before the passage of Chapala and Megh could also be the reason for their rapid intensification.
Tropical storms and category 1 cyclones
Although a number of cyclones occurred over NIO during the 1997–2019 period, the phytoplankton bloom happened mostly for the storms and lower category cyclones. For instance, there were five cyclones that appeared over BoB in the pre-monsoon seasons that triggered phytoplankton bloom, four (80%) of them were either tropical storms or category 1 cyclones (Table 1). Similarly, out of 25 cyclones that made Chl-a blooms in BoB during the post-monsoon seasons, 20 of them were (80%) either tropical storms or category 1 cyclones. An analogues occurrence of cyclone-induced phytoplankton bloom is observed for the lower category cyclones in AS, where eight cyclones out of 18 (44.4%) were either tropical storms or category 1 cyclones. These analyses suggest that the slow-moving storms stay more time over the oceans and impart high momentum to upwell the subsurface nutrient-rich water, leading to the phytoplankton blooms in the open oceans with a time lag of 4–12 days, as illustrated in Fig. 7 (blue coloured bar chart). The bloom is as higher as about 20–500% with respect to the pre-cyclone Chl-a levels, and is even up to 1385% as for the case of BoB01 in 2003 and 3758% for the cyclone Gonu in AS (Supplementary Fig. 6); demonstrating the impact and scale of cyclone-induced primary productivity in the open oceans. This slow-moving cyclone-induced primary productivity is very important in the context of climate change, as there is a global slowdown in the translational speed of tropical cyclones.
Fig. 7: The change in cyclone-induced bloom and time lag.

Left: with respect to pre-cyclone Chl-a values (blue), 0.5 mg/m3 (dark blue) and 0.2 mg/m3 (magenta) for the post-monsoon cyclones in Bay of Bengal (BoB). Right: The time lag in days in cyclone-induced Chl-a bloom with respect to the pre-cyclone Chl-a (blue histograms, left) for the post-monsoon cyclones of BoB.

Full size image

To test robustness of the estimates of cyclone-induced change in Chl-a (i.e. Fig. 7), we also considered two other background Chl-a values (i.e. 0.2 and 0.5 mg/m3), which were also taken as the background Chl-a of the ocean basins and the Chl-a threshold for bloom detection. Since the value extracted from the 1° × 1° latitude-longitude region at the track (i.e. blue diagrams) is different from the basin average and bloom threshold values, there are significant differences in the amplitude of blooms, as displayed in Fig. 7. It shows that the pre-cyclone Chl-a values are between 0.5 and 0.2 mg/m3. Therefore, the change in Chl-a is higher with the estimates based on 0.2 mg/m3 (magenta histogram) and about 10 cyclones show a change in Chl-a of about 400%. The highest bloom of about 800–850% is found for Madi (2013), Hudhud (2014) and BOB05 (1999). On the other hand, the change in Chl-a with respect to 0.5 mg/m3 (dark blue histogram) is lower than that with the pre-cyclone estimates, and the change is mostly within 250%, although few cyclones show around 400%. The highest bloom is observed for the cyclone Madi (2013), about 300%. In AS, the Chl-a bloom is mostly between 300 and 1000%, except for Gonu in 2007. The change in Chl-a is about 6000% with respect to the basin average of 0.2 mg/m3, and about 3500% based on the bloom threshold for the cyclone Gonu. The assessment confirm that the cyclone-induced bloom (change in percent) in AS is about five times higher than that of BoB.
The impact of ENSO and IOD on Chl-a blooms
Several studies have examined the relationship between El Niño and Southern Oscillation (ENSO) and cyclone activity across different oceanic basins12,55,56. The influence of ENSO on tropical cyclone activity in BoB during the period 1997–2010 is also investigated by Girishkumar et al.19. We have chosen the dates after the cyclone passage, and considered the Niño and IOD indices to classify the cyclones occurred in the El Niño, La Niña, normal, positive IOD (PIOD) and negative IOD (NIOD) years, as listed in Table 1, Table 2 and Supplementary Table 2. In addition, we have prepared the composites of Chl-a and SSHA for 10 days before and after the passage of each cyclone to assess the inter-annual variability in Chl-a and SSHA with respect to the El Niño, La Niña, normal, PIOD and NIOD years (Fig. 8, for BoB). Out of the 25 cyclones, three of them occurred in El Niño, fourteen in La Niña, nine in normal, four in PIOD and five in NIOD years. The cyclones those occurred in PIOD or NIOD years also happened to be in the El Niño/La Niña years and therefore included in both analyses, and are shown in the figure. The number of cyclones are more in the La Niña years, which were mostly followed by the normal, PIOD and NIOD years. The magnitude of phytoplankton bloom is higher in the PIOD years than that in the NIOD years. In the El Niño years, the magnitude of bloom is comparatively smaller and the bloom in normal years is around 0.5 mg/m3.
Fig. 8: The differences in Chl-a during El Niño, La Niña, normal, PIOD and NIOD years.

The SSHA a 10-day before and b during the passage of cyclone) and Chl-a composite maps with respect to El Niño, La Niña, normal, PIOD and NIOD years in the BoB. The respective cyclone tracks are also shown.

Full size image

Supplementary Fig. 9 shows the composite of SSHA and Chl-a with respect to El Niño, La Niña, normal, PIOD and NIOD years in AS. Here, more number of cyclones occurred in the El Niño years as compared to that in the La Niña and normal years. In contrast, there are more number of cyclones in the PIOD years than the NIOD years in AS, but the amplitude of bloom is higher for the NIOD years. These are also the reasons for the differences in phytoplankton bloom in AS and BoB, as the impact of ENSO and IOD events is different in both basins. The response of cyclones in IOD years are similar to those in the La Niña years. Although the spatial extent of bloom is larger in the El Niño years owing to the higher number of cyclone occurrences, the magnitude of bloom is higher for the cyclones occurred in La Niña and NIOD years. The normal years exhibit bloom similar to that of the El Niño years. In BoB, the analysis of SSHA composite for the El Niño, La Niña, IOD and normal years is dominated by negative SSHA (suggesting the presence of cold-core eddies), but the normal years are more influenced by warm-core eddies (e.g. Fig. 8). In AS, on the other hand, the normal, La Niña and IOD years are dominated by cold-core eddies, whereas the El Niño years are overwhelmed by warm-core eddies (e.g. Supplementary Fig. 9). The influence of IOD is higher than that of ENSO, which is one of the reasons for the inter-annual variability of phytoplankton blooms. The characteristics of phytoplankton blooms in AS and BoB are in contrast with the differences in SST in IOD years in both basins, and this feature is also found with the cyclone-induced blooms. There are noticeable difference in Chl-a concentrations among the normal and El Niño, La Niña, PIOD or NIOD years, and are exhibited in Supplementary Figs. 7 and 8. More

1.
Ceballos, G. et al. Accelerated modern human–induced species losses: Entering the sixth mass extinction. Sci. Adv. 1, e1400253 (2015).
ADS  PubMed  PubMed Central  Article  Google Scholar 
2.
Cardillo, M. et al. Human population density and extinction risk in the world ’ s carnivores. PLoS Biol. 2, 909–914 (2004).
CAS  Article  Google Scholar 

3.
Pimm, S. L. et al. Can we defy nature’s end ?. Science 293, 2207–2208 (2001).
CAS  PubMed  Article  Google Scholar 

4.
Ripple, W. J. et al. Bushmeat hunting and extinction risk to the world’s mammals. R. Soc. Open Sci. 3, 160498 (2016).
ADS  PubMed  PubMed Central  Article  Google Scholar 

5.
Chapron, G. et al. Recovery of large carnivores in Europe’s modern human-dominated landscapes. Science 346, 1517–1519 (2014).
ADS  CAS  PubMed  Article  Google Scholar 

6.
Owen-Smith, N. Pleistocene extinctions: the pivotal role of megaherbivores. Paleobiology 13, 351–362 (1987).
Article  Google Scholar 

7.
Pradhan, N. M. B. & Wegge, P. Dry season habitat selection by a recolonizing population of Asian elephants Elephas maximus in lowland Nepal. Acta Theriol. (Warsz) 52, 205–214 (2007).
Article  Google Scholar 

8.
Hayward, M. W. et al. The reintroduction of large carnivores to the Eastern Cape South Africa: an assessment. Oryx 41(205), 214 (2007).
Google Scholar 

9.
Owen-Smith, N. Megaherbivores: the influence of very large body size on ecology. Trends Ecol. Evol. https://doi.org/10.1111/j.1523-1739.1989.tb00246.x (1989).
Article  Google Scholar 

10.
Karki, J. B., Jhala, Y. V. & Khanna, P. P. Grazing lawns in Terai Grasslands, Royal Bardia National Park, Nepal1. Biotropica 32, 423–429 (2000).
Article  Google Scholar 

11.
McNaughton, S. J. Serengeti migratory wildebeest: facilitation of energy flow by grazing. Science 191, 92–94 (1976).
ADS  CAS  PubMed  Article  Google Scholar 

12.
Skarpe, C. et al. The return of the giants: ecological effects of an increasing elephant population. Ambio 33, 276–282 (2004).
PubMed  Article  Google Scholar 

13.
Foose, T. J. & Van Strien, N. Asian Rhinos—Status Survey and Conservation Action Plan. Vol. 32 (1997).

14.
Subedi, N. et al. Population status, structure and distribution of the greater one-horned rhinoceros Rhinoceros unicornis in Nepal. Oryx 47, 352–360 (2013).
Article  Google Scholar 

15.
Talukdar, B. R. R. C. Asian Rhino specialist group report. Pachyderm 53, 25–27 (2013).
Google Scholar 

16.
Hedges, S., Sagar Baral, H., Timmins, R. & Duckworth, J. Bubalus arnee. IUCNRed List Threat. Species 2008 (2008).

17.
Leader-Williams, N. Fate riding on their horns—and genes?. ORYX 47, 311–312 (2013).
Article  Google Scholar 

18.
Amin, R., Thomas, K., Emslie, R. H., Foose, T. J. & VanStrien, N. An overview of the conservation status of and threats to rhinoceros species in the wild. Int. Zoo Yearb. 40, 96–117 (2006).
Article  Google Scholar 

19.
IUCN Red List Categories and Criteria, Version 3.1, second edition | IUCN Library System. IUCN (2012).

20.
Phillips, S. J. & Dudík, M. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography 31, 161–175 (2008).
Article  Google Scholar 

21.
Lacy, R. C. Vortex: A computer simulation model for population viability analysis. Wildl. Res. 20, 1–13 (1993).
Article  Google Scholar 

22.
IUCN/SSC. Guidelines for Reintroductions and Other Conservation Translocations. Version 1.0. Gland, Switzerland: IUCN Species Survival Commission. Ecologial Applications (2013).

23.
Guillera-Arroita, G. et al. Is my species distribution model fit for purpose? Matching data and models to applications. Glob. Ecol. Biogeogr. 24, 276–292 (2015).
Article  Google Scholar 

24.
Elith, J. & Leathwick, J. R. Species distribution models: ecological explanation and prediction across space and time. Annu. Rev. Ecol. Evol. Syst. 40, 677–697 (2009).
Article  Google Scholar 

25.
Veloz, S. D. Spatially autocorrelated sampling falsely inflates measures of accuracy for presence-only niche models. J. Biogeogr. 36, 2290–2299 (2009).
Article  Google Scholar 

26.
Kramer-Schadt, S. et al. The importance of correcting for sampling bias in MaxEnt species distribution models. Divers. Distrib. 19, 1366–1379 (2013).
Article  Google Scholar 

27.
Jhala, Y., Qureshi, Q. & Gopal, R. Status of Tigers, Copredators and Prey in India, 2014 (2015).

28.
Graham-Rowe, D. Biodiversity: endangered and in demand. Nature 480, S101–S103 (2011).
ADS  CAS  PubMed  Article  Google Scholar 

29.
Frankham, R. Genetic considerations in reintroduction programmes for top-order, terrestrial predators. In Reintroduction of Top-Order Predators 371–387 (Wiley-Blackwell, 2009). https://doi.org/10.1002/9781444312034.ch17

30.
Martin, E., Kumar Talukdar, B. & Vigne, L. Rhino poaching in Assam: challenges and opportunities. Pachyderm (2009).

31.
Rawat, G. S., Goyal, S. P. & Johnsingh, A. J. T. Ecological observations on the grasslands of Corbett Tiger Reserve India. Indian Forester 123(958), 963 (1997).
Google Scholar 

32.
Dinerstein, E. & Schaller, G. The Return of the Unicorns: Natural History and Conservation of Greater—One Horned Rhinoceros (2003). https://doi.org/10.7312/dine08450.

33.
N Subedi 2012 Effect of Mikania micrantha on thrdemography, Habitat use, and Nutrition of Greater one horned rhinoceros in Chitwan National Park Dr. Philos Thesis

34.
Mathur, V.B., Gopal, R., Yadav, S.P., P. R. S. Management Effectiveness Evaluation (MEE) of Tiger Reserves in India: Process and Outcomes 97 (2011).

35.
Amin, R., Thomas, K., Emslie, R. H., Foose, T. J. & Van Strien, N. V. An overview of the conservation status of and threats to rhinoceros species in the wild. Int. Zoo Yearb. 40, 96–117 (2006).
Article  Google Scholar 

36.
Singh, S. P., Sharma, A. & Talukdar, B. K. Translocation of Rhinos within Assam : a successful third round of the second phase of translocations under Indian Rhino Vision (IRV) 2020, 1–6 (2012).
Google Scholar 

37.
Jhala, Y. ., Qureshi, Q., Gopal, R. & Sinha, P. . Status of the Tigers, Co-predators, and Prey in India, 2010. (2011).

38.
Rai, S. After 250 years, rhinos set to make comeback in west UP|Meerut News—Times of India. Times of India (2016).

39.
Subedi, N., Lamichhane, B. R., Amin, R., Jnawali, S. R. & Jhala, Y. V. Demography and viability of the largest population of greater one-horned rhinoceros in Nepal. Glob. Ecol. Conserv. 12, 241–252 (2017).
Article  Google Scholar 

40.
IUCN Standards and Petitions Committee. IUCN Standards and Petitions Committee. Stand. Petitions Comm. 1, 1–60 (2019).
Google Scholar 

41.
Cockrill, W. R. The water buffalo: a review. Br. Vet. J. 137, 8–10 (1981).
CAS  PubMed  Article  Google Scholar 

42.
Kumar, S. et al. Mitochondrial DNA analyses of Indian water buffalo support a distinct genetic origin of river and swamp buffalo. Anim. Genet. 38, 227–232 (2007).
CAS  PubMed  Article  Google Scholar 

43.
Aiyadurai, A. Wildlife hunting and conservation in Northeast India—a need for an interdisciplinary understanding.pdf. Int. J. Gall. Conserv. 2, 61–73 (2011).
Google Scholar 

44.
Jhala, H. Y., Pokheral, C. P. & Subedi, N. Well being and conservation awareness of communities around Chitwan National Park, Nepal|Jhala|Indian Forester. Indian For. 145, 114–120 (2019).
Google Scholar 

45.
Heinen, J. T. & Kandel, R. Threats to a small population: a census and conservation recommendations for wild buffalo Bubalus arnee in Nepal. Oryx 40, 324–330 (2006).
Article  Google Scholar 

46.
Choudhury, A. The decline of the wild water buffalo in north-east India. Oryx 28, 70–73 (1994).
Article  Google Scholar 

47.
Magioncalda, W. A modern insurgency: India’s evolving naxalite problem. S. Asia Monitor 140 (2010).

48.
Melles, S. J., Fortin, M.-J., Lindsay, K. & Badzinki, D. Expanding northward: influence of climate change, forest connectivity, and population processes on a threatened species’ range shift. Glob. Chang. Biol. 17, 17–31 (2011).
ADS  Article  Google Scholar 

49.
Challender, D. W. S. & MacMillan, D. C. Poaching is more than an enforcement problem. Conserv. Lett. 7, 484–494 (2014).
Article  Google Scholar 

50.
Hirzel, A. H., Hausser, J., Chessel, D. & Perrin, N. Ecological-niche factor analysis: How to compute habitat-suitability maps without absence data?. Ecology 83, 2027–2036 (2002).
Article  Google Scholar 

51.
Stockman, A. K., Beamer, D. A. & Bond, J. E. An evaluation of a GARP model as an approach to predicting the spatial distribution of non-vagile invertebrate species. Divers. Distrib. 12, 81–89 (2006).
Article  Google Scholar 

52.
Merow, C., Smith, M. J. & Silander, J. A. A practical guide to MaxEnt for modeling species’ distributions: what it does, and why inputs and settings matter. Ecography (Cop.) 36, 1058–1069 (2013).
Article  Google Scholar 

53.
Phillips, S. J. & Dudı, M. Modeling of species distributions with Maxent : new extensions and a comprehensive evaluation. Ecograohy 31, 161–175 (2008).
Article  Google Scholar 

54.
Elith, J. et al. A statistical explanation of MaxEnt for ecologists. Divers. Distrib. 17, 43–57 (2011).
Article  Google Scholar 

55.
Phillips, S. J., Anderson, R. P., Dudík, M., Schapire, R. E. & Blair, M. E. Opening the black box: an open-source release of Maxent. Ecography (Cop.) 40, 887–893 (2017).
Article  Google Scholar 

56.
Congalton, R. G. A review of assessing the accuracy of classifications of remotely sensed data. Remote Sens. Environ. 37, 35–46 (1991).
ADS  Article  Google Scholar 

57.
Roy, D. P. et al. Landsat-8: Science and product vision for terrestrial global change research. Remote Sens. Environ. 145, 154–172 (2014).
ADS  Article  Google Scholar 

58.
Dinerstein, E. & Price, L. Demography and habitat use by greater one-horned rhinoceros in Nepal. J. Wildl. Manage. 55, 401 (1991).
Article  Google Scholar 

59.
Bontemps, S. et al. GLOBCOVER 2009 Products Description and Validation Report.

60.
Barsi, J., Lee, K., Kvaran, G., Markham, B. & Pedelty, J. The spectral response of the landsat-8 operational land imager. Remote Sens. 6, 10232–10251 (2014).
ADS  Article  Google Scholar 

61.
Laurie, A. Behavioural ecology of the Greater one-horned rhinoceros (Rhinoceros unicornis). J. Zool. 196, 307–341 (2009).
Article  Google Scholar 

62.
Dnerstein, E. & Price, L. Demography and habitat use by greater one-horned rhinoceros in Nepal author(s): Eric Dinerstein and Lori Price Source. J. Wildl. Manage. 55, 401–411 (1991).
Article  Google Scholar 

63.
Pettorelli, N. et al. The Normalized Difference Vegetation Index (NDVI): unforeseen successes in animal ecology. Climate Res. 46, 15–27 (2011).
ADS  Article  Google Scholar 

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

65.
Norberg, A. et al. A comprehensive evaluation of predictive performance of 33 species distribution models at species and community levels. Ecol. Monogr. https://doi.org/10.1002/ecm.1370 (2019).
Article  Google Scholar 

66.
Peterson, A. T. & Nakazawa, Y. Environmental data sets matter in ecological niche modelling: an example with Solenopsis invicta and Solenopsis richteri. Glob. Ecol. Biogeogr. https://doi.org/10.1111/j.1466-8238.2007.00347.x (2007).
Article  Google Scholar 

67.
Laurie, A. Behavioural ecology of the greater one-horned rhinoceros Rhinoceros unicornis. J. Zool. 196, 307–341 (1982).
Article  Google Scholar 

68.
Kriticos, D. J. et al. CliMond: Global high-resolution historical and future scenario climate surfaces for bioclimatic modelling. Methods Ecol. Evol. 3, 53–64 (2012).
Article  Google Scholar 

69.
Rodríguez, E. et al. An Assessment of the SRTM Topographic Products.

70.
Watson, J. E. M., Whittaker, R. J. & Dawson, T. P. Habitat structure and proximity to forest edge affect the abundance and distribution of forest-dependent birds in tropical coastal forests of southeastern Madagascar. Biol. Conserv. 120, 311–327 (2004).
Article  Google Scholar 

71.
Society WC. Last of the Wild Project, Version 2, 2005 (LWP-2): Global Human Influence Index (HII) Dataset (Geographic) Society. Conserv Wildl https://doi.org/10.7927/H4BP00QC (2005).
Article  Google Scholar 

72.
Phillips, S. J. et al. Sample selection bias and presence-only distribution models: Implications for background and pseudo-absence data. Ecol. Appl. 19, 181–197 (2009).
PubMed  Article  Google Scholar 

73.
Jiménez-Valverde, A. Insights into the area under the receiver operating characteristic curve (AUC) as a discrimination measure in species distribution modelling. Glob. Ecol. Biogeogr. 21, 498–507 (2011).
Article  Google Scholar 

74.
Allouche, O., Tsoar, A. & Kadmon, R. Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). J. Appl. Ecol. 43, 1223–1232 (2006).
Article  Google Scholar 

75.
Shabani, F., Kumar, L. & Ahmadi, M. Assessing accuracy methods of species distribution models: AUC, specificity, sensitivity and the true skill statistic. 18, (2018).

76.
Warren, D. L. & Seifert, S. Ecological niche modeling in Maxent: the importance of model complexity and the performance of model selection criteria. Ecol. Soc. Am. 21, 335–342 (2011).
Google Scholar 

77.
Wiltshire, K. H. & Tanner, J. E. Comparing maximum entropy modelling methods to inform aquaculture site selection for novel seaweed species. Ecol. Modell. 429, 109071 (2020).
Article  Google Scholar 

78.
Smeraldo, S. et al. Modelling risks posed by wind turbines and power lines to soaring birds: the black stork (Ciconia nigra) in Italy as a case study. Biodivers. Conserv. 29, 1959–1976 (2020).
Article  Google Scholar 

79.
Puri, K. & Joshi, R. A case study of greater one horned rhinoceros (Rhinoceros unicornis) in Kaziranga National Park of Assam India. NeBIO 9, 307–309 (2018).
Google Scholar 

80.
DNPWC 2015. Rhino Count Report 2015 DNPWC.

81.
Mukherjee, T., Sharma, L. K., Saha, G. K., Thakur, M. & Chandra, K. Past, present and future: combining habitat suitability and future landcover simulation for long-term conservation management of Indian rhino. Sci. Rep. https://doi.org/10.1038/s41598-020-57547-0 (2020).
Article  PubMed  PubMed Central  Google Scholar 

82.
CCF West Bengal. Estimation of Indian Rhinoceros (Rhinoceros unicornis) 2019 West Bengal (2019)

83.
Boyce, M. Population viability analysis: adaptive management for threatened and endangered species. 226–238 (1997).

84.
Khatri, T. B., Shah, D. N. & Mishra, N. Wild Water Buffalo Bubalus arnee in Koshi Tappu Wildlife Reserve, Nepal: status, population and conservation importance. J. Threat. Taxa 04, 3294–3301 (2012).
Article  Google Scholar 

85.
Phillips, S., Anderson, R., Dudík, M., Schapire, R. & Blair, M. Opening the black box: an open-source release of Maxent. Ecography 40, 887–893 (2017). More

1.
Martin, E. A. et al. The interplay of landscape composition and configuration: New pathways to manage functional biodiversity and agroecosystem services across Europe. Ecol. Lett. 22, 1083–1094. https://doi.org/10.1111/ele.13265 (2019).
Article  PubMed  Google Scholar 
2.
Daane, K. M., Vincent, C., Isaacs, R. & Ioriatti, C. Entomological opportunities and challenges for sustainable viticulture in a global market. Annu. Rev. Entomol. 63, 193–214. https://doi.org/10.1146/annurev-ento-010715-023547 (2018).
CAS  Article  PubMed  Google Scholar 

3.
Viers, J. H. et al. Vinecology: Pairing wine with nature. Conserv. Lett. 6, 287–299. https://doi.org/10.1111/conl.12011 (2013).
Article  Google Scholar 

4.
Pertot, I. et al. A critical review of plant protection tools for reducing pesticide use on grapevine and new perspectives for the implementation of IPM in viticulture. Crop Prot. 97, 70–84. https://doi.org/10.1016/j.cropro.2016.11.025 (2017).
ADS  CAS  Article  Google Scholar 

5.
Roehrich, R. & Boller, E. Tortricids in Vineyards. In Tortricid Pests. Their Biology, Natural Enemies and Control (eds Van der Geest, L. P. S. & Evenhuis, H. H.) 507–514 (Elsevier Science Publishers, 1991).

6.
Ioriatti, C., Lucchi, A. & Varela, L. G. Grape berry moths in Western European vineyards and their recent movement into the New World. In Arthropod Management in Vineyards: Pests, Approaches, and Future Directions (eds Bostanian, N. J., Vincent, C. & Isaacs, R.) 339–359 (Springer, 2012).

7.
Stellwaag, F. Die Weinbauinsekten der Kulturländer (Paul Parey, Berlin, 1928).
Google Scholar 

8.
Pavan, F., Girolami, V. & Sacilotto, G. Second generation of grape berry moths, Lobesia botrana (Den. & Schiff.) (Lep., Tortricidae) and Eupoecilia ambiguella (Hb.) (Lep., Cochylidae): Spatial and frequency distributions of larvae, weight loss and economic injury level. J. Appl. Entomol. 122, 361–368. https://doi.org/10.1111/j.1439-0418.1998.tb01513.x (1998).
Article  Google Scholar 

9.
Svobodova, E. et al. Determination of areas with the most significant shift in persistence of pests in Europe under climate change. Pest Manag. Sci. 70, 708–715. https://doi.org/10.1002/ps.3622 (2014).
CAS  Article  PubMed  Google Scholar 

10.
Reineke, A. & Thiéry, D. Grapevine insect pests and their natural enemies in the age of global warming. J. Pest Sci. 89, 313–328. https://doi.org/10.1007/s10340-016-0761-8 (2016).
Article  Google Scholar 

11.
Gutierrez, A. P., Ponti, L., Gilioli, G. & Baumgartner, J. Climate warming effects on grape and grapevine moth (Lobesia botrana) in the Palearctic region. Agric. For. Entomol. 20, 255–271. https://doi.org/10.1111/afe.12256 (2018).
Article  Google Scholar 

12.
Schartel, T. E. et al. Reconstructing the European Grapevine Moth (Lepidoptera: Tortricidae), Invasion in California: Insights from a successful eradication. Ann. Entomol. Soc. Am. 112, 107–117. https://doi.org/10.1093/aesa/say056 (2019).
Article  Google Scholar 

13.
Gutierrez, A. P. et al. Prospective analysis of the invasive potential of the European grapevine moth Lobesia botrana (Den. & Schiff.) in California. Agric. For. Entomol. 14, 225–238. https://doi.org/10.1111/j.1461-9563.2011.00566.x (2012).
Article  Google Scholar 

14.
Gilligan, T. M. et al. Discovery of Lobesia botrana ([Denis and Schiffermüller]) in California: An invasive species new to North America (Lepidoptera: Tortricidae). Proc. Entomol. Soc. Wash. 113, 14–30. https://doi.org/10.4289/0013-8797.113.1.14 (2011).
Article  Google Scholar 

15.
Selkoe, K. A. & Toonen, R. J. Microsatellites for ecologists: A practical guide to using and evaluating microsatellite markers. Ecol. Lett. 9, 615–629 (2006).
Article  Google Scholar 

16.
Reineke, A., Karlovsky, P. & Zebitz, C. P. W. Preparation and purification of DNA from insects for AFLP-analysis. Insect Mol. Biol. 7, 95–99. https://doi.org/10.1046/j.1365-2583.1998.71048.x (1998).
CAS  Article  PubMed  Google Scholar 

17.
Reineke, A. et al. A novel set of microsatellite markers for the European Grapevine Moth Lobesia botrana isolated using next-generation sequencing and their utility for genetic characterization of populations from Europe and the Middle East. Bull. Entomol. Res. 105, 408–416. https://doi.org/10.1017/S0007485315000267 (2015).
CAS  Article  PubMed  Google Scholar 

18.
Schuelke, M. An economic method for the fluorescent labeling of PCR fragments. Nat. Biotechnol. 18, 233–234 (2000).
CAS  Article  Google Scholar 

19.
Kalinowski, S. T., Taper, M. L. & Marshall, T. C. Revising how the computer program CERVUS accommodates genotyping error increases success in paternity assignment. Mol Ecol 16, 1099–1106. https://doi.org/10.1111/j.1365-294X.2007.03089.x (2007).
Article  PubMed  PubMed Central  Google Scholar 

20.
Excoffier, L. & Lischer, H. E. L. Arlequin suite ver. 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567. https://doi.org/10.1111/j.1755-0998.2010.02847.x (2010).
Article  PubMed  PubMed Central  Google Scholar 

21.
Raymond, M. & Rousset, F. GENEPOP (version 1.2): Population genetics software for exact tests and ecumenicism. J. Hered. 86, 248–249 (1995).
Article  Google Scholar 

22.
Goudet, J. FSTAT (Version 1.2): A computer program to calculate F-statistics. J. Hered. 86, 485–486. https://doi.org/10.1093/oxfordjournals.jhered.a111627 (1995).
Article  Google Scholar 

23.
Rice, W. R. Analyzing tables of statistical tests. Evolution 43, 223–225 (1989).
Article  Google Scholar 

24.
Chapuis, M.-P. & Estoup, A. Microsatellite null alleles and estimation of population differentiation. Mol. Biol. Evol. 24, 621–631. https://doi.org/10.1093/molbev/msl191 (2007).
CAS  Article  PubMed  Google Scholar 

25.
Pritchard, J. K., Stephens, M. & Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 155, 945–959 (2000).
CAS  PubMed  PubMed Central  Google Scholar 

26.
Falush, D., Stephens, M. & Pritchard, J. K. Inference of population structure using multilocus genotype data: Linked loci and correlated allele frequencies. Genetics 164, 1567–1587 (2003).
CAS  PubMed  PubMed Central  Google Scholar 

27.
Evanno, G., Regnaut, S. & Goudet, J. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol. Ecol. 14, 2611–2620. https://doi.org/10.1111/j.1365-294X.2005.02553.x (2005).
CAS  Article  PubMed  PubMed Central  Google Scholar 

28.
Earl, D. A. & vonHoldt, B. M. STRUCTURE HARVESTER: A website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361. https://doi.org/10.1007/s12686-011-9548-7 (2012).
Article  Google Scholar 

29.
Littell, R. C., Milliken, G. A., Stroup, W. W., Wolfinger, R. D. & Schabenberger, O. SAS for Mixed Models (SAS Institute, Cary, 2006).
Google Scholar 

30.
Saour, G. Flight ability and dispersal of European grapevine moth gamma-irradiated males (Lepidoptera: Tortricidae). Fla. Entomol. 99, 73–78. https://doi.org/10.1653/024.099.sp110 (2016).
Article  Google Scholar 

31.
Schmitz, V., Roehrich, R. & Stockel, J. Dispersal of marked and released Lobesia botrana in a small isolated vineyard and the effect of synthetic sex pheromone on moth movements. J. Int. Sci. Vigne Vin 30, 67–72 (1996).
CAS  Google Scholar 

32.
Sciarretta, A., Zinni, A., Mazzocchetti, A. & Trematerra, P. Spatial analysis of Lobesia botrana (Lepidoptera: Tortricidae) male population in a Mediterranean agricultural landscape in Central Italy. Environ. Entomol. 37, 382–390. https://doi.org/10.1093/ee/37.2.382 (2008).
CAS  Article  PubMed  Google Scholar 

33.
De Meeûs, T. Revisiting FIS, FST, wahlund effects, and null alleles. J. Hered. 109, 446–456. https://doi.org/10.1093/jhered/esx106 (2018).
Article  PubMed  Google Scholar 

34.
Delbac, L. & Thiéry, D. Damage to grape flowers and berries by Lobesia botrana larvae (Denis & Schiffernüller) (Lepidoptera: Tortricidae), and relation to larval age. Aust. J. Grape Wine R. 22, 256–261. https://doi.org/10.1111/ajgw.12204 (2016).
Article  Google Scholar 

35.
Fuentes-Contreras, E. et al. Measuring local genetic variability in populations of codling moth (Lepidoptera: Tortricidae) across an unmanaged and commercial orchard interface. Environ. Entomol. 43, 520–527. https://doi.org/10.1603/en13131 (2014).
Article  PubMed  Google Scholar 

36.
Dong, Z. K., Li, Y. F. & Zhang, Z. Y. Genetic diversity of melon aphids Aphis gossypii associated with landscape features. Ecol. Evol. 8, 6308–6316. https://doi.org/10.1002/ece3.4181 (2018).
Article  PubMed  PubMed Central  Google Scholar 

37.
Vialatte, A., Dedryver, C. A., Simon, J. C., Galman, M. & Plantegenest, M. Limited genetic exchanges between populations of an insect pest living on uncultivated and related cultivated host plants. P. Roy. Soc. B-Biol. Sci. 272, 1075–1082. https://doi.org/10.1098/rspb.2004.3033 (2005).
Article  Google Scholar 

38.
Manel, S. & Holderegger, R. T. years of landscape genetics. Trends Ecol. Evol. 28, 614–621. https://doi.org/10.1016/j.tree.2013.05.012 (2013).
Article  PubMed  Google Scholar 

39.
Tscharntke, T. et al. When natural habitat fails to enhance biological pest control—Five hypotheses. Biol. Conserv. 204(Part B), 449–458. https://doi.org/10.1016/j.biocon.2016.10.001 (2016).
Article  Google Scholar  More

Jaw shape variation and diet in small mammals
Using 2D geometric morphometrics (Fig. 2a), we found that jaw shape is a good proxy for diet among small extant mammals. In Fig. 3, taxa with negative PC1 scores have shorter jaws, and taxa with positive PC1 scores have longer jaws; taxa with positive PC2 scores have taller ascending rami and taxa with negative PC2 scores have shorter ascending rami. Among extant mammals, most dietary categories (excluding omnivores) can be distinguished along PC1 (Fig. 3a): herbivores plot at the negative end of PC1, insectivores towards the positive end, and carnivores in between. These categories are also statistically different from each other (Table 2), showing that jaw shape can distinguish between most major dietary types. However, our data cannot distinguish between carnivores and omnivores.
Fig. 2: Data acquired from the jaws of Mesozoic and extant small mammals.

a Jaw landmarking regimen used in this study. Modified from ref. 12. In orange: six fixed landmarks; in blue: 58 sliding semi landmarks. b Moment arm measurements taken in this study. Modified from ref. 19.

Full size image

Fig. 3: Scatter plots of the principal component analysis (PCA) results (PC1 vs. PC2).

a Extant taxa, b extinct taxa. Convex hulls shown for extant insectivores (yellow), carnivores (red), omnivores (purple) and herbivores (blue). Icon colors indicate known dietary categories of extant mammals and suggested dietary categories for Mesozoic mammals (obtained from the literature). See Table 1 for taxon names.

Full size image

Table 2 Summary of the Procrustes ANOVA (Type II, Conditional SS) performed for jaw shape data as a function of dietary group.
Full size table

Data on the jaw shape of Mesozoic mammals were projected onto the extant taxa morphospace (Fig. 3b). In order to determine whether jaw shape could be used as a dietary proxy in Mesozoic mammals, we obtained previous independent determinations of likely diets, which variously employed dental morphology, tooth wear facets and body size (e.g., see refs. 1,7,12,14,24,25,26,27,28,29,30,31,32). We saw a very good correspondence between previous proposed diets for Mesozoic mammals and their position on the morphospace. See Supplementary Fig. 6 for a principal components analysis scatter plot which includes multituberculates and haramiyidans; these taxa were excluded from our study because the vast majority of them have jaw shapes dissimilar to the other extinct and extant mammals in our sample (i.e., allotherians have shorter jaws and thus more negative PC1 scores).
Stem mammals
Most stem mammals plot within the morphospace of extant insectivores and have positive PC1 scores. One exception is Sinoconodon (taxon #2, Fig. 3), which plots within the morphospace of extant carnivores; Sinoconodon is considered a carnivore based on dental morphology5. Haramiyavia (#1) is thought to have been a plant-dominated omnivore23 based on dental morphology, but here it plots within the morphospace of extant insectivores. Both morganucodontans in this study, Morganucodon (#3) and Dinnetherium (#4), have similar PC1 scores to extant insectivores, echoing the findings of Gill et al.14.
Molar morphology indicates omnivorous or faunivorous diets for docodontans; here they mostly plot within the morphospace of extant insectivores, with the exception of Haldanodon (#6) and Docofossor (#7). Agilodocodon (#9) was previously considered a plant-dominated omnivore, with exudativorous dental features which indicated a diet mainly composed of plant sap33; more recently, Wible and Burrows34 challenged this hypothesis and suggested that the teeth of Agilodocodon most closely resemble those of extant insectivores. Here, Agilodocodon plots firmly within the morphospace of extant insectivores, close to the insectivorous dusky antechinus (Antechinus swainsonii, #61) and the elephant shrews (Elephantulus rufescens [#114] and E. brachyrhynchus [#115]), which are insect-dominated omnivores.
According to Ji et al.28 the swimming docodontan, Castorocauda (#5), has dental features indicative of feeding on aquatic invertebrates and small vertebrates, like fish. Castorocauda is often depicted as being carnivorous and, particularly, piscivorous7,28,33. The jaw shape of Castorocauda is similar to that of modern day insectivores; this docodontan might have been feeding on “soft” aquatic invertebrates (Fig. 3). The other Mesozoic mammal in our sample proposed to have been semi-aquatic, Teinolophos (#13), plots in a similar area of the morphospace to Castorocauda. Our extant sample also includes a semi-aquatic carnivore, the water opossum (Chironectes minimus, #69), which plots in the middle of the carnivore morphospace, far away from Castorocauda and Teinolophos.
Docofossor (#7) has skeletal features indicative of a fossorial lifestyle and a dentition similar to those of extant mammals foraging underground, such as moles, solenodons, and tenrecs35. This docodontan has previously been considered an insectivore7. Here, Docofossor plots within the morphospace of extant carnivores; however, it plots close to the burrowing Hispaniolan solenodon (Solenodon paradoxus, #109), which has an insectivorous diet. Among the extant insectivores in our sample, the burrowing vermivores (e.g., the hairy-tailed mole, Parascalops breweri [#108], and the Hispaniolan solenodon) have more negative PC1 scores than other insectivores (similar to that of Docofossor), and their PC1 values are more similar to those of carnivores.
The dental morphology of Haldanodon (#6) is indicative of an insectivorous diet. Here, it plots within the carnivore morphospace (very near extant herbivores), because of its tall coronoid process and comparatively shorter jaw. Docodon (#8) likely ate insects and other small invertebrates27 and, based on its diminutive size36, Microdocodon (#10) was probably insectivorous. Both of these docodontans plot within the insectivore morphospace.
Non-therian crown mammals
The jaw shape of non-therian crown mammals varies widely, plotting mostly within the morphospace of insectivores and carnivores. Fruitafossor (#11), a fossorial mammal with teeth similar to extant armadillos, has been considered an omnivore eating insects, small invertebrates and some plants26. Here, it plots within the insectivore morphospace, closely to the insectivorous and fossorial hairy-tailed mole (Parascalops breweri, #108), and shares similar PC1 scores with other fossorial taxa, such as Docofossor (#7) and the Hispaniolan solenodon (#109).
Extant monotremes eat insects and other small invertebrates. It has been proposed that the Early Cretaceous monotreme Teinolophos (#13) had a semiaquatic lifestyle (on the basis of its enlarged mandibular canal37) and ate in a similar manner to the insectivorous Kuehneotherium38. Here Teinolophos, and the australosphenidan Henosferus (#12), have PC1 scores similar to insectivores and omnivores.
The eutriconodontans are a very diverse group of insectivores and carnivores which had a wide range of body sizes, including some of the largest Mesozoic mammals known1. Here all eutriconodontans fall within or very close to the extant carnivore morphospace. In particular, Triconodon (#16) and Argentoconodon (#19) plot within the carnivore morphospace, Trioracodon (#17) and Volaticotherium (#18) plot between the carnivore and insectivore morphospaces, and Yanoconodon (#15) plots within the insectivore morphospace. Both gobiconodontids, Gobiconodon (#20) and Repenomamus (#21), have more negative PC1 scores and plot closer to the herbivore morphospace, but still remain within or close to the carnivore morphospace. Triconodon, Trioracodon, Gobiconodon, and Repenomamus are all considered carnivores based on craniodental morphology and body size1,7,31; additionally, there is direct evidence for the carnivorous diet of Repenomamus from fossilized stomach contents4. Yanoconodon and Volaticotherium are considered insectivores7.
“Symmetrodontans” like Spalacotherium (#22), Zhangheotherium (#24) and Maotherium (#25) have often been considered insectivores based on their craniodental morphology1,7 (note “symmetrodontans” likely do not represent a monophyletic group, but are often grouped together based on their tooth morphology1). Here, all “symmetrodontans” plot within the insectivore morphospace. Dryolestids are also commonly considered insectivorous1,29. Here, Crusafontia (#26) plots between the morphospace of extant carnivores and insectivores, while Amblotherium (#27) plots within the insectivore morphospace.
Vincelestes (#29) has previously been considered a carnivore on the basis of jaw shape12. Here, it plots near the morphospaces of both omnivores and herbivores. Bonaparte24 considered the incisor wear of Vincelestes reminiscent of Cenozoic carnivores, and Rougier25 considered its jaw morphology indicative of a forceful bite enabling the incorporation of tough plant matter into a primarily carnivorous/insectivorous diet.
Therian crown mammals
Extant marsupials have a large diversity of diets, including herbivory, but the extinct metatherians in our sample are considered to have been limited in diet to insectivory and carnivory (note that there are some putatively herbivorous/omnivorous extinct metatherians, like Glasbius and polydolopimorphians39,40). Their jaw shape is very similar to that of extant carnivores and insectivores (Fig. 3). Dental morphology indicates that Eodelphis (#32) and Deltatheridium (#30) were carnivorous, Didelphodon (#31) durophagous or molluscivorous31,32, and Alphadon (#33) is considered to have been insectivorous, on the basis of its jaw shape and body size12. Dental microwear indicates a broad diet consisting of vertebrates, plants, and hard-shelled invertebrates for Didelphodon; biomechanical analyses of its skull and jaw points towards a durophagous diet15,16. Biomechanical analyses of the resistance to bending and torsion of Eodelphis jaws, points to a durophagous diet in Eodelphis cutleri and non-durophagous faunivory for Eodelphis browni16. Here, Eodelphis, Deltatheridium and Didelphodon plot closely to the extant carnivores, while Alphadon plots closely to the extant insectivores.
Extant placentals also have a wide range of diets, but many of the extinct eutherians in this study (i.e., Sinodelphys [#34], Juramaia [#35], Eomaia [#36], Kennalestes [#40], Barunlestes [#44], and Kulbeckia [#43]) are considered insectivorous7,12. Here, we corroborate this hypothesis (Fig. 3): all extinct eutherians plot within the insectivore morphospace, with the exception of Asioryctes (#38) which plots in the insectivore/carnivore morphospace, and Juramaia and Sinodelphys, which plot just outside the insectivore morphospace.
Using jaw shape to infer diet in Mesozoic mammals
We performed a phylogenetic flexible discriminant analysis (phylo FDA) following Motani and Schmitz41 to determine the posterior probability of the Mesozoic taxa belonging to one of three dietary categories: insectivore, carnivore, or herbivore (we omitted omnivores as they are not well discriminated in Fig. 3). We used the first seven PC scores (of the PCA of Procrustes coordinates of jaw shape), which together accounted for 81.39% of the variance. The results of the analysis can be seen in Fig. 4 and the posterior probability values can be seen in Supplementary Data 1. We used the extant taxa of known diets as the training dataset for the discriminant analysis: these taxa were classified correctly 89.19% of the time. For the most part, we see a good separation between dietary groups among extant mammals (Fig. 4a), with some exceptions: the primarily herbivorous olingo (Bassaricyon gabbii, #94) plots with the carnivores (although mainly frugivorous, it can consume small vertebrates), and a couple of insectivores plot very near the carnivores (i.e., the little brown bat [Myotis lucifugus, #104] and the Hispaniolan solenodon [Solenodon paradoxus, #109]). These three taxa, alongside the carnivorous greater bulldog bat (Noctilio leporinus, #101), were the only extant taxa misclassified by the discriminant analysis.
Fig. 4: Phylogenetic flexible discriminant analysis results, showing discriminant axis 1 (DA1) and two (DA2), of all taxa in this study.

Extinct taxa are color coded based on their posterior probability of belonging to one of the established dietary categories. Convex hulls show the position of the extant taxa in the plot and are color coded based on their dietary categories.

Full size image

The Mesozoic mammals included in our sample have largely been considered faunivorous and the results of the phylo FDA (Fig. 4b) corroborate this hypothesis. The majority of them are classified as insectivorous, including most stem mammals, australophenidans, “symmetrodontans” and eutherians, among others. Among the eutriconodontans, Argentoconodon, Gobiconodon, Repenomamus, and Trioracodon, are classified as carnivores, Triconodon and Yanoconodon are classified as insectivores, but with moderate support (posterior probabilities: 48% and 52%, respectively), and Phascolotherium and Volaticotherium are more confidently classified as insectivores (posterior probabilities: 60% and 73%, respectively). Among the metatherians, Didelphodon and Eodelphis are classified as carnivores, while Alphadon and Deltatheridium are classified as insectivores with moderate support (posterior probabilities: 54% and 52%, respectively). The stem mammals, Haramiyavia, Sinoconodon, and Docofossor are all confidently classified as carnivores (posterior probabilities over 80%), and the crown mammals Crusafontia and Kennalestes are also classified as carnivores, but with moderate support (posterior probabilities: 54% and 52%, respectively). Two taxa in the analysis are classified as herbivores, because of their relatively tall ascending rami: Vincelestes (#29) and Haldanodon (#6). The dental morphology of Vincelestes points to a primarily faunivorous diet24, but it has been previously noted that its jaw morphology is indicative of a forceful bite; Rougier25 suggested that this jaw morphology might have enabled Vincelestes to incorporate tough plant matter into its diet, but it might also be indicative of durophagy. The dental morphology27 and body size of Haldanodon point towards an insectivorous diet; in this analysis, the posterior probability of Haldanodon being a herbivore is not high (only 40.3%). The evidence thus far suggests Haldanodon had a faunivorous diet; its jaw morphology might be indicative of the incorporation of tougher food sources into its diet.
Mechanical advantage of the jaws of small mammals
We obtained mechanical advantage (MA) data to test whether extant mammals of different dietary groups have distinct MA values (Table 3). The mechanical advantage measurements were standardized across all jaws to account for differences in jaw morphology (e.g., presence or absence of the angular process) (Fig. 2b); the outlever was measured at the anterior end of the jaw and at the first lower molar (m1). When measuring mechanical advantage at the jaw tip and considering extant taxa only, we find statistically significant differences in the mechanical advantage of the masseter (MAM) values in all pairwise dietary combinations except for carnivore-insectivore (Table 3). The mechanical advantage of the temporalis (MAT) is statistically distinct only between herbivores and insectivores, and carnivores and insectivores (Table 3). Herbivores and carnivores do not have statistically distinct MAT values. This may differ in a sample of larger ( > 5 kg) therians. When measuring the outlever at the m1, we find statistically significant differences in all pairwise comparisons of MAM between dietary groups, except for herbivore–omnivore and carnivore–insectivore. When considering MAT, we only find significant differences between omnivores and carnivores, insectivores and herbivores, and insectivores and carnivores.
Table 3 Pairwise p values (uncorrected significance) of one-way PERMANOVAs of the mechanical advantage values of the masseter (MAM) and temporalis (MAT) obtained in this study on extant taxa of known dietary preferences only (permutation N = 9999).
Full size table

Figure 5 shows the mechanical advantage of the masseter (left) and temporalis (right), measured at the jaw tip, in a phylogenetic context (see also Supplementary Fig. 7 for individual taxon names). Phylogeny appears to have a large influence on the mechanical advantage and diet of the jaws of small mammals. Most Mesozoic taxa have low (blue) to intermediate (green) MAM values. Most stem mammals have intermediate (green) to high (red) MAM values and non-therian crown mammals have low MAM values, with the exception of Fruitafossor and Vincelestes (the latter has the highest MAM value of all taxa, both extinct and extant). Most eutherians, both extinct and extant, have intermediate to low MAM values, with the exception of the relatively high values (yellow to orange) seen in elephant shrews (order Macroscelidea) and the four-toed hedgehog (order Eulipotyphla, Atelerix albiventris). Some members of the orders Carnivora (including canids and euplerids) and Afrosoricida have some of the lowest MAM values. Metatherians have MAM values ranging from low to intermediate (in the orders Dasyuromorphia and Didelphimorphia, as well as in the Mesozoic metatherians) to some of the highest in the order Diprotodontia (e.g., the sugar glider [Petaurus breviceps], the woylie [Bettongia penicillata], the cuscus [Phalanger orientalis]).
Fig. 5: Mechanical advantage values of the masseter and temporalis when biting at the jaw tip visualized in the context of the phylogeny used in this study.

See Supplementary Fig. 7 for individual taxon names.

Full size image

Most taxa have intermediate MAT values (Fig. 5 and Supplementary Fig. 7). Very low MAT values are seen in the extinct non-therian crown mammals Teinolophos and Zhangheotherium and a few extant taxa, including marsupials like the Western barred bandicoot (Perameles bougainville) and the numbat (Myrmecobius fasciatus), and placentals such as the striped treeshrew (Tupaia dorsalis) and the short-snouted elephant shrew (Elephantulus brachyrhynchus). The highest MAT values belong to members of the order Carnivora, including skunks (Mephitis macroura and Conepatus humboldtii), the least weasel (Mustela nivalis) and the tayra (Eira barbara). Some diprotodontians like the common ringtail possum (Pseudocheirus peregrinus) and the sugar glider (Petaurus breviceps) also have relatively high MAT values. Some extinct taxa also have relatively high MAT values, including the stem mammal Docofossor, and the non-therian crown mammals, Triconodon and Vincelestes.
Figures 6 and 7 present a visualisation of the mechanical advantage of the masseter and the temporalis (x axis, outlever measured at the jaw tip) and the PC1 scores of Fig. 3 (y axis, mainly describes jaw length) because, as previously mentioned, this is the axis in which dietary categories among extant mammals are best discriminated. In the y axis of Figs. 6a and 7a, herbivores have short jaws, carnivores have short to intermediate-length jaws and insectivores have intermediate-length to long jaws. In Fig. 6a, insectivores and carnivores have low mechanical advantage values of the masseter (i.e., when biting: less forcefulness, more speed), and herbivores have higher mechanical advantage values (i.e., when biting: more forcefulness, less speed). In Fig. 7a, insectivores have lower mechanical advantage values of the temporalis, while carnivores and herbivores have higher mechanical advantage values. Note that most carnivores have intermediate MAT values, but some mustelids (i.e., the least weasel [Mustela nivalis, #99], the American badger [Taxidea taxus, #96], and the North American river otter [Lontra canadensis, #98]), have the highest MAT values among extant mammals. Also note that, among insectivores, those with the highest MAT values are burrowing vermivores (i.e., the short-tailed shrew tenrec [Microgale brevicaudata, #111], the hairy-tailed mole [Parascalops breweri, #108], and the Hispaniolan solenodon [Solenodon paradoxus, #109]). By using a combination of their MAM and MAT values (as well as their jaw length), we can distinguish dietary categories among extant mammals. We decided to omit omnivores from these figures because, as seen in Fig. 3, they cannot be distinguished from other dietary groups on the basis of jaw shape.
Fig. 6: Scatter plot of the mechanical advantage of the masseter (x axis) vs. PC1 scores from Fig. 3 (y axis), which mainly describes jaw length.

a Extant taxa, b extinct taxa. Colors indicate known dietary categories of extant mammals and suggested dietary categories for Mesozoic mammals (obtained from the literature). Ovals indicate where extant taxa of known dietary categories plot, as in part a.

Full size image

Fig. 7: Scatter plot of the mechanical advantage of the temporalis (x axis) vs. PC1 scores from Fig. 3 (y axis), which mainly describes jaw length.

a Extant taxa, b extinct taxa. Colors indicate known dietary categories of extant mammals and suggested dietary categories for Mesozoic mammals (obtained from the literature). Ovals indicate where extant taxa of known dietary categories plot, as in part a.

Full size image

We also obtained additional mechanical advantage measurements, in which the outlever was measured at the first lower molar (m1), rather than the jaw tip (Supplementary Figs. 8, 10, 11, and 13). We made this alternative measurement because Grossnickle17 found that the length of the posterior portion of the jaw (measured from the jaw joint to the m1) is a strong predictor of diet in mammals. Compared to the mechanical advantage (MA) measurements at the jaw tip (Figs. 6a and 7a), we see a less clear distinction between dietary groups among extant mammals. There is considerable overlap between dietary groups in Supplementary Fig. 10 (jaw length~MAM). In Supplementary Fig. 11 (jaw length~MAT), there is a better separation between dietary groups.
Based on previous likely determinations of diet of Mesozoic mammals (see Supplementary Data 1 for the full list of sources), most taxa plot where it is “expected” of them, with some exceptions (Figs. 6b and 7b): 1) about half of the stem mammals (i.e., Haramiyavia, Sinoconodon, Morganucodon, Haldanodon, and Docofossor), most of which are thought to have been faunivorous, have higher MAM values than modern insectivores and carnivores, and 2) the docodontan Castorocauda has MAM and MAT values consistent with an insectivorous diet, as opposed to the carnivorous diet proposed for this taxon7,28,33. Most Mesozoic mammals have mechanical advantage values similar to modern insectivores, a few taxa are similar to carnivores (e.g., Sinoconodon, Triconodon, Trioracodon, Argentoconodon, Gobiconodon, Repenomamus, Deltatheridium, Didelphodon, and Eodelphis), and some are more similar to herbivores (e.g., Vincelestes and Fruitafossor). More

Model equations
We introduce a spatiotemporal model to describe the bacteria-phage interaction in the upper part of the soil with the depth H (we consider (H=1) m) in a typical agricultural field. Here we consider a 1D model where all abiotic and biotic components depend on time t and vertical coordinate h. The biotic component of the model consists of 4 compartments: phage-free bacteria (S) susceptible to infection by the phage, bacteria infected by the phage in its lysogenic ((I_1)) and lytic ((I_2)) states, and free phages (P). The total density of the host bacterial populations N is defined as (N = S + I_1 + I_2). The schematic diagram illustrating bacteria-phage interactions is similar to that of Egilmez and co-authors16. The local species interactions are described based on the classical modelling approach6,19. Our spatiotemporal model is of reaction-diffusion type and is described by the following equations

$$begin{aligned} begin{aligned} frac{partial S(t,h)}{partial t}&= D_b frac{partial ^2 S(t,h)}{partial h^2} +alpha (T) S(t,h) Big [1-frac{N(t,h)}{C(h)}Big ] – K_S S(t,h)P(t,h), \ frac{partial I_1(t,h)}{partial t}&= D_b frac{partial ^2 I_1(t,h)}{partial h^2} + {overline{alpha }}(T) I_1(t,h) Big [1- frac{N(t,h)}{C(h)}Big ] + K_1(T) S(t,h) P(t,h) – lambda _1(T) I_1(t,h), \ frac{partial I_2(t,h)}{partial t}&= D_b frac{partial ^2 I_2(t,h)}{partial h^2} + K_2(T) S(t,h) P(t,h) + lambda _1(T) I_1(t,h) – lambda _2 I_2(t,h), \ frac{partial P(t,h)}{partial t}&= D_P frac{partial ^2 P(t,h)}{partial h^2} -K N(t,h) P(t,h) – mu P(t,h) + b lambda _2 I_2(t,h). end{aligned} end{aligned}$$
(1)

In the above model, we parameterise the growth of susceptible bacteria via a standard logistic growth function6, where (alpha) is the maximal per capita growth rate and C is the carrying capacity of the environment; we assume that C(h) varies with depth. Infection of S by phages P at low temperatures results in lysogeny which is described by a mass action term (K_s S(t,h) P(t,h)). The growth of lysogenic bacteria (I_1) is described by a logistic function as in the case of S; however, with a different maximal growth rate ({overline{alpha }} (T)) as detailed in the next subsection. At high temperatures, the transition from the lysogenic to the lytic cycle of infection occurs: this is described by the term (lambda _1 (T) I_1(t,h)). Infection by the phage via the lytic cycle is modelled by the term (K_2 (T)S(t)P(t)). The death rate of infected bacteria due to lysis is modelled by (lambda _2 (T) I_2). The lysis of a bacterium results in the release of b new phages, the the burst size6. In the equation for P, KN(t)P(t) stands for the loss of phage due to binding to any type of bacteria (for simplicity, we assume that there is no saturation in binding). Finally, (mu P(t,h)) is the natural mortality or deactivation of phages.
According to this framework, the vertical displacement of the phage and bacteria are modelled by a diffusion term (first term in each equation), where (D_b) and (D_P) are the diffusion coefficients of bacteria and phage, respectively. The variation of the temperature T across the soil is described by the heat equation

$$begin{aligned} frac{partial T(t,h)}{partial t} = D_h frac{partial ^2 T(t,h)}{partial h^2}, end{aligned}$$
(2)

where (D_h) is the diffusion coefficient of heat transfer (see more detail in the next section). Models (1)–(2) should be supplied with appropriate boundary conditions. We assume that the model has the zero-flux boundary condition for all biotic components (bacteria and phage) at (h=0) and (h=H). For the temperature, we consider Dirichlet boundary conditions such that (T(t,0)= T_s (t)) and (T(t,H)= T_H), where (T_s (t)) is the surface temperature and (T_H) is a constant temperature in deeper soil layers.
Parameterisation of equation terms
Next we describe the functional forms of the dependence of model parameters on the temperature and the depth. Following the previous study16, we assume that the maximal bacterial growth rates (alpha (T)) and ({overline{alpha }}(T)) are described by

$$begin{aligned} alpha (T)= & {} exp left (-frac{(T-T_0)^2}{2sigma ^2}right )alpha _{text {max}}, end{aligned}$$
(3)

$$begin{aligned} {overline{alpha }}(T)= & {} alpha (T) left [1-frac{T^n}{T_1^n + T^n}right ] = alpha _{text {max}}exp {left (-frac{(T-T_0)^2}{2sigma ^2}right )} left [1-frac{T^n}{T_1^n + T^n}right ], end{aligned}$$
(4)

where (T_0=38.2 ^circ text {C}) is the optimal temperature; (T_1=34.8 ^circ text {C}) is the temperature corresponding to the switch between the lytic and the lysogenic cycles; (alpha _{text{max}}=23 text {day}^{-1}) is the maximal possible growth, (sigma =9.1 ^circ text {C}) describes the decay of growth with temperature T16,20.
In the equation for ({overline{alpha }}(T)), we assume that at a high temperature normal cell division of (I_1) stops since there is a transition to a lytic state in bacteria. In the soil bacteria grow anaerobically or microaerophillically, and the growth rates of B. pseudomallei under such conditions are yet to be studied. For simplicity they are assumed to be the same as under aerobic conditions. Realistic values of the above parameters are listed in Table 1. Note that in the model both (alpha (T)) and ({overline{alpha }}(T)) are in fact effective growth rates of the bacterial populations, i.e. they incorporate the replication of cells and as well as their mortality.
Table 1 Parameters used in the model along with their units and ranges.
Full size table

The overall adsorption rate of the phage K is estimated as (2 times 10^{-7} text {ml}^{-1} text { day}^{-1}) from Egilmez et al.16. The adsorption constants (K_1 (T)), (K_2 (T)) and the transition rate from lysogenic to lytic cycle (lambda _1(T)) depend on temperature as follows16:

$$begin{aligned} K_1(T)&= left(1-frac{T^n}{T_1^n + T^n}right) K_S, nonumber \ K_2(T)&= frac{T^n}{T_1^n + T^n} K_S, nonumber \ lambda _1(T)&= frac{T^n}{T_1^n + T^n} {lambda _1}_{text {max}} , end{aligned}$$
(5)

where (K_S) is the maximal phage adsorption constant ((K_S=epsilon K) where (epsilon =0.3) is the adsorption efficiency) and (lambda _{1_text {max}}=23 text {day}^{-1}) is the maximal transition rate which is assumed to be equal to the maximal growth rate of the bacteria16. The lysis rate of bacteria (lambda _2=20 text {day}^{-1}) (depending on 50 min latency time13) and the burst size (b = 100) in the model are assumed to be constant16. The temperature dependence of (alpha (T)), ({overline{alpha }}(T)), (K_1 (T)), (K_2 (T)) and (lambda _1(T)) are shown in Fig. 2. The mortality rate of phages (mu) is high near the surface due to ultraviolet radiation, but the role of ultraviolet radiation becomes negligible starting from a depth of a few centimetres because sunlight cannot penetrate the soil. For the above reason, we can assume (mu =3 text {day}^{-1}) to be constant.
Figure 2

(a) Temperature dependence of the adsorption constants (K_i) ((i=1,2)) of the phage (measured in (text {ml}^{-1} text {day}^{-1})). (b) Growth rates of susceptible (alpha (T)) and lysogenic ({overline{alpha }}(T)) bacteria and the transition rate (lambda _1(T)) from the lysogenic cycle to the lytic cycle (measured in (text {day}^{-1})). The corresponding analytical expressions for the temperature dependence are given by (3)–(5).

Full size image

The carrying capacity C of the bacteria varies with the depth of the soil, according to empirical observations21,22,23. This can be explained by the fact that the humus, oxygen, nitrogen contents, or/and water content in the soil generally decrease with depth24. We use a combined approach to parameterise C(h) based on the available empirical data. We assume that in the absence of phages, the bacteria achieve numbers close to the carrying capacity at a given depth. Firstly, we parameterise the dependence of the overall bacterial load on depth in paddy soils in Southern Asia using the existing data22. Then we re-scale the obtained curve based on the available observations of B. pseudomallei at a depth (h=30 text {cm})25,26. We approximate C(h) using the following simple Gaussian-type curve

$$begin{aligned} C(h)=(C_text {surf} -C_0)exp (-B h^2)+C_0, end{aligned}$$
(6)

where (C_text {surf}) gives the maximal number of bacteria near the surface (h), B determines how fast the bacterial abundance decreases with depth, (C_0) is background bacterial density which takes into account the fact that bacteria can survive even at large depths (e.g. (h=100 text {cm})). Based on our estimates (see supplementary material SM1 for more detail), we will use the following parameter values as defaults: (C_text {surf} = 1 times 10^6) (text {cell/ml}), (B=7.5 times 10^{-4}) (1/{text{cm}}^2), (C_0=10^4) (text {cell/ml}) . One can easily see that C(h) has a maximum at the surface and monotonically decreases with depth. We assume that the carrying capacity of the environment is not influenced by seasonal variation.
The coefficient (D_h) in the equation for the temperature distribution can be estimated as follows. Generally, (D_h) is related to (rho _s), (C_{rho s}) and (k_s) which are the bulk density, specific heat and thermal conductivity in soil, respectively, i.e. (D_h=k_s/(rho _s C_{rho s)}). We use the estimates for (rho _s), (C_{rho s}) and (k_s) from27 which gives (rho _s=110.52 text {kg}/text{m}^3), (C_{rho s} = 1130) (text {J/kg K}) and (k_s = 0.0967) (text {W/m K}) and, for the diffusion coefficient (D_h=7.7 times 10^{-8}) (text {m}^2 text{s}^{-1}). The variation of (T_s)—the surface temperature—is obtained from the historical weather report for the surface16. The bottom boundary temperature (T_H) at (h=H=1 text {m}) is considered to be (22 ^circ text{C}). The initial value of the temperature distribution (T_s (0)) is assumed to be linear, but this assumption does not affect long-term temperature dynamics.
The paddy fields in which we model the bacteria-phage interactions are flooded lands, where the soil is either mud or muddy water. Many factors can affect vertical dispersal of bacteria and phages in such soil. For instance, rain water can carry bacteria and phage up or down in the soil, which can be mathematically modelled by adding an advection term; however, for simplicity we ignore such effects in this paper. We also assume the phage and bacteria vertical diffusion coefficients to be constant; however, it is rather hard to provide accurate estimates for (D_p) and (D_b). In water, the diffusion coefficient of bacteria and phages can be estimated as (3.6times 10 ^{-10} text {m}^2 text{s}^{-1}= 0.3 text {cm}^2 text{day}^{-1}) and (2.8 times 10^{-12} text {m}^2 text{s}^{-1}= 0.002 text {cm}^2 text{day}^{-1}), respectively28, but the diffusivity in soil should be smaller than this. As such, these values should be considered as upper limits for (D_P) and (D_b), with the actual coefficients being orders of magnitude smaller. We undertook simulations with different combinations of diffusion coefficients in this range, and found that the patterns of vertical distribution do not largely depend on the diffusion coefficients provided (D_P< 10^{-3} text {cm}^2 text{day}^{-1}) and (D_b < 10^{-2} text {cm}^2 text{day}^{-1}), due to the strong external forcing of the system by temperature (see “Results” section for details). In our numerical simulations, we use both explicit and implicit numerical schemes. We take a 0.1 cm spatial step size to get a proper resolution. We separately compute the heat equation to define T(t) with a smaller time resolution and then apply the temperature obtained to model bacteria-phage interactions for a larger time resolution (for example, (Delta t cong 7 times 10^{-5}) day). We compute the average densities of the species (both in terms of spatial and temporal averaging) using a numerical right Riemann sum. The accuracy of our numerical simulation was verified by reducing both time and space steps and comparing the results obtained. We use daily and seasonal variation of temperatures (for the period of 2013–2016) in the provinces of Nakhon Phanom and Sa Kaeo in Thailand to parameterise the model (http://www.worldweatheronline.com). The unit of the densities of bacteria and phages are cells/ml. The summary of model parameters as well their values are provided in Table 1. More

Elements of TA
Identifying biogeographic trends within volatile data required several innovations in the MASTIF model20, building from multivariate state-space methods in previous applications41,52. Standard modeling options, such as generalized linear models and their derivatives, do not accommodate key features of the masting processes. First, multiple data types are not independent. Maturation status is binary with detection error, CCs are non-negative integers, also with detection error, and STs require a transport model (dispersal) linking traps to trees, and identification error in seed identification. Of course, a tree observed to bear seed, now or in the past, is known to be mature now. However, failure to observe seed does not mean that an individual is immature because there are detection errors and failed crop years41,64.
Second, seed production is quasiperiodic within an individual (serial dependence), quasi-synchronous between individuals (“mast years”), [and] there is dependence on environmental variation, and massive variation within and between trees41,53,65. Autoregressive error structures (AR(p) for p lag terms) impose a rigid assumption of dependence that is not consistent with quasiperiodic variation that can drift between dominant cycles within the same individual over time43. It does not allow for individual differences in mast periodicity.
Third, climate variables that affect fecundity operate both through interannual anomalies over time and as [a] geographic variation. The masting literature deals almost exclusively with the former, but our application must identify the latter: the potentially smooth variation of climate effects across regions must be extracted from the many individual time series, each dominated by local “noise.”
Finally, model fitting is controlled by the size classes that dominate a given site and thus is insensitive to size classes that are poorly represented. Large trees are relatively rare in eastern forests, making it hard to identify potential declines in large, old individuals41,53. Conversely, the shade-intolerant species that dominate second-growth forests often lack the smaller size classes needed to estimate maturation and early stages where fecundity may be increasing rapidly.
Several of the foregoing challenges are resolved in the MASTIF model by introducing latent states for individual maturation status and tree-year seed production. The dependent data types (maturation status, CCs, STs) become conditionally independent in the hierarchical MASTIF model (e.g., ref. 66). The serial dependence is handled as a conditional hidden Markov process for maturation that combines with CCs and STs by way of stochastic (latent) conditional fecundity. Maturation status and conditional fecundity must be estimated jointly, that is, not with separate models. The latent maturation/fecundity treatment avoids imposing a specific AR(p) structure. In the MASTIF model there is a posterior covariance in maturation/fecundity across all tree-year estimates that need not adhere to any specific assumption20. The dependence across individuals and years is automatic and available from the posterior distribution.
Separating the spatial from temporal components of climate effects is possible here, not only because the entire network is analyzed together but also because predictors in the model include both climate norms for the individual sites and interannual anomalies across sites35,52. TA depends on both of these components.
Extracting the trends from volatile data further benefits from random individual effects for each tree and the combination of climate anomalies and year effects over time. A substantial literature focuses on specific combinations of climate variables that best explain year-to-year fecundity variation, including combinations of temperature, moisture, and water balance during specific seasons over current and previous years19,25,41. Results vary for each study, presumably due to the differences in sites, species, size classes, duration, data type, and modeling assumptions. For TA, the goal is to accommodate the local interannual variation to optimize identification of trends in space and time. Thus, we include a small selection of important climate anomalies (spring minimum T of the current year, summer T of the current and previous year, and moisture D of the current and previous year). The climate anomalies considered here do not include every variable combination that could be important for all size classes of every species on every site. For this reason, we combine climate anomalies with year effects. Year effects in the model are fixed effects within an ecoregion and random between ecoregions (ecoregions are shown in Fig. 2 and listed in Supplementary Data 2). They are fixed within an ecoregion because they are not interpreted as exchangeable and drawn at random from a large population of possible years. They are random between ecoregions due to the uneven distribution of sites (Supplementary Data 1)20.
To optimize inference on size effects, the sampling of coefficients in posterior simulation is implemented as a weighted regression. This means that the contribution of tree diameter to fecundity is inversely proportional to the abundance of that size class in the data. This approach has the effect of balancing the contributions of abundant and rare sizes. Identifying size effects further benefits from the introduction of opportunistic field sampling, which can target the large individuals that are typically absent from field study plots.
MASTIF data network
Data included in the analysis come from published and unpublished sources and offer one or both of the two data types, CCs and STs (Supplementary Data 1). Both data types inform tree-year fecundity; they are plotted by year in Fig. 6.
Fig. 6: Distribution of observation trees by year in the North American region of the MASTIF network.

Sites are listed by ecoregion in the Supplementary Data 2.

Full size image

CCs in the MASTIF network are obtained by one of three methods. Most common are counts with binoculars that are recorded with an estimate of the fraction of the crop that was observed. A second CC method makes use of seeds collected per ground surface area relative to the crown area. This method is used where conspecific crowns are isolated and wind dispersal is limited. The crop fraction is the ratio of ground area for traps relative to the projected crown area. Examples include HNHR67 and BCEF68.
A third CC method is based on evidence for past cone production that is preserved on trees. This has been used for Abies balsamea at western Quebec sites69, Pinus ponderosa in the Rocky Mountains70, and for Pinus edulis at SW sites27.
ST data include observations on individual trees that combine with seed counts from traps. Because individual studies can report different subcategories of seeds, and few conduct rigorous tests of viability, we had to combine them using the closest description to the concept of “viable”. For example, we do not include empty conifer seeds. A dispersion model provides estimates of seeds derived from each tree. ST and CC studies are listed in Supplementary Data 1. The likelihoods for CCs and STs are detailed in ref. 20. Individually and in combination, the two data types provide estimates, with full uncertainty, for fecundity across all tree-years.
Fitted species had multiple years of observations from multiple sites, which included 211,146 trees and 2,566,594 tree-years from 123 species. Sites are shown in Fig. 2 of the main text by ecoregion, they are named in Fig. 1 and summarized in Supplementary Data 1. For TA the fits were applied to 7,723,671 trees on inventory plots. Mean estimates for the genus were used for inventory trees belonging to species for which there were not confident fits in the MASTIF model, which amounted to 7.2% of inventory trees. Detailed site information is available at the website MASTIF.
Covariates
Covariates in the model include as main effects tree diameter, tree canopy class (shading), and the climate variables in Fig. 1 of the main text and described in Table 1. A quadratic diameter term in the MASTIF model allows for changes in diameter response with size52. Shade classes follow the USDA Forest Inventory and Analysis (FIA)/National Ecological Observation Network (NEON) scheme that ranges from a fully exposed canopy that does not interact with canopies of other trees to fully shaded in the understory. Shading provides information on competition that has proved highly significant for fecundity in previous analyses41,52.
Table 1 Predictors in the model, not all of which are important for all species.
Full size table

To distinguish between the effects of spatial variation versus interannual variability, spring T and moisture D are included in the model as site means and site anomalies35. Spring minimum T affect phenology and frost risk during flowering and early fruit initiation. Summer mean T (June–August) is included both as a linear and quadratic term. Mean summer T is linked to thermal energy availability during the growing season, with the quadratic term allowing for potential suppression due to extreme heat. Moisture D (cumulative monthly PET-P (potential evapotranspiration[-] minus precipitation) for January–August) is included as a site mean and an annual anomaly. Moisture D is important for carbon assimilation and fruit development during summer in the eastern continent and, additionally, from the preceding winter in the western continent. For species that develop over spring and summer, anomalies incorporate the current and previous year. We did not include longer lags in covariates. For species that disperse seed in spring (Ulmus spp. and some members of Acer), only the previous year was used. Temperature anomalies were included for spring, but not summer, simply to reduce the number of times that temperature variables enter the model, and these two variables tended to be correlated at many sites.
Climate covariates were derived from gridded climate products and combined with local climate monitoring where it is available. Terraclimate71 provides monthly resolution, but it is spatially coarse. For both norms and trends, we used the period from 1990 to 2019 because global temperatures have been increasing consistently since the 1980s, and this span broadly overlaps with fecundity data (Fig. 6). CHELSA72 data are downscaled to a 1 km grid, but it does not extend to 2019. Our three-component climate scaling used regression to project CHELSA forward using Terraclimate, followed by downscaling to 1 km with CHELSA, with further downscaling to local climate data. Even where local climate data exist, they often do not span the full duration of field studies, making the link to gridded climate data important. Local climate data were especially important for mountainous sites in the Appalachians, Rockies, Sierra Nevada, and Cascades.
Of the full list of variables, a subset was retained, depending on species (some have narrow geographic ranges) and deviance information criteria of the fitted model (Supplementary Data 2). Year effects in the model allow for the interannual variation that is not absorbed by anomalies20.
Model fitting and TA
As mentioned above, model fitting applied the hierarchical Bayes model of ref. 20 to the combination of time series and opportunistic observations summarized in Fig. 1. Posterior simulation was completed with Markov chain Monte Carlo based on direct sampling, Metropolis, and Hamiltonian Markov chain. Model fitting used 211,146 trees and 2,566,594 tree-years from 123 species (Supplementary Data 2). Only species with multiple observation years were included.
The climate variable referenced as C in Eq. (1) of the main text is, in fact, a vector of climate variables described in the previous section, spring minimum T, summer mean T, and moisture D (Table 1). The anomalies and year effects in the fitted model contribute to the trends not explained by biogeographic variation as γ in Eq. (1). For main effects in the model, the partial derivatives are fitted coefficients, an example being the response to spring minimum temperature (partial f/partial {T}_{mathrm{sp}}={beta }_{{T}_{mathrm{sp}}}). For predictors involved in interactions, the partial derivatives are combinations of fitted coefficients and variables. For example, the response to moisture D, which interacts with tree size, is (partial [F], f/partial {D}={beta }_{{D}} + beta_{GD}G). The response to diameter G, which is quadratic and interacts with D, is (partial f/partial G={beta }_{G}+2{beta }_{{G}^{2}}G ,+{beta }_{GD}D).
Trend decomposition applied the fitted model to every tree in forest inventories from the USDA FIA program, the Canada’s National Forest Inventory, the NEON, and our MASTIF collaboration. Each tree in these inventories has a species and diameter. For trees that lack a canopy class, regression was used to predict it from distances and tree diameters based on inventories that include both location and canopy class, including NEON, FIA, and the MASTIF network. Although inventories differ in the minimum diameter they record, few trees are reproductive at diameters below the lower diameter limits in these surveys, so the effect on fecundity estimates is negligible. For the indirect effects of climate coming through tree growth rates, the same covariates were fitted to growth as previously defined for fecundity, using the change in diameter observed over multiple inventories. A Tobit model was used to accommodate the fact that a second measurement can be smaller than an earlier measurement. The Tobit thus treats negative growth as censored at zero. TA to inventory plots used 7,717,677 trees. Because not all species in the inventory data are included in the MASTIF network, mean fecundity parameters for the genus were used for unfitted species. Species fitted in the MASTIF network accounted for >90% of trees in inventory plots (Supplementary Data 2).
From the predictive distributions for every tree in the inventory data, we evaluated predictive mean trends aggregated to species and plot in Fig. 2b. We extracted specific terms that comprise the components in Fig. 4 and aggregated them too to the plot averages.
General form for TA
Equation 1 simplifies the model to highlight direct and indirect effects. Again, climate variables and tree size represent only a subset of the predictors in the model that are collected in a design vector ({{bf{x}}}_{t}=[{x}_{1,t},ldots ,{x}_{Q,t}]^{prime}), where the q = 1, …, Q predictors include shading from local competition, individual size, and climate and habitat variables (Table 1). On the proportionate scale, Eq. (1) can be written in terms of all predictors, including main effects and interactions, as

$$frac{{mathrm{d}}f}{{mathrm{d}}t}=mathop{sum }limits_{q=1}^{Q}left(frac{partial f}{partial {x}_{q}}+sum _{q^{prime} in {I}_{q}}frac{partial f}{partial ({x}_{q}{x}_{q^{prime} })}{x}_{q^{prime} }right)frac{{mathrm{d}}{x}_{q}}{{mathrm{d}}t}+gamma$$
(2)

where Iq are variables that interact with xq. In this application, interactions include tree diameter with moisture deficit and diameter squared. Each term in the summation consists of a main effect of xq and interactions that are multiplied by the rate of change in variable xq. For the simple case of only two predictors, Eq. (2) is recognizable as Eq. (1) of the main text, where x1, x2 have been substituted for variables G and C. In our application, predictors include additional climate and shading (Table 1).
Recognizing that environmental variables affect not only fecundity but also growth rate, we extract the size effect, that is, the xq that is G, and incorporate these indirect effects (through growth) by expanding g = dG/dt in Eq. (1) of the main text as

$$g=mathop{sum }limits_{q=1}^{Q}left(frac{partial g}{partial {x}_{q}}+mathop{sum}limits _{q^{prime} in {I}_{q}}frac{partial g}{partial ({x}_{q}{x}_{q^{prime} })}{x}_{q^{prime} }right){x}_{q}+nu$$
(3)

where ν is the component of growth that is not accommodated by other terms. This expression allows us to evaluate the full effect of climate variables, including those coming indirectly through growth.
Connecting fitted coefficients in MASTIF to TA
This section connects the continuous, deterministic Eq. (1) to the MASTIF model20 with the interpretation of responses, direct effects, and full effects of Fig. 5. To summarize key elements of the fitted model20, consider a tree i at site j that grows to reproductive maturity and then produces seed depending on its size, local competitive environment, and climate. We wish to estimate the effects of its changing environment and condition on fecundity using a model that includes spatial variation in predictors that are tracked longitudinally over years t. Fecundity changes through maturation probability ρij(t), which increases as trees increase in size, and through conditional fecundity ψij(t), the annual seed production of a mature tree. Let zij(t) = 1 be the event that a randomly selected tree i is mature in year t. Then, ρij(t) is the corresponding probability that the tree is mature, E[zij(t)] = ρij(t)(ρ is not to be confused with the probability that a tree that is now immature will make the transition to the mature state in an interval dt = 1. That is a different quantity detailed in the Supplement to ref. 41). Fecundity has expected value Fij(t) = ρij(t)ψij(t). On a proportionate (log) scale,

$${f}_{ij}(t)={mathrm{log}},{F}_{ij}(t)={mathrm{log}},{rho }_{ij}(t)+{mathrm{log}},{psi }_{it}(t)$$
(4)

The corresponding rate equation is

$$frac{{mathrm{d}}f}{{mathrm{d}}t}=frac{{mathrm{d}},{mathrm{log}},rho }{{mathrm{d}}t}+frac{{mathrm{d}},{mathrm{log}},psi }{{mathrm{d}}t}$$
(5)

The discretized and stochasticized version of Eq. (1) is

$$frac{{mathrm{d}}{F}_{ij}}{{mathrm{d}}t} = , frac{{F}_{ij,t+{mathrm{d}}t}-{F}_{ij,t}}{{mathrm{d}}t}+{epsilon }_{ij,t}\ = , {{Delta }}{F}_{ij,t}+{epsilon }_{ij,t}$$
(6)

where dt = 1 and ϵij,t is the integration error. When applied to a dynamic process model, this term further absorbs process error (see above), which is critical here to allow for conditional independence where observations are serially dependent. In simplest terms, ϵ is model miss-specification that allows for dependence in data.
The MASTIF model that provides estimates for TA is detailed in ref. 20. Elements of central interest for TA are

$${F}_{ij,t} = , {z}_{ij,t}{psi }_{ij,t}\ left[{z}_{ij,t}=1right] sim , {{Bernoulli}}left({rho }_{ij,t}right)\ {rho }_{ij,t} = , {{Phi }}({{boldsymbol{mu }}}_{ij,t})\ mathrm{log},{psi }_{ij,t} = ,{{bf{x}}}_{ij,t}^{prime}{boldsymbol{beta }}+{h}_{t}left(Tright)+{epsilon }_{ij,t}$$

where μij,t = α0 + αGGij,t describes the increase in maturation probability with size, Φ(⋅) is the standard normal distribution function (a probit), ϵij,t ~ N(0, σ2), and ht(T) can include year effects, h(T) = κt, or lagged effects, (h(T)=mathop{sum }nolimits_{k = 1}^{p}{kappa }_{k}{psi }_{ij,t-k}), that contribute to γ in Eq. (1) of the main text. If year effects are used, then γ includes the trend in year effects. (The generative version of this model writes individual states at t conditional on t − 1 and is given in the Supplement to ref. 20.). If an AR(p) model is used, then γ = κ1 (provided data are not detrended). Random individual effects in the fitted model are marginalized for prediction of trees that were not fitted, meaning that σ2 is the sum of model residual and random-effects variance. Again, the length-Q design vector xij,t includes individual attributes (e.g., diameter Gij,t), local competitive environment, and climate (Table 1). There is a corresponding coefficient vector β.
Moving to a difference equation (rate of change) for conditional log fecundity,

$${{Delta }}{f}_{ij,t}={{Delta }}mathrm{log},{rho }_{ij,t}+{{Delta }}mathrm{log},{psi }_{ij,t}$$

where

$${{Delta }}mathrm{log},{psi }_{ij,t} =mathrm{log},{psi }_{ij,t+1}-mathrm{log},{psi }_{ij,t}\ ={{Delta }}{{bf{x}}}_{ij,t}^{prime}{boldsymbol{beta }}+{gamma }_{ij,t}+{nu }_{ij,t}\ {{Delta }}{{bf{x}}}_{ij,t} ={{bf{x}}}_{i,t}-{{bf{x}}}_{ij,t-1}\ {nu }_{ij,t} sim N(0,2{sigma }^{2})$$

The variance in the last line is the variance of the difference Δϵij,t.
Elements of basic theory in Eq. (1) of the main text are linked to data through the modeling framework as

$${{Delta }}{f}_{ij,t}= +{beta }_{{T}_{sp}}{{Delta }}{T}_{sp,j}\ +left({beta }_{T}+2{beta }_{{T}^{2}}{T}_{j}right){{Delta }}{T}_{j}\ +left({beta }_{D}+{beta }_{GD}{G}_{ij,t}right){{Delta }}{D}_{j}\ +left({alpha }_{G}frac{phi ({{boldsymbol{mu }}}_{ij,t})}{{{Phi }}({{boldsymbol{mu }}}_{ij,t})}+{beta }_{G}+2{beta }_{{G}^{2}}{G}_{ij,t}+{beta }_{GD}{D}_{j}right){{Delta }}{G}_{ij}\ +{gamma }_{ij,t}+{nu }_{ij,t}$$
(7)

where ϕ(⋅) is the standard normal density function that comes from the rate of progress toward maturation. Again, the anomalies do not appear in this expression for trends because trends in the anomalies and year effects enter through γ.
The first four lines in Eq. (7) are, respectively, the effects of trends in spring minimum temperatures ΔTsp,j, summer mean temperature ΔTj, moisture deficit ΔDj, and size ΔGij, where the latter comes from growth on inventory plots. The contribution of maturation to change in fecundity is the first term in the fourth line, αGϕ(μij,t)/Φ(μij,t). A map of this term in Fig. 7b shows the strong contribution to fecundity in the East due to the young (Fig. 7a) and/or small (Fig. 4b) trees there. The sum of these terms dominates the patterns in Fig. 3c.
Fig. 7: Size and maturation effects on fecundity.

a Stand age variable in FIA data and b positive effect of maturation for increasing fecundity in the eastern continent. In the West, maturation does not contribute to rising fecundity because large trees are predominantly [mature] larger.

Full size image

All terms in Eq. (7) have units of mean change in proportionate fecundity, and these are mapped in figures of the main text. We focus on proportionate fecundity because it reflects the full effect of climate as opposed to total fecundity, which would often be dominated by one or a few trees of a single species. However, from proportionate fecundity we can obtain change in fecundity as ΔFij,t = Fij,t × Δfij. Stand-level effects on fecundity change at site j can be obtained from individual change as

$${{Delta }}{F}_{j}=mathop{sum }limits_{i=1}^{{n}_{j}}{{Delta }}{F}_{ij}=mathop{sum }limits_{i=1}^{{n}_{j}}{F}_{ij}{{Delta }}{f}_{ij,t}$$

Again, maps in Fig. 5 show mean proportionate effects for all trees on an inventory plot.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More