Philippa Kaur

1.
Wiseman, B. R. Plant-resistance to insects in integrated pest-management. Plant Dis. 78, 927–932. https://doi.org/10.1094/pd-78-0927 (1994).
Article  Google Scholar 
2.
Painter, R. H. Insect resistance in crop plants. Soil Sci. 72 (1951).

3.
Orr, D. B. & Boethel, D. J. Influence of plant antibiosis through four trophic levels. Oecologia 70, 242–249. https://doi.org/10.1007/BF00379247 (1986).
ADS  CAS  Article  PubMed  Google Scholar 

4.
Smith, C. M. & Clement, S. L. Molecular Bases of Plant Resistance to Arthropods. Annu. Rev. Entomol. 57, 309–328. https://doi.org/10.1146/annurev-ento-120710-100642 (2011).
CAS  Article  PubMed  Google Scholar 

5.
Radcliffe, R. H. in Radcliffe’s IPM world textbook Vol. https://ipmworld.umn.edu/ratcliffe-hessian-fly (eds Radcliffe E.B. & Hutchison W.D.) (University of Minnesota, 1997).

6.
Kosma, D. K., Nemacheck, J. A., Jenks, M. A. & Williams, C. E. Changes in properties of wheat leaf cuticle during interactions with Hessian fly. Plant J 63, 31–43 (2010).
CAS  PubMed  Google Scholar 

7.
Smiley, R. W., Gourlie, J. A., Whittaker, R. G., Easley, S. A. & Kidwell, K. K. Economic impact of Hessian fly (Diptera: Cecidomyiidae) on spring wheat in Oregon and additive yield losses with Fusarium crown rot and lesion nematode. J Econ Entomol 97, 397–408 (2004).
Article  Google Scholar 

8.
Harris, M. O., Sandanayaka, M. & Griffin, A. Oviposition preferences of the Hessian fly and their consequences for the survival and reproductive potential of offspring. Ecol. Entomol. 26, 473–486. https://doi.org/10.1046/j.1365-2311.2001.00344.x (2001).
Article  Google Scholar 

9.
Ganehiarachchi, G. A. S. M., Anderson, K. M., Harmon, J. & Harris, M. O. Why oviposit there? Fitness consequences of a gall midge choosing the plant’s youngest leaf. Environ Entomol 42, 123–130 (2013).
CAS  Article  Google Scholar 

10.
Kanno, H. & Harris, M. O. Physical features of grass leaves influence the placement of eggs within the plant by the Hessian fly. Entomol. Exp. Appl. 96, 69–80. https://doi.org/10.1046/j.1570-7458.2000.00680.x (2000).
Article  Google Scholar 

11.
Harris, M. O. & Rose, S. Chemical, color, and tactile cues influencing oviposition behavior of the Hessian fly (Diptera, Cecidomyiidae). Environ. Entomol. 19, 303–308. https://doi.org/10.1093/ee/19.2.303 (1990).
Article  Google Scholar 

12.
Kanno, H. & Harris, M. O. Leaf physical and chemical features influence selection of plant genotypes by hessian fly. J. Chem. Ecol. 26, 2335–2354 (2000).
CAS  Article  Google Scholar 

13.
Cervantes, D. E., Eigenbrode, S. D., Ding, H. J. & Bosque-Perez, N. A. Oviposition responses by Hessian fly, Mayetiola destructor, to wheats varying in surface waxes. J. Chem. Ecol. 28, 193–210 (2002).
CAS  Article  Google Scholar 

14.
Morris, B. D., Foster, S. P. & Harris, M. O. Identification of 1-octacosanal and 6-methoxy-2-benzoxazolinone from wheat as ovipositional stimulants for Hessian fly, Mayetiola destructor. J. Chem. Ecol. 26, 859–873 (2000).
CAS  Article  Google Scholar 

15.
Harris, M. O., Rose, S. & Malsch, P. The role of vision in the host plant-finding behavior of the Hessian fly. Physiol. Entomol. 18, 31–42. https://doi.org/10.1111/j.1365-3032.1993.tb00446.x (1993).
Article  Google Scholar 

16.
Rohfritsch, O. A fungus associated gall midge, Lasioptera arundinis (Schiner), on Phragmites australis (Cav) Trin. Bull. Soc. Bot. France Lett. Bot. 139, 45–59. https://doi.org/10.1080/01811797.1992.10824942 (1992).
Article  Google Scholar 

17.
Schmid, R. B., Knutson, A., Giles, K. L. & McCornack, B. P. Hessian fly (Diptera: Cecidomyiidae) biology and management in wheat. J. Integr. Pest Manag. 9, 12. https://doi.org/10.1093/jipm/pmy008 (2018).
Article  Google Scholar 

18.
Gagné, R. J. & Hatchett, J. H. Instars of the Hessian Fly (Diptera: Cecidomyiidae). Ann. Entomol. Soc. Am. 82, 73–79. https://doi.org/10.1093/aesa/82.1.73 (1989).
Article  Google Scholar 

19.
Lidell, M. C. & Schuster, M. F. Distribution of the Hessian fly and its control in Texas. Southwestern Entomologist 15, 133–145 (1990).
Google Scholar 

20.
Morgan, G., Sansone, C. & Knutson, A. Hessian fly in Texas wheat. E-350 (Texas A&M, 2005).

21.
Flanders, K. L., Reisig, D. D., Buntin, G. D., Herbert, J. D. A. & Johnson, D. W. Biology and management of Hessian fly in the Southeast. ANR1069 (Alabama Cooperative Extension System, 2013).

22.
Wellso, S. G. Aestivation and Phenology of the Hessian Fly (Diptera: Cecidomyiidae) in Indiana. Environ. Entomol. 20, 795–801. https://doi.org/10.1093/ee/20.3.795 (1991).
Article  Google Scholar 

23.
Boyd, M. L. & Bailey, W. C. Hessian fly management on wheat. G7180 (Missouri Extension, University of Missouri-Columbia, 2000).

24.
Ando, K. et al. Genome-wide associations for multiple pest resistances in a Northwestern United States elite spring wheat panel. PLoS One 13, e0191305/0191301-e0191305/0191325. https://doi.org/10.1371/journal.pone.0191305 (2018).

25.
Anderson, K. M. & Harris, M. O. Susceptibility of North Dakota Hessian Fly (Diptera: Cecidomyiidae) to 31 H Genes Mediating Wheat Resistance. J. Econ. Entomol. 112, 2398–2406 (2019).
Article  Google Scholar 

26.
Sardesai, N., Nemacheck, J. A., Subramanyam, S. & Williams, C. E. Identification and mapping of H32, a new wheat gene conferring resistance to Hessian fly. Theor. Appl. Genet. 111, 1167–1173 (2005).
CAS  Article  Google Scholar 

27.
Zhu, L., Liu, X. & Chen, M.-S. Differential accumulation of phytohormones in wheat seedlings attacked by avirulent and virulent Hessian fly (Diptera: Cecidomyiidae) larvae. J. Econ. Entomol. 103, 178–185 (2010).
CAS  Article  Google Scholar 

28.
Mithöfer, A. & Boland, W. Recognition of Herbivory-Associated Molecular Patterns. Plant Physiol. 146, 825. https://doi.org/10.1104/pp.107.113118 (2008).
CAS  Article  PubMed  PubMed Central  Google Scholar 

29.
Stuart, J. J., Chen, M.-S., Shukle, R. & Harris, M. O. Gall midges (Hessian flies) as plant pathogens. Annu. Rev. Phytopathol. 50, 339–357 (2012).
CAS  Article  Google Scholar 

30.
Liu, X. et al. Gene expression of different wheat genotypes during attack by virulent and avirulent Hessian fly (Mayetiola destructor) larvae. J. Chem. Ecol. 33, 2171–2194 (2007).
CAS  Article  Google Scholar 

31.
Subramanyam, S. et al. Expression of two wheat defense-response genes, Hfr-1 and Wci-1, under biotic and abiotic stresses. Plant Sci. 170, 90–103. https://doi.org/10.1016/j.plantsci.2005.08.006 (2006).
CAS  Article  Google Scholar 

32.
Wu, J. et al. Differential responses of wheat inhibitor-like genes to Hessian fly, Mayetiola destructor, attacks during compatible and incompatible interactions. J. Chem. Ecol. 34, 1005–1012 (2008).
CAS  Article  Google Scholar 

33.
Giovanini, M. P. et al. A novel wheat gene encoding a putative chitin-binding lectin is associated with resistance against Hessian fly. Mol. Plant Pathol. 8, 69–82 (2007).
CAS  Article  Google Scholar 

34.
Liu, X. et al. Reactive oxygen species are involved in plant defense against a gall midge. Plant Physiol. 152, 985. https://doi.org/10.1104/pp.109.150656 (2010).
CAS  Article  PubMed  PubMed Central  Google Scholar 

35.
Bari, R. & Jones, J. D. G. Role of plant hormones in plant defence responses. Plant Mol. Biol. 69, 473–488. https://doi.org/10.1007/s11103-008-9435-0 (2009).
CAS  Article  PubMed  Google Scholar 

36.
Denancé, N., Sánchez-Vallet, A., Goffner, D. & Molina, A. Disease resistance or growth: the role of plant hormones in balancing immune responses and fitness costs. Frontiers in Plant Science 4. https://doi.org/10.3389/fpls.2013.00155 (2013)

37.
Dinh, S. T., Baldwin, I. T. & Galis, I. The HERBIVORE ELICITOR-REGULATED1 gene enhances abscisic acid levels and defenses against herbivores in Nicotiana attenuate plants. Plant Physiol. 162, 2106–2124 (2013).
CAS  Article  Google Scholar 

38.
War, A. R., Paulraj, M. G., War, M. Y. & Ignacimuthu, S. Role of salicylic acid in induction of plant defense system in chickpea (Cicer arietinum L.). Plant signaling & behavior 6, 1787–1792. https://doi.org/10.4161/psb.6.11.17685 (2011).

39.
Nguyen, D., Rieu, I., Mariani, C. & van Dam, N. M. How plants handle multiple stresses: hormonal interactions underlying responses to abiotic stress and insect herbivory. Plant Mol. Biol. 91, 727–740. https://doi.org/10.1007/s11103-016-0481-8 (2016).
CAS  Article  PubMed  PubMed Central  Google Scholar 

40.
Lee, A. et al. Inverse correlation between jasmonic acid and salicylic acid during early wound response in rice. Biochem. Biophys. Res. Commun. 318, 734–738. https://doi.org/10.1016/j.bbrc.2004.04.095 (2004).
CAS  Article  PubMed  Google Scholar 

41.
Kunkel, B. N. & Brooks, D. M. Cross talk between signaling pathways in pathogen defense. Curr. Opin. Plant Biol. 5, 325–331. https://doi.org/10.1016/S1369-5266(02)00275-3 (2002).
CAS  Article  PubMed  Google Scholar 

42.
Farmer, E. E., Alméras, E. & Krishnamurthy, V. Jasmonates and related oxylipins in plant responses to pathogenesis and herbivory. Curr. Opin. Plant Biol. 6, 372–378. https://doi.org/10.1016/S1369-5266(03)00045-1 (2003).
CAS  Article  PubMed  Google Scholar 

43.
Loake, G. & Grant, M. Salicylic acid in plant defence—the players and protagonists. Curr. Opin. Plant Biol. 10, 466–472. https://doi.org/10.1016/j.pbi.2007.08.008 (2007).
CAS  Article  PubMed  Google Scholar 

44.
Felton, G. W., Bi, J. L., Summers, C. B., Mueller, A. J. & Duffey, S. S. Potential role of lipoxygenases in defense against insect herbivory. J. Chem. Ecol. 20, 651–666. https://doi.org/10.1007/BF02059605 (1994).
CAS  Article  PubMed  Google Scholar 

45.
Audenaert, K., De Meyer, G. B. & Höfte, M. M. Abscisic Acid Determines Basal Susceptibility of Tomato to Botrytis cinerea and Suppresses Salicylic Acid-Dependent Signaling Mechanisms. Plant Physiol. 128, 491. https://doi.org/10.1104/pp.010605 (2002).
CAS  Article  PubMed  PubMed Central  Google Scholar 

46.
Mohr, P. G. & Cahill, D. M. Suppression by ABA of salicylic acid and lignin accumulation and the expression of multiple genes, in Arabidopsis infected with Pseudomonas syringae pv. tomato. Functional & Integrative Genomics 7, 181–191, https://doi.org/10.1007/s10142-006-0041-4 (2007).

47.
Harris, M. O., Dando, J. L., Griffin, W. & Madie, C. Susceptibility of cereal and non-cereal grasses to attack by Hessian fly (Mayetiola destructor (Say)). N. Zeal. J. Crop Hortic. Scie.ce 24, 229–238. https://doi.org/10.1080/01140671.1996.9513957 (1996).
Article  Google Scholar 

48.
Gitelson, A. A. & Merzlyak, M. N. Signature analysis of leaf reflectance spectra: algorithm development for remote sensing of chlorophyll. J. Plant Physiol. 148, 494–500. https://doi.org/10.1016/S0176-1617(96)80284-7 (1996).
CAS  Article  Google Scholar 

49.
Foster, S. P. & Harris, M. O. Foliar chemicals of wheat and related grasses influencing oviposition by Hessian fly, Mayetiola destructor (Say) (Diptera: Cecidomyiidae). J Chem Ecol 18, 1965–1980 (1992).
CAS  Article  Google Scholar 

50.
Gagne, R. J., Hatchett, J. H., Lhaloui, S. & El Bouhssini, M. Hessian fly and barley stem gall midge, two different species of mayetiola (Diptera: Cecidomyiidae) in Morocco. Ann. Entomol. Soc. Am. 84, 436–443. https://doi.org/10.1093/aesa/84.4.436 (1991).
Article  Google Scholar 

51.
Cherif, A., Kinoshita, N., Taylor, D. & Mediouni Ben Jemâa, J. Molecular characterization and phylogenetic comparisons of three Mayetiola species (Diptera: Cecidomyiidae) infesting cereals in Tunisia. Applied Entomology and Zoology 52, 543–551, https://doi.org/10.1007/s13355-017-0507-y (2017).

52.
Gould, F. Simulation models for predicting durability of insect-resistant germ plasm: hessian fly (Diptera: Cecidomyiidae)-resistant Winter Wheat. Environ. Entomol. 15, 11–23. https://doi.org/10.1093/ee/15.1.11 (1986).
Article  Google Scholar 

53.
Chen, M.-S., Liu, X., Wang, H. & El-Bouhssini, M. Hessian fly (Diptera: Cecidomyiidae) interactions with barley, rice, and wheat seedlings. J Econ Entomol 102, 1663–1672 (2009).
Article  Google Scholar 

