Ecology

1.Nobel, P. S. Physicochemical and Environmental Plant Physiology (Academic Press, 2005).
Google Scholar 
2.Hsiao, T. C. Plant responses to water stress. Annu. Rev. Plant Physiol. 24, 519–570 (1973).CAS 
Article 

Google Scholar 
3.Schönherr, J. Resistance of plant surfaces to water loss : transport properties of cutin, suberin and associated lipids. In Encyclopedia Plant Physiology, NS Vol. 12B (eds Lange, O. L. et al.) 154–179 (Springer, 1982).
Google Scholar 
4.Lendzian, K. J. Gas permeability of plant cuticles: oxygen permeability. Planta 155, 310–315 (1982).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
5.Langenfeld-Heyser, R. Physiological functions of lenticels. In Trees—Contributions to Modern Tree Physiology (eds Rennenberg, H. et al.) 43–56 (Backhuys, 1997).
Google Scholar 
6.Riederer, M. & Schreiber, L. Protecting against water loss: analysis of the barrier properties of plant cuticles. J. Exp. Bot. 52, 2023–2032 (2001).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
7.Kerstiens, G. Parameterization, comparison, and validation of models quantifying relative change of cuticular permeability with physicochemical properties of diffusants. J. Exp. Bot. 57, 2525–2533 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Schönherr, J. Characterization of aqueous pores in plant cuticles and permeation of ionic solutes. J. Exp. Bot. 57, 2471–2491 (2006).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
9.Groh, B., Hübner, C. & Lendzian, K. J. Water and oxygen permeance of phellems isolated from trees: the role of waxes and lenticels. Planta 215, 794–801 (2002).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
10.Lendzian, K. J. Survival strategies of plants during secondary growth: barrier properties of phellems and lenticels towards water, oxygen, and carbon dioxide. J. Exp. Bot. 57, 2535–2546 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Renault, H. et al. (2017) A phenol-enriched cuticle is ancestral to lignin evolution in land plants. Nat. Commun. 8, 14713 (2017).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
12.Haas, K. Phytochemische und rasterelektronenmikroskopische Untersuchungen zum Oberflächenwachs von Laubmoosen (Bryatae) (Grauer, 1999).
Google Scholar 
13.Clémençon, H., Emmett, V. & Emmett, E. E. Cytology and Plectology of the Hymenomycetes (J Cramer, 2012).
Google Scholar 
14.Moore, D., Gange, A. C., Gange, E. G. & Boddy, L. Fruit bodies: their production and development in relation to environment. In Ecology of Saprotrophic Basidiomycetes (eds Boddy, L. et al.) 79–103 (Elsevier Academic Press, 2008).
Google Scholar 
15.Halbwachs, H., Simmel, J. & Bässler, C. Tales and mysteries of fungal fruiting: how morphological and physiological traits affect a pileate lifestyle. Fungal Biol. Rev. 30, 36–61 (2016).Article 

Google Scholar 
16.Sakamoto, Y. Influences of environmental factors on fruiting body induction, development and maturation in mushroom-forming fungi. Fungal Biol. Rev. 32, 236–248 (2018).Article 

Google Scholar 
17.Straatsma, G., Ayer, F. & Egli, S. Species richness, abundance, and phenology of fungal fruit bodies over 21 years in a Swiss forest plot. Mycol. Res. 105(5), 515–523 (2001).Article 

Google Scholar 
18.Kües, U. & Liu, Y. Fruiting body production in basidiomycetes. Appl. Microbiol. Biotechnol. 54, 141–152 (2000).PubMed 
Article 
PubMed Central 

Google Scholar 
19.Beluhan, S. & Ranogajec, A. Chemical composition and non-volatile components of Croatian wild edible mushrooms. Food Chem. 124, 1076–1082 (2011).CAS 
Article 

Google Scholar 
20.Beecher, T. M., Magan, N. & Burton, K. S. Water potentials and soluble carbohydrate concentrations in tissues of freshly harvested and stored mushrooms (Agaricusbisporus). Postharvest Biol. Technol. 22, 121–131 (2001).CAS 
Article 

Google Scholar 
21.Bonnier G, Mangin L (1884) Recherches sur la respiration et la transpiration des champignons. Ann. Sc. Natur., sér. VI, t. XVII:210–30522.Moser, M. Transpirationsschutz bei höheren Pilzen. Schweizerische Zeitschrift für Pilzkunde 42(4), 50–54 (1964).
Google Scholar 
23.Pieschel, E. Über die Transpiration und die Wasserversorgung der Hymenomyceten. Bot. Archiv. VIII, 64–104 (1924).
Google Scholar 
24.Seybold, A. Weitere Beiträge zur Transpirationsanalyse. IV. Über die Transpiration der Hutpilze. Planta 16, 518–525 (1932).Article 

Google Scholar 
25.Becker, M., Kerstiens, G. & Schönherr, J. Water permeability of plant cuticles: permeance, diffusion and partition coefficients. Trees 1, 54–60 (1986).CAS 
Article 

Google Scholar 
26.Schreiber, L. & Schönherr, J. Water and Solute Permeability of Plant Cuticles. Measurement and Data Analysis (Springer, 2009).
Google Scholar 
27.Schönherr, J. & Lendzian, K. J. A simple and inexpensive method of measuring water permeability of isolated plant cuticular membranes. Z Pflanzenphysiol 102, 321–327 (1981).Article 

Google Scholar 
28.Weast, R. C. CRC Handbook of Chemistry and Physics: Humidity Constant (CRC Press, 1983).
Google Scholar 
29.Riederer, M. & Schneider, G. Comparative study of the composition of waxes extracted from isolated leaf cuticules and from whole leaves of Citrus: evidence for selective extraction. Physiol. Plant 77, 373–384 (1989).CAS 
Article 

Google Scholar 
30.Lendzian, K. J. & Kerstiens, G. Sorption and transport of gases and vapors in plant cuticles. Rev. Environ. Cont. Tox. 121, 65–128 (1991).CAS 

Google Scholar 
31.Kerstiens, G., Federholzner, R. & Lendzian, K. J. Dry deposition and cuticular uptake of pollutant gases. Agric. Ecosyst. Environ. 42, 239–253 (1992).CAS 
Article 

Google Scholar 
32.Metzler, H. & Krause, B. Angewandte Statistik (Dt Verlag Wiss, 1983).
Google Scholar 
33.Baur, P. Lognormal distribution of water permeability and organic solute mobility in plant cuticles. Plant Cell Environ. 20, 167–177 (1997).ADS 
CAS 
Article 

Google Scholar 
34.Stamets, P. Growing Gourmet and Medicinal Mushrooms (Ten Speed Press, 1993).
Google Scholar 
35.Moser, M. Fungal growth and fructification under stress conditions. Ukrainian Botanical J. 50, 5–12 (1993).
Google Scholar 
36.Pinna, S., Gevry, M. F., Côté, M. & Sirois, L. Factors influencing fructification phenology of edible mushrooms in a boreal mixed forest of Eastern Canada. For. Ecol. Manag. 260(3), 294–301 (2010).Article 

Google Scholar 
37.Buller, A. H. R. Researches on Fungi. II. Further Investigations Upon the Production and Liberation of Spores in Hymenomyctes (Hafner Publishing Co, 1922).
Google Scholar 
38.Kües, U. Life history and developmental processes in the basidiomycete Coprinus cinereus. Microbiol. Mol. Biol. Rev. 64, 316–353 (2000).PubMed 
PubMed Central 
Article 

Google Scholar 
39.Money, N. More g’s than the space shuttle: ballistospore discharge. Mycologia 90, 547–558 (1998).Article 

Google Scholar 
40.Husher, J. et al. Evaporative cooling of mushrooms. Mycologia 91, 351–352 (1999).Article 

Google Scholar 
41.Dressaire, E., Yamada, L., Song, B. & Roper, M. Mushrooms use convectively created airflows to disperse their spores. Proc. Natl. Acad. Sci. U. S. A. 113, 2833–2838 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
42.De Groot, P. W., Schaap, P. J., Sonnenberg, A. S., Visser, J. & Van Griensven, L. J. The Agaricus bisporus hypAgene encodes a hydrophobin and specifically accumulates in peel tissue of mushroom caps during fruit body development. J. Mol. Biol. 257, 1008–1018 (1996).PubMed 
Article 
PubMed Central 

Google Scholar 
43.Wösten, H. A. B. & Wessels, J. G. H. The emergence of fruiting bodies in basidiomycetes. In Growth, Differentiation and Sexuality. The Mycota (A Comprehensive Treatise on Fungi as Experimental Systems for Basic and Applied Research), Vol. 1 (eds. Kües, U. & Fischer, R.) (Springer, 2006).44.Itoh, Y. H., Sugai, A., Uda, I. & Itoh, T. The evolution of lipids. Adv. Space Res. 28, 719–724 (2001).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
45.Segré, D., Ben-Eli, D., Deamer, D. W. & Lancet, D. The lipid world. Origins Life Evol. Biosphere 31, 119–145 (2001).ADS 
Article 

Google Scholar 
46.Samson, R. A., Stalpers, J. A. & Verkerke, W. A simplified technique to prepare fungal specimens for scanning electronmicroscopy. Cytobios 24, 7–11 (1979).CAS 
PubMed 
PubMed Central 

Google Scholar  More

1.Whyte, I. Studying elephant movements, in studying elephants, in African Wildl. Found. Tech. Ser. 7. African Wildl. Found. (ed Kangwana, K.) 75–89 (1996).2.Rasmussen, L. E. L. & Krishnamurthy, V. How chemical signals integrate Asian elephant society: The known and the unknown. Zoo Biol. 19, 405–423 (2000).CAS 
Article 

Google Scholar 
3.Nair, S., Balakrishnan, R., Seelamantula, C. S. & Sukumar, R. Vocalizations of wild Asian elephants (Elephas maximus ): structural classification and social context. J. Acoust. Soc. Am. 126, 2768–2778 (2009).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
4.Stoeger, A. S. & Manger, P. Vocal learning in elephants: neural bases and adaptive context. Curr. Opin. Neurobiol. 28, 101–107 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Moss, C. J. & Poole, J. H. Relationships and social structure in African elephants. Primate Soc. Relationsh.: An Integr. Approach 315-325 (1983).
Google Scholar 
6.Foerder, P., Galloway, M., Barthel, T., Moore, D. E. & Reiss, D. Insightful problem solving in an asian elephant. PLoS ONE 6, e23251 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
7.Lee, P. C. Allomothering among African elephant. Animal Behaviour 35, 278-291 (1987).8.Byrne, R. W., Bates, L. & Moss, C. J. Comparative cognition & behavior reviews. Elephant Cogn. 4, 65–79 (2009).
Google Scholar 
9.Bates, L. A., Poole, J. H. & Byrne, R. W. Elephant cognition. Curr. Biol. 18, 544–546 (2008).Article 
CAS 

Google Scholar 
10.Vance, E. A., Archie, E. A. & Moss, C. J. Social networks in African elephants. Comput. Math. Organ. Theory 15, 273–293 (2009).Article 

Google Scholar 
11.de Silva, S. & Wittemyer, G. A comparison of social organization in Asian elephants and African savannah elephants. Int. J. Primatol. 33, 1125–1141 (2012).Article 

Google Scholar 
12.Shoshani, J. Understanding proboscidean evolution: a formidable task. Trends Ecol. Evol. 13, 480–487 (1998).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
13.Rohland, N. et al. Proboscidean mitogenomics: chronology and mode of elephant evolution using mastodon as outgroup. PLoS Biol. 5, 1663–1671 (2007).CAS 
Article 

Google Scholar 
14.de Flamingh, A. Genetic structure of the savannah elephant population (Loxodonta africana (Blumenbach 1797)) in the Kavango-Zambezi Transfrontier Conservation Area. ProQuest Diss. Theses 102 (2013).15.Grubb, P., Groves, C. P., Dudley, J. P. & Shoshani, J. Living African elephants belong to two species: Loxodonta africana (Blumenbach, 1797) and Loxodonta cyclotis (Matschie, 1900). Elephant 2, 1–4 (2000).Article 

Google Scholar 
16.Roca, A. L., Georgiadis, N., Pecon-Slattery, J. & O’Brien, S. J. Genetic evidence for two species of elephant in Africa. Science 293, 1473–1477 (2001).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
17.Roca, A. L. et al. Elephant natural history: a genomic perspective. Annu. Rev. Anim. Biosci. 3, 139–167 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
18.Wasser, S. K. et al. Assigning African elephant DNA to geographic region of origin: applications to the ivory trade. Proc. Natl. Acad. Sci. U. S. A. 101, 14847–14852 (2004).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
19.Ishida, Y. et al. Distinguishing forest and savanna African elephants using short nuclear DNA sequences. J. Hered. 102, 610–616 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
20.Comstock, K. E. et al. Patterns of molecular genetic variation among African elephant populations. Mol. Ecol. 11, 2489–2498 (2002).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
21.Palkopoulou, E. et al. A comprehensive genomic history of extinct and living elephants. Proc. Natl. Acad. Sci. U. S. A. 115, E2566–E2574 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
22.Shoshani, J. & Eisenberg, J. F. Elephas maximus. Mamm. Species 182, 1–8 (1982).Article 

Google Scholar 
23.Sukumar, R. The Living Elephants (Oxford University Press, 2003).
Google Scholar 
24.Sukumar, R. A brief review of the status, distribution and biology of wild Asian elephants Elephas maximus. Int. Zoo Yearb. 40, 1–8 (2006).Article 

Google Scholar 
25.Olivier, R. Distribution and status of the Asian elephant. Oryx 14(4), 379–424. https://doi.org/10.1017/S003060530001601X (1978).Article 

Google Scholar 
26.Santiapillai, C. The Asian elephant conservation: a global strategy. Gajah 18, 21–39 (1997).
Google Scholar 
27.Sukumar, R. Ecology of the Asian elephant in Southern India. i. movement and habitat utilization patterns. J. Trop. Ecol. 5, 1–18 (1989).Article 

Google Scholar 
28.Vidya, T. N. C., Fernando, P., Melnick, D. J. & Sukumar, R. Population genetic structure and conservation of Asian elephants (Elephas maximus) across India. Anim. Conserv. 8, 377–388 (2005).Article 

Google Scholar 
29.Fleischer, R. C., Perry, E. A., Muralidharan, K., Stevens, E. E. & Wemmer, C. M. Phylogeography of the Asian elephant (Elephas maximus) based on mitochondrial DNA. Evolution (N. Y.) 55, 1882–1892 (2001).CAS 

Google Scholar 
30.Fernando, P., Pfrender, M. E., Encalada, S. E. & Lande, R. Mitochondrial DNA variation, phylogeography and population structure of the Asian elephant. Heredity (Edinb). 84, 362–372 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
31.Fernando, P. Elephants in Sri Lanka: past, present, and future. Loris 22, 38–44 (2000).
Google Scholar 
32.Hendavitharana, W., Dissanayake, S. & de Silva, M. The survey of elephants in Sri Lanka. Gajah 12, 1–30 (1994).
Google Scholar 
33.Eggert, L. S., Rasner, C. A. & Woodruff, D. S. The evolution and phylogeography of the African elephant inferred from mitochondrial DNA sequence and nuclear microsatellite markers. Hungarian Q. 49, 1993–2006 (2008).
Google Scholar 
34.Ishida, Y., Georgiadis, N. J., Hondo, T. & Roca, A. L. Triangulating the provenance of African elephants using mitochondrial DNA. Evol. Appl. 6, 253–265 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
35.Liu, C. Z., Wang, L., Xia, X. J. & Jiang, J. Q. Characterization of the complete mitochondrial genome of cape elephant shrew, Elephantulus edwardii. Mitochondrial DNA Part B Resour. 3, 738–739 (2018).Article 

Google Scholar 
36.Fernando, P. et al. DNA analysis indicates that Asian elephants are native to Borneo and are therefore a high priority for conservation. PLoS Biol. 1, 110–115 (2003).CAS 
Article 

Google Scholar 
37.Ahlering, M. A. et al. Genetic diversity, social structure, and conservation value of the elephants of the Nakai Plateau, Lao PDR, based on non-invasive sampling. Conserv. Genet. 12, 413–422 (2011).Article 

Google Scholar 
38.Goossens, B. et al. Habitat fragmentation and genetic diversity in natural populations of the Bornean elephant: implications for conservation. BIOC 196, 80–92 (2016).
Google Scholar 
39.Shoshani, J., Golenberg, E. M. & Yang, H. Elephantidae phylogeny: Morphological versus molecular results. Acta Theriol. (Warsz) 43, 89–122 (1998).Article 

Google Scholar 
40.Vidya, T. N. C. & Sukumar, R. Amplification success and feasibility of using microsatellite loci amplified from dung to population genetic studies of the Asian elephant (Elephas maximus). Curr. Sci. 88, 489–492 (2005).CAS 

Google Scholar 
41.Vidya, T. N. C., Varma, S., Dang, N. X., Van Thanh, T. & Sukumar, R. Minimum population size, genetic diversity, and social structure of the Asian elephant in Cat Tien National Park and its adjoining areas, Vietnam, based on molecular genetic analyses. Conserv. Genet. 8, 1471–1478 (2007).Article 

Google Scholar 
42.Suwattana, D., Jirasupphachok, J., Kanchanapangka, S. & Koykul, W. Tetranucleotide microsatellite markers for molecular testing in Thai domestic elephants (Elephas maximus indicus). Thai J. Vet. Med. 40, 405–409 (2010).
Google Scholar 
43.Eggert, L. S. et al. Using genetic profiles of African forest elephants to infer population structure, movements, and habitat use in a conservation and development landscape in Gabon. Conserv. Biol. 28, 107–118 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
44.Kinuthia, J. et al. The selection of a standard STR panel for DNA profiling of the African elephant (Loxodonta africana) in Kenya. Conserv. Genet. Resour. 7, 305–307 (2015).Article 

