Ecology

1.Benelli, G. Green synthesized nanoparticles in the fight against mosquito-borne diseases and cancer—a brief review. Enzyme Microbial Technol 95, 58–68 (2016).CAS 
Article 

Google Scholar 
2.Dash, A. P., Valecha, N. & Anvikar, A. R. Malaria in India: challenges and opportunities. J. Biosci 33(4), 583–928 (2008).CAS 
PubMed 
Article 

Google Scholar 
3.World Malaria Report: Geneva: World Health Organization. Accessed 18th July 2017.4.Olotu, A. et al. Seven-year efficacy of RTS, S/AS01 malaria vaccine among young African children. N. Engl. J. Med 374, 2519–2529 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Solomona, S., Plattnerb, G. K., Knuttic, R. & Friedlingsteind, P. Irreversible climate change due to carbon dioxide emissions. Proc. Natl. Acad. Sci. U.S.A. 106, 1704–1709 (2009).ADS 
Article 

Google Scholar 
6.Shaalan, E. A. S., Canyonb, D., Younesc, M. W. F., Abdel-Wahaba, H. & Mansoura, A. H. A review of botanical phytochemicals with mosquitocidal potential. Environ. Int. 3, 1149–1166 (2005).Article 
CAS 

Google Scholar 
7.Sundukov, Y. N. First record of the ground beetle Trechoblemus postilenatus (Coleoptera, Carabidae) in Primorskii krai. Far East Entomol. 165, 16 (2006).
Google Scholar 
8.Soni, N. & Prakash, S. Green nanoparticles for mosquito control. Sci. World J. 214, 1–6 (2014).Article 

Google Scholar 
9.Abinaya, M. et al. Structural characterization of Bacillus licheniformis Dahb1 exopolysaccharide antimicrobial potential and larvicidal activity on malaria and Zika virus mosquito vectors. Environ. Sci. Pollut. Res 25, 5 (2018).Article 
CAS 

Google Scholar 
10.Shawkey, A. M., Rabeh, M. A., Abdulall, A. K. & Abdellatif, A. O. Green nanotechnology: anticancer activity of silver nanoparticles using Citrullus colocynthis aqueous extracts. Adv. Life Sci. Technol. 13, 60–70 (2013).
Google Scholar 
11.Thomas, S., Ravishankaran, S. & Johnson Amala Justin, N. A. Resting and feeding preferences of Anopheles stephensi in an urban setting, perennial for malaria. Malar. J. 16(11), 1–7 (2017).
Google Scholar 
12.Murugan, K. et al. Sargassum wightii-synthesized ZnO nanoparticles reduce the fitness and reproduction of the malaria vector Anopheles stephensi and cotton bollworm Helicoverpa armigera. Physiol. Mol. Plant Pathol. 101, 202–213 (2018).CAS 
Article 

Google Scholar 
13.Kalimuthu, K., Panneerselvam, C., Murugan, K. & Hwang, J. S. Green synthesis of silver nanoparticles using Cadaba indica Lam leaf extract and its larvicidal and pupicidal activity against Anopheles stephensi and Culex quinquefasciatus. J. Entomol. Acarol. Res. 45(2), e11 (2013).Article 

Google Scholar 
14.Patra, A., Raja, A. S. M. & Shah, N. Current developments in (Malaria) mosquito protective methods: a review paper. Int. J. Mosquito Res. 6(1), 38–45 (2019).
Google Scholar 
15.Wahab, R., Ahmad, J. & Ahmad, N. Application of multi-dimensional (0D, 1D, 2D) nanostructures for the cytological evaluation of cancer cells and their bacterial response. Colloids Surf. A Physicochem. Eng. Asp. 583, 123953 (2019).CAS 
Article 

Google Scholar 
16.Bhadra, J., Alkareem, A. & Al-Thani, N. A review of advances in the preparation and application of polyaniline based thermoset blends and composites. J. Polym. Res. 27(5), 1–20 (2020).Article 
CAS 