54.
Ratcliffe, R. H., Safranski, G. G., Patterson, F. L., Ohm, H. W. & Taylor, P. L. Biotype status of Hessian fly (Diptera, Cecidomyiidae) populations from the eastern United-States and their response to 14 Hessian fly resistance genes. J. Econ. Entomol. 87, 1113–1121. https://doi.org/10.1093/jee/87.4.1113 (1994).
Article  Google Scholar 

55.
Tooker, J. F. & Frank, S. D. Genotypically diverse cultivar mixtures for insect pest management and increased crop yields. J. Appl. Ecol. 49, 974–985. https://doi.org/10.1111/j.1365-2664.2012.02173.x (2012).
Article  Google Scholar 

56.
Erb, M., Meldau, S. & Howe, G. A. Role of phytohormones in insect-specific plant reactions. Trends Plant Sci. 17, 1–20 (2012).
Article  Google Scholar 

57.
Williams, C. E., Collier, C. C., Nemacheck, J. A., Liang, C. Z. & Cambron, S. E. A lectin-like wheat gene responds systemically to attempted feeding by avirulent first-instar Hessian fly larvae. J. Chem. Ecol. 28, 1411–1428. https://doi.org/10.1023/a:1016200619766 (2002).
CAS  Article  PubMed  Google Scholar 

58.
Herrera-Vasquez, A., Salinas, P. & Holuigue, L. Salicylic acid and reactive oxygen species interplay in the transcriptional control of defense genes expression (vol 6, 171, 2015). Frontiers in Plant Science 8, https://doi.org/10.3389/fpls.2017.00964 (2017).

59.
Hatchett, J. H., Kreitner, G. L. & Elzinga, R. J. Larval Mouthparts and Feeding Mechanism of the Hessian Fly (Diptera: Cecidomyiidae). Ann. Entomol. Soc. Am. 83, 1137–1147. https://doi.org/10.1093/aesa/83.6.1137 (1990).
Article  Google Scholar 

60.
Schotzko, D. J. & Bosque-Perez, N. A. Relationship between Hessian fly infestation density and early seedling growth of resistant and susceptible wheat. J. Agric. Urban Entomol. 19, 95–107 (2002).
Google Scholar 

61.
Ratcliffe, R. H. et al. Biotype composition of Hessian fly (Diptera: Cecidomyiidae) populations from the southeastern, midwestern, and northwestern United States and virulence to resistance genes in wheat. J Econ. Entomol. 93, 1319–1328 (2000).
CAS  Article  Google Scholar 

62.
Song, S., Gong, W., Zhu, B. & Huang, X. Wavelength selection and spectral discrimination for paddy rice, with laboratory measurements of hyperspectral leaf reflectance. ISPRS J. Photogram. Remote Sens. 66, 672–682. https://doi.org/10.1016/j.isprsjprs.2011.05.002 (2011).
ADS  Article  Google Scholar 

63.
Dechant, B., Cuntz, M., Vohland, M., Schulz, E. & Doktor, D. Estimation of photosynthesis traits from leaf reflectance spectra: correlation to nitrogen content as the dominant mechanism. Remote Sens. Environ. 196, 279–292. https://doi.org/10.1016/j.rse.2017.05.019 (2017).
ADS  Article  Google Scholar 

64.
Ollinger, S. V. Sources of variability in canopy reflectance and the convergent properties of plants. New Phytol. 189, 375–394. https://doi.org/10.1111/j.1469-8137.2010.03536.x (2011).
CAS  Article  PubMed  Google Scholar 

65.
Almeida Trapp, M., De Souza, G. D., Rodrigues-Filho, E., Boland, W. & Mithöfer, A. Validated method for phytohormone quantification in plants. Frontiers in Plant Science 5, https://doi.org/10.3389/fpls.2014.00417 (2014).

66.
Davis, T. S., Bosque-Pérez, N. A., Popova, I. & Eigenbrode, S. D. Evidence for additive effects of virus infection and water availability on phytohormone induction in a staple crop. Frontiers in Ecology and Evolution 3, https://doi.org/10.3389/fevo.2015.00114 (2015). More

Radiation experiments
We experimentally quantified radiolytic hydrogen (H2) production in (i) pure water, (ii) seawater, and (iii) seawater-saturated sediment. We irradiated these materials with α- or γ-radiation for fixed time intervals and then determined the concentrations of H2 produced. Sediment samples were slurried with natural seawater to achieve a slurry porosity (φ) of ~0.83, which is the average porosity of abyssal clay in the South Pacific Gyre34. The seawater source is described below. To avoid microbiological uptake of radiolytic H2 during the course of the experiment, seawater and marine sediment slurries were pre-treated with HgCl2 (0.05% solution) or NaN3 (0.1% wt/vol). To ensure that addition of these chemicals did not impact radiolytic H2 yields, irradiation experiments with pure water plus HgCl2 or NaN3 were also conducted. HgCl2 or NaN3 addition had no statistically significant impact on H2 yields5,6,10.
Experimental samples were irradiated in 250 mL borosilicate vials. A solid-angle 137Cs source (beam energy of 0.67 MeV) was used for the γ-irradiation experiments at the Rhode Island Nuclear Science Center (RINSC). The calculated dose rate for sediment slurries was 2.19E−02 Gy h−1 accounting for the (i) source activity, (ii) distance between the source and the samples, (iii) sample vial geometry, and (iv) attenuation coefficient of γ-radiation through air, borosilicate, and sediment slurry. 210Po (5.3 MeV decay−1)-plated silver strips with total activities of 250 μCi were used for the α-irradiation experiments. For α-irradiation of each sediment slurry, a 210Po-plated strip was placed inside the borosilicate vial and immersed in the slurry. Calculated total absorbed doses were 4 Gy and 3 kGy for γ-irradiation and α-irradiation experiments, respectively.
The settling time of sediment grains in the slurries (1 week) was long compared to the time span of each experiment (tens of minutes to an hour for α-experiments, hours to days for γ-experiments). Therefore, we assumed that the suspension was homogenous during the course of each experiment.
H2 concentrations were measured by quantitative headspace analysis via gas chromatography. For headspace analysis, 30 mL of N2 was first injected into the sample vial. To avoid over-pressurization of the sample during injection, an equivalent amount of water was allowed to escape through a separate needle. The vials were then vigorously shaken for 5 min to concentrate the H2 into the headspace. Finally, a 500-μL-headspace subsample was injected into a reduced gas analyzer (Peak Performer 1, PP1). The reduced gas analyzer was calibrated using a 1077 ppmv H2 primary standard (Scott-Marrin, Inc.). A gas mixer was used to dilute the H2 standard with N2 gas to obtain various H2 concentrations and produce a five-point linear calibration curve (0.7, 2, 5, 20, and 45 ppm). H2 concentrations of procedural blanks consisting of sample vials filled with non-irradiated deionized 18-MΩ water were also determined. The concentration detection limit obtained using this protocol was 0.8–1 nM H2. Relative error was less than 5%. Radiation experiments were performed at a minimum in triplicate.
Sample selection and experimental radiolytic H2 yields, G(H2)
Millipore Milli-Q system water was used for our pure-water experiments. For seawater experiments, we used bottom water collected in the Hudson Canyon (water depth, 2136 m) by RV Endeavor expedition EN534. Salinity of North Atlantic bottom water in the vicinity of the Hudson Canyon (34.96 g kg−1) is similar to that of mean open-ocean bottom water (34.70 g kg−1)44,45.
The 20 sediment samples used for the experiments were collected by scientific coring expeditions in three ocean basins (expedition KN223 to the North Atlantic46, expedition KN195-3 to the Equatorial and North Pacific47, International Ocean Discovery Program (IODP) Expedition 329 to the South Pacific Gyre34, MONA expedition to the Guaymas Basin48, expedition EN32 to the Gulf of Mexico49, and expedition EN20 to the Venezuela Basin50). To capture the dominant sediment types present in the global ocean, we selected samples typical of five common sediment types [abyssal clay (11 samples), nannofossil-bearing clay or calcareous marl (2 samples), clay-bearing diatom ooze (3 samples), calcareous ooze (2 samples), and lithogenous sediment (2 samples)]. The locations, lithological descriptions, and mineral compositions of the samples are given in Supplementary Tables 1,  2,  3, and Supplementary Fig. 1. Additional chemical and physical descriptions of the sediment samples used in the radiation experiments can be found in the expedition reports for the expeditions on which the samples were collected34,46.
Energy-normalized radiolytic H2 yields are commonly expressed as G(H2)-values (molecules H2 per 100 eV absorbed)1. As shown in Supplementary Fig. 2, for all irradiated samples (pure water, seawater, and marine sediment slurries), H2 production increased linearly with absorbed α- and γ-ray-dose. We calculated G(H2)-values for each sample and radiation type (α or γ) as the slope of the least-square regression line of radiolytic H2 concentration versus absorbed dose (Supplementary Fig. 2). The error on the yields is less than 10% for each sample. G(H2)-values for each sample and radiation type (α or γ) are reported in Supplementary Table 3.
Although radiolytic OH• is known to react with dissolved organic matter51, total organic content does not appear to significantly impact radiolytic H2 production, since the most organic-rich sediment (e.g., Guaymas Basin and Gulf of Mexico sediment) did not yield particularly high H2 (Supplementary Table 3).
Calculated radiolytic production rates of H2 and oxidants in the cored sediment of individual sites
We calculated radiolytic H2 production rates (PH2, in molecules H2 cm−3 yr−1) for the cored sediment column at nine sites with oxic subseafloor sediment in the North Pacific, South Pacific, and North Atlantic; and seven sites with anoxic subseafloor sediment in the Bering Sea, South Pacific, Equatorial Pacific, and Peru Margin (see Supplementary Fig. 3 for site locations). For these calculations, we used the following equation from Blair et al.2:

$$P_{{mathrm{H}}_2} = {sum} A _{{mathrm{m}},i}rho left( {1 – {{upvarphi}} } right)E_i{mathrm{G}}({mathrm{H}}_2)_i$$
(1)