Google Scholar 
45.Hedges, S. Monitoring elephant populations and assessing threats. Universities Press (India) Pvt. Ltd., Hyderabad, India 259–292 (2012).46.Eggert, L. S., Ramakrishnan, U., Mundy, N. I. & Woodruff, D. S. Polymorphic microsatellite DNA markers in the African elephant (Loxondonta africana) and their use in the Asian elephant (Elephas maximus). Mol. Ecol. 9, 2222–2224 (2000).Article 

Google Scholar 
47.Nyakaana, S., Arctander, P. & Siegismund, H. R. Population structure of the African savannah elephant inferred from mitochondrial control region sequences and nuclear microsatellite loci. Heredity (Edinb). 89, 90–98 (2002).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
48.Kongrit, C. et al. Isolation and characterization of dinucleotide microsatellite loci in the Asian elephant (Elephas maximus). Mol. Ecol. Resour. 8, 175–177 (2007).Article 
CAS 

Google Scholar 
49.Fernando, P., Vidya, T. N. C. & Melnick, D. J. Isolation and characterization of tri- and tetranucleotide microsatellite loci in the Asian elephant, Elephas maximus. Mol. Ecol. Resour. 8, 232–233 (2001).Article 

Google Scholar 
50.Archie, E. A., Moss, C. J. & Alberts, S. C. Characterization of tetranucleotide microsatellite loci in the African Savannah Elephant (Loxodonta africana africana). Mol. Ecol. Notes 3, 244–246 (2003).CAS 
Article 

Google Scholar 
51.Lieckfeldt, D., Schmidt, A. & Pitra, C. Isolation and characterization of microsatellite loci in the great bustard, Otis tarda. Mol. Ecol. Notes 1, 133–134 (2001).CAS 
Article 

Google Scholar 
52.Nyakaana, S., Okello, J. B. A., Muwanika, V. & Siegismund, H. R. Six new polymorphic microsatellite loci isolated and characterized from the African savannah elephant genome. Mol. Ecol. Notes 5, 223–225 (2005).CAS 
Article 

Google Scholar 
53.Okello, J. B. A. et al. Population genetic structure of savannah elephants in Kenya: conservation and management implications. J. Hered. 99, 443–452 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
54.Nyakaana, S. & Arctander, P. Isolation and characterization of microsatellite loci in the African elephant, Loxodonta africana. Mol. Ecol. 10, 1436–1437 (1998).
Google Scholar 
55.Comstock, K. E., Wasser, S. K. & Ostrander, E. A. Polymorphic microsatellite DNA loci identified in the African elephant (Loxodonta africana). Mol. Ecol. 9, 1004–1006 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
56.Hartl, G. B., Hartl, K. F., Hemmer, W. & Nadlinger, K. Electrophoretic and chromosomal variation in captive Asian elephants (Elephas maximus). Zoo Biol. 14, 87–95 (1995).Article 

Google Scholar 
57.Bourgeois, S. et al. Single-nucleotide polymorphism discovery and panel characterization in the African forest elephant. Ecol. Evol. 8, 2207–2217 (2018).PubMed 
PubMed Central 

Google Scholar 
58.Sharma, R. et al. Two different high throughput sequencing approaches identify thousands of De Novo genomic markers for the genetically depleted Bornean elephant. PLoS ONE 7, e49533 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
59.Reddy, P. C. et al. Comparative sequence analyses of genome and transcriptome reveal novel transcripts and variants in the Asian elephant Elephas maximus. J. Biosci. 40, 891–907 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
60.Mondol, S. et al. New evidence for hybrid zones of forest and savanna elephants in Central and West Africa. Mol. Ecol. 24, 6134–6147 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
61.Hou, Z. C. et al. Elephant transcriptome provides insights into the evolution of eutherian placentation. Genome Biol. Evol. 4, 713–725 (2012).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
62.Tollis, M. et al. Elephant Genomes Reveal Insights into Differences in Disease Defense Mechanisms between Species. bioRxiv 2020.05.29.124396 (2020).63.Rohland, N. et al. Genomic DNA sequences from mastodon and woolly mammoth reveal deep speciation of forest and savanna elephants. PLoS Biol. 8, 16–19 (2010).Article 
CAS 

Google Scholar 
64.Lynch, V. J. et al. Elephantid genomes reveal the molecular bases of woolly mammoth adaptations to the Arctic. Cell Rep. 12, 217–228 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
65.Yang, H., Golenberg, E. M. & Shoshani, J. Phylogenetic resolution within the elephantidae using fossil DNA sequence from the American mastodon (Mammut americanum) as an outgroup. Proc. Natl. Acad. Sci. U. S. A. 93, 1190–1194 (1996).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
66.Orlando, L., Hänni, C. & Douady, C. J. Mammoth and elephant phylogenetic relationships: Mammut americanum, the missing outgroup. Evol. Bioinforma. 3, 45–51 (2007).CAS 
Article 

Google Scholar 
67.Eggert, L. S., Eggert, J. A. & Woodruff, D. S. Estimating population sizes for elusive animals: the forest elephants of Kakum National Park, Ghana. Mol. Ecol. 12, 1389–1402 (2003).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
68.Vidya, T. N. C. Evolutionary history and population genetic structure of Asian elephants in India. Indian J. Hist. Sci. 51, 391–405 (2016).
Google Scholar 
69.Schuttler, S. G., Whittaker, A., Jeffery, K. J. & Eggert, L. S. African forest elephant social networks: fission-fusion dynamics, but fewer associations. Endanger. Species Res. 25, 165–173 (2014).Article 

Google Scholar 
70.Ahlering, M. A. et al. Identifying source populations and genetic structure for savannah elephants in human-dominated landscapes and protected areas in the Kenya-Tanzania borderlands. PLoS ONE 7, e52288 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
71.Vidya, T. N. C., Fernando, P., Melnick, D. J. & Sukumar, R. Population differentiation within and among Asian elephant (Elephas maximus) populations in southern India. Heredity (Edinb). 94, 71–80 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
72.Sukumar, R., Ramakrishnan, U. & Santosh, J. A. Impact of poaching on an Asian elephant population in Periyar, southern India: a model of demography and tusk harvest. Anim. Conserv. 1, 281–291 (1998).Article 

Google Scholar 
73.Mondol, S., Mailand, C. R. & Wasser, S. K. Male biased sex ratio of poached elephants is negatively related to poaching intensity over time. Conserv. Genet. 15, 1259–1263 (2014).Article 

Google Scholar 
74.Breuer, T., Maisels, F. & Fishlock, V. The consequences of poaching and anthropogenic change for forest elephants. Conserv. Biol. 30, 1019–1026 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
75.Mailand, C. & Wasser, S. K. Isolation of DNA from small amounts of elephant ivory. Nat. Protoc. 2, 2228–2232 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
76.Lee, E. et al. The identification of elephant ivory evidences of illegal trade with mitochondrial cytochrome b gene and hypervariable D-loop region. J. Forensic Leg. Med. 20, 174–178 (2015).Article 

Google Scholar 
77.Chakraborty, S., Boominathan, D., Desai, A. A. & Vidya, T. N. C. Using genetic analysis to estimate population size, sex ratio, and social organization in an Asian elephant population in conflict with humans in Alur, southern India. Conserv. Genet. 15, 897–907 (2014).Article 

Google Scholar 
78.Fernando, P. & Pastorini, J. Range-wide status of Asian elephants. Gajah 35, 15–20 (2011).
Google Scholar 
79.Ishida, Y., Gugala, N. A., Georgiadis, N. J. & Roca, A. L. Evolutionary and demographic processes shaping geographic patterns of genetic diversity in a keystone species, the African forest elephant (Loxodonta cyclotis). Ecol. Evol. 8, 4919–4931 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
80.Kongrit, C. Genetic tools for the conservation of wild Asian elephants. Int. J. Biol. 9, 1 (2017).Article 

Google Scholar 
81.McComb, K., Shannon, G., Sayialel, K. N. & Moss, C. Elephants can determine ethnicity, gender, and age from acoustic cues in human voices. Proc. Natl. Acad. Sci. U. S. A. 111, 5433–5438 (2014).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
82.Prithiviraj, F. & Melnick, D. J. Molecular sexing eutherina mammals. Mol. Ecol. Notes 1, 350–353 (2001).Article 

Google Scholar 
83.Vandebona, H. et al. DNA fingerprints of the Asian elephant in Sri Lanka, Elephas maximus maximus, using multilocus probe 33.15 (Jeffreys). J. Natl. Sci. Found. Sri Lanka 32, 83–96 (2004).CAS 
Article 

Google Scholar 
84.Gugala, N. A., Ishida, Y., Georgiadis, N. J. & Roca, A. L. Development and characterization of microsatellite markers in the African forest elephant (Loxodonta cyclotis). BMC Res. Notes 9, 4–9 (2016).Article 
CAS 

Google Scholar 
85.Zhang, L. et al. Asian elephants in China: estimating population size and evaluating habitat suitability. PLoS ONE 10, 1–13 (2015).
Google Scholar 
86.Vartia, S. et al. A novel method of microsatellite genotyping-by-sequencing using individual combinatorial barcoding. R. Soc. Open Sci. 3, 150565 (2016).ADS 
MathSciNet 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
87.Tighe, A. J. et al. Testing PCR amplification from elephant dung using silica-dried swabs. Pachyderm 59, 56–65 (2018).
Google Scholar 
88.Bourgeois, S. et al. Improving cost-efficiency of faecal genotyping: new tools for elephant species. PLoS ONE 14, e0210811 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
89.Hedges, S., Johnson, A., Ahlering, M., Tyson, M. & Eggert, L. S. Accuracy, precision, and cost-effectiveness of conventional dung density and fecal DNA based survey methods to estimate Asian elephant (Elephas maximus) population size and structure. Biol. Conserv. 159, 101–108 (2013).Article 

Google Scholar 
90.Moßbrucker, A. M. et al. Non-invasive genotyping of Sumatran elephants : implications for conservation The Sumatran elephant (Elephas maximus sumatranus) is one of three currently recognized subspecies. Trop. Conserv. Sci. 8, 745–759 (2015).Article 

Google Scholar 
91.Ishida, Y. et al. Short amplicon microsatellite markers for low quality elephant DNA. Conserv. Genet. Resour. 4, 491–494 (2012).Article 

Google Scholar 
92.Thitaram, C. et al. Evaluation and selection of microsatellite markers for an identification and parentage test of Asian elephants (Elephas maximus). Conserv. Genet. 9, 921–925 (2008).CAS 
Article 

Google Scholar 
93.Lorenz, T. C. Polymerase chain reaction: basic protocol plus troubleshooting and optimization strategies. J. Vis. Exp. 2, 1–15. https://doi.org/10.3791/3998 (2012).CAS 
Article 

Google Scholar 
94.Litt, M. & Luty, J. A. Hypervariable amplification. Am. J. Hum. Genet. 44, 397–401 (1989).CAS 
PubMed 
PubMed Central 

Google Scholar 
95.Park, Y. J., Lee, J. K. & Kim, N. S. Simple sequence repeat polymorphisms (SSRPs) for evaluation of molecular diversity and germplasm classification of minor crops. Molecules 14, 4546–4569 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
96.Vieira, M. L. C., Santini, L., Diniz, A. L. & Munhoz, C. D. F. Microsatellite markers: what they mean and why they are so useful. Genet. Mol. Biol. 39, 312–328 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
97.Stafne, E. T., Clark, J. R., Weber, C. A., Graham, J. & Lewers, K. S. Simple sequence repeat (SSR) markers for genetic mapping of raspberry and blackberry. J. Am. Soc. Hortic. Sci. 130, 722–728 (2005).CAS 
Article 

Google Scholar 
98.Tommasini, L. et al. The development of multiplex simple sequence repeat (SSR) markers to complement distinctness, uniformity and stability testing of rape (Brassica napus L.) varieties. Theor. Appl. Genet. 106, 1091–1101 (2003).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
99.Norrgard, K. Forensics, DNA Fingerprinting, and CODIS. Nat. Educ. 1, 35 (2008).
Google Scholar 
100.Maroju, P. A. et al. Schrodinger’s scat: A critical review of the currently available tiger (Panthera Tigris) and leopard (Panthera pardus) specific primers in India, and a novel leopard specific primer. BMC Genet. 17, 1–6 (2016).Article 

Google Scholar 
101.Waits, L. P. & Pearkau, D. Noninvasive genetic sampling tools for wildlife biologists: a review of applications and recommendations for accurate data collection. J. Wildl. Manag. 69, 1419–1433 (2005).Article 

Google Scholar 
102.Baird, N. A. et al. Rapid SNP discovery and genetic mapping using sequenced RAD markers. PLoS ONE 3, 1–7 (2008).Article 
CAS 

Google Scholar 
103.Miller, M. R., Dunham, J. P., Amores, A., Cresko, W. A. & Johnson, E. A. Rapid and cost-effective polymorphism identification and genotyping using restriction site associated DNA (RAD) markers. Genome Res. 17, 240–248 (2007).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
104.Delord, C. et al. A cost-and-time effective procedure to develop SNP markers for multiple species: a support for community genetics. Methods Ecol. Evol. 9, 1959–1974 (2018).Article 

Google Scholar 
105.Magwanga, R. O. et al. GBS mapping and analysis of genes conserved between Gossypium tomentosum and Gossypium hirsutum cotton cultivars that respond to drought stress at the seedling stage of the BC2F2generation. Int. J. Mol. Sci. 19, 1614 (2018).PubMed Central 
Article 
CAS 

Google Scholar 
106.Elshire, R. J. et al. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS ONE 6, 1–10 (2011).Article 
CAS 