Google Scholar 
17.Jaganathana, A. et al. (+16), Earthworm-mediated synthesis of silver nanoparticles: a potent toolagainst hepatocellular carcinoma, Plasmodium falciparum parasites and malaria mosquitoes. Parasitol. Int. 65(2016), 276–284 (2016).Article 
CAS 

Google Scholar 
18.Abdelkhalek, A. & Al-Askar, A. A. Green synthesized ZnO nanoparticles mediated by Mentha spicata extract induce plant systemic resistance against Tobacco mosaic virus. Appl. Sci. 10, 15 (2020).Article 
CAS 

Google Scholar 
19.Ishwarya, R. et al. Facile green synthesis of zinc oxide nanoparticles using Ulva lactuca seaweed extract and evaluation of their photocatalytic, antibiofilm and insecticidal activity. J. Photochem. Photobiol. 2018(178), 249–258 (2018).Article 
CAS 

Google Scholar 
20.Murugan, K. et al. Nano-insecticides for the control of human and crop pests. In Short Views on Insect Genomics and Proteomics. Entomology in Focus (eds Raman, C. et al.) 229–251 (Springer, 2016).
Google Scholar 
21.Bauer, A. W., Kirby, W. M., Sherris, J. C. & Turck, M. Antibiotic susceptibility testing by a standardized single disk method. Am. J. Clin. Pathol. 45(4), 493–496 (1966).CAS 
PubMed 
Article 

Google Scholar 
22.Anitha, J. et al. Earthworm-mediated synthesis of silver nanoparticles: a potent tool against hepatocellular carcinoma, Plasmodium falciparum parasites and malaria mosquitoes. Parasitol. Int. 65, 276–284 (2016).Article 
CAS 

Google Scholar 
23.Wahab, R., Khan, F. & Al-Khedhairy, A. A. Hematite iron oxide nanoparticles: apoptosis of myoblast cancer cells and their arithmetical assessment. RSC Adv. 8(44), 24750–24759 (2018).ADS 
CAS 
Article 

Google Scholar 
24.Ashley, E. A. et al. Spread of artemisinin resistance in Plasmodium falciparum malaria. N. Engl. J. Med. 371, 411–423 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
25.Rajan, R., Chandran, K., Harper, S. L., Yun, S. I. & Kalaichelvan, P. T. Plant extract synthesized nanoparticles: an ongoing source of novel biocompatible materials. Ind. Crop Prod. 70, 356–373 (2015).CAS 
Article 

Google Scholar 
26.Suresh, U. et al. Tackling the growing threat of dengue: Phyllanthus niruri-mediated synthesis of silver nanoparticles and their mosquitocidal properties against the dengue vector Aedes aegypti (Diptera: Culicidae). Parasitol. Res. 114, 1551–1562 (2015).PubMed 
Article 

Google Scholar 
27.Natarajan, K., Selvaraj, S. & Murty, V. R. Microbial production of silver nanoparticle. Digest J. Nanomat. Biostruct. 5, 135–140 (2010).
Google Scholar 
28.Song, Y. J., Jang, H. K. & Kim, S. B. Biological synthesis of gold nanoparticles using Magnolia kobus and Diopyros kaki leaf extract. Process Biochem. 44, 1133–1138 (2009).CAS 
Article 

Google Scholar 
29.Krishnan, R. & Maru, G. B. Isolation and analysis of polymeric polyphenol fractions from black tea. Food Chem. 94, 331–340 (2006).CAS 
Article 

Google Scholar 
30.Shankar, S., Rai, A., Ahmad, A. & Sastry, M. Rapid synthesis of Au, Ag and bimetallic Au core-Ag shell nanoparticles using Neem (Azadirachta indica) leaf broth. J. Colloid Interface Sci. 275, 496–550 (2004).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
31.Chandran, S. P., Chaudhary, M., Pasricha, R., Ahmad, A. & Sastry, M. Synthesis of gold nanotriangles and silver nanoparticles using Aloe vera plant extract. Biotechnol. Prog. 22, 577–583 (2006).CAS 
PubMed 
Article 