where i is alpha, beta, or gamma radiation; Am is radioactivity per mass solid; φ is porosity; ρ is density solid; (E_i) is decay energy; and ({mathrm{G}}({mathrm{H}}_2)_i) is radiolytic yield.
We calculated radiolytic oxidant production rates for these sediment columns from the H2 production rates. Because H2 production and oxidant production are stochiometrically balanced in water radiolysis [2H2O → H2 + H2O2], the calculated radiolytic H2 production rates (in electron equivalents) are equal to radiolytic oxidant production rates (in electron equivalents).
The in situ γ- and α-radiation dosages in marine sediment are, respectively, 13 and 15 orders of magnitude lower than the dosage used in our experiments. Because the measured G(H2) for pure water in our γ-irradiation experiment (dose rate = 2.19E-02 Gy h−1) is statistically indistinguishable from previously published G(H2) values at much higher dose rates (ca. 1.00E+3 Gy h−1)5, we infer that the γ-irradiation G(H2) value is constant with dose rate over five orders of magnitude. Similarly, our experimental pure water H2 yields following α-particle irradiation from a 210Po-source (dose rate of 2.55E+03 Gy h−1) are indistinguishable from the yield obtained by Crumière et al.6 [G(H2) = 1.30 ± 0.13] for air-saturated deionized water exposed to a cyclotron-generated He2+ particle beam at higher dose rate (dose rate 1.62E+05 Gy h−1). The close similarity in H2 yields obtained in both experiments implies that (i) radiolytic H2 yield from α-particle irradiation is identical to that from cyclotron-generated He2+ particle irradiation, and (ii) this yield is constant over a two-orders-of-magnitude range dose rate. Therefore, we use our experimentally determined α- and γ-irradiation G(H2) values for the low radiation dose rate found in the subseafloor. Because the G(H2) of β irradiation has not been experimentally determined for water-saturated materials, we assume that the G(H2) of β-radiation matches the G(H2) of γ-radiation for the same sediment types. In pure water, their G(H2) values differ by only 17%1. Because β radiation, on average, contributes only 11% of the total radiolytic H2 production from the U, Th series and K decay in marine sediment, these estimates of total H2 production differ by only 2–5% relative to estimates where the G(H2) of β radiation is assumed equal to that for pure water or for α radiation of the same sediment types.
To calculate H2 production rates for the entire sediment column at seven South Pacific sites and two North Atlantic sites, we measured downcore sediment profiles of U, Th, and K (i) 187 sediment samples from IODP Expedition 329 Sites U1365, U1366, U1367, U1368, U1369, U1370, and U137134,52, and (ii) 40 samples from KN223 expedition Sites 11 and 12 (ref. 46). Total U and Th (ppm) and K2O (wt%) for these sites are reported in the EarthChem SedDB data repository. We measured U, Th, and K abundances using standard atomic emission and mass spectrometry techniques (i.e. ICP-ES and ICP-MS) in the Analytical Geochemistry Facilities at Boston University. Sample preparation, analytical protocol, and data are reported in Dunlea et al.52. The precision for each element is ~2% of the measured value, based on three separate digestions of a homogenized in-house standard of deep-sea sediment.
To calculate H2 production rates for the sediment columns at North Pacific coring Sites EQP10 and EQP11 (ref. 47), we used radioactive element content data from Kyte et al.53, who measured chemical concentrations at high resolution in bulk sediment in core LL44-GPC3. Because Site EQP11 was cored at the same location as LL44-GPC3 (ref. 53) and the sediment retrieved at all three sites is homogeneous abyssal clay, we assume the radioactive element abundances measured in core LL44-GPC3 to be representative of Sites EQP10 and EQP11 (ref. 47). Calculated radiolytic H2 production rates for South Pacific sites are listed in Supplementary Table 4 and for North Atlantic and North Pacific sites in Supplementary Table 5.
For Bering Sea Sites U1343 and U1345 (ref. 54), sedimentary U, Th, and K content measurements are unavailable. Since sediment recovered at these two sites is primarily siliciclastic with a varying amount of diatom-rich clay, we use U, Th, and K concentration values reported for upper continental crust by Li and Schoonmaker for these Bering Sea sites55. Finally, we calculate downhole radiolytic H2 production rates for ODP Leg 201 Sites 1225, 1226, 1227, and 1230 (ref. 35). Sediment compositions for these sites include nannofossil-rich calcareous ooze (Site 1225), alternation of nannofossil (calcareous) ooze and diatom ooze (Site 1226), and siliciclastic with diatom-rich clay intervals (Sites 1227 and 1230). Because sedimentary U, Th, and K measurements are not available for Leg 201 sites, we used average U, Th, and K concentration values measured in North Atlantic46 and South Pacific Sites34,52 with corresponding lithologies.
We use isotopic abundance values reported in Erlank et al.56 to calculate the abundance of 238U, 235U, 232Th, and 40K from the measured ICP-MS values of total U, Th, and K concentration. We then converted radionuclide concentrations to activities using Avogadro’s number and each isotope’s decay constant2. We refer to Blair et al. for a detailed explanation of activities and radiolytic yield calculations2.
Calculation of global radiolytic H2 and oxidant production rates in marine sediment
We calculated global radiolytic H2 production in ocean sediment by applying Eq. (1) (ref. 2) globally. As with the rates at individual sites, we calculated global radiolytic oxidant production (in electron equivalents) from global H2 production and the stochiometry of water radiolysis [2H2O → H2 + H2O2].
Our global radiolytic H2 production calculation spatially integrates calculations of sedimentary porewater radiolysis rates that are based on (i) our experimentally constrained radiolytic H2 yields for the principal marine sediment types, (ii) measured radioactive element content of sediment cores in three ocean basins (North Atlantic46, North Pacific53, and South Pacific34,52), and (iii) global distributions of sediment lithology57, sediment porosity58, and sediment column length59,60.
To generate the global map of radiolytic H2 production, we created global maps of seafloor U, Th, and K concentrations, density, G(H2)-α values, and G(H2)-γ-and-β by assigning each grid cell in our compiled seafloor lithology map (Supplementary Fig. 4) its lithology-specific set of input variables (Supplementary Table 6). Because our model assumes that lithology is constant with depth, U, Th, and K content, grain density, and G(H2)-values are constant with depth.
The G(H2)-values (α, β, and γ radiation), radioactive element content (sedimentary U, Th, and K concentration), density, porosity, and sediment thickness are determined as follows.
Radiolytic yield [G(H2)] for α,β-&-γ radiation
Radiolytic yields for the main seafloor lithologies are obtained by averaging experimentally derived yields for the respective lithologies (Supplementary Table 6). We assume that G(H2)-β values equal G(H2)-γ values.
Sediment lithology
For these calculations of radiolytic chemical production, we generally used seafloor lithologies and assumed that sediment type is constant with sediment depth. For seafloor lithology, the geographic database of global bottom sediment types57 was compiled into five lithologic categories: abyssal clay, calcareous ooze, siliceous ooze, calcareous marl, and lithogenous (Supplementary Fig. 4). Some areas of the seafloor are not described in the database57. These include (i) high-latitude regions (as the seafloor lithology database extends from 70°N to 50°S)57 and (ii) some discrete areas located along continental margins (e.g., Mediterranean Sea, Timor Sea, South China Sea, Supplementary Fig. 4). We used other data sources to identify seafloor lithologies for these regions. We added an opal belt (siliceous ooze) in the Southern Ocean between 57°S and 66°S61,62. The geographic extent of this opal belt was based on DeMaster62 and Dutkiewicz et al.61. We defined the areas of the seafloor from 50°S to 57°S, from 66°S to 90°S, and in the Arctic Ocean as mostly composed of lithogenous material, based on (i) drillsite lithologies in the Southern Ocean [ODP: Site 695 (ref. 63), Site 694 (ref. 63), Site 1165 (ref. 64), Site 739 (ref. 65)], the Bering Sea and Arctic Ocean [International Ocean Discovery Program (IODP): Sites U1343 and U1345 (ref. 54), Site M0002 (ref. 66), ODP: Site 910 (ref. 67), Site 645 (ref. 68)] and between 50°S and 57°S [Deep Sea Drilling Project (DSDP): Site 326 (ref. 69), Ocean Drilling Program (ODP): Site 1138 (ref. 70), Site 1121 (ref. 71)], and Dutkiewicz et al.61.
In the North and South Atlantic, sediment type can be very different at depth than at the seafloor. For these regions, we departed from our assumption that sediment lithology is the same at depth as at the seafloor. Subseafloor lithologies at ODP Sites [1063 (ref. 72), 951 (ref. 73), 925 (ref. 74), and 662 (ref. 75)] and IODP Sites [U1403 (ref. 76) and U1312 (ref. 77)] indicate that sediment in the Atlantic Ocean basin is generally 30–90% biogenic carbonate content and detrital clay78, even where the seafloor lithology is abyssal clay57. Therefore, regions in the Atlantic Ocean described as abyssal clay in the seafloor lithology database57 were characterized as calcareous marl for our calculations (Supplementary Fig. 4). Because abyssal clay catalyzes radiolytic H2 production at a higher rate than calcareous marl, this characterization may underestimate production of radiolytic H2 and radiolytic oxidants in these Atlantic regions.
Radioactive element content
For four of the five lithologic types in our global maps (abyssal clay, siliceous ooze, calcareous ooze, and calcareous marl), we average U, Th, and K concentrations from sites in the North Atlantic46, North Pacific53, and South Pacific34,52. The average U, Th, and K concentration values are consistent with data reported in Li and Schoonmaker55 for the characteristic U, Th, and K content found in abyssal clay and calcareous ooze. For lithogenous sediment, we use U, Th, and K concentration values reported for upper continental crust by Li and Schoonmaker55. Lithology-specific radioactive element values are given in Supplementary Table 6 and used to calculate Am,i in Eq. (1).
Density
Characteristic density values for calcite, quartz, terrigenous clay, and opal-rich sediment were extracted from the Proceedings of the Integrated Ocean Drilling Program Volume 320/321 and are assigned to calcareous ooze, lithogenous sediment, abyssal clay, and siliceous ooze, respectively79.
Global porosity
For global porosity, we use a seafloor porosity data set by Martin et al.58 and accounted for sediment compaction with depth by using separate sediment compaction length scales for continental-shelf (0–200 m water depth; c0 = 0.5 × 10−3), continental-margin (200–2500 m; c0 = 1.7 × 10−3), and abyssal sediment ( >3500 m; c0 = 0.85 × 10−3)80,81. Once the porosity was 0.1%, the depth integration was halted.
Global sediment thickness
We calculated global depth-integrated radiolytic H2 production by summing the seafloor production rates over sediment depth in one-meter intervals (Fig. 3 in main text). Sediment thickness is from Whittaker et al., supplemented with Laske and Masters where needed82,83.
Ocean depth
For porosity calculations, water depths were determined using the General Bathymetric Chart of the Oceans84, resampled to a 5-arc minute grid, i.e. the resolution of the Naval Oceanographic Office’s Bottom Sediment Type (BTS) database “Enhanced dataset”57.
Dissolved H2 concentration profiles
H2 concentrations from South Pacific Sites U1365, U1369, U1370, and U1371, and the measurement protocol, are described in ref. 1. H2 concentrations from North Atlantic KN223 Sites 11, 12, and 15, and North Pacific Site EQP11 were determined using the same protocol and are posted on SedDB (see “Data availability”). The detection limit for H2 ranged between 1 and 5 nM H2, depending on site, and is displayed as gray vertical lines in Fig. 2 of the main text. H2 concentrations for Equatorial Pacific Site 1225 and Peru Trench Site 1230 were measured by the “headspace equilibration technique”, which measures steady-state H2 levels reached following laboratory incubation of the sediment samples85,86.
For comparison to these measured H2 concentrations, we use diffusion-reaction calculations to quantify what in situ H2 concentrations would be in the absence of H2-consuming reactions. The results of these calculations are represented as solid circles (•) in Fig. 2 of the main text. Temporal changes in H2 concentration due to diffusive processes and radiolytic H2 production in situ are expressed by Eq. (2):

$$frac{{partial {mathrm{H}}_2(x,t)}}{{partial t}} = frac{D}{{varphi F}}frac{{partial ^2{mathrm{H}}_2(x,t)}}{{partial x^2}} + P(x)$$
(2)

with
D: the diffusion coefficient of H2(aq) at in situ temperature
(varphi): porosity
F: formation factor
x: depth
Z: sediment column thickness
({mathrm{H}}_2): hydrogen concentration
P: radiolytic H2 production rate
t : time.
With constant radiolytic H2 production, P(x) = P with depth,
and at steady-state,

$$frac{{partial ^2{mathrm{H}}_2(x)}}{{partial x^2}} = – frac{{Pvarphi F}}{D}.$$
(3)

We integrate Eq. (3) over the length x twice,

$${mathrm{H}}_2(x) = – frac{1}{2}frac{{Pvarphi F}}{D}x^2 + Ax + B$$
(4)

where A and B in Eq. (4) are constants of integration. We use two boundary conditions to derive the value of these constants.
Boundary condition 1: concentration of H2 at the sediment-water interface, x = 0, is zero due to diffusive loss to the overlying water column.
Boundary condition 2: concentration of H2 at the basement–sediment-water interface, x = Z, is zero due to diffusive loss to the underlying basement.
With these boundary conditions, (A = frac{1}{2}frac{{Pvarphi F}}{D}Z) and B = 0
and

$${mathrm{H}}_2(x) = frac{1}{2}frac{{Pvarphi F}}{D}(xZ – x^2).$$
(5)

In cases where we expect radiolytic H2 production rates to significantly vary with depth due to changes in lithology, we adapted the boundary conditions and applied a two-layer diffusion model to account for this variation.
Calculation of Gibbs Energies for the Knallgas reaction
For H2 concentrations above the detection limits at South Pacific IODP Expedition 329 sites (Supplementary Fig. 5)34, we quantified in situ Gibbs energies (ΔGr) of the Knallgas reaction (H2 + ½O2 → H2O). In situ ΔGr values depend on pressure (P), temperature (T), ionic strength, and chemical concentrations, all of which are explicitly accounted for in our calculations:

$$Delta G_{mathrm{r}} = Delta G^circ _{mathrm{r}}left( {T,P} right) + 2.3,RT,{mathrm{log}}_{10}Q$$
(6)

where:
ΔGr: in situ Gibbs energy of reaction (kJ mol H2−1)
ΔG°r(T,P): Gibbs energy of reaction under in situ T and P conditions (kJ mol H2−1)
R: gas constant (8.314 kJ−1 mol K−1)
Q: activity quotient of compounds involved in the reaction.
We use the measured composition of the sedimentary pore fluid to determine values of Q.
For a more complete overview of in situ Gibbs energy-of-reaction calculations in subseafloor sediment, see Wang et al.87.
Calculation of organic oxidation rates (net rates of O2 reduction and DIC production)
We calculated the vertical distribution of net O2 reduction rates at nine sites where the sediment is oxic from seafloor to basement and the vertical distribution of DIC production rates at seven sites where the subseafloor sediment is anoxic (see Supplementary Fig. 3 for site locations). Dissolved O2 concentrations are from Røy et al.47 and D’Hondt et al.88. DIC concentrations are from ODP Leg 201 (ref. 35), and the Proceedings of the IODP Expedition 323 (Sites U12343, U1345)54 and IODP Expedition 329 (Site U1371 (ref. 34)).
The net rates are calculated using the MatLab program and numerical procedures of Wang et al.89, modified by using an Akima spline, rather than a 5-point running mean, to generate a best-fit line to the chemical concentration data. Details of the calculation protocol for O2 production rates and DIC production rates are respectively described in the supplementary information of D’Hondt et al.88 and in Walsh et al.90. The DIC reaction rates and their first standard deviations calculated for the seven sites are given in Supplementary Table 7.
To facilitate comparisons of radiolytic chemical rates to net DIC production rates, rates are converted to electron equivalents (2 electrons per H2, 4 electrons per O2, 4 electrons per organic C oxidized).
Estimation of sediment ages
We estimated sediment ages for Sites U1343 and U1343 using the sediment-age model of Takahashi et al.54, which is based on biostratigraphic and magnetostratigraphic data. Because detailed chronostratigraphic data are not available for the remaining sites (Equatorial Pacific sites (1225 and 1226), Peru Trench Site 1230 and Peru Basin Site 1231, South Pacific sites U1365, U1366, U1367, U1369, U1370, and U1371, North Pacific sites EQP9 and EQP10, and North Atlantic sites KN223-11 and KN223-12), we used the mean sediment accumulation rate for each of these sites (Supplementary Fig. 3) to convert its sediment depth (in meters below seafloor) to sediment age (in millions of years, Ma). Mean sediment accumulation rate was calculated by dividing sediment thickness by basement age91 (Supplementary Table 8). For Sites 1225, 1226, 1230, 1231, U1365, U1366, U1367, U1369, U1370, and U1371, sediment thickness was determined by drilling to basement34,35. For Sites EQP9, EQP10, KN223-11, and KN223-12, sediment thicknesses were determined from acoustic basement reflection data. More

1.
Desiderà, E. et al. Acoustic fish communities: Sound diversity of rocky habitats reflects fish species diversity. Mar. Ecol. Prog. Ser. 608, 183–197 (2019).
ADS  Article  Google Scholar 
2.
Iorio, L. D., Gervaise, C., Lossent, J., Valentini-Poirier, C.-A. & Boissery, P. Benthic biophonic assemblages, their environmental divers, eco-acoustic scores at the level of the Western Mediterranean basin, and their implications for large-scale ecosystem monitoring. J. Acoust. Soc. Am. 144, 1692–1692 (2018).
ADS  Article  Google Scholar 

3.
Buscaino, G. et al. Temporal patterns in the soundscape of the shallow waters of a Mediterranean marine protected area. Sci. Rep. 6, 1–13 (2016).
Article  CAS  Google Scholar 

4.
McNett, G. D., Luan, L. H. & Cocroft, R. B. Wind-induced noise alters signaler and receiver behavior in vibrational communication. Behav. Ecol. Sociobiol. 64, 2043–2051 (2010).
Article  Google Scholar 

5.
Hildebrand, J. A. Anthropogenic and natural sources of ambient noise in the ocean. Mar. Ecol. Prog. Ser. 395, 5–20 (2009).
ADS  Article  Google Scholar 

6.
Popper, A. N. & Hawkins, A. D. An overview of fish bioacoustics and the impacts of anthropogenic sounds on fishes. J. Fish Biol. 94, 692–713 (2019).
PubMed  PubMed Central  Article  Google Scholar 

7.
Belkovich, V. M., Bibikov, N. G., Dubrovsky, N. A., Suhoruchenko, M. N. & Zhuravlev, V. A. Preliminary estimates of low-frequency sound effect on sea animals in the Eastern Arctic. (1994).

8.
Mate, B. R., Stafford, K. M. & Ljungblad, D. K. A change in sperm whale (Physeter macroephalus) distribution correlated to seismic surveys in the Gulf of Mexico. J. Acoust. Soc. Am. 96, 3268–3269 (1994).
ADS  Article  Google Scholar 

9.
Dunlop, R. A. et al. The behavioural response of humpback whales (Megaptera novaeangliae) to a 20 cubic inch air gun. Aquat. Mamm. 41, 412–433 (2015).
Article  Google Scholar 

10.
Kunc, H. P., McLaughlin, K. E. & Schmidt, R. Aquatic noise pollution: Implications for individuals, populations, and ecosystems. Proc. R. Soc. B Biol. Sci. 283, 20160839 (2016).
Article  Google Scholar 