Google Scholar 
107.Chandrasekara, C. H. W. M. R. B., Wijesundera, W. S. S., Perera, H. N., Chong, S. S. & Rajan-Babu, I. S. Cascade screening for fragile X syndrome/CGG repeat expansions in children attending special education in Sri Lanka. PLoS ONE 10, 1–10 (2015).
Google Scholar 
108.Felsenstein, J. 2002. {PHYLIP}(Phylogen. I. P. ver. 3. 6a3.—P. by the author. PHYLIP(Phylogeny Inference Package) ver. 3.6a3. (2002).109.Nei, M. & Li, W. H. Mathematical model for studying genetic variation in terms of restriction endonucleases. Proc. Natl. Acad. Sci. U. S. A. 76, 5269–5273 (1979).ADS 
CAS 
PubMed 
PubMed Central 
MATH 
Article 

Google Scholar 
110.Rambaut, A. FigTree ver.1. 3.1: tree figure drawing tool. http://tree.bio.ed.ac.uk/software/figtree. (2009). More

1.Stearns, S. C. The Evolution of Life Histories (Oxford University Press, 1992).
Google Scholar 
2.Caughley, G. & Sinclair, A. R. E. Wildlife Ecology and Management (Blackwell Science, 1994).
Google Scholar 
3.Servanty, S. et al. Influence of harvesting pressure on demographic tactics: Implications for wildlife management. J. Appl. Ecol. 48(4), 835–843 (2011).Article 

Google Scholar 
4.Marboutin, E., Bray, Y., Péroux, R., Mauvy, B. & Lartiges, A. Population dynamics in European hare: Breeding parameters and sustainable harvest rates. J. Appl. Ecol. 40(3), 580–591 (2003).Article 

Google Scholar 
5.Stoneberg, R. P. & Jonkel, C. L. Age determination of black bears by cementum layers. J. Wildlife Manage. 30(2), 411–414 (1966).Article 

Google Scholar 
6.Roth, V. L. & Shoshani, J. Dental identification and age determination in Elephas maximus. J. Zool. 214, 567–588 (1988).Article 

Google Scholar 
7.Dutta, S. & Sengupta, P. Men and mice: Relating their ages. Life Sci. 152, 244–248 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Dimmick, R. W. & Pelton, M. R. Criteria of sex and age. In Research and Management Techniques for Wildlife and Habitats 5th edn, (ed. Bookhout, T.
A.) 169–214 (The Wildlife Society, Bethesda, MA, US, 1994).
Google Scholar 
9.Morris, P. A review of mammalian age determination methods. Mamm. Rev. 2, 69–103 (1972).Article 

Google Scholar 
10.Augusteyn, R. C. On the relationship between rabbit age and lens dry weight: Improved determination of the age of rabbits in the wild. Mol. Vis. 13, 2030–2034 (2007).CAS 
PubMed 
PubMed Central 

Google Scholar 
11.Augusteyn, R. C. Growth of the lens: In vitro observations. Clin. Exp. Optom. 91(3), 226–239 (2008).PubMed 
Article 
PubMed Central 

Google Scholar 
12.Augusteyn, R. C. Growth of the eye lens: I. Weight accumulation in multiple species. Mol. Vis. 20, 410–426 (2014).PubMed 
PubMed Central 

Google Scholar 
13.Lord, D. R. The lens as an indicator of age in cottontail rabbits. J. Wildl. Manage. 23, 358–360 (1959).Article 

Google Scholar 
14.Forsyth, D. M., Garel, M. & McLeod, S. R. Estimating age and age class of harvested hog deer from eye lens mass using frequentist and Bayesian methods. Wildlife biol. 22(4), 137–143 (2016).Article 

Google Scholar 
15.Dudzinski, M. L. & Mykytowycz, R. The eye lens as an indicator of age in the wild rabbit in Australia. CSIRO Wildl. Res. 6, 156–159 (1961).Article 

Google Scholar 
16.Myers, K. & Gilbert, N. Determination of age of wild rabbits in Australia. J. Wildl. Manage. 32, 841–849 (1968).Article 

Google Scholar 
17.Wheeler, S. H. & King, D. R. The use of eye-lens weights for aging wild rabbits, Oryctolagus cuniculus (L.) in Australia. Aust. Wildl. Res. 7, 79–84 (1980).Article 

Google Scholar 
18.Tablado, Z., Revilla, E. & Palomares, F. Breeding like rabbits: Global patterns of variability and determinants of European wild rabbit reproduction. Ecography 32, 310–320. https://doi.org/10.1111/j.1600-0587.2008.05532.x (2009).Article 

Google Scholar 
19.Ferreira, C. et al. Biometrical analysis reveals major differences between the two subspecies of the European rabbit. Biol. J. Linn. Soc. 116, 106–116 (2015).Article 

Google Scholar 
20.Branco, M., Monnerot, M., Ferrand, N. & Templeton, A. R. Postglacial dispersal of the European rabbit (Oryctolagus cuniculus) on the Iberian Peninsula reconstructed from nested clade and mismatch analyses of mitochondrial DNA genetic variation. Evolution 56, 792–803. https://doi.org/10.1111/j.0014-3820.2002.tb01390.x (2002).Article 
PubMed 
PubMed Central 

Google Scholar 
21.Gómez, A. & Lunt, D. H. Refugia within refugia: Patterns of phylogeographic concordance in the Iberian Peninsula. In Phylogeography in Southern European Refugia (eds Weiss, S. & Ferrand, N.) 155–188 (Springer, 2006).
Google Scholar 
22.Geraldes, A. et al. Reduced introgression of the Y chromosome between subspecies of the European rabbit (Oryctolagus cuniculus) in the Iberian Peninsula. Mol. Ecol. 17, 4489–4499 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
23.Carneiro, M., Ferrand, N. & Nachman, M. W. Recombination and speciation: Loci near centromeres are more differentiated than loci near telomeres between subspecies of the European rabbit (Oryctolagus cuniculus). Genetics 181, 593–606 (2009).PubMed 
PubMed Central 
Article 

Google Scholar 
24.Rafati, N. et al. A genomic map of clinal variation across the European rabbit hybrid zone. Mol. Ecol. 27, 1457–1478. https://doi.org/10.1111/mec.14494 (2018).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
25.Vaquerizas, P. H. et al. The paradox of endangered European rabbits regarded as pests in the Iberian Peninsula: Subspecies differences in trends matter. Endang. Species Res. 43, 99–102 (2020).Article 

Google Scholar 
26.Arques, J. & Peiró, V. Estructura de Sexos y Edades de una población de Conejos (Oryctolagus cuniculus) del sudeste de España. Mediterránea. Serie de Estudios Biológicos 18, 1–33 (2005).
Google Scholar 
27.Trout, R. C. & Smith, G. C. The reproductive productivity of the wild rabbit (Oryctolagus cuniculus) in southern England on sites with different soils. J. Zool. 237(3), 411–422 (1995).Article 

Google Scholar 
28.Boussès, P., Arthur, C. & Chapuis, J. L. Rôle du facteur trophique sur la biologie des populations de lapins (Oryctolagus cuniculus L.) des Iles Kerguelen. Revue d’écologie 43, 329–343 (1988).
Google Scholar 
29.Bonino, N. & Donadio, E. Body parameters and sexual dimorphism in the European wild rabbit (Oryctolagus cuniculus) introduced in Argentina. Mastozool. Neotrop. 17(1), 123–127 (2010).
Google Scholar 
30.Carneiro, M. et al. Rabbit genome analysis reveals a polygenic basic for phenotypic change during domestication. Science 345(6200), 1074–1079 (2014).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
31.Myers, K. The rabbit in Australia. In Dynamics of Numbers in Populations (eds den Boer, P. J. & Gradwell, G. R.) 478–506 (Proceedings of the NATO Advanced Study Institute Oosterbeek, 1970).
Google Scholar 
32.Delibes-Mateos, M., Villafuerte, R., Cooke, B. & Alves, P. C. Oryctolagus cuniculus (Linnaeus, 1758). In Lagomorphs: Pikas, Rabbits and Hares of the World (eds Smith, A. T. et al.) 99–104 (John Hopkins University Press, 2018).
Google Scholar 
33.Carneiro, M. et al. The genomic architecture of population divergence between subspecies of the European rabbit. PLoS. Genet. 10(8), e1003519. https://doi.org/10.1371/journal.pgen.1003519 (2014).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
34.Bonino, N. & Soriguer, R. Genetic lineages of feral populations of the Oryctolagus cuniculus (Leporidae, Lagomorpha) in Argentina. Mammalia 72, 355–357 (2008).Article 

Google Scholar 
35.Branco, M. & Ferrand, N. Biochemical and population genetics of the rabbit, Oryctolagus cuniculus, carbonic anhydrases I and II, from the Iberian Peninsula and France. Biochem. Genet. 41, 391–404. https://doi.org/10.1023/B:BIGI.0000007774.39262.8e (2003).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
36.Geraldes, A., Ferrand, N. & Nachman, M. W. Contrasting patterns of introgression at X-linked loci across the hybrid zone between subspecies of the European rabbit (Oryctolagus cuniculus). Genetics 173, 919–933 (2006).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Lo Valvo, M., Scala, A. & Scalisi, M. Biometric characterization and taxonomic considerations of European rabbit Oryctolagus cuniculus (Linnaeus 1758) in Sicily (Italy). World Rabbit Sci. 22(3), 207–214. https://doi.org/10.4995/wrs.2014.1467 (2014).Article 

Google Scholar 
38.Miller, G. S. Catalogue of the Mammals of Western Europe in the Collection of the British Museum (Trustees of the British Museum, 1912).
Google Scholar 
39.Sharples, C. M., Fa, J. E. & Bell, D. J. Geographical variation in size in the European rabbit Oryctolagus cuniculus (Lagomorpha: Leporidae) in western Europe and North Africa. Zool. J. Linn. Soc-Lond. 117, 141–158. https://doi.org/10.1111/j.1096-3642.1996.tb02153.x (1996).Article 

Google Scholar 
40.Carro, F., Ortega, M. & Soriguer, R. C. Is restocking a useful tool for increasing rabbit densities?. Global Ecol. Conserv. 17, e00560. https://doi.org/10.1016/j.gecco.2019.e00560 (2019).Article 

Google Scholar 
41.Angulo, E. & Villafuerte, R. Modelling hunting strategies for the conservation of wild rabbit populations. Biol. Conserv. 115, 291–301 (2003).Article 

Google Scholar 
42.Delibes-Mateos, M., Delibes, M., Ferreras, P. & Villafuerte, R. Key role of European rabbits in the conservation of the Western Mediterranean Basin Hotspot. Conserv. Biol. 22, 1106–1117 (2008).PubMed 
Article 
PubMed Central 

Google Scholar 
43.Garrido, J. L., Ferreres, J. & Gortázar, C. Las especies cinegéticas españolas en el siglo XXI. (eds. Garrido, J. L., Ferreres, J. & Gortázar, C.)
(Independently published, Ciudad Real, Spain, 2019).
Google Scholar 
44.Ríos-Saldaña, C. et al. Control of the European rabbit in central Spain. Eur. J. Wildlife Res. 59, 573–580. https://doi.org/10.1007/s10344-013-0707-x (2013).Article 

Google Scholar 
45.Lees, A. C. & Bell, D. J. A conservation paradox for the 21st century: The European wild rabbit Oryctolagus cuniculus, an invasive alien and an endangered native species. Mammal Rev. 38, 304–320 (2008).Article 

Google Scholar 
46.Cooke, B. D. Rabbits: Manageable environmental pests or participants in new Australian ecosystems?. Wildlife Res. 39, 279–289 (2013).ADS 
Article 

Google Scholar 
47.Calvete, C., Angulo, E. & Estrada, R. Conservation of European wild rabbit populations when hunting is age and sex selective. Biol. Conserv. 121(4), 623–634 (2005).Article 

Google Scholar 
48.Delibes-Mateos, M., Ramírez, E., Ferreras, P. & Villafuerte, R. Translocations as a risk for the conservation of European wild rabbit Oryctolagus cuniculus lineages. Oryx 42(2), 259–264 (2008).Article 

Google Scholar 
49.Andersen, J. & Jensen, B. Studies on the European hare. XXVIII. The weight of the eye lens in the European hares of known age. Acta Theriol. 17, 87–92 (1972).Article 

Google Scholar 
50.Suchentrunck, F., Willing, R. & Hartl, G. B. On eye lens weights and other age criteria of the Brown hare (Lepus europaeus Pallas, 1778). Z. Säugetierkd. 56, 365–374 (1991).
Google Scholar 
51.Villafuerte, R. et al. Large-scale assessment of myxomatosis prevalence in European wild rabbits (Oryctolagus cuniculus) 60 years after first outbreak in Spain. Res. Vet. Sci. 114, 281–286 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
52.Rouco, C., Villafuerte, R., Castro, F. & Ferreras, P. Effect of artificial warren size on a restocked European wild rabbit population. Anim. Conserv. 14, 117–123 (2011).Article 

Google Scholar 
53.Southern, N. The ecology and population dynamics of the wild rabbit (Oryctolagus cuniculus). Ann. Appl. Biol. 27, 509–514 (1940).Article 

Google Scholar 
54.Dunnet, G. M. Growth rate of young rabbits, Oryctolagus cuniculus (L.). CSIRO Wildl. Res. 1, 66–67 (1956).Article 

Google Scholar 
55.Ferreira, A. & Ferreira, A. J. Post-weaning growth of endemic Iberian wild rabbit subspecies, Oryctolagus cuniculus algirus, kept in a semi-extensive enclosure: Implications for management and conservation. World Rabbit Sci. 22, 129–136. https://doi.org/10.4995/wrs.2014.1673 (2014).Article 

Google Scholar 
56.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/ (Vienna,
Austria, 2020).
Google Scholar 
57.du Sert, N. P. et al. Reporting animal research: Explanation and elaboration for the ARRIVE guidelines 2.0. PLoS Biol. 18(7), e3000411. https://doi.org/10.1371/journal.pbio.3000411 (2020).MathSciNet 
CAS 
Article 

Google Scholar 
58.Burnham, K. P. & Anderson, D. R. Monte Carlo insights and extended examples. In Model Selection and Multimodel Inference, (eds. Burnham K. P. &
Anderson D. R.) https://doi.org/10.1007/978-0-387-22456-5_5). (Springer, New York, NY, US, 2002).59.Pastore, M. Overlapping: A R package for estimating overlapping in empirical distributions. J. Open Source Softw. 32, 1023 (2018).ADS 
Article 

Google Scholar 
60.Pastore, M. & Calcagnì, A. Measuring distribution similarities between samples: A distribution-free overlapping Index. Front. Psychol. 10, 1089 https://doi.org/10.3389/fpsyg.2019.01089 (2019).Article 

Google Scholar 
61.Williams, C. & Moore, R. Phenotypic adaptation and natural selection in the wild rabbit, Oryctolagus cuniculus, Australia. J. Anim. Ecol. 58(2), 495–507. https://doi.org/10.2307/4844 (1989).Article 

Google Scholar  More

1.Fernández-Brime, S., Muggia, L., Maier, S., Grube, M. & Wedin, M. Bacterial communities in an optional lichen symbiosis are determined by substrate, not algal photobionts. FEMS Microbiol. Ecol. 95, fiz012 (2019).PubMed 
Article 
CAS 

Google Scholar 
2.Grube, M. et al. Exploring functional contexts of symbiotic sustain within lichen-associated bacteria by comparative omics. ISME J. 9, 412–424 (2015).CAS 
PubMed 
Article 

Google Scholar 
3.Spribille, T. et al. Basidiomycete yeasts in the cortex of ascomycete macrolichens. Science 353, 488–492 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
4.Peksa, O. & Škaloud, P. Do photobionts influence the ecology of lichens? A case study of environmental preferences in symbiotic green alga Asterochloris (Trebouxiophyceae). Mol. Ecol. 20, 3936–3948 (2011).PubMed 
Article 

Google Scholar 
5.Řídká, T., Peksa, O., Rai, H., Upreti, D. K. & Škaloud, P. Photobiont diversity in Indian Cladonialichens, with special emphasis on the geographical patterns. In Terricolous Lichens in India, Volume 1: Diversity Patterns and Distribution Ecology (eds Rai, H. & Upreti, D. K.) 53–71 (Springer, New York, 2014).
Google Scholar 
6.Rolshausen, G. et al. Expanding the mutualistic niche: Parallel symbiont turnover along climatic gradients. Proc. R. Soc. B 287, 1924 (2020).Article 

Google Scholar 
7.Kosecka, M. et al. Trentepohlialean algae (Trentepohliales, Ulvophyceae) show preference to selected mycobiont lineages in lichen symbioses. J. Phycol. 56, 979–993 (2020).CAS 
PubMed 
Article 

Google Scholar 
8.Muggia, L., Pérez-Ortega, S., Kopun, T., Zellnig, G. & Grube, M. Photobiont selectivity leads to ecological tolerance and evolutionary divergence in a polymorphic complex of lichenized fungi. Ann. Bot. 114, 463–475 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
9.Vančurová, L., Muggia, L., Peksa, O., Řídká, T. & Škaloud, P. The complexity of symbiotic interactions influences the ecological amplitude of the host: A case study in Stereocaulon (lichenized Ascomycota). Mol. Ecol. 27, 3016–3033 (2018).PubMed 
Article 
CAS 

Google Scholar 
10.Beck, A., Kasalicky, T. & Rambold, G. Myco-photobiontal selection in a Mediterranean cryptogam community with Fulgensia fulgida. New Phytol. 153, 317–326 (2002).Article 

Google Scholar 
11.Nelsen, M. P. & Gargas, A. Actin type intron sequences increase phylogenetic resolution: An example from Asterochloris (Chlorophyta: Trebouxiophyceae). Lichenologist 38, 435–440 (2006).Article 

Google Scholar 
12.Nelsen, M. P. & Gargas, A. Dissociation and horizontal transmission of co-dispersed lichen symbionts in the genus Lepraria (Lecanorales: Stereocaulaceae). New Phytol. 177, 264–275 (2008).CAS 
PubMed 
Article 

Google Scholar 
13.Škaloud, P. & Peksa, O. Evolutionary inferences based on ITS rDNA and actin sequences reveal extensive diversity of the common lichen alga Asterochloris. Mol. Phylogenet. Evol. 54, 36–46 (2010).PubMed 
Article 

Google Scholar 
14.Škaloud, P., Steinová, J., Řídká, T., Vančurová, L. & Peksa, O. Assembling the challenging puzzle of algal biodiversity: Species delimitation within the genus Asterochloris (Trebouxiophyceae, Chlorophyta). J. Phycol. 51, 507–527 (2015).PubMed 
Article 
CAS 

Google Scholar 
15.Steinová, J., Škaloud, P., Yahr, R., Bestová, H. & Muggia, L. Reproductive and dispersal strategies shape the diversity of mycobiont-photobiont association in Cladonia lichens. Mol. Phylogenet. Evol. 134, 226–237 (2019).PubMed 
Article 

Google Scholar 
16.Vančurová, L., Peksa, O., Němcová, Y. & Škaloud, P. Vulcanochloris (Trebouxiales, Trebouxiophyceae), a new genus of lichen photobiont from La Palma, Canary Islands, Spain. Phytotaxa 219, 118–132 (2015).Article 

Google Scholar 
17.Vančurová, L. et al. Symbiosis between river and dry lands: Phycobiont dynamics on river gravel bars. Algal Res. 51, 102062. https://doi.org/10.1016/j.algal.2020.102062 (2020).Article 

Google Scholar 
18.Yahr, R., Vilgalys, R. & Depriest, P. T. Strong fungal specificity and selectivity for algal symbionts in Florida scrub Cladonia lichens. Mol. Ecol. 13, 3367–3378 (2004).CAS 
PubMed 
Article 

Google Scholar 
19.Tschermak-Woess, E. Asterochloris phycobiontica gen. et spec. nov., der Phycobiont der Flechte Varicellaria carneonivea. Plant Syst. Evol. 135, 279–294 (1980).Article 

Google Scholar 
20.Moya, P. et al. Molecular phylogeny and ultrastructure of the lichen microalga Asterochloris mediterranea sp. Nov. from Mediterranean and Canary Islands ecosystems. Int. J. Syst. Evol. Microbiol. 65(6), 1838–1854 (2015).PubMed 
Article 
CAS 

Google Scholar 
21.Kim, J. I. et al. Asterochloris sejongensis sp. nov. (Trebouxiophyceae, Chlorophyta) from King George Island, Antarctica. Phytotaxa 295, 60–70 (2017).Article 

Google Scholar 
22.Kim, J. I. et al. Taxonomic study of three new Antarctic Asterochloris(Trebouxio-phyceae) based on morphological and molecular data. Korean Soc. Phycol. 35(1), 17–32 (2020).CAS 

Google Scholar 
23.Pino-Bodas, R. & Stenroos, S. Global biodiversity patterns of the photobionts associated with the genus Cladonia (Lecanorales, Ascomycota). Microb. Ecol. https://doi.org/10.1007/s00248-020-01633-3 (2020).Article 
PubMed 

Google Scholar 
24.Fernandez-Mendoza, F. et al. Population structure of mycobionts and photobionts of the widespread lichen Cetraria aculeata. Mol. Ecol. 20, 1208–1232 (2011).CAS 
PubMed 
Article 

Google Scholar 
25.Helms, G. Taxonomy and symbiosis in associations of Physciaceae and Trebouxia. Göttingen, Doctoral thesis. Germany: Georg-August Universität Göttingen (2003).26.Kroken, S. & Taylor, J. W. Phylogenetic species, reproductive mode, and specificity of the green alga Trebouxia forming lichens with the fungal genus Letharia. Bryologist 103, 645–660 (2000).CAS 
Article 

Google Scholar 
27.Leavitt, S. D. et al. Fungal specificity and selectivity for algae play a major role in determining lichen partnerships across diverse ecogeographic regions in the lichen- forming family Parmeliaceae (Ascomycota). Mol. Ecol. 24, 3779–3797 (2015).PubMed 
Article 

Google Scholar 
28.Magain, N., Miadlikowska, J., Goffinet, B., Sérusiaux, E. & Lutzoni, F. Macroevolution of Specificity in Cyanolichens of the Genus Peltigera Section Polydactylon (Lecanoromycetes, Ascomycota). Syst. Biol. 66, 74–99 (2016).
Google Scholar 
29.Mark, K. et al. Contrasting cooccurrence patterns of photobiont and cystobasidiomycete yeast associated with common epiphytic lichen species. New Phytol. 227, 1362–1375 (2020).CAS 
PubMed 
Article 

Google Scholar 
30.O’Brien, H. E., Miadlikowska, J. & Lutzoni, F. Assessing population structure and host specialization in lichenized cyanobacteria. New Phytol. 198, 557–566 (2013).PubMed 
Article 

Google Scholar 
31.Piercey-Normore, M. D. & DePriest, P. T. Algal switching among lichen symbionts. Am. J. Bot. 88, 1490–1498 (2001).CAS 
PubMed 
Article 

Google Scholar 
32.Yahr, R., Vilgalys, R. & DePriest, P. T. Geographic variation in algal partners of Cladonia subtenuis (Cladoniaceae) highlights the dynamic nature of a lichen symbiosis. New Phytol. 171, 847–860 (2006).CAS 
PubMed 
Article 

Google Scholar 
33.Muggia, L. et al. The symbiotic playground of lichen thalli—A highly flexible photobiont association in rock inhabiting lichens. FEMS Microbiol. Ecol. 85, 313–323 (2013).CAS 
PubMed 
Article 

Google Scholar 
34.Sadowska-Deś, A. D. et al. Integrating coalescent and phylogenetic approaches to delimit species in the lichen photobiont Trebouxia. Mol. Phylogenet. Evol. 76, 202–210 (2014).PubMed 
Article 

Google Scholar 
35.Wirtz, N. et al. Lichen fungi have low cyanobiont selectivity in maritime Antarctica. New Phytol. 160, 177–183 (2003).Article 

Google Scholar 
36.Ertz, D., Guzow-Krzemińska, B., Thor, G., Łubek, A. & Kukwa, M. Photobiont switching causes changes in the reproduction strategy and phenotypic dimorphism in the Arthoniomycetes. Sci Rep. 8, 4952 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
37.Marshall, W. A. & Chalmers, M. O. Airborne dispersal of Antarctic terrestrial algae and cyanobacteria. Ecography 20, 585–594 (1997).Article 

Google Scholar 
38.Printzen, C., Domaschke, S., Fernandez-Mendoza, F. & Perez-Ortega, S. Biogeography and ecology of Cetraria aculeata, a widely distributed lichen with a bipolar distribution. MycoKeys 6, 33–53 (2013).Article 

Google Scholar 
39.Pardo-De la Hoz, C. J. et al. Contrasting symbiotic patterns in two closely related lineages of trimembered lichens of the genus Peltigera. Front. Microbiol. 9, 2770–2770 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
40.Singh, G. et al. Fungal–algal association patterns in lichen symbiosis linked to macroclimate. New Phytol. 214, 317–329 (2017).PubMed 
Article 

Google Scholar 
41.Cordeiro, L. M. C. et al. Molecular studies of photobionts of selected lichens from the coastal vegetation of Brazil. FEMS Microbiol. Ecol. 54, 381–390 (2005).CAS 
PubMed 
Article 

Google Scholar 
42.Pullaiah, T. (ed.) (2019) Biodiversity in Bolivia in: Global Biodiversity Volume 4: Selected Countries in the Americas and Australia, Chapter 1. 1–574 (Apple Academic Press, 2001).43.Ibisch, P. L. & Mérida, G. Biodiversity: the richness of Bolivia. State of knowledge and conservation. 1–638 (Ministry of Sustainable Development, Editorial FAN, 2004).44.Werth, S. & Sork, V. L. Ecological specialization in Trebouxia (Trebouxiophyceae) photobionts of Ramalina menziesii (Ramalinaceae) across six range-covering ecoregions of western North America. Am. J. Bot. 101, 1127–1140 (2014).PubMed 
Article 

Google Scholar 
45.Rolshausen, G., Dal Grande, F., Sadowska-Deś, A. D., Otte, J. & Schmitt, I. Quantifying the climatic niche of symbiont partners in a lichen symbiosis indicates mutualistmediated niche expansions. Ecography 41, 1380–1392 (2018).Article 

Google Scholar 
46.Blaha, J., Baloch, E. & Grube, M. High photobiont diversity insymbioses of the euryoecious lichen Lecanora rupicola (Lecanoraceae, Ascomycota). Biol. J. Linn. Soc. 88, 283–293 (2006).Article 

Google Scholar 
47.Vargas Castillo, R. & Beck, A. Photobiont selectivity and specificity in Caloplaca species in a fog-induced community in the Atacama Desert, northern Chile. Fungal Biol. 116, 665–676 (2012).PubMed 
Article 

Google Scholar 
48.Bačkor, M., Klemová, K., Bačkorová, M. & Ivanova, V. Comparison of the phytotoxic effects of usnic acid on cultures of free-living alga Scenedesmus quadricauda and aposymbiotically grown lichen photobiont Trebouxia erici. J. Chem. Ecol. 36, 405–411 (2010).PubMed 
Article 
CAS 

Google Scholar 
49.Galloway, D. J. Lichen biogeography. In: Nash T. H. (ed.) Lichen Biology, 315–335 (Cambridge University Press, 2008).50.Dal Grande, F. et al. Environment and host identity structure communities of green algal symbionts in lichens. New Phytol. 217, 277–289 (2017).Article 
CAS 

Google Scholar 
51.Dal Grande, F., Widmer, I., Wagner, H. H. & Scheidegger, C. Vertical and horizontal photobiont transmission within populations of a lichen symbiosis. Mol. Ecol. 21, 3159–3172 (2012).Article 

Google Scholar 
52.Otálora, M. A. G. et al. Multiple origins of high reciprocal symbiotic specificity at an intercontinental spatial scale among gelatinous lichens (Collemataceae, Lecanoromycetes). Mol. Phylogenet. Evol. 56, 1089–1095 (2010).PubMed 
Article 
CAS 

Google Scholar 
53.Ruprecht, U., Fernández-Mendoza, F., Türk, R. & Fryday, A. High levels of endemism and local differentiation in the fungal and algal symbionts of saxicolous lecideoid lichens along a latitudinal gradient in southern South America. Lichenologist 52, 287–303 (2020).PubMed Central 
Article 
PubMed 

Google Scholar 
54.Ahti, T., Cladoniaceae in Flora Neotropica Monograph 78. 2–362 (New York Botanical Garden Press, 2000).55.Parnmen, S., Leavitt, S. D., Rangsiruji, A. & Lumbsch, H. T. Identification of species in the Cladia aggregata group using DNA barcoding (Ascomycota: Lecanorales). Phytotaxa 115, 1–14 (2013).Article 

Google Scholar 
56.Guzow-Krzemińska, B. et al. New species and records of lichens from Bolivia. Phytotaxa 397, 257–279 (2019).Article 

Google Scholar 
57.Sipman, H. J. M. Survey of Lepraria species with lobed thallus margins in the tropics. Herzogia 17, 23–35 (2004).
Google Scholar 
58.Saag, L., Saag, A. & Randlane, T. World survey of the genus Lepraria (Stereocaulaceae, lichenized Ascomycota). Lichenologist 41, 25–60 (2009).Article 

Google Scholar 
59.Guzow-Krzemińska, B. et al. Phylogenetic placement of Lepraria cryptovouauxii sp. nov. (Lecanorales, Lecanoromycetes, Ascomycota) with notes on other Lepraria species from South America. MycoKeys 53, 1–22 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
60.Orange, A., James, P. W. & White, F. J. Microchemical Methods for the Identification of Lichens. 1–101 (British Lichen Society, 2001).61.Cubero, O. F., Crespo, A., Fatehi, J. & Bridge, P. D. DNA extraction and PCR amplification method suitable for fresh, herbarium-stored, lichenized, and other fungi. Plant Syst. Evol. 216, 243–249 (1999).CAS 
Article 

Google Scholar 
62.White, T. J., Bruns, T., Lee, S. & Taylor, J. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics in PCR Protocols: A Guide to Methods and Applications (eds Innis, M. A., Gelfand, D. H., Sninsky, J. J. & White, T. J.) 345–322 (Academic Press, 1990).63.Sherwood, A. R., Garbary, D. J. & Sheath, R. G. Assessing the phylogenetic position of the Prasiolales (Chlorophyta) using rbcL and 18S rRNA gene sequence data. Phycologia 39, 139–146 (2000).Article 

Google Scholar 
64.Nelsen, M. P., Rivas Plata, E., Andrew, C. J., Lücking, R. & Lumbsch, H. T. Phylogenetic diversity of trentepohlialean algae associated with lichen-forming fungi. J. Phycol. 47, 282–290 (2011).PubMed 
Article 

Google Scholar 
65.Widmer, I., Dal Grande, F., Cornejo, C. & Scheidegger, C. Highly variable microsatellite markers for the fungal and algal symbionts of the lichen Lobaria pulmonaria and challenges in developing biont-specific molecular markers for fungal associations. Fungal Biol. 114, 538–544 (2010).CAS 
PubMed 
Article 

Google Scholar 
66.Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403–410 (1990).CAS 
Article 

Google Scholar 
67.Katoh, K., Misawa, K., Kuma, K. & Miyata, T. MAFFT: A novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 30, 3059–3066 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
68.Okonechnikov, K., Golosova, O. & Fursov, M. UGENE team. Unipro GENE: A unified bioinformatics toolkit. Bioinformatics 28, 1166–1167 (2012).CAS 
PubMed 
Article 

Google Scholar 
69.Huelsenbeck, J. P. & Ronquist, F. MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics 17, 754–755 (2001).CAS 
PubMed 
Article 

Google Scholar 
70.Stamatakis, A. RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
71.Chernomor, O., von Haeseler, A. & Minh, B. Q. Terrace aware data structure for phylogenomic inference from supermatrices. Syst. Biol. 65, 997–1008 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
72.Nguyen, L. T., Schmidt, H. A., von Haeseler, A. & Minh, B. Q. IQ-TREE: A fast and effective stochastic algorithm for estimating maximum likelihood phylogenies. Mol. Biol. Evol. 32, 268–274 (2015).CAS 
Article 

Google Scholar 
73.Miller, M. A., Pfeiffer, W. & Schwartz, T. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. Gateway Computing Environments Workshop (GCE): 1–8 (2010).74.Lanfear, R., Frandsen, P. B., Wright, A. M., Senfeld, T. & Calcott, B. PartitionFinder 2: New methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Mol. Biol. Evol. 34, 772–773 (2016).
Google Scholar 
75.Rambaut, A. FigTreev1.4.2. Retrieved from: http://tree.bio.ed.ac.uk/software/figtree/ (2006–2014).76.Oksanen, et al., vegan: Community ecology package manual. Retrieved from https://cran.rproject.org/package=vegan (2017).77.Fick, S. E. & Hijmans, R. J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
78.Borcard, D., Legendre, P., Avois-Jacquet, C. & Tuomisto, H. Dissecting the spatial structure of ecological data at multiple scales. Ecology 85, 1826–1832 (2004).Article 

Google Scholar 
79.Lefeuvre, P. BoSSA: A bunch of structure and sequence analysis manual. Retrieved from https://cran.r-project.org/package=BoSSA (2018).80.McArdle, B. H. & Anderson, M. J. Fitting multivariate models to community data: A comment on distance-based redundancy analysis. Ecology 82, 290–297 (2001).Article 

Google Scholar 
81.Stenroos, S., Pino-Bodas, R., Hyvönen, J., Lumbsch, T. H. & Ahti, T. Phylogeny of family Cladoniaceae (Lecanoromycetes, Ascomycota) based on sequences of multiple loci. Cladistics 35, 351–384 (2019).Article 

Google Scholar 
82.R Core Team. R: A Language and Environment for Statistical Computing.83.https://www.R-project.org/ (2017).84.RStudio Team. RStudio: Integrated Development for R. RStudio, Inc., Boston, MA URL http://www.rstudio.com/ (2018).85.Aktas, C. Manipulating DNA Sequences and Estimating Unambiguous Haplotype Network with Statistical Parsimony. Retrieved from https://CRAN.R project.org/package=haplotypes (2020).86.Piel, W. H., Chan, L., Dominus, M. J., Ruan, J., Vos, R. A., and V. Tannen. TreeBASE v. 2: A Database of Phylogenetic Knowledge. In: e-BioSphere (2009) More

1.Bâki Iz, H. Is the global sea surface temperature rise accelerating?. Geod. Geodyn. 9, 432–438 (2018).Article 

Google Scholar 
2.Bopp, L. et al. Multiple stressors of ocean ecosystems in the 21st century: Projections with CMIP5 models. Biogeosciences 10, 6225–6245 (2013).ADS 
Article 

Google Scholar 
3.Sarmiento, J. L. et al. Response of ocean ecosystems to climate warming. Global Biogeochem. Cycles 18 (2004).4.Oliver, E. C. J. et al. Longer and more frequent marine heatwaves over the past century. Nat. Commun. 9, 1–12 (2018).CAS 
Article 

Google Scholar 
5.Frölicher, T. L., Fischer, E. M. & Gruber, N. Marine heatwaves under global warming. Nature 560, 360–364 (2018).ADS 
PubMed 
Article 
CAS 

Google Scholar 
6.Beaugrand, G. & Kirby, R. R. How do marine pelagic species respond to climate change? Theories and observations. Ann. Rev. Mar. Sci. 10, 169–197 (2018).PubMed 
Article 

Google Scholar 
7.Verity, P. G., Smetacek, V. & Smayda, T. J. Status, trends and the future of the marine pelagic ecosystem. Environ. Conserv. 29, 207–237 (2002).Article 

Google Scholar 
8.Palko, B. J., Beardsley, G. L. & Richards, W. J. Synopsis of the biological data on dolphin-fishes, Coryphaena hippurus Linnaeus and Coryphaena equiselis Linnaeus. NOAA Tech. Rep. NMFS Circ. 443, 1–28 (1982).
Google Scholar 
9.Oxenford, H. A. Biology of the dolphinfish (Coryphaena hippurus) in the western central Atlantic: A review. Sci. Mar. 63, 277–301 (1999).Article 

Google Scholar 
10.Moltó, V. et al. A global review on the biology of the dolphinfish (Coryphaena hippurus) and its fishery in the Mediterranean Sea: advances in the last two decades. Rev. Fish. Sci. Aquac. (2020).11.FAO. Coryphaena hippurus (Linnaeus, 1758). Species fact sheets. http://www.fao.org/fishery/species/3130/en (2019).12.Morales-Nin, B., Cannizzaro, L., Massuti, E., Potoschi, A. & Andaloro, F. An overview of the FADs fishery in the Mediterranean Sea. Proc. Tuna Fish. Fish Aggreg. Dev. Symp. 184–207 (2000).13.Morales-Nin, B. Mediterranean FADs fishery: An overview. In Second International Symposium on Tuna Fisheries and Fish Aggregating Devices (2011).14.Giorgi, F. Climate change hot-spots. Geophys. Res. Lett. 33, L08707 (2006).ADS 
Article 

Google Scholar 
15.Durrieu de Madron, X. et al. Marine ecosystems’ responses to climatic and anthropogenic forcings in the Mediterranean. Prog. Oceanogr. 91, 97–166 (2011).16.Adloff, F. et al. Mediterranean sea response to climate change in an ensemble of twenty first century scenarios. Clim. Dyn. 45, 2775–2802 (2015).Article 

Google Scholar 
17.Darmaraki, S. et al. Future evolution of marine heatwaves in the mediterranean sea. Clim. Dyn. 53, 1371–1392 (2019).Article 

Google Scholar 
18.Bignami, S., Sponaugle, S. & Cowen, R. K. Effects of ocean acidification on the larvae of a high-value pelagic fisheries species, mahi-mahi Coryphaena hippurus. Aquat. Biol. 21, 249–260 (2014).Article 

Google Scholar 
19.Norton, J. G. Apparent habitat extensions of dolphinfish (Coryphaena hippurus) in response to climate transients in the California current*. Sci. Mar. 63, 239–260 (1999).Article 

Google Scholar 
20.Chang, S.-K. & Maunder, M. N. Aging material matters in the estimation of von Bertalanffy growth parameters for dolphinfish (Coryphaena hippurus). Fish. Res. 119–120, 147–153 (2012).Article 

Google Scholar 
21.Furukawa, S. et al. Age, growth, and reproductive characteristics of dolphinfish Coryphaena hippurus in the waters off west Kyushu, northern East China Sea. Fish. Sci. 78, 1153–1162 (2012).CAS 
Article 

Google Scholar 
22.Asch, R. G., Stock, C. A. & Sarmiento, J. L. Climate change impacts on mismatches between phytoplankton blooms and fish spawning phenology. Glob. Chang. Biol. 25, 2544–2559 (2019).ADS 
PubMed 
Article 

Google Scholar 
23.Shoji, J. et al. Possible effects of global warming on fish recruitment: shifts in spawning season and latitudinal distribution can alter growth of fish early life stages through changes in daylength. ICES J. Mar. Sci. 68, 1165–1169 (2011).Article 

Google Scholar 
24.R Core Team. R: A Language and Environment for Statistical Computing. Version 3.6.2. https://www.R-project.org/ (R Foundation for Satistical Computing, 2019).25.Wickham, H. ggplot2. Wiley Interdiscip. Rev. Comput. Stat. 3, 180–185 (2011).Article 

Google Scholar 
26.Morrongiello, J. R., Thresher, R. E. & Smith, D. C. Aquatic biochronologies and climate change. Nat. Clim. Chang. 2, 849–857 (2012).ADS 
Article 

Google Scholar 
27.Schismenou, E. et al. Seasonal changes in otolith increment width trajectories and the effect of temperature on the daily growth rate of young sardines. Fish. Oceanogr. 25, 362–372 (2016).Article 

Google Scholar 
28.Schismenou, E. et al. Disentangling the effects of inherent otolith growth and model-simulated ecosystem parameters on the daily growth rate of young anchovies. Mar. Ecol. Prog. Ser. 515, 227–237 (2014).ADS 
Article 

Google Scholar 
29.Catalán, I. A. et al. Daily otolith growth and ontogenetic geochemical signatures of age-0 anchovy (Engraulis encrasicolus) in the gulf of cádiz (SW Spain). Mediterr. Mar. Sci. 15, 781–789 (2014).Article 

Google Scholar 
30.Tanner, S. E. et al. Regional climate, primary productivity and fish biomass drive growth variation and population resilience in a small pelagic fish. Ecol. Indic. 103, 530–541 (2019).Article 

Google Scholar 
31.Ito, S., Okunishi, T., Kishi, M. J. & Wang, M. Modelling ecological responses of Pacific saury (Cololabis saira) to future climate change and its uncertainty. ICES J. Mar. Sci. 70, 980–990 (2013).Article 

Google Scholar 
32.Vinagre, C., Ferreira, T., Matos, L., Costa, M. J. & Cabral, H. N. Latitudinal gradients in growth and spawning of sea bass, Dicentrarchus labrax, and their relationship with temperature and photoperiod. Estuar. Coast. Shelf Sci. 81, 375–380 (2009).ADS 
Article 

Google Scholar 
33.Suthers, I. M. & Sundby, S. Role of the midnight sun: Comparative growth of pelagic juvenile cod (Gadus morhua) from the Arcto-Norwegian and a Nova Scotian stock. ICES J. Mar. Sci. 53, 827–836 (1996).Article 

Google Scholar 
34.Pepin, P. et al. Once upon a larva: Revisiting the relationship between feeding success and growth in fish larvae. ICES J. Mar. Sci. 72, 359–373 (2015).Article 

Google Scholar 
35.Fablet, R. et al. Shedding light on fish otolith biomineralization using a bioenergetic approach. PLoS ONE 6, e27055 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Lorenzen, K. Toward a new paradigm for growth modeling in fisheries stock assessments: Embracing plasticity and its consequences. Fish. Res. 180, 4–22 (2016).Article 

Google Scholar 
37.Campos-Candela, A., Palmer, M., Balle, S., Álvarez, A. & Alós, J. A mechanistic theory of personality-dependent movement behaviour based on dynamic energy budgets. Ecol. Lett. 22, 213–232 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
38.Payne, M. R. et al. Uncertainties in projecting climate-change impacts in marine ecosystems. ICES J. Mar. Sci. 73, 1272–1282 (2016).Article 

Google Scholar 
39.Fernandes, J. A. et al. Can we project changes in fish abundance and distribution in response to climate? Glob. Chang. Biol. (2020).40.Ramírez-Romero, E. et al. Assessment of the skill of coupled physical-biogeochemical models in the NW Mediterranean. Front. Mar. Sci. (2020).41.Rountrey, A. N., Coulson, P. G., Meeuwig, J. J. & Meekan, M. Water temperature and fish growth: Otoliths predict growth patterns of a marine fish in a changing climate. Glob. Chang. Biol. 20, 2450–2458 (2014).ADS 
PubMed 
Article 

Google Scholar 
42.Moltó, V., Ospina-Alvarez, A., Gatt, M., Palmer, M. & Catalán, I. A. A Bayesian approach to recover the theoretical temperature-dependent hatch date distribution from biased samples: The case of the common dolphinfish (Coryphaena hippurus). Preprint at: https://arxiv.org/abs/2004.01000 (2020).43.Catalán, I. A. et al. Critically examining the knowledge base required to mechanistically project climate impacts: A case study of Europe’s fish and shellfish. Fish Fish. 1–17 (2019).44.Morrongiello, J. R., Walsh, C. T., Gray, C. A., Stocks, J. R. & Crook, D. A. Environmental change drives long-term recruitment and growth variation in an estuarine fish. Glob. Chang. Biol. 20, 1844–1860 (2014).ADS 
PubMed 
Article 

Google Scholar 
45.Baudron, A. R., Needle, C. L., Rijnsdorp, A. D. & Tara Marshall, C. Warming temperatures and smaller body sizes: Synchronous changes in growth of North Sea fishes. Glob. Chang. Biol. 20, 1023–1031 (2014).46.Pauly, D. & Cheung, W. W. L. Sound physiological knowledge and principles in modeling shrinking of fishes under climate change. Glob. Chang. Biol. 24, e15–e26 (2018).ADS 
PubMed 
Article 

Google Scholar 
47.Wenger, A. S., Whinney, J., Taylor, B. & Kroon, F. The impact of individual and combined abiotic factors on daily otolith growth in a coral reef fish. Sci. Rep. 6, 1–10 (2016).Article 
CAS 

Google Scholar 
48.García, A. et al. Climate-induced environmental conditions influencing interannual variability of Mediterranean bluefin (Thunnus thynnus) larval growth. Fish. Oceanogr. 22, 273–287 (2013).Article 

Google Scholar 
49.Pimentel, M., Pegado, M., Repolho, T. & Rosa, R. Impact of ocean acidification in the metabolism and swimming behavior of the dolphinfish (Coryphaena hippurus) early larvae. Mar. Biol. 161, 725–729 (2014).CAS 
Article 

Google Scholar 
50.FAO-CopeMed II. Report of the CopeMed II-MedSudMed Workshop on the Status of Coryphaena hippurus Fisheries in the Western-Central Mediterranean, Cádiz, Spain, 8–9 October 2019. CopeMed Technical Documents No. 54 (GCP/INT/028SPA-GCP/INT/362/EC). 1–22 (2019).51.Tittensor, D. P. et al. A protocol for the intercomparison of marine fishery and ecosystem models: Fish-MIP v1. 0. Geosci. Model Dev. 11, 1421–1442 (2018).52.Massutí, E. & Morales-Nin, B. Reproductive biology of dolphin-fish (Coryphaena hippurus L.) off the island of Majorca (western Mediterranean). Fish. Res. 30, 57–65 (1997).53.Massutí, E. & Morales-Nin, B. Seasonality and reproduction of dolphin-fish (Coryphaena hippurus) in the Western Mediterranean*. Sci. Mar. 59, 357–364 (1995).
Google Scholar 
54.Potoschi, A., Reñones, O. & Cannizzaro, L. Sexual development, maturity and reproduction of dolphinfish (Coryphaena hippurus) in the western and central Mediterranean*. Sci. Mar. 63, 367–372 (1999).Article 

Google Scholar 
55.Alemany, F. et al. Influence of physical environmental factors on the composition and horizontal distribution of summer larval fish assemblages off Mallorca island (Balearic archipelago, western Mediterranean). J. Plankton Res. 28, 473–487 (2006).Article 

Google Scholar 
56.Torres, A. P. et al. Decapod crustacean larval communities in the Balearic Sea (western Mediterranean): Seasonal composition, horizontal and vertical distribution patterns. J. Mar. Syst. 138, 112–126 (2014).Article 

Google Scholar 
57.Massutí, E., Deudero, S., Sánchez, P. & Morales-Nin, B. Diet and Feeding of Dolphin (Coryphaena hippurus) in Western Mediterranean Waters. Bull. Mar. Sci. 63, 329–341 (1998).
Google Scholar 
58.Merten, W., Appeldoorn, R., Rivera, R. & Hammond, D. Diel vertical movements of adult male dolphinfish (Coryphaena hippurus) in the western central atlantic as determined by use of pop-up satellite archival transmitters. Mar. Biol. 161, 1823–1834 (2014).Article 

Google Scholar 
59.D’Ortenzio, F. & D’Alcalà, M. R. On the trophic regimes of the Mediterranean Sea: A satellite analysis. Biogeosciences 5, 2959–2983 (2008).
Google Scholar 
60.IPCC. Climate change 2014: Synthesis report. In Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. (Core Writing Team, Pachauri, R.K., Meyer, L.A. eds.). (IPCC, 2014).61.Grazzini, F. & Viterbo, P. Record-breaking warm sea surface temperature of the Mediterranean Sea. ECMWF Newsl. 98, 30–31 (2003).
Google Scholar 
62.Olita, A., Sorgente, R., Ribotti, A., Natale, S. & Gaberšek, S. Effects of the 2003 European heatwave on the Central Mediterranean Sea surface layer: a numerical simulation. Eur. Geosci. Union 3, 85–125 (2006).
Google Scholar 
63.Garrabou, J. et al. Mass mortality in Northwestern Mediterranean rocky benthic communities: Effects of the 2003 heat wave. Glob. Chang. Biol. 15, 1090–1103 (2009).ADS 
Article 

Google Scholar 
64.Hobday, A. J. et al. A hierarchical approach to defining marine heatwaves. Prog. Oceanogr. 141, 227–238 (2016).ADS 
Article 

Google Scholar 
65.Hobday, A. J. et al. Categorizing and naming marine heatwaves. Oceanography 31, 162–173 (2018).Article 

Google Scholar 
66.Schlegel, R. W. Marine Heatwave Tracker. http://www.marineheatwaves.org/tracker. https://doi.org/10.5281/zenodo.3787872 (2020).67.Ricker, W. E. Computation and interpretation of biological statistics of fish populations. Bull. Fish. Res. Bd. Can. 191, 1–382 (1975).68.Solano-Fernández, M., Montoya-Márquez, J. A., Gallardo-Cabello, M. & Espino-Barr, E. Age and growth of the Dolphinfish Coryphaena hippurus in the coast of Oaxaca and Chiapas, Mexico. Rev. Biol. Mar. Oceanogr. 50, 491–505 (2015).Article 

Google Scholar 
69.Höhne, L. et al. Environmental determinants of perch (Perca fluviatilis) growth in gravel pit lakes and the relative performance of simple versus complex ecological predictors. Ecol. Freshw. Fish 00, 1–17 (2020).
Google Scholar 
70.Kuhn, M. caret: Classification and Regression Training. R package. Version 6.0-86. https://CRAN.R-project.org/package=caret (2020).71.Su, Y.-S. & Yajima, M. R2jags: Using R to Run ‘JAGS’. R Package Version 0.5-7. https://CRAN.R-project.org/package=R2jags (2015).72.Plummer, M. rjags: Bayesian Graphical Models Using MCMC. R Package Version 4-10. https://CRAN.R-project.org/package=rjags (2015).73.Then, A. Y., Hoenig, J. M., Hall, N. G. & Hewitt, D. A. Evaluating the predictive performance of empirical estimators of natural mortality rate using information on over 200 fish species. ICES J. Mar. Sci. 72, 82–92 (2015).Article 

Google Scholar 
74.Massutí, E., Morales-Nin, B. & Moranta, J. Otolith microstructure, age, and growth patterns of dolphin, Coryphaena hippurus, in the western Mediterranean. Fish. Bull. 97, 891–899 (1999).
Google Scholar 
75.Copemed II. Report of the CopeMed II-MedSudMed Workshop on Stock Assessment of Coryphaena hippurus in the Western-Central Mediterranean. Málaga, Spain 13–15 September 2016. Copemed II Technical Documents No. 44 (GCP/INT/028/SPA – GCP/INT/006/EC). Málaga, 2016. 1–31. http://www.faocopemed.org/pdf/publications/CopeMedII_TD44.pdf (2016). More

Selection and initial metabolic profiling of organismsIn order to maximize the chance of obtaining communities with diverse taxonomic profiles from different environmental compositions, the organisms selected were drawn from a number of bacterial taxa known to employ varying metabolic strategies. In addition, given the growing relevance of synthetic microbial communities to industrial and biotechnological applications73,74,75,76, we chose to employ bacterial species that have previously been used as model organisms and have well-characterized metabolic capabilities. This criterion, paired with the availability of flux-balance models associated with a majority of these organisms, allows us to explore the metabolic mechanisms observed in our various experimental conditions with higher confidence. These selection principles resulted in a set of 15 candidate bacterial organisms (Acinetobacter baylyi, Bacillus licheniformis, Bacillus subtilis, Corynebacterium glutamicum, Escherichia coli, Lactococcus lactis, Methylobacterium extorquens, Pseudomonas aeruginosa, Pseudomonas fluorescens, Pseudomonas putida, Salmonella enterica, Streptomyces coelicolor, Shewanella oneidensis, Streptococcus thermophilus, and Vibrio natriegens) spanning three bacterial phyla (Actinobacteria, Firmicutes, and Proteobacteria, Supplementary Table 1, Supplementary Fig. 1a).A microtiter plate-based phenotypic assay was used to assess the metabolic capabilities of each of the 15 candidate organisms. Each organism, stored in glycerol at −80 °C, was initially grown in 3 mL of Miller’s LB broth (Sigma–Aldrich, St. Louis, MO) for 18 h with shaking at 300 rpm at each organism’s recommended culturing temperature (Supplementary Table 1). To maximize oxygenation of the cultures and prevent biofilm formation, culture tubes were angled at 45° during this initial growth phase. Candidate organism Streptococcus thermophilus was found to have produced too little biomass in this time period and was grown for an additional 8 h. Each culture was then separately washed three times by centrifuging at 6000 × g for 2 min, removing the supernatant, suspending the pellet in 1 mL of M9 minimal medium with no carbon source, and vortexing or triturating to homogenize. The cultures were then diluted to OD600 0.5 ± 0.1 as read by a microplate reader (BioTek Instruments, Winooski, VT) and distributed into each well of three PM1 Phenotype MicroArray Plates (Biolog Inc., Hayward, CA) per organism at final OD600 of 0.05 ± 0.01. The carbon sources in the PM1 plates (Supplementary Table 2, Supplementary Fig. 1b) were resuspended in 150 µl of M9 minimal media prepared from autoclaved M9 salts (BD, Franklin Lakes, NJ) and filter-sterilized MgSO4 and CaCl2 prior to inoculation. The cultures in each PM1 plate were incubated at each organism’s recommended culturing temperature with shaking at 300 rpm for 48 h. After this growing period, the OD600 of each culture was measured by a microplate reader to quantify growth. To account for evaporation in the outer wells of the plates, which could yield in inflated OD readings, three ‘evaporation control’ plates with no carbon source were inoculated with bacteria at a final OD600 of 0.05 and incubated at 30 °C for 48 h. The averaged OD600 readings of these plates were subtracted from the readings of the bacterial growth plates to correct for evaporation. A one-tailed t-test was performed using these corrected OD600 values to determine significance of growth above the value of the negative controls (p  More

Song description: terminologySong organisationBlue whale vocal sequences are traditionally referred to as ‘calls’19,20,21, however, as they meet the criterion of ‘song’ as used in the bioacoustic community22, in this study we use the term ‘song’ to refer to regularly-repeated whale vocalisations. The song is repeated in a sequence with regular intervals, defined as the Inter-Call Interval (ICI), measured as the time interval between the beginning of the one song and the beginning of the following song. Note that although we use the term ‘song’, we chose to keep the definition ‘ICI’ as this nomenclature is used traditionally in the whale literature, rather than ‘ISI’, which usually designates Inter-Series (or Sequence) Interval. Songs are composed of units and we used the term ‘unit’ to designate parts of the song that are separated by a silence (see reviewed criteria in23). Units were divided into subunits: subunits are defined as such when there is a sudden change in the sound structure for instance becoming harmonic or noisy.Sound typesA sound can be of different types: (1) the simpler one is the simple tone, which is either pure, with the same frequency all along, or showing frequency and/or amplitude modulations; (2) harmonic sounds are sounds with multiple tones at frequencies that are integer multiples of the frequency of the original wave, called the fundamental frequency ((F_{0})). When one of the harmonics has a greater amplitude than the others, it is called ‘resonance frequency’; (3) pulsed sounds are, as defined in24, the repetition of similar “pulses” or short signals with a constant pulse rate, often aurally perceived by humans as amplitude modulated sounds. On spectrogram representation, using a long analysis time window, these sounds are characterized by sidebands with regular spacing. The frequency difference ((Delta f)) between each sideband is the pulse rate of the sound. In their recent study, Patris et al. made the difference between what they defined as ‘tonal pulsed sounds’ and ‘non-tonal pulsed sounds’24. Following their criterium, the sidebands of the tonal pulsed sounds show a harmonic relationship, meaning that the frequency of each sideband divided by the pulsed rate is a positive integer. If it is not the case, then the sound is a non-tonal pulsed sound.Nonlinear phenomenaNonlinear phenomena are observed in a variety of birds25, anurans26 and mammals27,28, including marine mammals (e.g., manatee29) and more particularly cetaceans (right whales30,31, killer whales30,32 and humpback whales33). They have been well described by a variety of authors27,34 and include: (1) frequency jumps, that are characterized by sudden (F_{0}) changes which moves up or down abruptly and discontinuously, and is different from continuous, smooth modulation27; (2) subharmonics, that are additional spectral components and can suddenly appear at integer fractional values of an identifiable (F_{0}) (e.g., (F_{0}/2), (F_{0}/3, ldots)) and as harmonics of these values. On a spectrogram, it results as bands of energy evenly spaced below (F_{0}) and between its harmonics throughout the spectrum; (3) biphonation, that is the simultaneous occurrence of two independent fundamental frequencies (F_{0}) and (G_{0}). Biphonation can be visible on a spectrogram as two distinct frequency contours35. Alternatively, if one source ((F_{0})) vibrates at a much lower frequency than the other ((G_{0})), biphonation will appear as visible sidebands at linear combinations of (F_{0}) and (G_{0}) (m(G_{0}) ± n(F_{0}), where m and n are integers), because the airflow is then modulated by the frequency difference. This is equivalent to considering that the lower (F_{0}) amplitude-modulates the higher frequency (G_{0}) (carrier frequency)28; (4) finally, deterministic chaos are broadband, noise-like segments. These episodes of non-random noise appear via abrupt transitions and can also contain some periodic energy, which appears as banding in a spectrogram. In extreme cases there are no repeating periods at all27,34.Analysis of the Chagos song and comparison with the Indian Ocean pygmy blue whale song types and Omura’s whale song typesChagos songThe Chagos song was composed of 3 units (Fig. 3). The 3-unit song was repeated in stereotyped series with an ICI of (190.79 pm 1.49) s (Fig. 7b).The first unit of the Chagos song is divided into 3 subunits (Fig. 3): in 2017, subunit 1 was pulsed with a rate (Delta f_{u1su1}) = 3.22 ± 0.01 Hz. Using Patris et al. ’s criterion24, we concluded that this subunit is a non-tonal pulsed sound, since the sidebands do not have a harmonic relationship. The carrier frequency (where the peak of energy lies) was 35.74 ± 0.02 Hz for 73% of the measured songs, 32.47 ± 0.05 Hz for 23% of the songs and 38.9 ± 0.06 Hz for 3% of the measured songs. One song had a carrier frequency of 29.18 Hz. This subunit 1 lasted 3.02 ± 0.03 s in duration. Subunit 2 was often less obvious (likely due to propagation effects, lower source level or possibly to deterministic chaos) so that it could not be measured for all of the songs sampled; it is also a short (1.53 ± 0.05 s) non-tonal pulsed unit with a pulse rate ((Delta f_{u1su2})) of approximately 3 Hz and a slightly different carrier frequency, induced by a frequency jump. The carrier frequency was of 36.02 ± 0.03 Hz for 87% of the measurements, 39.16 ± 0.05 Hz for 8% of the measured songs and 32.97 ± 0.1 Hz for 5%. Finally, subunit 3 was a tonal unit showing a frequency modulation. The subunit started at 29.55 ± 0.02 Hz down to 29.35 ± 0.02 Hz over approximately 3.5 s, then down to 28.10 ± 0.09 Hz as a decrease to 27.62 ± 0.04 Hz over 3 s. The total duration of this subunit was 6.40 ± 0.07 s, and the total duration of the unit 1 was 11.36 ± 0.08 s.Unit 2 was a pure tone following after a silence of 3.06 ± 0.1 s. Its peak frequency was 22.34 ± 0.05 Hz and its duration was 3.24 ± 0.07 s. Finally, unit 3, also a pure tone, followed after a silence of 14.38 ± 0.23 s. It had a peak frequency of 17.44 ± 0.05 Hz and lasted 2.94 ± 0.15 s. The third unit was sometimes absent. This could be due to a variation in the song or due to propagation losses. When unit 3 was present, the total song duration was 34.38 ± 0.4 s.The frequency for the beginning of the third subunit of the unit 1 of the Chagos song (point 1 in Fig. 3a) decreased by approximatively 0.33 Hz/year across years (Fig. 4).This phenomenon will be examined in details in a further study.Figure 3Spectrogram (a) and waveforms (b) of a Chagos song recorded on the eastern side of the Chagos Archipelago (DGS) in August 2017. Detailed waveforms show the signal structure of the units within the song. Spectrogram parameters: Hamming window, 1024-point FFT length, 90% overlap. Note that the axes differ among plots. (c) Measurements (mean ± standard error (s.e.)) of the acoustic features. N is the number of measurements, (u_{i}su_{j}) stands for (unit_{i} subunit_{j}) where i and j are the unit and subunit numbers, (Delta f) designates the frequency difference between the sidebands, f and d are the frequency and duration of the feature indicated in subscript, and when present, the number in brackets refers to the point measured as indicated on the spectrogram. (F_{x}) or (G_{x}) designate the xth harmonic of a sound, and Cf designates the carrier frequency of a sound.Full size imageFigure 4The decline in frequency of the Chagos song from 2002 to 2017: spectrogram representation of five songs recorded at Diego Garcia in years 2002, 2005, 2012, 2015 and 2017. Spectrogram parameters: Hamming window, 1024-point FFT length, 90% overlap.Full size imageIndian Ocean pygmy blue whale songsThis section describes the structural, temporal and frequency features of the pygmy blue whale song-types commonly reported in the Indian Ocean. Note that as the frequency of at least parts of these songs are known to vary within and across years36,37,38,39,40,41, the frequency values obtained here are only valid for the years sampled.Madagascan pygmy blue whale The Madagascan pygmy blue whale song had 2 units (Fig. 5a). Unit 1 was divided into 2 subunits. In 2004, subunit 1 was a noisy pulsed sound, characteristic of deterministic chaos, with a pulse rate (Delta f_ {u1su1}) = 1.44 ± 0.01 Hz and of 4.76 ± 0.005 s duration. Subunit 2 was a tonal sound with harmonics. Its (F_ {0}), estimated as the mean frequency difference between the harmonics, was 7.04 ± 0.005 Hz. The maximum energy was in the (F_ {5}) (resonance frequency), which commenced at 35.31 ± 0.02 Hz and remained stable over 10.65 ± 0.13 s ((F_{5_{u1su2}}) in Fig. 5a). The frequency then remained stable over another 3.00 ± 0.16 s or in some songs increased to 35.91 ± 0.05 Hz [range = 34.84–37.05 Hz]. The total duration of subunit 2 was 13.65 ± 0.12 s, and unit 1 was 18.41 ± 0.15 s.Unit 2 followed after 27.74 ± 0.13 s. It had 2 subunits. Subunit 1 was a noisy pulsed sound, identified as deterministic chaos, it had a pulsed rate of (Delta f_ {u1su1}) = 1.25 ± 0.017 Hz, and a duration of 3.30 ± 0.05 s. Subunit 2 was a complex harmonic-like signal, with sidebands spaced by (Delta f_ {u2su2}) = 1.39 ± 0.003 Hz. Calculations of the ratio of the sideband frequencies over (Delta f) show that these 1.39 Hz-spaced bands do not have a harmonic relationship. However, relatively higher energy lies in frequency bands that have a harmonic relationship, where the band with the greatest energy started at 25.11 ± 0.02 Hz and ended at 24.33 ± 0.02 Hz ((G_{3_{u2su2}}) on Fig. 5). On the low signal-to-noise ratio (SNR) songs, only the harmonic bands were visible, this explains why this unit has been described previously as a harmonic signal when it is not7. The complex structure of subunit 2 can be explained by a phenomenon of biphonation, where there are two concurrent frequencies, with a lower fundamental frequency ((F_{0})) of 1.39 Hz, a higher fundamental frequency ((G_{0})) of 8.37 Hz (resonance frequency (G_{3}) starting at 25.11 Hz), and the sidebands at m(G_{0}) ± n(F_{0}) consistent with the amplitude modulation of (G_{0}) by (F_{0}). This biphonation event lasted for 16.04 ± 0.19 s. Finally, subunit 2 ended in a tonal sound with the harmonics ((G_{0}) = 7.94 Hz ± 0.003 Hz), that decreased in frequency from 24.30 ± 0.02 Hz to 23.05 ± 0.04 Hz over 4.88 ± 0.11 s (measured for the harmonic where there is the greatest energy ((G_{3_{u2su2}}))). Unit 2 was 23.63 ± 0.73 s in duration. The total duration of the Madagascan pygmy blue whale song was 68.68 ± 0.34 s.Sri Lankan pygmy blue whale The Sri Lankan pygmy blue whale song had 3 units (Fig. 5b). In 2009, unit 1 was a pulsed, non-tonal sound of a duration of 22.25 ± 0.11 s. The pulse rate was (Delta f_{u1}) = 3.28 ± 0.09 Hz. The carrier frequency of unit 1 started at 29.87 ± 0.09 Hz (‘Cf’ on Fig. 5b), and slightly down swept to 29.68 ± 0.09 Hz over 4.57 ± 0.06 s, then the frequency decreased to 25.85 ± 0.09 Hz over 17.68 ± 0.09 s.Unit 2 followed after 16.45 ± 0.12 s of silence. Unit 2 was a tonal sound with harmonics spaced by 12.21 ± 0.08 Hz. The maximum of energy was in the (F_{5_{u2}}) and started at 56.55 ± 0.12 Hz, increased to 60.63 ± 0.03 Hz over 4.87 ± 0.09 s, then increased to 60.80 ± 0.03 Hz over 8.80 ±0.09 s, and finally increased sharply to 70.13 ± 0.18 Hz overe 0.92 ± 0.06 s. Unit 2 was 14.60 ± 0.07 s in duration.Unit 3 followed after 2.20 ± 0.06 s of silence. It started as a non-tonal pulsed sound lasting 4.46 ± 0.01 s, with a pulse rate (Delta f_{u3}) = 3.29 ± 0.12 Hz and a carrier frequency starting at 103.47 ± 0.05 Hz and slightly decreasing to 102.91 ± 0.03 Hz. It then continued as a pure tone starting at 102.63 ± 0.05 Hz down to 102.41 ± 0.04 Hz during 24.19 ± 0.14 s and then suddenly peaked to 108.08 ± 0.06 Hz. Unit 3 lasted 29.25 ± 0.10 s in total, and the entire song was 84.76 ± 0.16 s in duration.Australian pygmy blue whale The Australian pygmy blue whale song is the most complex of the pygmy blue whale songs. It is traditionally described as a 3-unit signal, although multiple variations in the unit order (or syntax) are found42. The song variants change the order and repetition of the unit types. Here, for simplicity, we selected and thus described only the common traditional 3-unit song (Fig. 5c).Unit 1 was 48.83 ± 0.20 s in duration. It had 2 subunits: subunit 1 was a pulsed sound, with a pulse rate (Delta f^{s}_{u1su1}) = 1.21 ± 0.01 Hz at the beginning of the subunit, pulsing accelerated to reach (Delta f^{e}_{u1su1}) = 1.71 ± 0.01 Hz at the end of the unit. Following the ratio “band frequency/pulse rate” criterion, this unit is a non-tonal pulsed sound. However, it is a biphonation sound, as higher energy bands, which do have a harmonic relationship and are spaced by approximately 9 Hz, are obvious on the spectrogram (grey arrows on Fig. 5). The higher fundamental frequency (G_{0}) was at (sim) 9.10 Hz. The resonance frequency of this harmonic sound was the (G_{1_{u1su1}}). It started at 18.20 ± 0.02 Hz and ended at 18.47 ± 0.02  Hz, and was 23.85 ± 0.16 s in duration. Subunit 2 is also a biphonation sound, with a (F_{0}) at 2.80 ± 0.03 Hz at the beginning of the unit ((Delta f^{s}_{u1su2}) in Fig. 5c), decreasing to 1.78 ± 0.01 Hz at the end of the subunit ((Delta f^{e}_{u1su2})), which gives an impression of a decreasing pulse rate when listening to the song. This change in (F_{0}) frequency creates the complicated pattern of intersecting sidebands toward the end of unit 2. The harmonic bands are spaced by approximately 20 Hz (= (G_{0}), precise measurements are given below). Subunit 2 had two variations: subunit 2 was continuous in 42.9% of the sampled songs, but was interrupted by a short gap in 57.1%. In the continuous subunit case (N = 48), the fundamental frequency ((G_{0_{u1su2}})), which is here the band with the most energy, started at 20.22 ± 0.03 Hz and ended at 20.71 ± 0.02 Hz. The subunit lasted 23.26 ± 0.2 s. In the interrupted subunit case (N = 64), the fundamental frequency ((G_{0_{u1su2}})) started at 20.12 ± 0.03 Hz and slightly increased to 20.44 ± 0.02 Hz over 15.27 ± 0.21 s. Then, there was a silence of 3.32 ± 0.08 s followed by the resumption of the subunit at 20.29 ± 0.03 Hz increasing to 20.48 ± 0.17 Hz over 5.71 ± 0.17 s. In this case, the total duration of the subunit (gap included) was 24.31 ± 0.14 s.Unit 2 followed after 7.30 ± 0.09 s. It started as a slightly noisy pulsed sound (possibly deterministic chaos) with a rate (Delta f_{u2}) = 2.77 ± 0.06 Hz during 4.54 ± 0.07 s, then continued as a tonal sound with harmonics. The (F_{0_{u2}}) started at 20.11 ± 0.06 Hz, increased to 22.61 ± 0.02 Hz over 5.14 ± 0.10 s, and then slowly increased to 23.84 ± 0.02 Hz over 23.84 ± 0.02 s. Unit 2 was 23.12 ± 0.12 s in duration.Unit 3 followed after 24.28 ± 0.09 s of silence. It started as a tonal sound with harmonics spaced by 8.93 ± 0.05 Hz. The resonance frequency ((F_{1_{u3}})) started at 7.59 ± 0.02 Hz then increased to 18.26 ± 0.01 Hz over 3.76 ± 0.05 s, with the appearance of sidebands with non-harmonic relationship, spaced by (Delta f_{u3}) = 3.19 ± 0.09 Hz. These non-tonal pulses stopped approximately 3.5 s before the end of the unit, which ends on the harmonic sound, slightly down swept to 18.05 ± 0.02 Hz. These sidebands could be subharmonics, ((F_{0}/3, 2F_{0}/3), etc). Alternatively, they could suggest a biphonation sound. This third unit lasted 18.82 ± 0.12 s in duration, and the whole 3-unit song was 123.54 ± 0.29 s in duration.Figure 5Spectrograms (upper panels) and waveforms (middle panels) of the song of the Madagascan, Sri Lankan and Australian pygmy blue whales, including detailed waveforms to show the internal signal structure. The Madagascan song was recorded off Crozet Island (CTBTO records, site H04S1) in April 2004, the Sri Lankan song was recorded at DGN (CTBTO records, site H08N1) in April 2009 and the Australian song was recorded at Perth Canyon in March 2008 (IMOS records). (Spectrogram parameters: Hamming window, 1024-point FFT length, 90% overlap. Note that the axes differ among plots.) And measurements (mean ± s.e., lower panels) of the acoustic features of the different song types. N is the number of measurements, (u_{i}su_{j}) stands for (unit_{i} subunit_{j}) where i and j are the unit and subunit numbers, (Delta f) designates the frequency difference between the sidebands, f and d are the frequency and duration of the feature indicated in subscript, and when present, the number in brackets refers to the point measured as indicated on the corresponding spectrogram. (F_{x}) or (G_{x}) designate the xth harmonic of a sound, and Cf designates the carrier frequency of a sound.Full size imageOmura’s whale songsAll Omura’s whale songs showed energy between 15 and 55 Hz and peaks of energy around 20 and 40–45 Hz (Fig. 6 lower panels).Ascension Island Omura’s whale Omura’s whale songs recorded in 2005 off Ascension Island started as a tonal sound at 19.84 ± 0.03 Hz. This tone was 3.21 ± 0.08 s in duration but less than 1 s after its beginning, it was overlapped by a noisy pulsed sound, typical of deterministic chaos. The pulse rate was estimated at (Delta f) = 1.44 ± 0.05 Hz. This deterministic chaos lasted for 5.20 ± 0.07 s. Finally, 2.65 ± 0.06 s after the beginning of the song, three tonal components appeared at harmonically independent frequencies, characteristic of triphonation: two tones starting simultaneously, one at 20.88 ± 0.02 Hz and the other at 21.85 ± 0.03 Hz, lasting respectively 4.08 ± 0.23  s and 3.65 ± 0.16 s, and a third tone starting a bit later, 4.48 ± 0.07 s after the beginning of the song, at a frequency of 47.22 ± 0.03 Hz and lasting 3.33 ± 0.09 s. The duration of the total component was 7.64 ± 0.11 s (Fig. 6a).Madagascan Omura’s whale song The following description of the Madagascan Omura’s whale song uses the description provided by Moreira et al.18 and observation from the spectrogram (Fig. 6b). In 2015, Cerchio et al. described the Madagascan Omura’s whale song recorded in 2013–2014 as a single-unit amplitude-modulated low frequency vocalization, with a 15–50 Hz bandwidth15. More recently, Moreira et al. reported a 2-unit song, with the first unit commencing as an amplitude-modulated component with bimodal energy at 20.75 Hz and 40.04 Hz, followed by a harmonic component with a low harmonic at 20.0 Hz and an upper harmonic at 41.0 Hz, as well as an additional tone at (sim) 30 Hz. Unit 1 was characterized as sometimes followed by a tonal unit at 16 Hz18. The ICI was 189.7 s (s.d. 16.47 s, measured from 118 series with (ge) 20 consecutive songs) and ranged from 145.5 to 237.6 s43.Based on the song example recorded in December 2015 in Nosy Be, Madagascar, and provided by S. Cerchio, we observed a 2-unit song (Fig. 6b). The first unit started as chaotic, with no visible sidebands. After (sim) 3 s the signal had a bi- or triphonation event (whilst the deterministic chaos still continues), with first a tone at 40.04 Hz, another tone with a harmonic relationship at 20.02 Hz but starting circa 2.6 s later and a third one at 27.8 Hz starting 4.4 s after the beginning of the first tone, whilst the chaotic sound ends (the chaotic sound lasted circa 9.3 s). The tones of the bi- or triphonic sound all ended at the same time, 11.7 s after the beginning of the song. The second unit seems to be optional15,17,18,44. It followed after 2.8 s silence. It was a tonal sound of 4.9 s in duration with a peak frequency of 16.6 Hz. (Note that the observations here are purely qualitative since only based on 1 song).Diego Garcia Omura’s whale song (DGC) The ‘Diego Garcia Croak’—DGC—recently attributed to the Omura’s whales17 was comprised of one unit (Fig. 6c), although sometimes a second unit was present. The first unit was tonal at the start, with a frequency of 17.91 ± 0.03 Hz, quickly becoming a noisy pulsed sound, characteristic of deterministic chaos, with a pulsed rate of 2.09 ± 0.07 Hz estimated on 41 songs. This chaotic component was 2.76 ± 0.06 s in duration to then became pulsed, although still slightly noisy, with a pulse rate of 2.21 ± 0.005 Hz. This part showed a peak of energy around 19.46 ± 0.08 Hz, and another one around 43.51 ± 0.11 Hz (Fig. 6c, lower panel), and lasted 4.07 ± 0.05 s. Finally, the unit ended as a tonal sound at 17.62 ± 0.04 Hz lasting 5.29 ± 0.12 s. This whole unit had a duration of 10.56 ± 0.14 s. In some occurrences (N = 12), a second tonal unit was present after a silence of 39.89 ± 0.5 s. Unit 2 started at 13.51 ± 0.06 to 13.46 ± 0.04 Hz and lasted for 3.81 ± 0.20 s. When the second unit was present, the entire song was 54.74 ± 0.19 s in duration. Note that in our study, out of the 80 songs measured only 12 had unit 2.Australian Omura’s whale song The Omura’s whale song recorded in 2013 off western Australia had two units (Fig. 6d). Unit 1 was a noisy pulsed sound with a pulse rate of 1.65 ± 0.06 Hz with deterministic chaos, and a duration of 6.28 ± 0.06 s.The peak in energy was at 25.32 ± 0.14 Hz followed by a gap of 2.53 ± 0.04 s, and then a second noisy pulsed unit, with a pulsed rate of 1.80 ± 0.02 Hz estimated on 83 songs. This unit lasted 4.08 ± 0.03 s and had a peak of energy at 25.25 ± 0.18 Hz and another one at 41.20 ± 0.18 Hz (Fig. 6d, lower panel). During the last third of unit 2, the song transitioned to a tonal sound, starting at 25.15 ± 0.02 Hz and swept down to 25.07 ± 0.02 Hz over 3.28 ± 0.04 s, then abruptly decreased to 19.8 ± 0.02 Hz and became tonal for 4.90 ± 0.07 s, forming a z-shape on the spectrogram representation. The whole song was 16.39 ± 0.08 s in duration.Figure 6Spectrograms (a–d), waveforms (e–h), acoustic measurements (mean ± (s.e.)—i–l), and Power Spectral Density (PSD—m–p) of the songs of the Omura’s whales from Ascension Island, Madagascar, Diego Garcia and Australia. The stars on the PSD (m–p) outline the peaks of energy. The Ascension Island song was recorded off Ascension Island (CTBTO records, site H10N1) in November 2005, the Madagascar song was recorded off Madagascar in December 2015 and provided by S. Cerchio, the Diego Garcia DGC song was recorded at DGN (CTBTO records, site H08N1) in October 2003 and the Australian song was recorded at Kimberley site in March 2013 (IMOS records). For the panels (a–d) and (i–j): N is the number of measurements, (u_{i}su_{j}) stands for (unit_{i} subunit_{j}) where i and j are the unit and subunit numbers, (Delta f) designates the frequency difference between the sidebands, f and d are the frequency and duration of the feature indicated in subscript, and when present, the number in brackets refers to the point measured as indicated on the corresponding spectrogram. (Spectrogram parameters: Hamming window, 1024-point FFT length, 90% overlap. Note that the axes differ among plots).Full size imageDeterministic chaosWe classified deterministic chaos as: ‘slight’, where sidebands were easily distinguished but the sound was noisy; ‘moderate’, where the sidebands were visible but difficult to measure; and ‘strong’, where the sound had no discernible structure. Where deterministic chaos was present, we identified its persistence, defined as the proportion of deterministic chaos over the duration of a song31.It was difficult to characterize the presence of deterministic chaos where the song (sub)unit was short and the pulse rate was low, as it is difficult to ascertain if the noisy structure (i.e., lack of structure) was part of the whale’s song (i.e., deterministic chaos) or whether it was due to an artefact, such as a sound propagation issue. This was the condition for the subunit 2 of unit 1 of the Chagos song. If this subunit had indeed a chaotic structure, this chaos was slight, and represented 4.5% of the entire duration of the song (Fig. 7a).Pygmy blue whale songs had only slight deterministic chaos, and of the entire song, it represented: 11.7% of the duration of the Madagascan song; 3.7% of the Australian song; and it was not present in the Sri Lankan pygmy blue whale song (Fig. 7a). In the Madagascan pygmy blue whale songs, slight deterministic chaos was in subunits 1 of both units 1 and 2, and in the Australian pygmy blue whale songs, deterministic chaos was present in subunit 1 of unit 2.In contrast, deterministic chaos was a significant proportion of all Omura’s whale songs (Figs. 6a–d and 7a). For the song of the Ascension Island Omura’s whale, moderate deterministic chaos was present across 68% of the duration of their song. For the Australian Omura’s whales, deterministic chaos was present across 63.2% of their song, it was moderate-to-strong in the first unit and slight in the second unit. The Madagascan Omura’s whales had strong deterministic chaos across 72% of their song, which excludes the tonal unit as the tonal part was not always present. The Diego Garcia DGC Omura’s whale song had a total chaos persistence of 65.2% (Fig. 7a), with a moderate deterministic chaos present in the first 2.7 s of the song, which represents 26.3% of the song duration (Fig. 7a medium grey section). The song then evolved to a more clearly pulsed sound, with a slightly noisy structure, classified as slight deterministic chaos. Here again, it was difficult to ascertain whether this lack of structure was a characteristic of the song or an artefact of the propagation. Yet, the slight lack of structure was consistently observed across the sampled songs.Inter-call-intervalsWhilst the Madagascan pygmy blue whale had a shorter ICI, all the other acoustic groups studied here had a similar ICI duration (Fig. 7b). Thus, ICI is not a key parameter in the distinction among species and cannot be used to determine whether Chagos-whales are a blue or an Omura’s whale.Figure 7(a) Proportion of deterministic chaos (i.e., chaos persistence) in the Chagos song compared with the three Indian Ocean pygmy blue whale song types (Madagascan, Sri Lankan and Australian) and the four Omura’s whale song types (Madagascan, Diego-Garcia DGC, Australian and Ascension Island). Chaos persistence is defined as the proportion of deterministic chaos across the entire song duration (given as a percentage). Shades of grey indicate the strength of the chaos: slight (light grey), moderate (medium grey) and strong (dark grey)). (b) Boxplot representation of the Inter-call Intervals (ICI expressed in s) for the different song-types measured in this study. On each box, the central mark is the median, the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points considered to be not outliers, and the outliers are plotted individually.Full size imageGeographic distributionChagos song was detected at 5 of our 6 recording sites at disparate locations across the Indian Ocean, from: the northern Indian Ocean, off Sri Lanka; on both sides of the central Indian Ocean, off the Chagos Archipelago; and in the far eastern Indian Ocean, off northern Western Australia (Fig. 2). The Chagos song was recorded off Sri Lanka (i.e., Trincomalee) in April. Blue whales were observed at the time the recordings were made, and the songs of the Sri Lankan pygmy blue whale were also recorded at the time. The acoustic recording had become degraded as they were made nearly forty years before, on 19 April 1984, and only six distinct Chagos songs were found. Unfortunately, these recordings were of poor SNR which prevented detailed acoustic measurement. The songs, however, had the distinct structure of the Chagos song (Fig. 3) and an ICI of (simeq) 200 s (range 200 to 209 s), consistent with the ICI rate measured for the Chagos song off the Chagos Archipelago (Fig. 7b). Further south in the northern Indian Ocean, 6,984 Chagos songs were detected in 2013 (from January to early December) at our recording site RAMA, but no songs were detected at this site in 2012, although recording had been made over a shorter period, from May to December, in that year. In the central Indian Ocean, a total of 486,316 Chagos songs were detected from January 2002 to March 2014 at DGN, and 737,089 Chagos songs from January 2002 to August 2018 at DGS. In the far eastern Indian Ocean, off Kimberley, northern Western Australia, low SNR Chagos songs were manually detected from January to May, in 41 out of the 331 recording days in 2012–2013. In the south-central Indian Ocean, at our recording site RTJ, no Chagos songs were detected in 2018.Figure 8 shows the average number of Chagos songs detected per day for each year of data at the sites located on: (a) the western (DGN); and (b) eastern (DGS) sides of the Chagos Archipelago; as well as (c) further north-east, at RAMA site. The number of songs varied over the years; fewer songs were recorded at both DGN and DGS sites in 2008. In comparison with the Chagos Archipelago sites, the number of songs detected at north-eastern RAMA was low in 2013, with an average of only 20 songs/day.Figure 8Average number of Chagos songs per day detected in each year of data on the (a) western (DGN) and (b) eastern (DGS) sides of the Chagos Archipelago, and (c) further north-east, at RAMA site.Full size imageSeasonalityFigure 9b shows the average seasonality of Chagos song occurrence on both sides of the Chagos Archipelago. On the western side of the central Indian Ocean (DGN site), Chagos songs were heard predominantly from September to January, with detections peaking in December and January. On the eastern side of the central Indian Ocean (DGS site), songs were detected from June to November, with detection peaks in August to October, depending on the year. In 2013, at the RAMA site (further north-east of the Chagos Archipelago), Chagos songs were detected from January to June (with peaks in May), and in November (Fig. 9a). Off Kimberley, in the north of Western Australia, low SNR Chagos songs were found from the 22 January 2012 to the 20 May 2012, with a peak in March (Fig. 9c).Figure 9(a) Seasonality of Chagos songs at RAMA in 2013, presented as a percentage of songs per month (i.e. monthly number of songs divided by total number of songs detected in the year); (b) Seasonality of Chagos song averaged over the years (±SE) on the western (DGN—gray) and eastern (DGS—orange) sides of the Chagos Archipelago. This average seasonality is calculated as such: the monthly number of songs is divided by the total number of songs detected in the corresponding year, and averaged over the years. Note that due to the low number of recording days at DGN in 2007 and 2014, and in 2007 at DGS, these years were removed from the averaging (DGN: 11 years and DGS: 15.5 years); (c) Hourly presence of Chagos songs in Kimberley (Western Australia) in 2012–2013. Note that the metric and thus the graphic representation used for this site is different from that for RAMA and DGN/DGS: in the Kimberley data set, Chagos songs were logged upon visual inspection of the spectrograms, and a metric of hourly presence/absence of the song per day was used (see the Methods section for details).Full size imageWe found strong evidence at both Chagos Archipelago sites (DGN and DGS) that the number of Chagos songs changes not only across months (Table 1; (p=0.02417), Table 2; (p < 0.001)) and years (Table 2; (p < 0.001), Table 2; (p < 0.001)), but also that there is an interaction between months and years (Table 1; (p < 0.001), Table 2; (p < 0.001); Fig. 10). This provides evidence to suggest that there is variation in the pattern of whale songs across years at both sites. Although Chagos songs were detected throughout the year, there were more songs detected at restricted times (Fig. 10). The timing of peaks in song detection was different between the sites. At DGN most songs were detected in 2 to 3 months, whereas at DGS songs were detected over a longer period, from 2 to 6 months. At DGN, where the Chagos song distribution in most years shows clear peaks towards December and January, in a few years, peaks were outside this time (e.g., 2005 in March and September, 2006 in September and 2008 in July and August; Fig. 10). Conversely, in DGS most songs were observed between June and November, although there were inter-annual differences (Fig. 10).Table 1 Assessing the likelihood of an effect on the number of Chagos songs per day at site DGN (n = 3917 days).Full size tableTable 2 Assessing the likelihood of an effect on the number of Chagos songs per day at site DGS (n = 5557 days).Full size tableFigure 10Number of Chagos songs per month for each year at DGN and DGS. Note that the scale of the y-axis differs among years to highlight the seasonal patterns. Months without data are indicated by ‘No Data’, and months with more than 50% of missing days are indicated by a black dot.Full size image More

California has a vast area and spans ten latitudes, and its internal geographical conditions and climate conditions vary widely13. Therefore, the California wildfires in history differed greatly in their frequency, size, intensity and extent of damage8. As the California wildfires are growing fiercer, they have become a hot topic worldwide. However, there is still a long way to go before the general conclusions from the wildfire literature can be applied in practice. For example, how the analyses of which types of wildfires are increasing the fastest can be used to guide the amendment of wildfire management policies? and how to guide fire fighting methods based on the results of the wildfire dominant factor model? To provide some practical reference for wildfire management work, we grouped the wildfires according to size (large fires, small fires) and ignition cause (natural fires and human-caused fires), and discussed their distribution characteristics separately using the administrative units from CAL FIRE and the weather division of California from the National Oceanic and Atmospheric Administration (NOAA) as the base map. While focusing on wildfires in the past two decades, the distribution of wildfires from 1920 to 1999 was used as prior information for comparison.Wildfire size distributionThe burned area of wildfires is an important indicator of their destructive power. Several studies have shown that 1% of large and extreme wildfires are responsible for 90% of the total damage caused by wildfires14,15. Besides, the Probability density distribution of wildfire burned area, that is the wildfire size, has an obvious heavy tail feature. Research from Strauss et al.16 and Holmes et al.17 indicate that the wildfire size distribution fits the Pareto distribution well. Based on their conclusions, five common heavy-tailed distributions were selected (which are Gamma, Lognormal, Pareto, Truncated Pareto and Weibull distribution) to fit the wildfire size distribution throughout California within the eighty years before year 2000 and twenty years after year 2000s, seeking the best description of the California wildfire size distribution. The estimated parameters and the goodness of fit test results are shown in Tables 1 and 2. The empirical wildfire size distribution and the fitting curve are shown in Figure. 1. Fig. 1 shows that the wildfire size distribution did not change much from the last century to the present. Also, all these fitting curves can capture the main feature of the empirical distribution. Table 1 lists the estimated shape and scale parameters for each distribution. It can be found that the shape parameter of current wildfire size distribution ((alpha)) decrease compared to the historical wildfires. The value of shape determines the thickness of the tail. A smaller shape value means a thicker tail. In the context of wildfires, it means the probability density of large wildfires increase. Table 2 shows the goodness of fit for each distribution by Akaike Information Criterion (AIC), Kolmogorov-Smirnov (K-S) and Cramer-VonMises (CvM) test score. For all the tests, the smaller the value of the test score, the better the fit. Among these five fits, the lognormal distribution is the best for wildfire size description in 1920–1999, following by the Pareto distribution; while the best fitting distribution in 2000–2019 changes to the truncated Pareto, the second-best fitting result is still from the Pareto distribution. Therefore, Pareto is appropriate to summarize the general feature of wildfire size distribution in California.
Table 1 Heavy-tailed distribution fitting results of wildfire size distribution.Full size tableTable 2 Goodness-of-fit test results of Akaike Information Criterion (AIC), Kolmogorov–Smirnov (K–S) test, and Cramer-Von Mises (CvM) test for heavy-tailed distribution fitting.Full size tableTo further explore the variation of wildfire size distribution within the entire state of California, the probability density of the logarithm of wildfire size was plotted for 1920 to 1999 and 2000 to 2019. As shown in Fig. 2, wildfires in 1920–1999 were mostly about 100–1000 acre (0.40–4.05 km(^2)) in size; while during 2000–2019, the number of small fires increased significantly,Figure 1The empirical histogram of wildfire size and the typical heavy tailed distribution fitting curves for wildfires in (a) 1920–1999 and (b) 2000–2019. The wildfire sizes are in acres (1000 acre = 4.05 km(^2)). The curve with different colors represent different types of distribution, the black, yellow, red, green, and blue curves represent the fitting result of Gamma, Log-normal, Pareto, Truncated Pareto, and Weibull distribution, separately. The tail of the distribution was truncated from the burned area of 2000 acres to show the fitting difference between different distributions.Full size image the majority of wildfire sizes were in the range of 10–100 acres (0.04–0.40 km(^2)). Wildfires were also divided into natural wildfires and human-caused wildfires based on their ignition causes. The red, green and blue dashed lines in the figure delineate the fitting results of Gamma, Lognormal and Weibull distribution separately, which capture the distribution characteristics for each type of the wildfires. The fitting parameters and the goodness of test results were attached in the supplementary information (Table S1). Figure 2b,e show that although the overall shape of the distribution of natural wildfires in 1920–1999 and 2000–2019 are similar, the proportion of extreme wildfires larger than 10,000 acres (40.47 km(^2)) has increased significantly in the last two decades. From Fig. 2c,f, it can be found that the shape of the fire size distribution of human-caused wildfires differs greatly, which is the result of the rapid increase of the proportion of small fires. Although human activity directly or indirectly ignited 44(%) of wildfires in the United States18 and 39(%) of wildfires in California (as shown in the statistical summary in Table 3), they are generally easily contained in the initial attack19. The rapidly growing population in California has led to increased human activities and community coverage, which has increased the incidence of human-caused wildfires20. However, the expansion of human land has reduced the continuity, which is essential for the spread of wildfires21. Also, the improvement of wildfire monitoring and fire fighting ability has made most of the small human-caused wildfires able to be extinguished during the first 24 h after discovery19. Together, these reasons lead to the rapid increase in the frequency of small human-caused fires in the past two decades.Figure 2Logarithm of California wildfire size empirical distribution in 1920–1999 and 2000–2019. The Gamma, Lognormal and Weibull distribution fitting results are indicated by the red, green and blue dash lines. The wildfire sizes are in acres (1000 acre = 4.05 km(^2)). (a–c) are the historical wildfires from 1920 to 1999, (d–f) are the wildfires from 2000 to 2019; (a,d) are the distribution of all wildfires, (b,e) are the distribution of natural wildfires, (c,f) are the distribution of human-caused wildfires.Full size imageTable 3 Statistical summary of wildfire ignition causes in CA from 2000 to 2019.Full size tableLarge and small fires are not only very different in the probability density distribution characteristics but also in prevention measures, response methods, and resources needed to be invested in fire fighting22,23. In order to discuss the spatiotemporal distribution of large and small wildfires, it is critical to determine the threshold of large wildfires. Therefore, the mean excess plot shown in Fig. 3 was used to determine the threshold of the large fire. The linear part’s starting point is the threshold of the extreme value in the original distribution17,24. As shown in Fig. 3, 500 acres (2.02 km(^2)) would be appropriate to separate the large fires and small fires for the entire California. Also, as shown in Fig. 1, 500 acres is an appropriate starting point of the heavy tail. Based on the historical record from CAL FIRE, the frequency of large wildfires accounted for 19.68 (%) of the total (1247 out of 6336 wildfires), while the burned area of large wildfires accounted for 97.04 (%) of the total burned area (13,089.68 out of 13,488.19 thousand acres, that is 52,972.05 out of 54,584.77 km(^2)) in the past two decades. According to the size class of fire defined by national wildfire coordinating group (NWCG), the large fire in this study refers to the wildfires of or larger than class E.Figure 3Mean excess plot for wildfires burned areas.Full size imageTemporal variation of wildfires in CA from 1920–1999 and 2000 to 2019Based on the wildfire history records provided by the CAL FIRE Fire Perimeter database, the frequency and burned area of wildfires in CA from 1920 to 2019 were extracted, and separated into two time periods: 1920–1999 and 2000–2019. California has seen an average of 317 wildfires a year over the past 20 years, which were included in the Fire Perimeter database, burning an average of 674,410 acres (2,729.24 km(^2)). Figure 4 shows the changes in the annual wildfire frequency (a–e) and burned area (f–j) over time. The red lines represent the segmented linear regression trend in 1920–1999 and 2000–2019, separately. The grey areas depicted the 95(%) confidence interval. Comparing the slope of the fitting line, it is apparent that in most cases, the frequency and burned area growth of wildfires in the past two decades are much higher than that during the 80 years in history, if the breakpoint is fixed to the year 2000. Also, the 95(%) confidence intervals of the regression lines over the past two decades are generally larger than that between 1920 and 1999. Although the sample size in these two time periods is different, it can be seen from the spread of data points that the uncertainty of wildfire frequency and burned area have increased significantly in the past two decades. From the view of fire frequency, the rapid increase in the number of small fires brings greater uncertainty than that of large fires, and the uncertainty of natural fires is higher than that of human-caused fires. In terms of the burned area, the uncertainty comes mainly from large wildfires and natural wildfires. When it comes to the increase rate, Fig. 4b,c,g,h show that in the large and small wildfire group, the accelerated increase of wildfire frequency was mainly contributed by the small fires, while the accelerated increase of burned area was from the large fires. The frequency of large wildfires and the burned area of small wildfires in the recent 20 years even have the trend of decrease. This trend suggests that it would be efficient for the fire management department to pay more attention to the regions with the potential risk of extreme fires and prevent small fires from burning continuously and becoming large fires. Figure 4d,e,i,j display the trend for the natural and human-caused wildfires. The increase of the human-caused wildfire frequency is much faster than that of the natural wildfires in both time periods. However, the increases in the burned area due to the increasing frequency of wildfires with different causes are similar. It shows that the human-caused small wildfires have the strongest growth trend in the recent twenty years. In the view of wildfire management, while human activities increase the likelihood of wildfires ignition, large natural fires are more threatening in terms of size and destruction.Figure 4Temporal distribution of wildfire frequency and burned area from 1920 to 2019. The red line indicates the segmented linear regression results for 1920–1999 and 2000–2019. The gray areas indicate the 95(%) confidence interval. (R^2) represents the coefficient of determination and p represents the p-value. (a–e) are the temporal distribution of wildfire frequency, (f–j) are the temporal distribution of the burned area of wildfires; (a,f) are the distribution for all wildfires; (b,g) are plots of large fires, which have the burned area larger than 500 acres (2.02 km(^2)), while (c,h) are plots of small fires, which have the burned area in the range of 10 acres (0.04 km(^2)) to 500 acres (2.02 km(^2)); (d,i,e,j) divided wildfires into natural fires and human-caused fires. The small plot in (h) zooms in to the burned area of 0–50 thousand acres.Full size imageCalifornia’s Mediterranean climate is characterized by hot and dry summers, which leads to a high wildfire ignition risk25,26. Also, the hot and dry Santa Ana wind events have accelerated the spread of wildfires each fall27. The precipitation in California was concentrated in the winter, and the temperature was moderate28, allowing wildland vegetation to grow fast and storing fuel for next year. However, the significant climate change after the year 2000 has affected the seasonal distribution of wildfires.Figure 5 compiles box plots of the seasonal variation of wildfire frequency and burned area distribution in 1920–1999 and 2000–2019, which were divided into different groups by size and ignition cause as well. The boxes and points in the plots represent the wildfire frequency or total burned area in this month each year. In general, the peak season for wildfires was late summer and early autumn. In terms of the frequency, from 1920 to 1999, the wildfire season started in June, and the most frequent occurrence was observed in August. In most years, the number of wildfires in July and August were similar, followed by June and September. However, from 2000 to 2019, the frequency of wildfires in July increased significantly and became much more considerable than in other months. Meanwhile, the start of the wildfire season has also advanced to May, and the duration has extended. From May to September, the overall fire frequency of all wildfires, large wildfires, and small wildfires increased each month. The number of natural fires also increased between June to September. The frequency of human-caused wildfires, on the other hand, increased each month. Similar to the previous discussions, the increase of wildfire frequency in July in the past two decades mainly came from small fires and human-caused wildfires. It is worth noting that there has been a major increase in the natural wildfires in July in the past two decades. In terms of the burned area, the month with the largest total burned area of wildfires in 2000–2019 has been advanced to July, compared to August in 1920–1999. Natural wildfires and human-caused wildfires contributed similarly to the burned area growth. There is no noticeable change in the total burned area in months other than the wildfire season.Figure 5Seasonal variation of wildfire frequency and burned area from 1920 to 2019. The threshold of large and small wildfires is 500 acre (2.02 km(^2)). (a–j show the seasonal variation of fire frequency, (k–t) show the seasonal variation of burned area; (a,b,k,l) are plots for all CA wildfires, (c–f) and (m–p) divided fires into large and small fire size group, (g–j) and (q–t) divided fires into natural and human-caused wildfire groups. The small plots in (o) and (p) zoom in to the burned area of 0–10 thousand acres.Full size imageSpatial distribution of wildfires in CA from 2000 to 2019CAL FIRE has 21 operational units throughout the state that are designated to address fire suppression over a certain geographic area and six ‘Contract Counties’ (Kern, Los Angeles, Marin, Orange, Santa Barbara and Ventura) for fire protection services. Due to the complex environmental and terrain conditions in California, the risk of wildfires varies significantly from region to region, and the causes of extreme wildfires are also completely different. In order to provide fire managers with more effective fire suppression measures, this study used kernel density estimation (KDE) to analyze hot spot regions of all the wildfires, natural fires and human-caused fires from 2000 to 2019, the KDE for wildfires in 1920–1999 were also added for comparison. The resolution of KDE analyses was 500 m. The results are shown in Figs. 6 and 7. Figure 6 treated all the fires equally, and shows the spatial density of wildfire numbers; while Fig. 7 weighted the wildfires with their burned area, and represents the burned area-weighted spatial density of wildfire occurrence.Comparing the spatial density distribution of all wildfires in different time periods in this study, as shown in Fig. 6a,d, it is evident that the coverage of wildfire occurrence has increased significantly. From 1920 to 1999, the only hot spot with a very high wildfire density was Los Angeles County (LAC). In the past two decades, not only did the hot spot of LAC expand to Ventura county (VNC) but also the wildfire density in the southwest corner of Riverside Unit (RRU) and San Diego Unit (MVU) on the south coast and the southwest corner of San Bernardino Unit (BDU) have grown to a very high level. In the eastern part of the San Joaquin Drainage under the central California climate division, namely the Sierra Nevada Mountains (identified in Fig. 10), wildfire density has increased from very low to very high. Among them, Nevada-Yuba-Placer Unit (NEU) and Tuolumne-Calaveras Unit (TCU) are the newly emerged high-density wildfire regions. Moreover, the spatial density distributions were grouped by causes, and Fig. 6b,e represent the natural wildfires, and c,f represent the human-caused wildfires. It can be found that while the high-density areas of natural wildfires have not shifted in both time periods, the density has increased. In contrast, the density of human-caused wildfires has increased notably in western and central California in the past two decades. Before the year 2000, there were almost no human-caused wildfires along the west coastline, but almost every county along the west-coast is characterized by an increase of human-caused wildfires in the past two decades. San Benito-Monterey Unit (BEU) and San Luis Obispo Unit (SLU) even became the new hot spots. Meanwhile, the coverage area of the original human-caused wildfire hot spots on the south coast has been further expanded. From 1920 to 1999, the density of human-caused wildfires in the Sierra Nevada Mountain was very low in central California. Still, in the past two decades, it has become a new wildfire ignition hot spot. The counties in northern California, such as Siskiyou Unit (SKU), Shasta-Trinity Unit (SHU), Tehama-Glenn Unit (TGU), etc., have been almost no human-caused wildfires from 1920 to 1999, but widespread human-caused wildfires have emerged in the past two decades.After inducing the wildfire burned area into the KDE calculation, the spatial density distribution has changed significantly. In general, as shown in 7a,d, the regions where large wildfires are concentrated are SKU and Sonoma-Lake-Napa Unit (LNU) in Northern California and MVU in the South Coast. Although the number of wildfires in the central Sierra Nevada Mountains has increased significantly, the total burned area did not significantly change. Thereafter, the wildfires with different causes were separated, and it can be found from 7b,e that natural wildfires with large burned areas were concentrated in northern California. In the past two decades, the region with a very high-density of wildfire occurrence in the northernmost SKU has expanded significantly, and a new hot spot of wildfires has also appeared in Lassen-Modoc Unit (LMU). However, the high-density wildfire area between Tuolumne-Calaveras Unit (TCU) and Madera-Mariposa-Merced Unit (MMU) did not arise in the past two decades. In the distribution of human-caused wildfires, as shown in 7c,f, the density of wildfires in MVU in the southernmost part of California has surpassed that of historical hot spots, VNC and LAC. Meanwhile, the density of wildfires at the junction of TCU and MMU in the central region has also increased.Comparing 6 and 7, it is obvious that the spatial distribution of wildfire density and burned area-weighted wildfire density are not entirely consistent. CAL FIRE Units along the South Coast, which are in the climate division of South Coast Drainage, are prominent in both densities, and are mainly composed of human-caused wildfires. The SKU and LMU units in the northernmost part of North Coast Drainage are the areas where natural wildfires were concentrated, and the distribution of SKU wildfires is relatively wider. The Units adjacent to the Sierra Nevada Mountains in central California, which are the units in the northeast of San Joaquin Drainage, show a low wildfire density when the burned area was added to the calculation, even though the number of wildfires has increased rapidly in the past two decades. This distribution is related to the vegetation cover and land use in California. In northern California, the evergreen and deciduous forests are the dominant vegetation, the forests are dense and less developed by human, and the population density is relatively low28,29. Wildfires are difficult to be detected early-on in these remote areas, and there is enough fuel to keep them burning and spreading. On the other hand, shrubs are the dominant vegetation in southern California. Also, most of the southern CA areas have been developed and associated with a higher level of human activity, leading to wildfires in southern California has a greater social and economic impact on human lives and society30.Figure 6Kernel density distribution of wildfire occurrence in CA during 1920–1999 (a–c), and 2000–2019 (d–f). (a–f) are wildfire density distribution maps for all wildfires, natural wildfires and human-caused wildfires in CA, separately.Full size imageFigure 7Kernel density distribution of burned area weighted wildfire occurrence in CA during 1920–1999 (a–c), and 2000–2019 (d–f). (a–f) are wildfire density distribution maps for all wildfires, natural wildfires and human-caused wildfires in CA, separately.Full size imageFrom the discussion above, it can be found that while the frequency and spatial density distribution of human-caused wildfires have changed significantly in the past two decades, the changes in burned area were relatively small because of the high proportion of small wildfires. Also, unlike natural fires, human-caused fires can be prevented or controlled in the early stage by taking effective measures19. Therefore, the human-caused wildfires were further classified to generate a more detailed spatial density distribution map. The anthropogenic causes were subdivided by CAL FIRE into 15 types. The spatial distribution of wildfires with different causes are shown in the supplementary figures (Supplementary Fig. 1). In this study, human-caused wildfires were classified into three categories: transportation (railroad, vehicle, aircraft), human activity (equipment use, smoking, campfire, debris, arson, playing with fire, firefighter training, non-firefighter training, escaped prescribed fire, illegal alien campfire) and construction (powerline, structure). As shown in Fig. 8, hot spots for all three broad types of wildfires include areas along the Sierra Nevada Range and along the southern coast. However they differ in the density level and coverage. Among them, the number and coverage of wildfires caused by human subjective behavior are larger than those caused by traffic and construction. Besides, the wildfires caused by human activities also led to the emergence of a unique hot spot in the northernmost edge of CA, which is the SKU county. Therefore, for the wildfire management purpose, it would be proactive to provide wildfire education to residents in regions with high wildfire risk, update the wildfire risk map in time, and issue early warnings of wildfire risk to the public during the fire season, to increase the public’s awareness of wildfire prevention.Figure 8KDE Analysis of human-caused wildfires in CA from 2000 to 2019. (a) Transportation (railroad, vehicle, aircraft); (b) Human Activity (equipment use, smoking, campfire, debris, arson, playing with fire, firefighter training, non-firefighter training, escaped prescribed fire, illegal alien campfire); (c) Human Construction (power line, structure).Full size imageMultivariate analysis of California wildfiresThe occurrence and spread of wildfires are related to human activities and environmental variables. In order to formulate effective suppression and control policies for wildfire management, it is essential to understand the relationship between the spatial distribution of wildfires and various variables. From the KDE analysis, the spatial distributions of the wildfire density calculated with and without burned area were obtained, which also shows the areas with high wildfire risk from 2000 to 2019. According to the research from Faivre et al.7, 12 variables that have potential correlations with wildfires, involving human-related variables, geographic conditions, fuel, and climate variables were selected to conduct the subsequent analyses.Table 4 calculated the spatial correlation between the burned area-weighted wildfire density and potential anthropogenic and environmental variables within the wildfire perimeters, as well as the interrelation between each variable. It can be derived from the first column that among the human-related variables, except for the distance to the road, other variables are positively correlated with the wildfire occurrence density. It means that in areas where wildfires have occurred in the last two decades, the farther away from the power line, the higher the wildfire density; the closer to the road, the higher the wildfire density; and the greater the density of houses and population, the higher the density of wildfires. Among environmental variables such as topography, vegetation cover, and climate, only elevation is negatively correlated with wildfire density. That is, the higher the elevation, the lower the wildfire density. From the correlations among various variables, it can be found that there is a strong correlation between the distance from the wildfire perimeter to the road and power line, population, and house density, as well as elevation and two climate variables. For further analyses, one variable would be removed between the two variables whose correlation is greater than 0.5. Therefore, the distance to power line, population density and elevation were removed in the multivariate analysis.Table 4 Spatial Correlation Analysis between 12 selected variables wildfire occurrence density: distance to power line (DP), distance to road (DR), housing density (DH), population density (DP), elevation, aspect, slope, tree, shrub, grass, maximum temperature (Tmax), maximum vapor pressure deficit (VPDmax).Full size tableThe principal component analysis (PCA) was implemented on the remaining variables and the two types of wildfire spatial densities obtained from KDE, to classify the variables and evaluate their relationships. The eigenvalue matrix was attached in the supplement information (Supplementary Table S3.). Both PCA results require five principal components to explain at least 80(%) of the data variance. The interrelations of the variables and the fire occurrence density decomposed by PC1 and PC2 are shown in Fig. 9. There is a strong and similar interrelationship between the two types of fire densities and the driver variables. The length and orientation of the variables indicate that the wildfire densities have the strongest correlation with the grass cover and the other two variables of vegetation cover (shrub and tree), namely fuel cover in general. Meanwhile, the correlation between the climate variables and the wildfire densities is also significant, especially for the maximum vapor pressure deficit (VPDmax). Besides, the human-related variables are moderately correlated with the wildfire densities, while topographic variables are almost orthogonal with the wildfire densities, which means their correlations are weak.Figure 9PCA loading plots with (a) fire occurrence density, (b) burned area weighted fire occurrence density. The variables include distance to road (DR), housing density (DH), aspect, slope, tree, shrub, grass, maximum temperature (Tmax), maximum vapor pressure deficit (VPDmax), wildfire density (FOD) and burned area weighted wildfire density ((FOD_A)).Full size imageBased on the analyses above, the Logistic Regression (LR) was implemented on the selected nine variables to further determine their relationship with wildfire occurrence. The coefficient, standard error and the significance level for each variable were shown in Table 5. The positive and negative sign of the coefficient represents the positive or negative correlation with the wildfire occurrence, and the p-value indicates whether the correlation is significant. The results reveal that the climate variables are the most critical in whether the wildfires can be ignited or not, followed by the variables of distance to road, and the cover of grass. The sign of the coefficient of the human-related variables is negative, which means that in general, wildfires ignited far from the human communities. Similarly, the areas where trees are dominant vegetation cover have fewer wildfire ignitions. Overall, logistic regression results show that the areas with high temperature, high VPD, grass as the dominant vegetation cover, and away from human communities have a higher risk of wildfire ignition.Table 5 Logistic regression results of uncorrelated explanatory variables for California wildfires occurrence (2000–2019).Full size table More