Google Scholar 
32.Benelli, G. Plant-synthesized nanoparticles: an eco-friendly tool against mosquito vectors? In Nanoparticles in the Fight Against Parasites Parasitology Research Monographs (ed. Mehlhorn, H.) 155–172 (Springer, 2015).
Google Scholar 
33.Sadraei, R. A simple method for preparation of nano-sized ZnO. Res. Rev. J. Chem. 5(2), 45–49 (2016).CAS 

Google Scholar 
34.Priyadarshini, K. A. et al. Biolarvicidal and pupicidal potential of silver nanoparticles synthesized using Euphorbia hirta against Anopheles stephensi Liston (Diptera: Culicidae). Parasitol. Res. 111(3), 997–1006 (2012).PubMed 
Article 

Google Scholar 
35.Satheeshkumar, K. & Kathireswari, P. Biological synthesis of Silver nanoparticles (Ag-NPS) by Lawsonia inermis (Henna) plant aqueous extract and its antimicrobial activity against human pathogens. Int. J. Curr. Microbiol. Appl. Sci. 5, 926–937 (2016).
Google Scholar 
36.Nareshkumar, G. et al. Electron channeling contrast imaging for III-nitride thin film structures. Mat. Sci. Semicon. Proc. 2016(47), 44–50 (2016).Article 
CAS 

Google Scholar 
37.Gandhi, S. & Madhusudhan, N. Retrieval of exoplanet emission spectra with HyDRA. Mon. Not. R. Astron. Soc. 47, 1–20 (2017).
Google Scholar 
38.Murugan, K. et al. Mosquitocidal and antiplasmodial activity of Senna occidentalis (Cassiae) and Ocimum basilicum (Lamiaceae) from Maruthamalai hills against Anopheles stephensi and Plasmodium falciparum. Parasitol. Res. 114, 3657–3664 (2015).PubMed 
Article 

Google Scholar 
39.Dinesh, D. et al. Mosquitocidal and antibacterial activity of green-synthesized silver nanoparticles from Aloe vera extracts: towards an effective tool against the malaria vector Anopheles stephensi?. Parasitol. Res. 114, 1519–1529 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
40.Pati, F. et al. Printing three-dimensional tissue analogues with decellularized extracellular matrix bioink. Nat. Commun. 5, 3935 (2014).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
41.Baxter, J. B. & Aydil, E. S. Nanowire based dye sensitized solar cells. Appl. Phys. Lett. 86, 53114 (2005).ADS 
Article 
CAS 

Google Scholar 
42.Reddy, K. M. et al. Selective toxicity of zinc oxide nanoparticles to prokaryotic and eukaryotic systems. Appl. Phys. Lett. 90(21), 213902–213903 (2007).ADS 
PubMed Central 
Article 
CAS 

Google Scholar 
43.Chwalibog, A. et al. Visualization of interaction between inorganic nano-particles and bacteria or fungi. Int. J. Nanomedicine. 2010(5), 1085–1094 (2010).Article 
CAS 

Google Scholar 
44.Saha, S., Dhanasekaran, D., Chandraleka, S. & Panneerselvam, C. A Synthesis, characterization and antimicrobial activity of cobalt metal complex against multi drug resistant bacterial and fungal pathogen Facta universitatis series. Phys. Chem. Technol. 7(1), 73–80 (2009).CAS 

Google Scholar 
45.Vivek, M., Kumar, P. S., Steffi, S. & Sudha, S. Biogenic silver nanoparticles by Gelidiella acerosa extract and their antifungal effects Avicenna. J. Med. Biotechnol. 3(3), 143 (2011).CAS 