11.
Rako-Gospić, N. & Picciulin, M. Underwater noise: Sources and effects on marine life. in World Seas: An Environmental Evaluation Volume III: Ecological Issues and Environmental Impacts (ed. Sheppard, C.) vol. 3 367–389 (Elsevier, 2019).

12.
Duarte, C. M., Chapuis, L., Collin, S. P., Costa, D. P., Devassy, R. P., Eguiluz, V. M., Erbe, C., Gordon, T. A. C., Halpern, B. S., Harding, H. R., Havlik, M. N., Meekan, M., Merchant, N. D., Miksis-Olds, J. L., Parsons, M., Predragovic, M., Radford, A. N., Radford, C. A., Simpson, S. D., Slabbekoorn, H., Staaterman, E., Van Opzeeland, I. C., Winderen, J., Zhang, X. & Juanes, F. The soundscape of the Anthropocene ocean. Science 371(6529), eaba4658 (2021).
CAS  PubMed  Article  Google Scholar 

13.
Cox, K., Brennan, L. P., Gerwing, T. G., Dudas, S. E. & Juanes, F. Sound the alarm: A meta-analysis on the effect of aquatic noise on fish behavior and physiology. Glob. Chang. Biol. 24, 3105–3116 (2018).
ADS  PubMed  Article  Google Scholar 

14.
La Manna, G., Manghi, M., Perretti, F. & Sarà, G. Behavioral response of brown meagre (Sciaena umbra) to boat noise. Mar. Pollut. Bull. 110, 324–334 (2016).
PubMed  Article  CAS  Google Scholar 

15.
Nedelec, S. L. et al. Motorboat noise impacts parental behaviour and offspring survival in a reef fish. Proc. R. Soc. B Biol. Sci. 284, 20170143 (2017).
Article  Google Scholar 

16.
Maxwell, R. J. et al. Does motor noise from recreational boats alter parental care behaviour of a nesting freshwater fish?. Aquat. Conserv. Mar. Freshw. Ecosyst. 28, 969–978 (2018).
Article  Google Scholar 

17.
Picciulin, M., Sebastianutto, L., Codarin, A., Calcagno, G. & Ferrero, E. A. Brown meagre vocalization rate increases during repetitive boat noise exposures: A possible case of vocal compensation. J. Acoust. Soc. Am. 132, 3118–3124 (2012).
ADS  PubMed  Article  Google Scholar 

18.
de Jong, K., Amorim, M. C. P., Fonseca, P. J., Klein, A. & Heubel, K. U. Noise affects acoustic courtship behavior similarly in two species of gobies. Proc. Meet. Acoust. 27, 010018 (2016).
Article  Google Scholar 

19.
Nedelec, S. L. et al. Motorboat noise disrupts co-operative interspecific interactions. Sci. Rep. 7, 1–8 (2017).
ADS  CAS  Article  Google Scholar 

20.
Codarin, A., Wysocki, L. E., Ladich, F. & Picciulin, M. Effects of ambient and boat noise on hearing and communication in three fish species living in a marine protected area (Miramare, Italy). Mar. Pollut. Bull. 58, 1880–1887 (2009).
CAS  PubMed  Article  Google Scholar 

21.
Holles, S., Simpson, S. D., Radford, A. N., Berten, L. & Lecchini, D. Boat noise disrupts orientation behaviour in a coral reef fish. Mar. Ecol. Prog. Ser. 485, 295–300 (2013).
ADS  Article  Google Scholar 

22.
Pine, M. K., Jeffs, A. G., Wang, D. & Radford, C. A. The potential for vessel noise to mask biologically important sounds within ecologically significant embayments. Ocean Coast. Manag. 127, 63–73 (2016).
Article  Google Scholar 

23.
de Jong, K., Amorim, M. C. P., Fonseca, P. J., Fox, C. J. & Heubel, K. U. Noise can affect acoustic communication and subsequent spawning success in fish. Environ. Pollut. 237, 814–823 (2018).
PubMed  Article  CAS  Google Scholar 

24.
Montgomery, J. C., Jeffs, A., Simpson, S. D., Meekan, M. & Tindle, C. Sound as an orientation cue for the pelagic larvae of reef fishes and decapod crustaceans. Adv. Mar. Biol. 51, 143–196 (2006).
PubMed  Article  Google Scholar 

25.
Hussey, N. E. et al. Aquatic animal telemetry: A panoramic window into the underwater world. Science 348, 1255642 (2015).
PubMed  Article  CAS  Google Scholar 

26.
Simpfendorfer, C. A., Heupel, M. R. & Collins, A. B. Variation in the performance of acoustic receivers and its implication for positioning algorithms in a riverine setting. Can. J. Fish. Aquat. Sci. 65, 482–492 (2008).
Article  Google Scholar 

27.
Abecasis, D. et al. A review of acoustic telemetry in Europe and the need for a regional aquatic telemetry network. Anim. Biotelemetry 6, 12 (2018).
Article  Google Scholar 

28.
Kessel, S. T. et al. A review of detection range testing in aquatic passive acoustic telemetry studies. Rev. Fish Biol. Fish. 24, 199–218 (2014).
Article  Google Scholar 

29.
Patterson, T. A., Thomas, L., Wilcox, C., Ovaskainen, O. & Matthiopoulos, J. State-space models of individual animal movement. Trends Ecol. Evol. 23, 87–94 (2008).
PubMed  Article  Google Scholar 

30.
Alós, J., Palmer, M., Balle, S. & Arlinghaus, R. Bayesian state-space modelling of conventional acoustic tracking provides accurate descriptors of home range behavior in a small-bodied coastal fish species. PLoS ONE 11, 1–23 (2016).
Article  CAS  Google Scholar 

31.
Jonsen, I. D. et al. State-space models for bio-loggers: A methodological road map. Deep Sea Res. Top. Stud. Oceanogr. 88, 34–46 (2013).
ADS  Article  Google Scholar 

32.
Bauchot, M. . Serranidae. in Guide FAO d’identification des espèces pour les besoins de la pêche (eds. Fisher, W., Bauchot, M. L. & Scheneider, M.) 1301–1329 (Organisation des Nations Unies pour l’Alimentation et l’Agriculture, 1987).

33.
Börger, L., Dalziel, B. D. & Fryxell, J. M. Are there general mechanisms of animal home range behaviour? A review and prospects for future research. Ecol. Lett. 11, 637–650 (2008).
PubMed  Article  Google Scholar 

34.
Alós, J. et al. Selective exploitation of spatially structured coastal fish populations by recreational anglers may lead to evolutionary downsizing of adults. Mar. Ecol. Prog. Ser. 503, 219–233 (2014).
ADS  Article  Google Scholar 

35.
Palmer, M., Balle, S., March, D., Alós, J. & Linde, M. Size estimation of circular home range from fish mark-release-(single)-recapture data: Case study of a small labrid targeted by recreational fishing. Mar. Ecol. Prog. Ser. 430, 87–97 (2011).
ADS  Article  Google Scholar 

36.
Mills, S. C. et al. Hormonal and behavioural effects of motorboat noise on wild coral reef fish. Environ. Pollut. 262, 114250 (2020).
CAS  PubMed  Article  Google Scholar 

37.
Hitt, S., Pittman, S. J. & Nemeth, R. S. Diel movements of fishes linked to benthic seascape structure in a Caribbean coral reef ecosystem. Mar. Ecol. Prog. Ser. 427, 275–291 (2011).
ADS  Article  Google Scholar 

38.
March, D., Palmer, M., Alós, J., Grau, A. & Cardona, F. Short-term residence, home range size and diel patterns of the painted comber Serranus scriba in a temperate marine reserve. Mar. Ecol. Prog. Ser. 400, 195–206 (2010).
ADS  Article  Google Scholar 

39.
Codling, E. A., Plank, M. J. & Benhamou, S. Random walk models in biology. J. R. Soc. Interface 5, 813–834 (2008).
PubMed  PubMed Central  Article  Google Scholar 

40.
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  Google Scholar 

41.
Gardiner, C. W. Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences (Springer, Berlin, 1990).
Google Scholar 

42.
Follana-Berná, G., Palmer, M., Lekanda-Guarrotxena, A., Grau, A. & Arechavala-Lopez, P. Fish density estimation using unbaited cameras: Accounting for environmental-dependent detectability. J. Exp. Mar. Bio. Ecol. 527, 151376 (2020).
Article  Google Scholar 

43.
Klimley, A. P., Voegeli, F., Beavers, S. C. & Le Boeuf, B. J. Automated listening stations for tagged marine fishes. Mar. Technol. Soc. J. 32, 94–101 (1998).
Google Scholar 

44.
Hedger, R. D. et al. The optimized interpolation of fish positions and speeds in an array of fixed acoustic receivers. ICES J. Mar. Sci. 65, 1248–1259 (2008).
Article  Google Scholar 

45.
Merchant, N. D., Blondel, P., Dakin, D. T. & Dorocicz, J. Averaging underwater noise levels for environmental assessment of shipping. J. Acoust. Soc. Am. 132, 343–349 (2012).
ADS  Article  Google Scholar 

46.
Plummer, M. rjags: Bayesian graphical models using MCMC. (2016).

47.
Su, Y. S. & Yajima, M. R2jags: Using R to run ‘JAGS’. (2015).

48.
R Development Core Team. R: A language and environment for statistical computing. (2017).

49.
Plummer, M., Best, N., Cowles, K. & Vines, K. CODA: convergence diagnosis and output analysis for MCMC. R News 6, 7–11 (2006).
Google Scholar 

50.
Gelman, A. et al. Bayesian data analysis. (Taylor & Francis Group, 2013).

51.
Fagan, W. F. & Calabrese, J. M. The correlated random walk and the rise of movement ecology. Bull. Ecol. Soc. Am. 95, 204–206 (2014).
Article  Google Scholar 

52.
Payne, N. L., Gillanders, B. M., Webber, D. M. & Semmens, J. M. Interpreting diel activity patterns from acoustic telemetry: The need for controls. Mar. Ecol. Prog. Ser. 419, 295–301 (2010).
ADS  Article  Google Scholar 

53.
Burger, J. & Gochfeld, M. On developing bioindicators for human and ecological health. Environ. Monit. Assess. 66, 23–46 (2001).
CAS  PubMed  Article  Google Scholar 

54.
Alós, J., March, D., Palmer, M., Grau, A. & Morales-Nin, B. Spatial and temporal patterns in Serranus cabrilla habitat use in the NW Mediterranean revealed by acoustic telemetry. Mar. Ecol. Prog. Ser. 427, 173–186 (2011).
ADS  Article  Google Scholar 

55.
Alós, J., Palmer, M., Rosselló, R. & Arlinghaus, R. Fast and behavior-selective exploitation of a marine fish targeted by anglers. Sci. Rep. 6, 1–13 (2016).
Article  CAS  Google Scholar 

56.
Campioni, L. et al. Individual and spatio-temporal variations in the home range behaviour of a long-lived, territorial species. Oecologia 172, 371–385 (2013).
ADS  PubMed  Article  Google Scholar 

57.
Topping, D. T. & Szedlmayer, S. T. Home range and movement patterns of red snapper (Lutjanus campechanus) on artificial reefs. Fish. Res. 112, 77–84 (2011).
Article  Google Scholar 

58.
Jorgensen, S. J. et al. Limited movement in blue rockfish Sebastes mystinus: Internal structure of home range. Mar. Ecol. Prog. Ser. 327, 157–170 (2006).
ADS  Article  Google Scholar 

59.
Harding, H. R. et al. Fish in habitats with higher motorboat disturbance show reduced sensitivity to motorboat noise. Biol. Lett. 14, 20180441 (2018).
PubMed  PubMed Central  Article  Google Scholar 

60.
Holmes, L. J., McWilliam, J., Ferrari, M. C. O. & McCormick, M. I. Juvenile damselfish are affected but desensitize to small motor boat noise. J. Exp. Mar. Bio. Ecol. 494, 63–68 (2017).
Article  Google Scholar 

61.
Huntingford, F. A. et al. Current issues in fish welfare. J. Fish Biol. 68, 332–372 (2006).
Article  Google Scholar 

62.
Sierra-Flores, R., Atack, T., Migaud, H. & Davie, A. Stress response to anthropogenic noise in Atlantic cod Gadus morhua L. Aquac. Eng. 67, 67–76 (2015).
Article  Google Scholar 

63.
de Jong, K. et al. Predicting the effects of anthropogenic noise on fish reproduction. Rev. Fish Biol. Fish. 30, 1–24 (2020).
Article  Google Scholar 

64.
McCormick, M. I., Fakan, E. P., Nedelec, S. L. & Allan, B. J. M. Effects of boat noise on fish fast-start escape response depend on engine type. Sci. Rep. 9, 1–10 (2019).
Article  CAS  Google Scholar 

65.
Hawkins, A. D. & Popper, A. N. A sound approach to assessing the impact of underwater noise on marine fishes and invertebrates. ICES J. Mar. Sci. 74, 635–651 (2017).
Article  Google Scholar 

66.
Soudijn, F. H., van Kooten, T., Slabbekoorn, H. & de Roos, A. M. Population-level effects of acoustic disturbance in Atlantic cod: a size-structured analysis based on energy budgets. Proc. Royal Soc. B. Sci. 287(1929), 20200490 (2020).
Article  Google Scholar  More

1.
Bogdanov, S., Ruoff, K. & Persano Oddo, L. Physico-chemical methods for the characterisation of unifloral honeys: a review. Apidologie 35, S4–S17 (2004).
Article  Google Scholar 
2.
Kwakman, P. H. S., te Velde, A. A., de Boer, L., Vandenbroucke-Grauls, C. M. J. E. & Zaat, S. A. J. Two major medicinal honeys have different mechanisms of bactericidal activity. PLoS ONE 6, e17709 (2011).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

3.
Lu, J. et al. The effect of New Zealand Kanuka, Manuka and Clover Honeys on bacterial growth dynamics and cellular morphology varies according to the species. PLoS ONE 8, e55898 (2013).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

4.
Salonen, A., Ollikka, T., Grönlund, E., Ruottinen, L. & Julkunen-Tiitto, R. Pollen analyses of honey from Finland. Grana 48, 281–289 (2009).
Article  Google Scholar 

5.
Balkanska, R., Stefanova, K. & Stoikova-Grigorova, R. Main honey botanical components and techniques for identification: a review. J. Apic. Res. https://doi.org/10.1080/00218839.2020.1765481 (2020).
Article  Google Scholar 

6.
Soares, S., Amaral, J. S., Oliveira, M. B. P. P. & Mafra, I. A comprehensive review on the main honey authentication issues: production and origin. Compr. Rev. Food Sci. Food Saf. 16, 1072–1100 (2017).
CAS  PubMed  Article  Google Scholar 

7.
Beckmann, K., Beckh, G., Luellmann, C. & Speer, K. Characterization of filtered honey by electrophoresis of enzyme fractions. Apidologie 42, 59–66 (2011).
CAS  Article  Google Scholar 

8.
Anklam, E. A review of the analytical methods to determine the geographical and botanical origin of honey. Food Chem. 63, 549–562 (1998).
CAS  Article  Google Scholar 

9.
Von Der Ohe, W., Persano Oddo, L., Piana, M. L., Morlot, M. & Martin, P. Harmonized methods of melissopalynology. Apidologie 35, 18–25 (2004).
Article  Google Scholar 

10.
Bell, K. L. et al. Pollen DNA barcoding: Current applications and future prospects. Genome 59, 629–640 (2016).
PubMed  Article  PubMed Central  Google Scholar 

11.
Guertler, P., Eicheldinger, A., Muschler, P., Goerlich, O. & Busch, U. Automated DNA extraction from pollen in honey. Food Chem. 149, 302–306 (2014).
CAS  PubMed  Article  Google Scholar 

12.
Hawkins, J. et al. Using DNA metabarcoding to identify the floral composition of honey: A new tool for investigating honey bee foraging preferences. PLoS ONE 10, e0134735 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

13.
Valentini, A., Miquel, C. & Taberlet, P. DNA barcoding for honey biodiversity. Diversity 2, 610–617 (2010).
CAS  Article  Google Scholar 

14.
Prosser, S. W. J. & Hebert, P. D. N. Rapid identification of the botanical and entomological sources of honey using DNA metabarcoding. Food Chem. 214, 183–191 (2017).
CAS  PubMed  Article  Google Scholar 

15.
Olivieri, C., Marota, I., Rollo, F. & Luciani, S. Tracking plant, fungal, and bacterial DNA in honey specimens. J. Forensic Sci. 57, 222–227 (2012).
CAS  PubMed  Article  Google Scholar 

16.
Snowdon, J. A. & Cliver, D. O. Microorganisms in honey. Int. J. Food Microbiol. 31, 1–26 (1996).
CAS  PubMed  Article  Google Scholar 

17.
Manirajan, B. A. et al. Diversity, specificity, co-occurrence and hub taxa of the bacterial-fungal pollen microbiome. FEMS Microbiol. Ecol. 94, 1–11 (2018).
Article  CAS  Google Scholar 

18.
Anderson, K. E. et al. Microbial ecology of the hive and pollination landscape: Bacterial associates from floral nectar, the alimentary tract and stored food of honey bees (Apis mellifera). PLoS ONE 8, e83125 (2013).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

19.
Aizenberg-Gershtein, Y., Izhaki, I. & Halpern, M. Do honeybees shape the bacterial community composition in floral nectar?. PLoS ONE 8, e83125 (2013).
ADS  Article  CAS  Google Scholar 

20.
Fridman, S., Izhaki, I., Gerchman, Y. & Halpern, M. Bacterial communities in floral nectar. Environ. Microbiol. Rep. 4, 97–104 (2012).
PubMed  Article  PubMed Central  Google Scholar 

21.
Nevas, M. et al. High prevalence of Clostridium botulinum types A and B in honey samples detected by polymerase chain reaction. Int. J. Food Microbiol. 72, 45–52 (2002).
CAS  PubMed  Article  PubMed Central  Google Scholar 

22.
Bonilla-Rosso, G. & Engel, P. Functional roles and metabolic niches in the honey bee gut microbiota. Curr. Opin. Microbiol. 43, 69–76 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

23.
Engel, P., Martinson, V. G. & Moran, N. A. Functional diversity within the simple gut microbiota of the honey bee. Proc. Natl. Acad. Sci. U. S. A. 109, 11002–11007 (2012).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

24.
Oksanen, J. et al. Package ‘vegan’ Title Community Ecology Package Version 2.5-6. (2019).

25.
Larsson, J. Area-Proportional Euler and Venn Diagrams with Ellipses [R package eulerr version 6.1.0].

26.
Warton, D. I. et al. So many variables: Joint modeling in community ecology. Trends Ecol. Evol. 30, 766–779 (2015).
PubMed  Article  Google Scholar 

27.
Tjur, T. Coefficients of determination in logistic regression models—A new proposal: The coefficient of discrimination. Am. Stat. 63, 366–372 (2009).
MathSciNet  MATH  Article  Google Scholar 

28.
Fielding, A. H. & Bell, J. F. A review of methods for the assessment of prediction errors in conservation presence/absence models (2020).https://doi.org/10.1017/S0376892997000088

29.
Guisan, A. et al. Measuring model accuracy: Which metrics to use? in Habitat Suitability and Distribution Models 241–269 (Cambridge University Press, 2017). doi:https://doi.org/10.1017/9781139028271.022.

30.
Tikhonov, G. et al. Joint species distribution modelling with the r-package HMSC. Methods Ecol. Evol. 11, 442–447 (2020).
PubMed  PubMed Central  Article  Google Scholar 

31.
Moran, N. A., Hansen, A. K., Powell, J. E. & Sabree, Z. L. Distinctive gut microbiota of honey bees assessed using deep sampling from individual worker bees. PLoS ONE 7, 1–10 (2012).
Google Scholar 

32.
Fünfhaus, A., Ebeling, J. & Genersch, E. Bacterial pathogens of bees. Curr. Opin. Insect Sci. 26, 89–96 (2018).
PubMed  Article  Google Scholar 

33.
Fries, I. Nosema ceranae in European honey bees (Apis mellifera). J. Invertebr. Pathol. 103, (2010).

34.
Balvočiute, M. & Huson, D. H. SILVA, RDP, Greengenes, NCBI and OTT—how do these taxonomies compare?. BMC Genomics 18, 114 (2017).
PubMed  PubMed Central  Article  Google Scholar 

35.
Meiklejohn, K. A., Damaso, N. & Robertson, J. M. Assessment of BOLD and GenBank—their accuracy and reliability for the identification of biological materials. PLoS One 14 (2019).

36.
Nilsson, R. H. et al. The UNITE database for molecular identification of fungi: handling dark taxa and parallel taxonomic classifications. Nucleic Acids Res. 47, 259–264 (2018).
Article  CAS  Google Scholar 

37.
Sickel, W. et al. Increased efficiency in identifying mixed pollen samples by meta-barcoding with a dual-indexing approach. BMC Ecol. 15, 1–9 (2015).
Article  CAS  Google Scholar 

38.
Cole, J. R. et al. The Ribosomal Database Project: Improved alignments and new tools for rRNA analysis. Nucleic Acids Res. 37, 141–145 (2008).
Article  CAS  Google Scholar 

39.
Bell, K. L., Loeffler, V. M. & Brosi, B. J. An rbcL reference library to aid in the identification of plant species mixtures by DNA metabarcoding. Appl. Plant Sci. 5, 1600110 (2017).
Article  Google Scholar 

40.
Edgar, R. C. & Flyvbjerg, H. Error filtering, pair assembly and error correction for next-generation sequencing reads. Bioinformatics 31, 3476–3482 (2015).
CAS  PubMed  Article  Google Scholar 

41.
Callahan, B. J., McMurdie, P. J. & Holmes, S. P. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. ISME J. 11, 2639–2643 (2017).
PubMed  PubMed Central  Article  Google Scholar 

42.
Rognes, T., Flouri, T., Nichols, B., Quince, C. & Mahé, F. VSEARCH: a versatile open source tool for metagenomics. PeerJ https://doi.org/10.7717/peerj.2584 (2016).
Article  PubMed  PubMed Central  Google Scholar 

43.
Vesterinen, E. J., Kaunisto, K. M. & Lilley, T. M. A global class reunion with multiple groups feasting on the declining insect smorgasbord. Sci. Rep. 10, 16595 (2020).
CAS  PubMed  PubMed Central  Article  Google Scholar 

44.
Barcaccia, G., Lucchin, M. & Cassandro, M. DNA barcoding as a molecular tool to track down mislabeling and food piracy. Diversity 8, 2 (2015).
Article  CAS  Google Scholar 

45.
Zábrodská, B. & Vorlová, L. Adulteration of honey and available methods for detection—a review. Acta Vet. Brno 83, S85–S102 (2014).
Article  Google Scholar 

46.
Hebert, P. D. N., Cywinska, A., Ball, S. L. & DeWaard, J. R. Biological identifications through DNA barcodes. Proc. R. Soc. B Biol. Sci. 270, 313–321 (2003).
CAS  Article  Google Scholar 

47.
DeSalle, R. & Goldstein, P. Review and Interpretation of Trends in DNA Barcoding. Front. Ecol. Evol. 7, 302 (2019).
Article  Google Scholar 

48.
Hawkins, J. et al. Using DNA metabarcoding to identify the floral composition of honey: A new tool for investigating honey bee foraging preferences. PLoS One 10 (2015).

49.
De Vere, N. et al. Using DNA metabarcoding to investigate honey bee foraging reveals limited flower use despite high floral availability. Sci. Rep. 7, 1–10 (2017).
Article  CAS  Google Scholar 

50.
Lucek, K. et al. Metabarcoding of honey to assess differences in plant-pollinator interactions between urban and non-urban sites. https://doi.org/10.1007/s13592-019-00646-3.

51.
Bruni, I. et al. A DNA barcoding approach to identify plant species in multiflower honey. Food Chem. 170, 308–315 (2015).
CAS  PubMed  Article  Google Scholar 

52.
Laha, R. C. et al. Meta-barcoding in combination with palynological inference is a potent diagnostic marker for honey floral composition. AMB Express 7, 132 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

53.
Utzeri, V. J., Ribani, A. & Fontanesi, L. Authentication of honey based on a DNA method to differentiate Apis mellifera subspecies: Application to Sicilian honey bee (A. m. siciliana) and Iberian honey bee (A. m. iberiensis) honeys. Food Control 91, 294–301 (2018).
CAS  Article  Google Scholar 

54.
Bovo, S. et al. Shotgun metagenomics of honey DNA: Evaluation of a methodological approach to describe a multi-kingdom honey bee derived environmental DNA signature. PLoS ONE 13, e0205575 (2018).
PubMed  PubMed Central  Article  CAS  Google Scholar 

55.
Bovo, S., Utzeri, V. J., Ribani, A. & Cabbri, R. Shotgun sequencing of honey DnA can describe honey bee derived environmental signatures and the honey bee hologenome complexity. https://doi.org/10.1038/s41598-020-66127-1.

56.
Vesterinen, E. J., Puisto, A. I. E., Blomberg, A. S. & Lilley, T. M. Table for five, please: Dietary partitioning in boreal bats. Ecol. Evol. 8, 10914–10937 (2018).
PubMed  PubMed Central  Article  Google Scholar 

57.
Vesterinen, E. J. et al. What you need is what you eat? Prey selection by the bat Myotis daubentonii. Mol. Ecol. 25, 1581–1594 (2016).
CAS  PubMed  Article  Google Scholar 

58.
Functional Genomics Unit, University of Helsinki, Finland. www.helsinki.fi/en/infrastructures/genome-analysis/biomedicum-functional-genomics-unit.

59.
Zhang, J., Kobert, K., Flouri, T. & Stamatakis, A. PEAR: A fast and accurate Illumina Paired-End reAd mergeR. Bioinformatics 30, 614–620 (2014).
CAS  PubMed  Article  Google Scholar 

60.
Schmieder, R. & Edwards, R. Quality control and preprocessing of metagenomic datasets. Bioinformatics 27, 863–864 (2011).
CAS  PubMed  PubMed Central  Article  Google Scholar 

61.
Edgar, R. C. UPARSE: Highly accurate OTU sequences from microbial amplicon reads. Nat. Methods 10, 996–998 (2013).
CAS  PubMed  PubMed Central  Article  Google Scholar 

62.
Wang, Q., Garrity, G. M., Tiedje, J. M. & Cole, J. R. Naïve Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl. Environ. Microbiol. 73, 5261–5267 (2007).
CAS  PubMed  PubMed Central  Article  Google Scholar 

63.
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 

64.
Benson, D. A., Karsch-Mizrachi, I., Lipman, D. J., Ostell, J. & Sayers, E. W. GenBank. Nucleic Acids Res. 39, D28–D31 (2011).
Article  CAS  Google Scholar 

65.
Huson, D. H., Auch, A. F., Qi, J. & Schuster, S. C. MEGAN analysis of metagenomic data. Genome Res. 17, 377–386 (2007).
CAS  PubMed  PubMed Central  Article  Google Scholar 

66.
Lee, T., Alemseged, Y. & Mitchell, A. Dropping Hints: Estimating the diets of livestock in rangelands using DNA metabarcoding of faeces. Metabarcoding Metagenomics 2, e22467 (2018).
Article  Google Scholar 

67.
Alberdi, A., Garin, I., Aizpurua, O. & Aihartza, J. The foraging ecology of the Mountain Long-eared bat Plecotus macrobullaris revealed with DNA mini-barcodes. PLoS One 7, (2012).

68.
Bolger, A. M., Lohse, M. & Usadel, B. Genome analysis Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).
CAS  PubMed  PubMed Central  Article  Google Scholar 

69.
Wood, D. E., Lu, J. & Langmead, B. Improved metagenomic analysis with Kraken 2. Genome Biol. 20, 1–13 (2019).
Article  CAS  Google Scholar 

70.
National Center for Biotechnology Information (NCBI); Bethesda (MD): National Library of Medicine (US). https://www.ncbi.nlm.nih.gov/ (1988).

71.
Lu, J., Breitwieser, F. P., Thielen, P. & Salzberg, S. L. Bracken: Estimating species abundance in metagenomics data. PeerJ Comput. Sci. 2017, e104 (2017).
Article  Google Scholar 

72.
Breitwieser, F. P. & Salzberg, S. L. Pavian: Interactive analysis of metagenomics data for microbiome studies and pathogen identification. Bioinformatics 36, 1303–1304 (2020).
CAS  PubMed  Article  Google Scholar 

73.
DIN (Deutsches Institut für Normung),. Untersuchung von Honig – Bestimmung der relativen Pollenhäufigkeit. DIN 10760, 2002–2005 (2002).
Google Scholar 

74.
Persano Oddo, L. et al. Main European unifloral honeys: descriptive sheets 1. Apidologie 35, 38–81 (2004).
Article  Google Scholar 

75.
Piper, A. M. et al. Prospects and challenges of implementing DNA metabarcoding for high-throughput insect surveillance. Gigascience 8, 1–22 (2019).
Article  Google Scholar 