Google Scholar 
46.Chobu, M., Nkwengulila, G., Mahande, A. M., Mwangonde, B. J. & Kweka, E. J. Direct and indirect effect of predators on Anopheles gambiae sensu stricto. Acta Trop. 142, 131–137 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
47.Murugan, K. et al. Hydrothermal synthesis of titanium dioxide nanoparticles: mosquitocidal potential and anticancer activity on human breast cancer cells (MCF-7). Parasitol. Res. 115, 1085–1096 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
48.Subramaniam, J. et al. Eco-friendly control of malaria and arbovirus vectors using the mosquitofish Gambusia affinis and ultra-low dosages of Mimusops elengi-synthesized silver nanoparticles: towards an integrative approach?. Environ. Sci. Pollut. Res. Int. 22(24), 20067–20083 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
49.Murugan, K. et al. Predation by Asian bullfrog tadpoles, Hoplobatrachus tigerinus, against the dengue vector, Aedes aegypti, in an aquatic environment treated with mosquitocidal nanoparticles. Parasitol. Res. 114, 3601–3610 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
50.Mahesh Kumar, P. et al. Mosquitocidal activity of Solanum xanthocarpum fruit extract and copepod Mesocyclops thermocyclopoides for the control of dengue vector Aedes aegypti. Parasitol. Res. 111, 609–618 (2012).PubMed 
Article 

Google Scholar 
51.Khooshe-Bast, Z., Sahebzadeh, N., Ghaffari-Moghaddam, M. & Mirshekar, A. Insecticidal effects of zinc oxide nanoparticles and Beauveria bassiana TS11 on Trialeurodes vaporariorum (Westwood, 1856) (Hemiptera: Aleyrodidae). Acta Agric Slov. 107(2), 299 (2016).CAS 
Article 

Google Scholar 
52.Ahmad, J., Wahab, R., Siddiqui, M. A., Saquib, Q. & Al-Khedhairy, A. A. Cytotoxicity and cell death induced by engineered nanostructures (quantum dots and nanoparticles) in human cell lines. J. Biol. Inorg. Chem. 25(2), 325–338 (2020).CAS 
PubMed 
Article 

Google Scholar 
53.Wahab, R. et al. Gold quantum dots impair the tumorigenic potential of glioma stem-like cells via β-catenin downregulation in vitro. Int. J. Nanomed. 14, 1131–1148 (2019).CAS 
Article 

Google Scholar 
54.Wahab, R., Saquib, Q. & Faisal, M. Zinc oxide nanostructures: a motivated dynamism against cancer cells. Process Biochem. 98(June), 83–92 (2020).CAS 
Article 

Google Scholar 
55.Wahab, R. et al. Microwave plasma-assisted silicon nanoparticles: cytotoxic, molecular, and numerical responses against cancer cells. RSC Adv. 9(23), 13336–13347 (2019).ADS 
CAS 
Article 

Google Scholar 
56.Anitha, J., Selvakumar, R. & Murugan, K. Chitosan capped ZnO nanoparticles with cell specific apoptosis induction through P53 activation and G2/M arrest in breast cancer cells—In vitro approaches. Int. J. Biol. Macromol. 136, 686–696 (2019).CAS 
PubMed 
Article 

Google Scholar 
57.Wahab, R. et al. Zinc oxide quantum dots: Multifunctional candidates for arresting C2C12 cancer cells and their role towards caspase 3 and 7 genes. RSC Adv. 6(31), 26111–26120 (2016).ADS 
CAS 
Article 

Google Scholar 
58.Liu, J. & Wang, Z. Increased oxidative stress a selective anticancer therapy. Oxid. Med. Cell. Longev. 2015, 294303 (2015).PubMed 
PubMed Central 

Google Scholar 
59.Droese, S. & Brandt, U. Molecular mechanisms of superoxide production by the mitochondrial respiratory chain. Adv. Exp. Med. Biol. 748, 145–169 (2012).CAS 
Article 

Google Scholar 
60.Gupta, S. C. et al. Upsides and downsides of reactive oxygen species for cancer: the roles of reactive oxygen species in tumorigenesis, prevention, and therapy. Antioxid. Redox Signal. 16, 1295–1322 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar  More

Experimental designOver 3 years, six paired organic British Blacktail flocks with intact beaks (i.e. not beak-trimmed) were visited between 1 and 40 weeks of age. Within each pair, one flock was ramp reared (RR) and one flock was control reared without ramps (CR). All flocks were kept on one farm which possessed two rearing houses and six laying sheds of approximately 2000 birds per flock. The site was multi-age, meaning that of the six laying sheds there were three different ages on the site at one time.The availability of this commercial facility enabled us to design an experiment whereby we allocated two rearing treatments, one with ramps provided to access elevated structures and a control with elevated structures but no ramps and to alternate these treatments between the two rearing houses available to avoid treatment x house confounds. Each rearing flock was moved independently to a laying house with no mixing, so we were able to continue data collection and examine any long-term effects of the rearing treatment during the laying period. Rearing flocks were systematically allocated so that each laying house received one RR flock and one CR flock during the experiment.Observations were made in the mornings at three time points during the rearing period at 1, 3 and 15–16 weeks, and three in the laying period at 16–17, 24 and 40 weeks of age. See Table 5 for a summary of experimental design, flock and housing information.Table 5 Experimental design for each ramp reared and control reared flock for the 6 replicates. There were two rearing sheds used, Rear1 (R1) and Rear2 (R2), with 6 different laying sheds named A1, A2, B1, B2, C1 and C2.Full size tableThe rearing sheds were static with 142.7 m2 of floor space covered with wood shavings. Rearing sheds were both set up with feed tracks giving mini pellet feed up to 11 weeks of age then pellet grower feed and 7 nipple drinker lines. The lighting schedule was 23 h light in the first day reducing gradually over the rearing period to 10 h light at 7 weeks of age. A minimum light intensity of 10 lx is required, but with windows and pop-holes light intensity was higher in the houses. The temperature was maintained at 30 °C during the first few days then slowly reduced to match the temperature in the laying sheds. Shed heating was provided by gas spot lamps, whole shed heating through hot pipes running along the length of the shed and hot air fans run by a biomass boiler. All flocks had access to the outside range by 10 weeks of age through two pop holes (each L: 2 m by H: 0.4 m). Flocks were moved the short distance from the rearing to the laying house at between 15 to 16 weeks of age in one night using transport modules.All rearing flocks had access to six elevated structures (ES) (see Fig. 3) from four days of age when the chicks were released from the brooding circles. Each ES comprised nine metal perches (length 302 cm, width 3.5 cm), with three perches (25 cm apart) at three different heights (43 cm, 73 cm and 103 cm). Two plastic grids (width 60 cm, length 115 cm) were fixed within the ES to provide platforms at different heights (Fig. 3). In each replicate, the RR flock had one ramp attached to each ES. Three of the ES were fitted with plastic grid ramps (width 60 cm, length 74 cm, angle 35.5°) leading up to the low perch and three ES had ramps (width 60 cm, length 115 cm, angle 40°) leading up to the middle perch. The CR flock had six ES without ramps.Figure 3Elevated structure dimensions used in the ramp reared sheds, (a) shows the high ramp (b) shows the low ramp. The control sheds elevated structures were identical to these but without ramps.Full size imageThe six single-tier laying houses on-site were mobile organic units with approximately 345m2 of floor space. See Fig. 4 for a schematic plan of their layout. All had a raised area comprising plastic slats over supports (approx. 70 cm from the litter) and a ground-level litter area covered with wood chip. Four of the sheds (Fig. 4a) were set up with the slatted area spanning the whole width of the shed and halfway down the length. In two of the sheds (Fig. 4b) the litter area was either side of the elevated slatted area. Nest boxes ran down the centre of the slatted area, dividing this into two sections. Intermittent ramps were installed at the level change, resulting in 4 m of ramp access and 4 m without ramps in the shed with litter at the end and 8 m of ramp access and 13 m without ramps in the shed with litter at the sides. In sheds A1 and A2 the height of the slatted area resulted in a steeper ramp angle of 45° compared to 30° in the other sheds. There were four pop holes at ground level with two on each side of the house (L: 2.35 m by H: 0.4 m) leading to the range from the litter area on both sides of the sheds. All sheds had aerial perches at 1 m high with 18 cm of perching space per bird resulting in approximately 360 m of perch length running the length of the slatted area. Feed tracks and drinker lines matched those in the rearing sheds. The lighting schedule was 16 h of light and varied between summer and winter with the lights set to turn off at the same time as natural dusk. The birds were fed on organic mini pellets throughout lay. Enrichment was provided to the flocks in the form of pecking objects such as buckets and boots. Replicates 5 and 6 were provided with pecking blocks and alfalfa hay nets hung on the litter area.Figure 4Plan view of the laying house layout (a) for replicates 1, 2, 4 and 5 and (b) for replicates 3 and 6. Images not to scale.Full size imageAssessments of behaviourObservations were made at three time points during the rearing period at 1, 3 and 15–16 weeks. On the first visit at 1 week of age, the total number of chicks on each ES was counted once in spot counts in the morning. At 3 and 15–16 weeks, observations of the movements up and down the ES were made. Three of the 6 structures were chosen at random in each shed. The number of chicks present on the different parts of the ES was counted at the beginning and end of the recording period to allow a comparison with the 1-week counts. The recordings involved 5-min continuous sampling where all movements down the ES were recorded and the area the chicks moved down from was noted. This was then repeated for movements up the ES. Focal bird recordings were taken at 3 and 15–16 weeks of age. Records were made for each of 3 randomly selected ES. When 10 focal birds had been observed (approximately 30 birds per flock), or 10 min had passed recordings stopped. A focal bird was chosen if it was performing orientation behaviour, indicating a downwards or upwards transition. This was described as the bird rotating its head to look in the direction of movement. Behaviours performed after the orientation behaviour were tallied, thus recorded as counts per behaviour (see Table 6). Recordings were stopped if birds completed a transition or moved away from transitioning.Table 6 An ethogram of behaviours of focal birds during up and down movements.Full size tableAt 15–16 weeks of age, three types of interactions were recorded for feather pecking. These included severe feather pecking (SFP), gentle feather pecking (GFP) and aggressive pecks (AP)28. A quadrat area 2 m by 2 m was randomly selected, with the number of birds in each quadrat counted at the beginning and end of the recording period. The number of SFP, GFP and AP were recorded over three minutes of continuous recordings in three different areas of the house, selected randomly at each end and the middle of the shed. Feather pecks and aggressive pecks were recorded as bouts: a series of pecks not separated by more than 5 s28. Rates of pecking were calculated as the number of pecks per bird per second.In the laying shed around 16–17 and at 24 weeks of age 3-min continuous sampling and focal bird recordings were taken for transitions between the slats and litter. Four recordings were made at 2-m lengths along the elevated slatted area: two areas with ramps (RA) and two areas without ramps (NRA) were selected. Separate recordings were taken for upwards and downwards movements and the number of birds in the recording area were counted at the start and end of the scans. At 16–17 and 24 weeks of age, feather pecking observations were taken using the same procedure as for the 15–16-week observations during rear.Welfare assessments and production data rearing phaseFeather scores of 20 birds per flock were recorded at 16–17 weeks of age by walking in a straight line down the centre of the shed, selecting a bird at random then counting two birds to the left of this and visually feather scoring that bird. Birds were not handled to minimise disturbance and plumage was scored using the method from Bright et al.33. The neck, back, rump, tail and wings were scored using a four-point scale 0 (best) to 4 (worst). Data were obtained from the farm records for percentage cumulative mortality and body weight.Welfare assessments and production data laying phaseAt 16–17 and 24 weeks of age, the attitude of the flocks was assessed using the approach distance and reactions to novel objects methodology developed by Whay et al.34. Distance to approach birds before they moved away was recorded by walking through the house selecting a bird at random and counting two birds to the left. The bird had to be standing up and facing the researcher, who approached the bird at a steady pace and recorded the distance before the bird moved away. This was repeated on 20 birds in each flock. Reactions to a novel object (blue folder at 17 weeks of age and a white and blue tub at 24 weeks of age) were assessed by placing a novel object on the ground and recording the time taken for the first bird to interact with it and then how many birds were within a 30 cm radius after 60 s. The novel object test was repeated in 4 areas per flock. Range use was recorded by counting the number of birds near to the house (5 m) in the middle range (5–20 m) and far (the rest of the range). Feather scores of 20 birds per flock were recorded at 17 and 24 weeks of age using the same procedure as for the 16-week assessment for birds at rear.At 40 weeks of age, feather cover and keel bone fractures were scored. Up to 100 birds per shed were caught from four different locations (25 litter, 25 slats, 25 perches, 25 nest boxes). In four sheds only 50 birds were caught as the birds were fearful and showed signs of distress. Feather cover was scored by picking the bird up and scoring the body and flight feathers separately using a the AssureWel three-point scale 0 (best) to 2 (worst)35. The keel damage was then scored using a 0 (no damage) to 2 scale based on the technique used by Wilkins et al.36. Validation for keel bone palpations was conducted. A score of 94% matched scores compared to an experienced gold standard assessor and 85% match at dissection for scoring a break. At 24 weeks of age, the number of floor eggs were counted over 1 day.Data were collected from the farm records on laying house percentage of daily eggs, average egg weight (grams), average hen body weight and feed conversion ratio.During the 16 week recordings in the final rearing flocks, the lighting inside the shed was considerably reduced compared to previous flocks. This resulted in poor visibility for feather cover and feather pecking observations, so these were not taken during this visit. Data were not obtained on keel fractures and feather cover scores at 40 weeks for the first laying flocks visited as their sheds were destroyed by strong winds.Statistical analysisData were analysed using SPSS 24 (IMB) or MLwiN 3.0. The statistical package MLwiN was chosen as it is designed for multilevel modelling and can therefore accommodate data nested within levels with repeated measures. Such models account for dependence between responses caused by grouping of birds within sheds, and repeated measures taken from the same sheds on different visits within and between replicates. Including visit and replicate as nested effects ensures that dependences (e.g. due to differing times of year when data were collected) are accounted for. All residuals were checked for normal distributions using a Shapiro-Wilks test or plotted graphically and no transformations were needed to meet the assumptions of the tests. All results are reported in the format mean ± SD unless when stated as the percentage of birds performing a behaviour during transitions.Assessments of behaviourAt rear, from the counts of chicks on structures and counts of transitions up and down the structures, a normal model (generalised linear model) was used with a four-level hierarchy (bird within shed within visit within replicate). The same normal model and four level hierarchy were used for the counts of transitions in the laying shed.For the focal bird behaviours of birds transitioning at rear and lay, the data were presented as the percentage of each behaviour calculated for the birds in the recording session for the two rearing treatments. The direction (up or down) was analysed separately. For the focal birds at lay, all were included in the analysis for the pre-transitioning behaviours, only birds that attempted a transition were included for analysis of the transitioning behaviours. Pre-transition behaviours for birds that moved-away and did not transition were analysed separately. Owing to the low occurrence of behaviours during the focal recordings for transitions up and down the ramps, data were coded as yes or no, and a Binomial model was used for analysis for both the rear and lay focal transition data with four hierarchical levels (bird within shed within visit within replicate).Welfare and production dataFor the Novel object test, human approach, feather pecking and feather cover data a normal model was used in MLwiN with four-levels (Bird within Shed within visit within replicate). Floor eggs were analysed using a two-tailed t-test in SPSS, due to limited data. Ordinal data such a keel bone fracture scores and feather cover recorded at 40 weeks of age were converted to binomial data due to a lack of data for some scores, these were therefore analysed using a binomial model in MLwiN with two levels (Bird within shed).Production data at rear (body weight in grams) and lay (% eggs daily, egg weight in grams, body weight in grams and feed conversion ratio) were obtained from farm records and analysed in SPSS using a general linear model with treatment (CR and RR) as a fixed factor and age (3, 8 and 14 weeks at rear and 20, 30 and 70 weeks at lay) as a random factor to account for repeated results. Cumulative percentage mortality was analysed at 14 weeks of age using a t-test to compare the treatment groups.Ethical approvalEthical approval for this project was granted by the University of Bristol’s Animal Welfare and Ethical review body under UIN: UB/16/040 and all methods were conducted in accordance with the review body and UK legislation. 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