76.
Ovaskainen, O. & Abrego, N. Joint species distribution modelling joint species distribution modelling (Cambridge University Press, Cambridge, 2020). https://doi.org/10.1017/9781108591720.
Google Scholar  More

Figure 2

Schematic figure of the labeling process. Participants usually upload several images in a single report. The best photo is picked by the validator who first marks the harassing or non-appropriate photos as hidden. All the non-best photos are marked as not classified. In some rare events, two or three images are annotated from the same report. The mosquito images are classified into four different categories (Aedes albopictus, Aedes aegypti, other species or can not tell) and also the confidence of the label is marked as probable or confirmed. In this paper we excluded the not classified, the hidden and the can not tell images.

Full size image

Between 2014 and 2019, 7686 citizen-made mosquito photos were labeled through Mosquito Alert by entomology experts, with labels indicating whether Ae. albopictus appear in the photos. The photos were included in reports that Mosquito Alert participants uploaded, and each report could contain several photos, see Fig. 2. The entomology experts usually labeled the best photo of the report, but sometimes they labeled two (420 times) or three (49 times) for a single report, meaning that the dataset consisted of 7168 reports. For 6699 reports, only one image was labeled by the experts; for 420 reports two were labeled; for 49 reports three were labeled. Although these reports usually contain several photos, only the ones with expert labels were used in the analysis, as cannot be assumed that all of the photos in a report would have been given the same label.
The main goals of Mosquito Alert during this 6 year period were to monitor Ae. albopictus spreading and provide early detection of Ae. aegypti in Spain. Although people participate in Mosquito Alert all over the world, the majority of the participants and the majority of the photos are in Spain (see Fig. 1). As Ae. aegypti has not been reported in Spain in recent times, most Mosquito Alert participants lived in areas where Ae. aegypti is not present, so most of the photos are of Ae. albopictus. For the detailed yearly distribution of the photos, see Table 1.
Table 1 The collected and expert validated dataset for the period 2014–2019.
Full size table

A popular deep learning model, ResNet5026 was trained and evaluated on the collected dataset with yearly cross-validation. ResNet50 was used because of its wide popularity and its proven classification power in various datasets. As presenting infinitesimal increments of the classification power is not a goal of this paper, we do not report various ImageNet state-of-the-art model performances. Yearly cross-validation was used to rule out any possibility of information leakage (possibility of a user submitting multiple reports for the same mosquito).
The trained model is not only capable of generating highly accurate predictions, but it can also ease the human annotator workload by auto-marking the images where the neural network is confident and more accurate, leaving more uncertain cases for the entomology experts. Moreover, while visualizing the erroneous predictions a few re-occurring patterns were identified, which can serve as a proposal for how to make images that can be best processed by the model.
Several aspects of the dataset were explored as follows.
Classification
Since Mosquito Alert was centered around Ae. albopictus during the relevant time period (2014–2019), the collected dataset is biased towards this species (Table 1). We explored training classifiers on the Mosquito Alert dataset alone and also tied training on a balanced dataset, where 3896 negative samples were added from the IP10227 dataset of various non-mosquito insects as negative samples. From the IP102 dataset, images similar to mosquitoes, and images of striped insects were selected. Although the presented mosquito alert dataset is filtered to contain only mosquito images, in later use, non-mosquito images might be uploaded by the citizens. Training the CNN on a combination of mosquito and non-mosquito images can improve the model to make correct predictions, classifying non-tiger mosquitoes for those cases too. For testing, in each fold, only the Mosquito Alert dataset was used.
The trained classifiers achieved an extremely high area under the receiver operating characteristic curve (ROC AUC) score of 0.96 (see Fig. 3). The fact that the ROC AUC score for each fold was always over 0.95 proves the consistency of our classifier. Inspecting the confusion matrix shows us that the model tends to make more false positive predictions (assuming tiger mosquito is defined as the positive outcome) than false negatives, resulting in high sensitivity. The augmentation of the Mosquito Alert dataset with various insects from IP102 images to make it more balanced resulted in a slight performance boost and narrowed the gap between the number of false positive and false negative samples as expected, see Table 2.
Figure 3

Left: ROC curve calculated on the prediction of the 7686 images in the Mosquito Alert dataset with yearly cross-validation. The blue line shows the case when only the Mosquito Alert dataset was used for training, the orange when the training dataset was balanced out with the addition of non-tiger mosquito insect images from the IP102 dataset. Also a zoom into the part of the ROC curve, where the two methods differ the most is highlighted. Right: the confusion matrix was calculated on the same predictions when only the Mosquito Alert dataset was used for training. For both, a positive label means tiger mosquito is present.

Full size image

Table 2 Yearly cross-validation results with using the Mosquito Alert dataset alone and its IP102 augmented version.
Full size table

How to take a good picture?
Inspection of the weaknesses of a machine learning model is a fruitful way to gain a deeper understanding of the underlying problems and mechanisms. In our case, a careful review of the mispredicted images led us to useful insights into what makes a photo hard to classify for the deep learning model. On Fig. 8, a few selected examples are presented. Unlike humans, deep learning models rely more on textures than on shapes28. As a consequence, grid-like background patterns or striped objects may easily confuse the machine classifier. A larger rich training set can help to avoid these pitfalls, but we also have the option to advise the participants. If participants avoid confusing setups when taking photos, this can improve the accuracy of the automated classification. These guidelines can be added to the Mosquito Alert application to help participants make good images of mosquitoes.

Do not use striped structure (e.g. mosquito net or fly-flap) as a background.

Avoid complex backgrounds when possible. A few examples: patterned carpet, different nets, reflecting/shiny background, bumpy wallpaper.

Use clear, white background (e.g. a sheet of plain paper is perfect if possible) or hold the mosquito with finger pads.

Make sure that as much as possible the mosquito is in focus and covers a large area of the photo.

In general, it is desirable to have a clean white background with the mosquito centered, and with the image containing as little background as possible.
Dataset size impact on model performance
Modern deep CNNs tend to generate better predictions when trained on larger datasets. In this experiment, we trained a ResNet50 model on 10–20–(cdots )–90–100% of 6686 images and evaluated the model on the remaining 1000 images. The 1000 images were selected from the same year (2019) and all of them came from reports with only one photo. There were 709 tiger mosquitoes out of the 1000 test images. ROC AUC and accuracy were calculated with a 500 round bootstrapping of the 1000 test images.
Figure 4

Training a ResNet50 model on a subsampled training dataset. The model was tested against the same 1000 test images for all the steps and statistics of the test metric was calculated with a 500 round bootstrapping. The curve proves the diversity of the Mosquito Alert dataset and also suggests that in the future when the dataset will be even larger, the classification performance will increase.

Full size image

The mean and the standard deviation of the 500 rounds are shown in Fig. 4 for each training data size. From the figure, we can conclude that the predictive power of the model increases as more data are used. The shape of the curves also suggests that the dataset did not reach its plateau. In the upcoming years, as the dataset size increases, ROC AUC and accuracy enhancement is expected.
On measuring image quality
Through the examined period, Mosquito Alert outreach was promoting a mosquito-targeted data collection strategy. Participants were expected to report two mosquito species (Ae. aegypti and Ae. albopictus). By defining these species as positive samples and all the other potential species of mosquito as negative, the submission decision by participants becomes a binary classification problem. In the majority of cases, when participants submit an image we should expect them to think of having a positive sample. Later, based on entomological expert validation, the true label for the image was obtained.
The main goal of such a surveillance system is to keep the sensitivity of the users as high as possible while keeping their specificity at an acceptable level. Therefore, measuring the sensitivity and specificity of the users would be a plausible quality measure. Unfortunately, there is no available information regarding the non-submitted mosquitoes (the true negative and false negative ones), meaning it is impossible to measure sensitivity. The specificity can be measured only in a special case, when there are no false positive images submitted by the user, resulting in a specificity of 1. Based on the latter argument, focusing on metrics derived from the ratio of the submitted tiger mosquito images vs. all submitted images is not meaningful. Instead, the quality can be measured by the usefulness of the photos from the viewpoint of the expert validator or a CNN, as presented in the next chapter.
Quality evolution of the images through time and space
The Mosquito Alert dataset is a unique collection of mosquito images, because, among other things, it is built from 5 consecutive years (not counting 2014, where less than 100 reports were submitted) and it also provides geolocation tags. This uniqueness of the dataset provides potential identification of time and spatial evolution and dependence of the citizen-based mosquito image quality. To explore such an evolution, we performed two different experiments. Geolocation tags were converted to country, region, and city-level information via the geopy Python package. It was found, that the vast majority (95% of all) of the reports were coming from Spain so we performed the analysis only for the Spanish data.
Figure 5

Number of submitted reports and the fraction of their ratio where the entomology expert annotator could tell if tiger mosquito was presented on the photo or not. The charts are shown for the four cities, where Mosquito Alert was the most popular.

Full size image

First, we explored the fraction of the photos, where the entomology expert marked “can not tell”, because the photo was not descriptive enough to decide which species were presented. Figure 5 shows the ratio of the useful mosquito reports, when mosquito decision was possible, compared to all the mosquito reports. The chart shows the above-mentioned ratio for four Spanish cities, which have the most reports submitted (the same information is showed on Supplementary Fig. S1 as a heatmap over Spain). The Mann–Kendall test on the fraction of useful reports shows p-values of 0.09, 0.09, 0.81, 0.22 for Barcelona, Valencia, Málaga, and Girona, which does not justify the presence of a significant trend in image quality, although any conclusions drawn from five data points must be handled with a pinch of salt. It does not mean anything about the individual participants’ quality progression, because Mosquito Alert is highly open and dynamic, and active participants can constantly change. Of note, through these years, the tiger mosquitoes have widely spread from the east coast to the southern and western regions of Spain29. New (and naive) citizen scientists living in the newly colonized regions have been systematically called to action and participation, thus, limiting the overall learning rate of the Mosquito Alert participants’ population. Our results suggest, that either a dynamic balance exists between naive and experienced participants over the period of data recollection, or mosquito photographing skills are independent of the user experience level. The expectation would be that as the population in Spain became more aware of the presence of tiger mosquitoes and their associated public health risks, the system should experience an increase in the useful report ratio, at least for tiger mosquitoes, and most tiger mosquito photos maybe classified automatically.
Figure 6

1000 random samples were selected for each years data. Separated ResNet50 models were trained on each of the years and each model was tested on the rest of the years data. Metrics were calculated with a 500 round random sampling with replacement from the test data. Left: mean of the 500 round bootstrapped accuracy calculations. Right: mean of the 500 round bootstrapped ROC AUC calculations.

Full size image

Second, we subsampled randomly 1000 images from all years between 2015 and 2019. Then we trained a different ResNet50 on data from the different years and generated predictions for the rest of the data, for each year separately. This way we can explore if data from any year is a “better training material” than the others. The results see Fig. 6, shows that 2015 is the worst training material, providing 0.83–0.84 ROC AUC score for the test period, while the rest (period 2016–2019) is similar, ROC AUC varies between 0.90 and 0.93. The reason why the 2015 data found to be the least favourable for training is its class imbalance, meaning that data from 2015 is extremely biased towards tiger mosquitoes (94%), so when training on 2015 data, the model does not see enough non-tiger mosquito samples, while for the other years lower class imbalance was found (70–80%), see Table 1. In general, machine learning models for classification require a substantial amount of examples for each possible class, in our case tiger and non-tiger mosquitoes, therefore worse performance is expected when training on the 2015 data.
Other than the varying class imbalance, we can conclude that the Mosquito Alert dataset quality is consistent, we did not find any concerning difference between training and testing our model for any of the 2016–2017–2018–2019 data pairs.
Pre-filtering the images before expert validation
Generating human annotations for an image classification task is a labour-intensive and expensive part of any project especially if the annotation requires expert knowledge. Therefore, having a model that generates accurate predictions for a well-defined subset of the data saves a lot of time and cost. We assume that the trained classifier is more accurate when the prediction probability is whether high or low and more inaccurate when it is close to 0.5. With this assumption in mind one can tune the (p_{low}) and (p_{high}) probabilities, in a way that images with a prediction probability (p_{low}< p < p_{high}) are discarded and sent to human validation. Figure 7 Randomly selecting 100,000 (p_{low}) and (p_{high}) thresholds on the predictions which were created via yearly cross-validation. Each time only samples were kept where the predicted probability were out of the ([p_{low};p_{high}]) interval. Each point shows the kept data fraction and the prediction accuracy. Varying the lower and upper predicted probability almost 98% of the images are correctly predicted while keeping 80% of all the images. Full size image Varying (p_{low}) and (p_{high}) provides a trade-off between prediction accuracy and the portion of images sent to human validation. Based on Fig. 7 sending 20% of the images to human validation while having an almost 98% accurate prediction for 80% of the dataset is a fruitful way to combine human labour-power and machine learning together. More

Variation in wood density
The values of Chinese fir’s wood physical properties varied considerably among different geographic sources and Tukey-HSD testing showed that some of these differences were statistically significant (Fig. 1). The maximum value (HNYX-T) of wood all-dry density (WDD) was 62.70% higher than the minimum (FJYK-P). The WDD of each source was consistent with the classification and performance indexes of conifer trees in the timber strength grade for structural use, a standard in China’s forestry industry39: FJYK-P was at level S10 ( HNZJJ-P  > FJYK-P, for which the maximum 58.0% higher than the minimum value. According to the wood grading standards in the grain compression index, HNZJJ-P and FJYK-P were at level II (29.1–44.0 MPa) and the rest of geographic sources were at level III (44.1–59.0 MPa) (Table 3).
The compression strength perpendicular to the grain of total tensile (CPG.TT) among geographic sources was ranked as follows: HNYX-T  > JXCS-R  > HNYX-P  > HNZJJ-P  > FJYK-P (Table 4). Its maximum value (HNYX-T) was 29.3% higher than the minimum (FJYK-P). The ranking for compression strength perpendicular to the grain of total radial (CPG.TR) was slightly different: HNYX-T  > JXCS-R  > HNYX-P  > HNZJJ-P  > FJYK-P, for which the maximum was 42.1% higher than the minimum value. Compression strength perpendicular to the grain of part radial (CPG.PR) had the same rank order as CPG.TT, with a maximum value (HNYX-T) 35.0% higher than the minimum (FJYK-P). Finally, compression strength perpendicular to the grain of part tensile (CPG.PT) was ranked as HNYX-T  > JXCS-R  > HNZJJ-P  > HNYX-P  > FJYK-P for the five geographic sources of Chinese fir.
Table 4 The statistical analysis of wood mechanical properties of Chinese fir.
Full size table

Factors influencing wood physical properties
Climate factors effect on wood physical properties
The influence of precipitation on the three kinds of density was consistent. Pre in January, October, November, and December was positively related to wood density, while it was negatively correlated with density in others months, especially in May (r = − 0.39), June (r = − 0.59), and August (r = − 0.64). On a seasonal scale, Pre in summer was negatively correlated with density (r = − 0.77), but it was positively correlated with autumn (r = 0.22). MaxT was positively correlated with density during the whole year, except in May (r = − 0.34), and likewise with wood density but most strongly in summer (r = 0.75). MinT was positively correlated with density, especially in Jan (r  > 0.7), though it was not significantly so in February and October (r  0.45). Pre showed no significant correlation with TSR.LD, RSR.LD, DDS.LD, and DDS.RD, whose correlation coefficients were 0.1–0.3. But Pre was negatively correlated with VSR.LD most of the year (except July, October). AveT was negatively correlation with TSR.RD, RSR.RD, and VSR.RD in January, February, March, and winter; however, AveT showed no significant correlation with DDS.RD. AveT was negatively correlated with TSR.LD, RSR.LD, DDS.LD, and VSR.LD during the whole year. In general, MinT had a significant positive relationship to TSR.RD (r = 0.47), RSR.RD (r = 0.48), and VSR.RD (r = 0.52), except in October, and it was negatively correlated with DDS.RD. MinT was positively related to RSR.LD, VSR.LD, yet negative related to DDS.LD. MaxT was negatively correlated with TSR.RD, RSR.RD, VSR.RD in January, February, May, and December, and winter. MaxT showed no significant correlation with DDS.RD, RSR.LD, DDS.LD or VSR.LD (Fig. 2c).
Pre had significant negative correlations with all of the mechanical properties in May, June, August, and summer, as evince by Fig. 2b, which also showed positive correlations in October. As we can seen, the effects of Pre on wood density and mechanical properties have the same tendency. Pre in all other months was not significantly correlated with mechanical properties (r  0.75), while it was showed no significant correlation in Feb and Oct (r  1000. Through stepwise regression modeling, 14 variables without multicollinearity were retained (i.e., MOE, MOR, TSG, CSG, CPG.TT, CPG.TR, CPG.PT, CPG.PR, DDS.RD, WDD, DDS.LD, TSR.RD, RSR.RD, VSR.LD).
PCA was applied to the above 14 selected physical variables. These results showed that the physical properties of wood loaded strongly on the first axis of the PCA, explaining 51.8% of variation in the 14 tested properties, while the second axis explained 11.0% of it. MOE, MOR, TSR.RD, RSR.RD, and VSR.LD loaded on the positive axis of PC1 and PC2. Both DDS.LD and DDS.RD loaded on the negative axis of PC1 and PC2, while TSG, CSG, CPG.TT, CPG.TR, CPG.PT, CPG.PR, and WDD loaded on the positive axis of PC1 and the negative axis of PC2 (Fig. 3). For a comprehensive evaluation of Chinese fir’s wood physical properties, we calculated the comprehensive scores of five geographic sources via the PCA. In this respect, significant differences were detected among the five geographic sources. Among them, the comprehensive score of HNYX-T was the highest whereas that of FJYK-P was the lowest (Fig. 4).
Figure 3

Sequence diagram plot of PCA analysis showing the relationship among physical properties of wood.

Full size image

Figure 4

Mean comprehensive score of PCA plot with 95% CI. Different letters (a, b, c, d, e) mean significant difference at 0.05 level.

Full size image More

1.
Brown, J. H. Why are there so many species in the tropics? J. Biogeogr. 41, 8–22 (2014).
PubMed  Article  PubMed Central  Google Scholar 
2.
Classen, A. et al. Temperature versus resource constraints: Which factors determine bee diversity on Mount Kilimanjaro, Tanzania? Glob. Ecol. Biogeogr. 24, 642–652 (2015).
Article  Google Scholar 

3.
Rahbek, C. et al. Humboldt’s enigma: What causes global patterns of mountain biodiversity? Science 365, 1108–1113 (2019).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

4.
Toranza, C. & Arim, M. Cross-taxon congruence and environmental conditions. BMC Ecol. 10, 18 (2010).
PubMed  PubMed Central  Article  Google Scholar 

5.
Gioria, M., Bacaro, G. & Feehan, J. Evaluating and interpreting cross-taxon congruence: Potential pitfalls and solutions. Acta Oecol. 37, 187–194 (2011).
ADS  Article  Google Scholar 

6.
Graham, C. H. et al. The origin and maintenance of montane diversity: Integrating evolutionary and ecological processes. Ecography (Cop.) 37, 711–719 (2014).
Article  Google Scholar 

7.
Westgate, M. J., Tulloch, A. I. T., Barton, P. S., Pierson, J. C. & Lindenmayer, D. B. Optimal taxonomic groups for biodiversity assessment: A meta-analytic approach. Ecography (Cop.) 40, 539–548 (2017).
Article  Google Scholar 

8.
Lomolino, M. V. Elevation gradients of species-density: Historical and prospective views. Glob. Ecol. Biogeogr. 8, 1–2 (2001).
Google Scholar 

9.
McCain, C. M. Global analysis of bird elevational diversity. Glob. Ecol. Biogeogr. 18, 346–360 (2009).
Article  Google Scholar 

10.
Peters, M. K. et al. Predictors of elevational biodiversity gradients change from single taxa to the multi-taxa community level. Nat. Commun. 7, 13736 (2016).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

11.
Sundqvist, M. K., Sanders, N. J. & Wardle, D. A. Community and ecosystem responses to elevational gradients: Processes, mechanisms, and insights for global change. Annu. Rev. Ecol. Evol. Syst. 44, 261–280 (2013).
Article  Google Scholar 

12.
Ruggiero, A. & Hawkins, B. A. Why do mountains support so many species of birds? Ecography (Cop.) 31, 306–315 (2008).
Article  Google Scholar 

13.
Mccain, C. M. & Colwell, R. K. Assessing the threat to montane biodiversity from discordant shifts in temperature and precipitation in a changing climate. Ecol. Lett. 14, 1236–1245 (2011).
PubMed  Article  PubMed Central  Google Scholar 

14.
Hawkins, B. A. et al. Energy, water, and broad-scale geographic patterns of species richness. Ecology 84, 3105–3117 (2003).
Article  Google Scholar 

15.
Currie, D. J. Energy and large-scale patterns of animal and plant species richness. Am. Nat. 137, 27–49 (1991).
Article  Google Scholar 

16.
Stein, A., Gerstner, K. & Kreft, H. Environmental heterogeneity as a universal driver of species richness across taxa, biomes and spatial scales. Ecol. Lett. 17, 866–880 (2014).
PubMed  Article  PubMed Central  Google Scholar 

17.
Costanza, J. K., Moody, A. & Peet, R. K. Multi-scale environmental heterogeneity as a predictor of plant species richness. Landsc. Ecol. 26, 851–864 (2011).
Article  Google Scholar 

18.
Vetaas, O. R., Paudel, K. P. & Christensen, M. Principal factors controlling biodiversity along an elevation gradient: Water, energy and their interaction. J. Biogeogr. https://doi.org/10.1111/jbi.13564 (2019).
Article  Google Scholar 

19.
Lande, R. Risks of population extinction from demographic and environmental stochasticity and random catastrophes. Am. Nat. 142, 911–927 (1993).
PubMed  Article  PubMed Central  Google Scholar 

20.
Kaspari, M., Alonso, L. & O’Donnell, S. Three energy variables predict ant abundance at a geographical scale. Proc. R. Soc. B Biol. Sci. 267, 485–489 (2000).
CAS  Article  Google Scholar 

21.
Pianka, E. R. Latitudinal gradients in species diversity: A review of concepts. Am. Nat. 100, 33–46 (1966).
Article  Google Scholar 

22.
Werenkraut, V. & Ruggiero, A. The richness and abundance of epigaeic mountain beetles in north-western Patagonia, Argentina: Assessment of patterns and environmental correlates. J. Biogeogr. 41, 561–573 (2014).
Article  Google Scholar 

23.
R Core Team. R version 3.6.2 ‘Dark and Stormy Night’ (2019). (Accessed 12 December 2019). https://www.r-project.org. 

24.
Blanchet, F. G., Cazelles, K. & Gravel, D. Co-occurrence is not evidence of ecological interactions. Ecol. Lett. 23, 1050–1063 (2020).
PubMed  Article  Google Scholar 

25.
Hodkinson, I. D. Terrestrial insects along elevation gradients: Species and community responses to altitude. Biol. Rev. 80, 489–513 (2005).
PubMed  Article  Google Scholar 

26.
Kampmann, D. et al. Mountain grassland biodiversity: Impact of site conditions versus management type. J. Nat. Conserv. 16, 12–25 (2008).
Article  Google Scholar 

27.
Janzen, D. H. et al. Changes in the arthropod community along an elevational transect in the Venezuelan Andes. Biotropica 8, 193–203 (1976).
Article  Google Scholar 

28.
Sirin, D., Eren, O. & Ciplak, B. Grasshopper diversity and abundance in relation to elevation and vegetation from a snapshot in Mediterranean Anatolia: Role of latitudinal position in altitudinal differences. J. Nat. Hist. 44, 1343–1363 (2010).
Article  Google Scholar 

29.
Alexander, G. & Hilliard, J. R. Altitudinal and seasonal distribution of Orthoptera in the Rocky Mountains of northern Colorado. Ecol. Monogr. 39, 385–432 (1969).
Article  Google Scholar 

30.
Mojica, A. S. & Fagua, G. Estructura de las comunidades de orthoptera (insecta) en un gradiente altitudinal de un bosque andino. Rev. Colomb. Entomol. 32, 200–213 (2006).
Google Scholar 

31.
Grytnes, J. A. Species-richness patterns of vascular plants along seven altitudinal transects in Norway. Ecography (Cop.) 26, 291–300 (2003).
Article  Google Scholar 

32.
McCain, C. M. & Grytnes, J.-A. Elevational gradients in species richness. Encyl. Life Sci. https://doi.org/10.1002/9780470015902.a0022548 (2010).
Article  Google Scholar 

33.
Xu, X. et al. Altitudinal patterns of plant species richness in the Honghe region of China. Pak. J. Bot. 49, 1039–1048 (2017).
Google Scholar 

34.
Kerr, J. T. & Packer, L. Habitat heterogeneity as a determinant of mammal species richness in high-energy regions. Nature 385, 252–254 (1997).
ADS  CAS  Article  Google Scholar 

35.
Röder, J. et al. Heterogeneous patterns of abundance of epigeic arthropod taxa along a major elevation gradient. Biotropica 49, 217–228 (2017).
Article  Google Scholar 

36.
Evans, K. L., Warren, P. H. & Gaston, K. J. Species-energy relationships at the macroecological scale: A review of the mechanisms. Biol. Rev. Camb. Philos. Soc. 80, 1–25 (2005).
PubMed  Article  Google Scholar 

37.
Kissling, W. D., Rahbek, C. & Böhning-Gaese, K. Food plant diversity as broad-scale determinant of avian frugivore richness. Proc. R. Soc. Biol. Sci. 274, 799–808 (2007).
Article  Google Scholar 

38.
Kissling, W. D., Field, R. & Böhning-Gaese, K. Spatial patterns of woody plant and bird diversity: Functional relationships or environmental effects? Glob. Ecol. Biogeogr. 17, 327–339 (2008).
Article  Google Scholar 

39.
Chown, S. L. & Gaston, K. J. Exploring links between physiology and ecology at macro scales: The role of respiratory metabolism in insects. Biol. Rev. 74, 87–120 (1999).
Article  Google Scholar 

40.
de Araújo, W. S. Different relationships between galling and non-galling herbivore richness and plant species richness: A meta-analysis. Arthropod. Plant. Interact. 7, 373–377 (2013).
Article  Google Scholar 

41.
Qian, H. & Kissling, W. D. Spatial scale and cross-taxon congruence of terrestrial vertebrate and vascular plant species richness in China. Ecology 91, 1172–1183 (2010).
PubMed  Article  PubMed Central  Google Scholar 

42.
Burrascano, S. et al. Congruence across taxa and spatial scales: Are we asking too much of species data? Glob. Ecol. Biogeogr. 27, 980–990 (2018).
Article  Google Scholar 

43.
Field, R. et al. Spatial species-richness gradients across scales: A meta-analysis. J. Biogeogr. 36, 132–147 (2009).
Article  Google Scholar 

44.
Giorgis, M. A. et al. Composición florística del Bosque Chaqueño Serrano de la provincia de Córdoba, Argentina. Kurtziana 36, 9–43 (2011).
Google Scholar 

45.
Cabido, M., Funes, G., Pucheta, E., Vendramani, F. & Díaz, S. A chorological analysis of the mountains from Central Argentina. Is all what we call Sierra Chaco really Chaco? Contribution to the study of the flora and vegetation of the Chaco: 12. Candollea 53, 321–331 (1998).
Google Scholar 

46.
Giorgis, M. A. et al. Changes in floristic composition and physiognomy are decoupled along elevation gradients in central Argentina. Appl. Veg. Sci. 20, 553–571 (2017).
Article  Google Scholar 

47.
Cabrera, A. L. Fitogeografia de la República Argentina. In Enciclopedia Argentina de Agricultura y Jardinería Vol. 14 (ed. Kugler, W. F.) 1–42 (ACME, New York, 1976).
Google Scholar 

48.
Giorgis, M. A. et al. Diferencias en la estructura de la vegetación del sotobosque entre una plantación de Pinus taedaL. (Pinaceae) y un matorral serrano (Cuesta Blanca, Córdoba). Kurtziana 31, 39–49 (2005).
Google Scholar 

49.
Martínez, G. A., Arana, M. D., Oggero, A. J. & Natale, E. S. Biogeographical relationships and new regionalisation of high-altitude grasslands and woodlands of the central Pampean Ranges (Argentina), based on vascular plants and vertebrates. Aust. Syst. Bot. 29, 473–488 (2016).
Article  Google Scholar 

50.
QGIS Development Team. QGIS Geographic Information System (Accessed 19 April 2019). (2019).

51.
Kent, M. The description of vegetation in the field. In Vegetation Description and Data Analysis: A Practical Approach (ed. Kent, M.) 65–116 (Wiley-Blackwell, Hoboken, 2012).
Google Scholar 

52.
Catálogo de las plantas vasculares del Cono Sur : (Argentina, Sur de Brasil, Chile, Paraguay y Uruguay). (Missouri Botanical Garden Press, 2008).

53.
Haddad, N., Tilman, D., Haarstad, J., Ritchie, M. & Knops, J. M. N. Contrasting effects of plant richness and composition on insect communities: a field experiment. Am. Nat. 158, 17–35 (2001).
CAS  PubMed  Article  PubMed Central  Google Scholar 

54.
Braun, H. & Zubarán, G. Tettigoniidae (Orthoptera) Species from Argentina and Uruguay (2019).

55.
Carbonell, C. S., Cigliano, M. M. & Lange, C. E. Acridomorph (Orthoptera) Species of Argentina and Uruguay. Version II [2019]. https://biodar.unlp.edu.ar/acridomorph/.

56.
Cigliano, M. M., Braun, H., Eades, D. C. & Otte, D. Orthoptera Species File. Version 5.0/5.0 (2018). http://orthoptera.speciesfile.org.

57.
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 

58.
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).
Article  Google Scholar 

59.
Carrara, R., Silvestro, V. A., Cheli, G. H., Campón, F. F. & Flores, G. E. Disentangling the effect of climate and human influence on distribution patterns of the darkling beetle Scotobius pilularius Germar, 1823 (Coleoptera: Tenebrionidae). Ann. Zool. 66, 693–701 (2016).
Article  Google Scholar 

60.
Aisen, S., Werenkraut, V., Márquez, M. E. G., Ramírez, M. J. & Ruggiero, A. Environmental heterogeneity, not distance, structures montane epigaeic spider assemblages in north-western Patagonia (Argentina). J. Insect Conserv. 21, 1–12 (2017).
Article  Google Scholar 

61.
Bilskie, J. Soil Water Status: Content and Potential (Campbell Scientific Inc., Logan, 2001).
Google Scholar 

62.
Tucker, C. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sens. Environ. 8, 127–150 (1979).
ADS  Article  Google Scholar 

63.
Wang, J., Rich, P. M., Price, K. P. & Dean-Kettle, W. Relations between NDVI, grassland production, and crop yield in the central great plains. Geocarto Int. 20, 5–11 (2005).
Article  Google Scholar 

64.
Oindo, B. O., de By, R. A. & Skidmore, A. K. Interannual variability of NDVI and bird species diversity in Kenya. Int. J. Appl. Earth Obs. Geoinf. 2, 172–180 (2000).
ADS  Article  Google Scholar 

65.
IGN. Modelo Digital de Elevaciones de la República Argentina. (Instituto Geográfico Nacional—Dirección General de Servicios Geográficos—Dirección de Geodesia, 2016).

66.
Riley, S. J., DeGloria, S. D. & Elliot, R. A terrain ruggedness index that quantifies topographic heterogeneity. Intermt. J. Sci. 5, 23–27 (1999).
Google Scholar 

67.
Stein, A. & Kreft, H. Terminology and quantification of environmental heterogeneity in species-richness research. Biol. Rev. 90, 815–836 (2015).
PubMed  Article  Google Scholar 

68.
Tilman, D. & Pacala, S. W. The maintenance of species richness in plant communities. In Species Diversity in Ecological Communities (eds Ricklefs, R. E. & Schulter, D.) 13–25 (University of Chicago Press, Chicago, 1993).
Google Scholar 

69.
Cleveland, W. S., Grosse, E. & Shyu, W. M. Local regresion models. In Statistical Models in S (eds Chambers, J. M. & Hastie, T. J.) 227 (Chapman and Hall, London, 1993).
Google Scholar 

70.
Szewczyk, T. & Mccain, C. M. A systematic review of global drivers of ant elevational diversity. PLoS ONE 11, e0155404 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

71.
Beck, J. et al. Elevational species richness gradients in a hyperdiverse insect taxon: A global meta-study on geometrid moths. Glob. Ecol. Biogeogr. 26, 412–424 (2017).
Article  Google Scholar 

72.
Bolker, B. M. Ecological Statistics: Contemporary Theory and Application (Oxford University Press Inc., Oxford, 2015).
Google Scholar 

73.
Grace, J. B. Structural Equation Modeling and Natural Systems (Cambridge University Press, Cambridge, 2006).
Google Scholar 

74.
Lefcheck, J. S. piecewiseSEM: Piecewise structural equation modelling in r for ecology, evolution, and systematics. Methods Ecol. Evol. 7, 573–579 (2016).
Article  Google Scholar 

75.
Jiménez-Alfaro, B., Chytrý, M., Mucina, L., Grace, J. B. & Rejmánek, M. Disentangling vegetation diversity from climate-energy and habitat heterogeneity for explaining animal geographic patterns. Ecol. Evol. 6, 1515–1526 (2016).
PubMed  PubMed Central  Article  Google Scholar  More

We have introduced and investigated a spatially explicit evolutionary food web model that allows us to explore the distribution of species in space and time, as well as the waxing and waning of species ranges with time. This is the first model that makes it possible to explore these features in the context of trophic networks, showing how they are influenced by competition as well as by predators and prey. Indeed, our model produces empirically well-known patterns in space and time, such as lifetime distributions, species-area relationships, distance decay of similarity, and temporal change of geographic range. On top, we obtain a variety of additional result and gain insights into the mechanisms that generate these patterns. While one might argue that some patterns must emerge trivially in a model like ours, the multitude of patterns emerging together is remarkable.
We find that most species only appear for short times and have small ranges. For species that conquer larger portions of the web we analysed the shape of the range evolution and found that basal species show the empirically observed “hat” pattern more often than species in higher trophic levels. This indicates that the trophic position of a species plays a major role for its range expansion success as well as the shape of its range expansion trajectory over time. To our knowledge, this has not been discussed in the literature so far.
Recently Zliobaite et al.11 analyzed which factors are more correlated to the rise and fall of the range expansion trajectory in fossil data sets of basal mammals. The range expansion follows the so called “hat pattern” that consists of five phases in a species lifetime: origination, expansion, peak, decline and extinction. They found that the temporal location of the peak of the hat pattern is more impacted by competition while the brims are more influenced by abiotic environmental factors. They also compared the range curves of different random walk models with the shape of empirical data and found that a random walk model with competition and environmental factor provides the most realistic looking curves.
Our results for basal species provide an illustration of the mechanism that might lead to these findings. There is one major difference in assumptions: the “environment” in our case is the trophic environment (network structure and abundance distribution). We do not model an abiotic surrounding, yet basal range curves look strikingly hat shaped. A species needs to fit into this trophic environment (network) to first establish a viable population (origination). To successfully increase its range it needs to disperse to neighbouring habitats and be a viable competitor there as well. This leads to the extinction of another species as the dispersing competitor takes its place. As neighbouring basal communities are similar in our systems, the chances are high that the species can spread on a large portion on the grid replacing other species (expansion). This continues until the species has reached the maximum range (peak). It is only a matter of time then until this process is repeated with the species having become the inferior competitor. The species is then successively replaced by a better adapted species (decline). The species thus ages as the network structure changes. At some point the species has vanished on all habitats (extinction). This means that we observe the same dynamics as suggested by Zliobaite et al., but with the trophic and not the abiotic environment as the main driver of the initial increase and later decrease of the range. The truth is probably that both the trophic and abiotic environment are important, as both play an important role in real ecosystems.
Regarding higher trophic species in our system the range expansion curves look more diverse and often do not resemble the hat shape. These species depend on the composition on the layer below. As this layer changes the fate of the higher species changes as well. The emergence of a new prey species that spreads over the network can save a consumer species from extinction. This cannot happen for basal species as these are more or less completely controlled by competition. The studies that we know often deal with basal species, so we are not convinced that the hat pattern is ubiquitous for all species. Future empirical work could focus on predator range expansion and try to find a case where a predator species could regain its range after the emergence of a new prey.
A qualitative difference between the basic trophic layer and higher layers occurred in our data also with respect to the similarity of networks in nearby habitats. The similarity index decays particularly slowly for the basal layer. Theory on distance decay suggests that spatial heterogeneity is a main driver in community turnover in two ways: (1) competitive species sorting along environmental gradients and (2) topological influences that let species with different dispersal abilities experience different landscapes6. As we use a homogeneous landscape we expect the first driver to be non-existent for the basal layer, as all habitats hold the same type and amount of resource. As species do not fundamentally differ in their dispersal abilities and we do not have a heterogeneous spatial topology, the second point is only weakly relevant for the basal species. They have slightly different chances of being chosen for dispersal as we choose the next disperser depending on the biomass density. As we have seen, biomass densities are quite similar for basal species. What remains is a temporal aspect: species that are older can reside on more habitats and have thus a higher chance of being chosen to disperse. To put the cart before the horse, this confirms the theory on distance decay: we expect a much faster decay in a heterogeneous environment, and this is exactly what we observe for higher trophic layers, which experience heterogeneity due to the spatial turnover in basal species composition. A trend to faster decay rates in higher trophic levels was also found in a meta-study7.
A model is always a simplification of reality. Some of the assumptions underlying our model are worth discussing. Species in nature are not restricted to the one-dimensional niche space that we assume. In fact, the original Webworld model characterised species by a large vector of traits31. We, in contrast, characterized species by three traits, all of which are based on body mass. We think that this is the reason why we do not observe super abundant species, but species densities are all of a similar order of magnitude. With a higher-dimensional trait space there must exist more diverse species types and probably also super abundant species, which have a globally optimal trait vector. Nevertheless, our simplification leads to an overall shape of the rank abundance curves that is realistic, as one would expect from a niche apportionment model40. Harpole and Tilman showed that diversity and evenness of grassland communities decreased when niche dimensionality was reduced by adding limiting nutrients to plot experiments41. In turn this indicates that rank abundance curves for less dimensional communities will be flatter. This is in line with the shape of our rank abundance curves.
The probably most interesting trait to add to each species would be its dispersal rate, and to let the dispersal rate evolve. We would expect that this would lead to higher-level species dispersing faster than lower-level species, thus making the differences in range between the different trophic levels smaller. In fact, a previous predator-prey model showed that the predator’s dispersal ability evolved in accordance to spatio-temporal fluctuations of the prey; with higher dispersal rates evolving for larger fluctuations in prey42. However, in the context of food webs the evolution of dispersal is poorly understood and a model like ours would be a good starting point.
Our choice of parameters is guided by the aim to make the model feasible. To be able to perform computer simulations on a large number of habitats we chose a relatively small value for the amount of resource R, so that the number of species of a local food web remained below 25. As we wanted to simulate food webs and not only basal communities we needed to choose a value of the efficiency (lambda ) that allows for the emergence of several trophic layers. In the original Webworld Model, (lambda ) was identified with the proportion of biomass that is passed from one trophic layer to the next. The value of 0.65 that we use here is much larger than the empirically established value of 0.143. However, in the original Webworld model no distinction was made between resident biomass and biomass fluxes. Therefore the variable B was identified with biomass, while its occurrence in Eq. (4) in fact suggests that it represents biomass flux. A further reason for the difference between the model value and the empirical value of (lambda ) is that the model does not take into account energy input due to the below-ground ecological processes. Due to all these simplifications of the model, we do not consider it important to provide an empirical justification of the precise value of the parameter (lambda ), but base its value on the condition that the model yields food webs with several trophic levels.
In contrast to other models, the model used here does not rely on an extrinsic extinction rate that randomly extirpates species that might be well adapted to the network. All extinction events are driven by the trophic dynamics, yet we observe an ongoing species turn over. We thus study the pure food web dynamics without a heterogeneous or fluctuating environment and still observe ecological reasonable species distributions. This indicates that incorporating abiotic environments and their fluctuations is not necessarily needed to study food web dynamics.
Rogge et al. analysed lifetime distributions and SAR curves in a model that is simpler than ours as it does not include population sizes22. The lifetime distributions that we find are considerably steeper (slope (-2.4)) compared to their value around (-1.7). It is also larger than the values reported for empirical findings, which lie between 1 and 2; Newman and Palmer pin them down to (1.7pm 0.3)10. Data for contemporary lifetime distributions show a power-law like shape23,44 with exponents that are in agreement with the exponent of paleological data10. It is noteworthy that there is no consensus whether lifetime distributions follow a power-law or an exponential law, as data often allow for both types of fit, due to (large) uncertainties in fossil data10,45. Exponentials, of course, have a changing slope in a double-logarithmic plot and can thus also be compatible with the exponent observed by us.
Curiously, our model shows steep lifetime distributions even though there is no external random extinction implemented as in other models22. One implication of our value (alpha > 2) is that our distribution has a well-defined mean. This features is shared with exponential distributions.
McPeek argues that lifetime distributions depend on the number and survival time of “transient” species, i.e. species that are on their way to extinction25. He reasons that the time to extinction is elongated for species that are similar, because the inferior competitor holds out longer when it competes with more similar species. If this applies to our type of model, this indicates that species in our system are, despite the one-dimensional niche axis, not as similar as species in the model of Rogge et al.22 that uses the same niche axis, as we observe shorter lifetimes. The difference is that interaction links in22 are binary (presence versus absence), whilst we use Gaussian feeding kernels. The fact that this difference affects the lifetime distributions emphasises the importance of considering details of the trophic interactions. The SAR curves on the contrary are flatter in our model than in the model by Rogge et al.  and in better agreement with empirical data.
O’Sullivan et al.18 found in a competitive metacommunity assembly models a similar collection of macroecological patterns (SAD, range size distribution RSD and SAR) as we did, when regional diversity was near equilibrium. They refer to the work of McGill16 who analysed the assumptions underlying models of macroecological patterns and found that three key ingredients seem to be sufficient for such patterns to emerge. Those are a left skewed SAD, clumping of populations in space, and species distributions in space that are uncorrelated from other species spatial distributions. O’Sullivan et al.18 report that all three ingredients occur in their model and are shaped by regional diversity equilibrium. The closer the system to regional equilibrium the stronger are the observed key patterns (SAD, RSD, Spatial non-correlation). They relate their finding with the theory of ecological structural stability, which revolves around the dynamics on a regional scale. Our communities, in contrast, are trophic communities, operate always near local and regional species equilibrium, i.e. in the regime where O’Sullivan and coauthors18 find the most prominent form of the basic patterns. Comparing the patterns we observe, we also see SADs that are left skewed, and a local clumping of species. We did not analyse the spatial correlation between species. As we have trophic layers of species there will be some correlation between predators and their prey as they can only persist in a habitat if prey is present. In addition to the results obtained by O’Sullivan et al.18, we also derive liefetime distributions, i.e., a paleoecological pattern that also seem to be connected to the metacommunity dynamics. This might indicate that spatial non-correlation is not the most important factor in the mechanisms producing macroecological patterns.
To conclude, our evolutionary food web model produces empirically well studied ecological and paleological patterns. We thus are armed with a valuable tool to broaden our understanding of the mechanisms behind those patterns. Our findings that trophic position influences geographic range and lifetime of a species might motivate further work regarding the interplay of abiotic and trophic factors on range expansion on evolutionary time scales.
More generally, evolutionary models can assist us in forming a deeper knowledge of the processes that lead to what is remnant in fossils. As recently pointed out by Marshall46 in his fifth law of paleobiology, extinction erases information. It is a strength of evolutionary food web models that they allow us to study processes whose extent eludes direct observations. More