Ecology

1.Barrett, S. C. & Hough, J. Sexual dimorphism in flowering plants. J. Exp. Bot. 64, 67–82 (2012).PubMed 
Article 
CAS 

Google Scholar 
2.Retuerto, R., Sánchez Vilas, J. & Varga, S. Sexual dimorphism in response to stress. Environ. Exp. Bot. 146, 1–4 (2018).Article 

Google Scholar 
3.Poissant, J., Wilson, A. J. & Coltman, D. W. Sex-specific genetic variance and the evolution of sexual dimorphism: a systematic review of cross-sex genetic correlations. Evolution 64, 97–107 (2010).PubMed 
Article 

Google Scholar 
4.Bonduriansky, R. & Chenoweth, S. F. Intralocus sexual conflict. Trends Ecol. Evol. 24, 280–288 (2009).PubMed 
Article 

Google Scholar 
5.Pennell, T. M., de Haas, F. J., Morrow, E. H. & van Doorn, G. S. Contrasting effects of intralocus sexual conflict on sexually antagonistic coevolution. Proc. Natl Acad. Sci. USA 113, E978–E986 (2016).CAS 
PubMed 
Article 

Google Scholar 
6.Charlesworth, B. & Charlesworth, D. A model for the evolution of dioecy and gynodioecy. Am. Nat. 112, 975–997 (1978).Article 

Google Scholar 
7.Lloyd, D. G. & Webb, C. Secondary sex characters in plants. Bot. Rev. 43, 177–216 (1977).Article 

Google Scholar 
8.Torimaru, T. & Tomaru, N. Relationships between flowering phenology, plant size, and female reproductive output in a dioecious shrub, Ilex leucoclada (Aquifoliaceae). Botany 84, 1860–1869 (2006).
Google Scholar 
9.Delph, L. F. & Meagher, T. R. Sexual dimorphism masks life history trade-offs in the dioecious plant Silene latifolia. Ecology 76, 775–785 (1995).Article 

Google Scholar 
10.Carroll, S. B. & Delph, L. F. The effects of gender and plant architecture on allocation to flowers in dioecious Silene latifolia (Caryophyllaceae). Int. J. Plant Sci. 157, 493–500 (1996).Article 

Google Scholar 
11.Delph, L. F., Gehring, J. L., Arntz, A. M., Levri, M. & Frey, F. M. Genetic correlations with floral display lead to sexual dimorphism in the cost of reproduction. Am. Nat. 166, S31–S41 (2005).PubMed 
Article 
PubMed Central 

Google Scholar 
12.Barrett, S. C., Yakimowski, S. B., Field, D. L. & Pickup, M. Ecological genetics of sex ratios in plant populations. Philos. Trans. R. Soc. Lond. B Biol. Sci. 365, 2549–2557 (2010).PubMed 
PubMed Central 
Article 

Google Scholar 
13.Muyle, A., Shearn, R. & Marais, G. A. The evolution of sex chromosomes and dosage compensation in plants. Genome Biol. Evol. 9, 627–645 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
14.Connallon, T. & Knowles, L. L. Intergenomic conflict revealed by patterns of sex-biased gene expression. Trends Genet. 21, 495–499 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
15.Ellegren, H. & Parsch, J. The evolution of sex-biased genes and sex-biased gene expression. Nat. Rev. Genet. 8, 689–698 (2007).CAS 
PubMed 
Article 

Google Scholar 
16.Rice, W. R. Sex chromosomes and the evolution of sexual dimorphism. Evolution 38, 735–742 (1984).PubMed 
Article 

Google Scholar 
17.Charlesworth, B., Jordan, C. Y. & Charlesworth, D. The evolutionary dynamics of sexually antagonistic mutations in pseudoautosomal regions of sex chromosomes. Evolution 68, 1339–1350 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Mank, J. E. The transcriptional architecture of phenotypic dimorphism. Nat. Ecol. Evol. 1, 1–7 (2017).Article 

Google Scholar 
19.Zemp, N. et al. Evolution of sex-biased gene expression in a dioecious plant. Nat. Plants 2, 16168 (2016).CAS 
PubMed 
Article 

Google Scholar 
20.Sanderson, B. J., Wang, L., Tiffin, P., Wu, Z. & Olson, M. S. Sex-biased gene expression in flowers, but not leaves, reveals secondary sexual dimorphism in Populus balsamifera. New Phytol. 221, 527–539.21.Delph, L. F. & Herlihy, C. R. Sexual, fecundity, and viability selection on flower size and number in a sexually dimorphic plant. Evolution: Int. J. Org. Evolution 66, 1154–1166 (2012).Article 

Google Scholar 
22.Golonka, A. M., Sakai, A. K. & Weller, S. G. Wind pollination, sexual dimorphism, and changes in floral traits of Schiedea (Caryophyllaceae). Am. J. Bot. 92, 1492–1502 (2005).PubMed 
Article 
PubMed Central 

Google Scholar 
23.Aloni, R., Aloni, E., Langhans, M. & Ullrich, C. I. Role of auxin in regulating Arabidopsis flower development. Planta 223, 315–328 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
24.Rocheta, M. et al. Comparative transcriptomic analysis of male and female flowers of monoecious Quercus suber. Front. Plant Sci. 5, 599 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
25.Zhao, D. & Tao, J. Recent advances on the development and regulation of flower color in ornamental plants. Front. Plant Sci. 6, 261 (2015).PubMed 
PubMed Central 

Google Scholar 
26.Moreau, C. et al. The b gene of pea encodes a defective flavonoid 3′, 5′-hydroxylase, and confers pink flower color. Plant Physiol. 159, 759–768 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
27.Hao, Z., Liu, S., Hu, L., Shi, J. & Chen, J. Transcriptome analysis and metabolic profiling reveal the key role of carotenoids in the petal coloration of Liriodendron tulipifera. Hortic. Res. 7, 1–16 (2020).Article 
CAS 

Google Scholar 
28.Hormaza, J. & Polito, V. Pistillate and staminate flower development in dioecious Pistacia vera (Anacardiaceae). Am. J. Bot. 83, 759–766 (1996).Article 

Google Scholar 
29.Boucher, L. D., Manchester, S. R. & Judd, W. S. An extinct genus of Salicaceae based on twigs with attached flowers, fruits, and foliage from the Eocene Green River Formation of Utah and Colorado, USA. Am. J. Bot. 90, 1389–1399 (2003).PubMed 
Article 
PubMed Central 

Google Scholar 
30.Manchester, S. R., Judd, W. S. & Handley, B. Foliage and fruits of early poplars (Salicaceae: Populus) from the Eocene of Utah, Colorado, and Wyoming. Int. J. Plant Sci. 167, 897–908 (2006).Article 

Google Scholar 
31.Wu, J. et al. Phylogeny of Salix subgenus Salix sl (Salicaceae): delimitation, biogeography, and reticulate evolution. BMC Evol. Biol. 15, 31 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
32.Liao, J., Cai, Z., Song, H. & Zhang, S. Poplar males and willow females exhibit superior adaptation to nocturnal warming than the opposite sex. Sci. Total Environ. 717, 137179 (2020).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
33.Dawson, T. E. & Bliss, L. Patterns of water use and the tissue water relations in the dioecious shrub, Salix arctica: the physiological basis for habitat partitioning between the sexes. Oecologia 79, 332–343 (1989).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
34.Lei, Y., Chen, K., Jiang, H., Yu, L. & Duan, B. Contrasting responses in the growth and energy utilization properties of sympatric Populus and Salix to different altitudes: implications for sexual dimorphism in Salicaceae. Physiol. Plant 159, 30–41 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
35.Ueno, N., Suyama, Y. & Seiwa, K. What makes the sex ratio female-biased in the dioecious tree Salix sachalinensis? J. Ecol. 95, 951–959 (2007).Article 

Google Scholar 
36.Jiang, H., Zhang, S., Lei, Y., Xu, G. & Zhang, D. Alternative growth and defensive strategies reveal potential and gender specific trade-offs in dioecious plants Salix paraplesia to nutrient availability. Front. Plant Sci. 7, 1064 (2016).PubMed 
PubMed Central 

Google Scholar 
37.Liao, J., Song, H., Tang, D. & Zhang, S. Sexually differential tolerance to water deficiency of Salix paraplesia-A female-biased alpine willow. Ecol. Evol. 9, 8450–8464 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
38.Saska, M. M. & Kuzovkina, Y. A. Phenological stages of willow (Salix). Ann. Appl. Biol. 156, 431–437 (2010).Article 

Google Scholar 
39.Thomas, R., Sheard, R. & Moyer, J. Comparison of conventional and automated procedures for nitrogen, phosphorus, and potassium analysis of plant material using a single digestion 1. Agron. J. 59, 240–243 (1967).CAS 
Article 

Google Scholar 
40.Arnon, D. I. Copper enzymes in isolated chloroplasts. Polyphenoloxidase in Beta vulgaris. Plant Physiol. 24, 1 (1949).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
41.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, Austria, 2020).42.Wickham, H. ggplot2: elegant graphics for data analysis (springer, 2016).43.Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
44.Dai, X. et al. The willow genome and divergent evolution from poplar after the common genome duplication. Cell Res. 24, 1274–1277 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
45.Wei, S., Yang, Y. & Yin, T. The chromosome-scale assembly of the willow genome provides insight into Salicaceae genome evolution. Hortic. Res. 7, 1–12 (2020).Article 
CAS 

Google Scholar 
46.Yu, G., Wang, L. G., Han, Y. & He, Q. Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
47.Thimm, O. et al. MAPMAN: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J. 37, 914–939 (2004).CAS 
PubMed 
Article 

Google Scholar 
48.Lohse, M. et al. M ercator: a fast and simple web server for genome scale functional annotation of plant sequence data. Plant Cell Environ. 37, 1250–1258 (2014).CAS 
PubMed 
Article 

Google Scholar 
49.Oliveros, J. C. Venny. An interactive tool for comparing lists with Venn’s diagrams. 2007–2015 http://bioinfogp.cnb.csic.es/tools/venny/index.html (2016).50.Kolde, R. Pheatmap: Pretty Heatmaps. R Package Version 1.0.12. https://CRANR-project.org/package=pheatmap (2019).51.Livak, K. J. & Schmittgen, T. D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods 25, 402–408 (2001).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
52.López-Ibáñez, J., Pazos, F. & Chagoyen, M. MBROLE 2.0-functional enrichment of chemical compounds. Nucleic Acids Res. 44, W201–W204 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
53.Peleg, Z. & Blumwald, E. Hormone balance and abiotic stress tolerance in crop plants. Curr. Opin. Plant Biol. 14, 290–295 (2011).CAS 
PubMed 
Article 

Google Scholar 
54.Kay, P., Groszmann, M., Ross, J., Parish, R. & Swain, S. Modifications of a conserved regulatory network involving INDEHISCENT controls multiple aspects of reproductive tissue development in Arabidopsis. N. Phytol. 197, 73–87 (2013).CAS 
Article 

Google Scholar 
55.Ditengou, F. A. et al. Characterization of auxin transporter PIN 6 plasma membrane targeting reveals a function for PIN 6 in plant bolting. N. Phytol. 217, 1610–1624 (2018).CAS 
Article 

Google Scholar 
56.Ogawa, M. et al. Gibberellin biosynthesis and response during Arabidopsis seed germination. Plant Cell 15, 1591–1604 (2003).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
57.Hedden, P. & Thomas, S. G. Gibberellin biosynthesis and its regulation. Biochem. J. 444, 11–25 (2012).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
58.Tyler, L. et al. DELLA proteins and gibberellin-regulated seed germination and floral development in Arabidopsis. Plant Physiol. 135, 1008–1019 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
59.Silverstone, A. L., Ciampaglio, C. N. & Sun, T. P. The Arabidopsis RGA gene encodes a transcriptional regulator repressing the gibberellin signal transduction pathway. Plant Cell 10, 155–169 (1998).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
60.Olszewski, N., Sun, T. P. & Gubler, F. Gibberellin signaling: biosynthesis, catabolism, and response pathways. Plant Cell 14, S61–S80 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
61.Middleton, A. M. et al. Mathematical modeling elucidates the role of transcriptional feedback in gibberellin signaling. Proc. Natl Acad. Sci. USA 109, 7571–7576 (2012).CAS 
PubMed 
Article 

Google Scholar 
62.Bévort, M. & Leffers, H. Down regulation of ribosomal protein mRNAs during neuronal differentiation of human NTERA2 cells. Differentiation 66, 81–92 (2000).PubMed 
Article 

Google Scholar 
63.Brothers, M. & Rine, J. Mutations in the PCNA DNA polymerase clamp of Saccharomyces cerevisiae reveal complexities of the cell cycle and ploidy on heterochromatin assembly. Genetics 213, 449–463 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
64.Li, C., Potuschak, T., Colón-Carmona, A., Gutiérrez, R. A. & Doerner, P. Arabidopsis TCP20 links regulation of growth and cell division control pathways. Proc. Natl Acad. Sci. USA 102, 12978–12983 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
65.Ray, S. & Pollard, J. W. KLF15 negatively regulates estrogen-induced epithelial cell proliferation by inhibition of DNA replication licensing. Proc. Natl Acad. Sci. USA 109, E1334–E1343 (2012).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
66.Halim, V., Vess, A., Scheel, D. & Rosahl, S. The role of salicylic acid and jasmonic acid in pathogen defence. Plant Biol. 8, 307–313 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
67.Zhu, F. et al. Salicylic acid and jasmonic acid are essential for systemic resistance against tobacco mosaic virus in Nicotiana benthamiana. Mol. Plant. Microbe Interact. 27, 567–577 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
68.Caarls, L., Pieterse, C. M., & Van Wees, S. How salicylic acid takes transcriptional control over jasmonic acid signaling. Front. Plant Sci. 6, 170 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
69.Wang, D., Weaver, N. D., Kesarwani, M. & Dong, X. Induction of protein secretory pathway is required for systemic acquired resistance. Science 308, 1036–1040 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
70.Checker, V. G., Kushwaha, H. R., Kumari, P. & Yadav, S. Role of phytohormones in plant defense: signaling and cross talk in Molecular aspects of plant-pathogen interaction (eds Singh, A. & Singh, I.) 159–184 (Springer, 2018).71.Niki, T., Mitsuhara, I., Seo, S., Ohtsubo, N. & Ohashi, Y. Antagonistic effect of salicylic acid and jasmonic acid on the expression of pathogenesis-related (PR) protein genes in wounded mature tobacco leaves. Plant Cell Physiol. 39, 500–507 (1998).CAS 
Article 

Google Scholar 
72.Shim, J. S. et al. AtMYB44 regulates WRKY70 expression and modulates antagonistic interaction between salicylic acid and jasmonic acid signaling. Plant J. 73, 483–495 (2013).CAS 
PubMed 
Article 

Google Scholar 
73.Romano, A. & Conway, T. Evolution of carbohydrate metabolic pathways. Res. Microbiol. 147, 448–455 (1996).CAS 
PubMed 
Article 

Google Scholar 
74.Akram, M. Citric acid cycle and role of its intermediates in metabolism. Cell Biochem. Biophys. 68, 475–478 (2014).CAS 
PubMed 
Article 

Google Scholar 
75.Plaxton, W. C. The organization and regulation of plant glycolysis. Annu. Rev. Plant Biol. 47, 185–214 (1996).CAS 
Article 

Google Scholar 
76.Montal, E. D. et al. PEPCK coordinates the regulation of central carbon metabolism to promote cancer cell growth. Mol. Cell 60, 571–583 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
77.Yang, J., Kalhan, S. C. & Hanson, R. W. What is the metabolic role of phosphoenolpyruvate carboxykinase? J. Biol. Chem. 284, 27025–27029 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
78.Huang, Y.-X. et al. Phosphoenolpyruvate carboxykinase (PEPCK) deficiency affects the germination, growth and fruit sugar content in tomato (Solanum lycopersicum L.). Plant Physiol. Biochem. 96, 417–425 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
79.Malone, S. et al. Phospho enol pyruvate carboxykinase in Arabidopsis: changes in gene expression, protein and activity during vegetative and reproductive development. Plant Cell Physiol. 48, 441–450 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
80.Murray, D. R. Nutritive role of seedcoats in developing legume seeds. Am. J. Bot. 74, 1122–1137 (1987).CAS 
Article 

Google Scholar 
81.Famiani, F. et al. Phosphoenolpyruvate carboxykinase and its potential role in the catabolism of organic acids in the flesh of soft fruit during ripening. J. Exp. Bot. 56, 2959–2969 (2005).CAS 
PubMed 
Article 

Google Scholar 
82.Osorio, S. et al. Alteration of the interconversion of pyruvate and malate in the plastid or cytosol of ripening tomato fruit invokes diverse consequences on sugar but similar effects on cellular organic acid, metabolism, and transitory starch accumulation. Plant Physiol. 161, 628–643 (2013).CAS 
PubMed 
Article 

Google Scholar 
83.Yuan, H., Zhang, J., Nageswaran, D. & Li, L. Carotenoid metabolism and regulation in horticultural crops. Hortic. Res. 2, 1–11 (2015).Article 
CAS 

Google Scholar 
84.Borghi, M. & Fernie, A. R. Floral metabolism of sugars and amino acids: implications for pollinators’ preferences and seed and fruit set. Plant Physiol. 175, 1510–1524 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
85.Sagawa, J. M. et al. An R2R3-MYB transcription factor regulates carotenoid pigmentation in Mimulus lewisii flowers. N. Phytol. 209, 1049–1057 (2016).CAS 
Article 

Google Scholar 
86.Tadmor, Y. et al. Genetics of flavonoid, carotenoid, and chlorophyll pigments in melon fruit rinds. J. Agric. Food Chem. 58, 10722–10728 (2010).CAS 
PubMed 
Article 

Google Scholar 
87.Chen, H. et al. A knockdown mutation of YELLOW-GREEN LEAF2 blocks chlorophyll biosynthesis in rice. Plant Cell Rep. 32, 1855–1867 (2013).CAS 
PubMed 
Article 

Google Scholar 
88.Grotewold, E. The genetics and biochemistry of floral pigments. Annu. Rev. Plant Biol. 57, 761–780 (2006).CAS 
PubMed 
Article 

Google Scholar 
89.Bennett, R. N. & Wallsgrove, R. M. Secondary metabolites in plant defence mechanisms. N. Phytol. 127, 617–633 (1994).CAS 
Article 

Google Scholar 
90.Erb, M. & Kliebenstein, D. J. Plant secondary metabolites as defenses, regulators, and primary metabolites: the blurred functional trichotomy. Plant Physiol. 184, 39–52 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
91.Falcone Ferreyra, M. L., Rius, S. & Casati, P. Flavonoids: biosynthesis, biological functions, and biotechnological applications. Front. Plant Sci. 3, 222 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
92.Hichri, I. et al. Recent advances in the transcriptional regulation of the flavonoid biosynthetic pathway. J. Exp. Bot. 62, 2465–2483 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
93.Jaakola, L. & Hohtola, A. Effect of latitude on flavonoid biosynthesis in plants. Plant Cell Environ. 33, 1239–1247 (2010).CAS 
PubMed 

Google Scholar 
94.Chomicki, G. et al. The velamen protects photosynthetic orchid roots against UV‐B damage, and a large dated phylogeny implies multiple gains and losses of this function during the Cenozoic. N. Phytol. 205, 1330–1341 (2015).CAS 
Article 

Google Scholar 
95.Zhang, Y., Feng, L., Jiang, H., Zhang, Y. & Zhang, S. Different proteome profiles between male and female Populus cathayana exposed to UV-B radiation. Front. Plant Sci. 8, 320 (2017).PubMed 
PubMed Central 

Google Scholar 
96.Hideg, É., Jansen, M. A. & Strid, Å. UV-B exposure, ROS, and stress: inseparable companions or loosely linked associates? Trends Plant Sci. 18, 107–115 (2013).CAS 
PubMed 
Article 

Google Scholar 
97.Kataria, S., Jajoo, A. & Guruprasad, K. N. Impact of increasing Ultraviolet-B (UV-B) radiation on photosynthetic processes. J. Photochem. Photobiol. B: Biol. 137, 55–66 (2014).CAS 
Article 

Google Scholar  More

Site informationThe study was conducted at the Kenauk Institute, an environmental research site, in western Quebec in July 2018, 2019 and 2020. All surveys occurred between 21h30 and 00h00 at night, with location and time of day randomized for each date of testing. Testing did not occur during inclement weather (rain or winds above 10 km/h). In 2018, during an initial field season, we surveyed bat populations using a traditional method (transect-based surveys) to determine which species were present. Six transects lasting 1.5 h each were laid out, and surveyed three times per season; three transects were located in open-canopy areas, and three were located in rugged, closed-canopy areas. Every 200 m, a flag marked a sampling point where we completed a 2-min static inventory using an Anabat SD2 (Titley Scientific, Columbia, MO). In this pilot study used to develop the main study, we observed all eight species known in Quebec, including the eastern red bat (Lasiurus borealis; 0.005 passes detected per minute in open-canopy habitat; 0.001 in closed-canopy), hoary bat (Lasiurus cinereus; 0.002 passes detected per minute in open-canopy habitat; 0.006 in closed-canopy) and tri-coloured bats (Perimyotis subflavus; 0.026 passes detected per minute in open-canopy habitat; none in closed-canopy). Species in the Eptesicus fuscus/Lasionycteris noctivagans acoustic complex were the most abundant (0.075 passes detected per minute in open-canopy habitat; 0.018 in closed-canopy) followed by Myotis species (Myotis leibii, Myotis septentrionalis, Myotis lucifugus; 0.075 passes detected per minute in open-canopy habitat; none in closed-canopy). Due to small sample sizes per species and because manual identification using spectrographic analyses can be unreliable for the differentiation of some bat species22, we pooled several bat species that had similar spectrograms into complexes. We pooled the big brown bat (Eptesicus fuscus) and the silver-haired bat (Lasionycteris noctivagans), and the Myotis species: little brown bat (Myotis lucifugus), northern long-eared Myotis (M. septentrionalis), and eastern small-footed bat (M. leibii)22. Therefore, these species are grouped together in analyses to minimize identification errors22. The big brown bat and silver-haired bat form the EPNO complex whereas the Myotis species form the MYSP complex. We identified to species the hoary bat (LACI), red bat (LABO), and tri-coloured bat (PESU)22. We identified bat passes visually using the output from the Anabat in the Anabat Insight software17,23.Detection efficiencyBecause total bat passes per minute were seven times higher in open-canopy habitats than in closed-canopy habitats, in 2019 we focused our surveying efforts in relatively open habitats. The Anabat (420 g) is too large to attach to a drone, thus in 2019 and 2020, we used Echometer Touch bat detectors (20 g; Wildlife Acoustics, Maynard, MA), commercially available and inexpensive detectors, attached to iPod 7 s (88 g; Apple Inc., Cupertino, CA). We do not directly compare between surveys done with the Anabat and the Echometer Touch, but merely used the 2018 Anabat surveys as a guide for expected bat species and distributions in 2019 and 2020. The UAV used was a commercially available Phantom 4 quadcopter from DJI (1.3 kg, DJI Technology Co. Inc., Shenzhen, China). To reduce sound interference from the drone, which could reduce the detection range of the instrument, we placed a 2-in. Sonoflat acoustic foam (Auralex, Indianapolis, IN) divider between the recorder and the drone, as recommended by past studies19,21 (Fig. 1).Figure 1Illustration of the three phases of the experiment design. A photograph of the UAV setup used in Phase 2 is presented in the top right corner. The setup consists of an Echometer Touch bat detector from Wildlife Acoustics and 2-inch Sonoflat acoustic foam from Auralex attached to a DJI Phantom 4 quadcopter using zip ties. (Images by Julian Herzog, Symbolon, FontAwesome retrieved from https://commons.wikimedia.org. Picture taken by the author).Full size imageIn both 2019 and 2020, we surveyed in three phases: (1) a 5-min recording from the ground without UAV; (2) a 5-min recording while the detector was attached to the UAV using zip ties and carabiners and while the UAV was manually flown in a 10–15 m diameter circle at canopy height (5––10 m above the pilot), depending on the survey site; and (3), identically to Phase 1, a 5-min recording taken from the ground without UAV (Fig. 1). The ground recorder, used sparsely in 2019 and consistently in 2020, was 1 m above the ground during phase 2. Based on surveys in 2018, seven sites were identified as having higher relative activity and were repeatedly monitored in 2019 and 2020 for bat activity. Of the seven study sites, five were located next to bodies of water and four were located near buildings; all were located in open areas. Open spaces and bodies of water are preferred hunting grounds for most bat species18, and make for an easier and safer drone flight. An additional bat detector (Echometer Touch 2, Wildlife Acoustics, Maynard USA) was used on the ground during Phase 2 to simultaneously monitor bat passes from the air and from the ground, to indicate whether bats were present but not detected due to UAV noise interference. In 2020, ten surveys were conducted with Echometer Touch 2 recorders on (1) the UAV, (2) on the ground, and (3) at a control site > 1 km from the current site. Control sites were only used in 2020. Because different bat detectors, as well as different classification software, detect and identify bats at different rates, we do not directly compare among different detectors or software24,25. In 2019, we used the Kaleidoscope software to identify bats automatically. We removed false identifications manually. In 2020, we used the Kaleidoscope software to identify all bats automatically. We also identified all passes visually and blind to the classification from Kaleidoscope. By classifying all bats using both software and visual identification, we aimed to determine whether our results were robust to identification technique.Data were collected beyond Phase 1 if the site had a bat density above three passes per 5 min (2019: N = 24 without ground detector; N = 5 with ground detector; 2020: N = 10 with ground detector; all sample sizes refer to experiments that included Phases 2 and 3). If insufficient bat activity was recorded at a given site after a 5-min period, data collection moved on to the next site, and data from that site was excluded from any analyses. Phase 1 was done to ensure there was an established bat presence, and to maximize sampling. The length of each phase was extended to 10 min if two passes were detected by the 5-min mark of Phase 1, allowing for the collection of more data, while maintaining the time proportions of each phase. While this process, necessary logistically to obtain a sufficient sample size, could lead to more bats detected during Phase 1, there should be no impact on Phase 3 compared to Phase 2, and thus, we used Tukey tests to examine Phase 3 relative to Phase 2, as well as Phase 1 compared with both other phases26.Each drone flight was performed by two field technicians: a pilot and an assistant. The UAV pilot held a basic operations pilot certificate for a small remotely-piloted aircraft system, visual line-of-sight (certificate number PC1917023611) in accordance with federal regulations enforced by Transport Canada. The assistant held the bat detector during Phases 1 and 3. During Phase 2, the assistant acted as the drone’s elevated launching and landing pad as the additional equipment obstructing the UAV’s landing gear. For take-off, they held the UAV upright above their head and gradually let go as the UAV gained altitude. For landing, the pilot gradually decreased the altitude of the drone until the landing gear was safely grasped by the assistant, who then held the UAV above their head until the propellers stopped moving. All methods were carried out in accordance with the guidelines of the Canadian Council for Animal Care. All experimental protocols were approved by McGill University animal care committee under protocol 2015-7599 and complied with the ARRIVE guidelines for animals.Statistical analyses were conducted using R 3.6.0 base package26. Generalized linear models (glm, Poisson distribution) were performed to determine the effect of phase (i.e., 1, 2, and 3) and detector location (detector on the UAV or on the ground) on the total number of bat passes. Tukey tests were then used to determine what phases and locations were significantly different from one another. To assess interspecific variation in detectability, the difference between the mean detection rate for Phase 1 and 3 and the detection rate in Phase 2 were calculated by species for each survey. A glm was then performed on the difference in detectability by species ([Average of Phases 1 and 3 − Average of Phase 2]–Species). Species were divided into four categories: MYSP (Myotis species complex), EPNO (big brown bat/silver-haired bat complex), LABO (eastern red bat), and LACI (hoary bat). No tri-coloured bats were detected, and are therefore absent from analyses. Detection phases were also divided into four categories in relation to the UAV flight: Phase 1 (pre-flight), Phase 2 from UAV-based detection (during flight), Phase 2 from ground-based detection (ground), and Phase 3 (post-flight).Detection capacityTo estimate the degree to which technological limitations affected the results gathered during the first experiment, a second experiment was conducted to estimate the impact of propeller-noise interference on the range of the bat detector. An Audio Generator SGA-8200 (Circuit-Test, Burnaby, Canada), connected to an Ultra Sound Advice S55/6 amplifier and loudspeaker (Ultra Sound Advice, London, UK) set to broadcast a 40 kHz sine wave at 40 dB SPLA @ 1 m, the highest dB setting, was used to replicate the high amplitude ultrasound reached by most bat species during their echolocation calls22. The Echometer Touch bat detector was moved away from the speaker along a measuring tape until the ultrasonic frequency could no longer be detected by the microphone. The procedure was then repeated with the detector attached to the flying UAV. As ambient sound perception cannot be evaluated when the microphone is attached to the UAV, the spectrogram on the Echometer Touch cellphone app (Wildlife Acoustics) connected to the detector was recorded with the screen video recording feature of the iPod 7 (Apple). These recordings were taken as the drone and bat detector were flown slowly along the ground to three distances (10 m, 15 m, 20 m) away from the ultrasound generator to better approximate the detection range. The videos were later visually assessed qualitatively by estimating the distance at which the signal from the speaker could no longer be distinguished from the noise interference of the drone.To quantify the spectral overlap of the drone with echolocation pulses, a spectral analysis of three 15 s recordings were performed using Avisoft SASLab Pro 4.40 (Avisoft Bioacoustics, Berlin Germany). These recordings included the drone flying, the drone motors running without propellers attached, and the ambient noise from the same location and time (control). Recordings were saved as 16 bit WAV files sampled at 256 kc/s and were normalized to 90% in SASLab Pro prior to parameterization. Spectrographs of those normalized recordings were generated using a Fast Fourier Transform length of 512 points, with a frame size of 100% and 75% overlap of Hann windows. This achieved a frequency resolution of 500 Hz and temporal resolution of 0.5 ms. Frequencies where noise was concentrated are evident from these spectrographs, but were confirmed by generating Logarithmic Power Spectra from each recording using Hann windowing achieving frequency resolution of 0.061 Hz. Noise is described at frequencies where the relative sound pressure level exceeded − 80 dB in those Power Spectra. More

1.Seebens, H. et al. No saturation in the accumulation of alien species worldwide. Nat. Commun. 8, 14435 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
2.Diagne, C. et al. InvaCost, a public database of the economic costs of biological invasions worldwide. Sci. Data 7, 277 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
3.Bradshaw, C. J. A. et al. Massive yet grossly underestimated global costs of invasive insects. Nat. Commun. 7, 12986 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
4.Meyerson, L. A. & Reaser, J. K. Biosecurity: moving toward a comprehensive approach. Bioscience 52, 593 (2002).Article 

Google Scholar 
5.Carvajal-Yepes, M. et al. A global surveillance system for crop diseases. Science 364, 1237–1239 (2019).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
6.Torres, A., David, M. & Bowman, Q. Risk management of international trade: emergency preparedness. Rev. Sci. Tech. Off. Int. Épizooties 21, 493–496 (2002).CAS 
Article 

Google Scholar 
7.Ricciardi, A. et al. Invasion science: a horizon scan of emerging challenges and opportunities. Trends Ecol. Evol. 32, 464–474 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
8.Giovani, B. et al. Science diplomacy for plant health. Nat. Plants 6, 902–905 (2020).PubMed 
Article 
PubMed Central 

Google Scholar 
9.Reaser, J. K. et al. The early detection of and rapid response (EDRR) to invasive species: a conceptual framework and federal capacities assessment. Biol. Invasions 22, 1–19 (2020).Article 

Google Scholar 
10.Delaney, D. G., Sperling, C. D., Adams, C. S. & Leung, B. Marine invasive species: validation of citizen science and implications for national monitoring networks. Biol. Invasions 10, 117–128 (2008).Article 

Google Scholar 
11.Crall, A. W. et al. Improving and integrating data on invasive species collected by citizen scientists. Biol. Invasions 12, 3419–3428 (2010).Article 

Google Scholar 
12.Maistrello, L. et al. Tracking the spread of sneaking aliens by integrating crowdsourcing and spatial modeling: the Italian invasion of halyomorpha halys. Bioscience https://doi.org/10.1093/biosci/biy112 (2018).Article 

Google Scholar 
13.Lepczyk, C. A., Boyle, O. D., Vargo, T. L. V. & Noss, R. F. Handbook of Citizen Science in Ecology and Conservation (University of California Press, Oakland, 2020).14.Devorshak, C. Plant pest risk analysis: concepts and applications. (CAB International, Wallingford, 2012).15.IPCC. Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part A: Global and Sectoral Aspects. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change Cambridge University Press, Cambridge and New York, 2014).16.Hijmans, R. J. & Graham, C. H. The ability of climate envelope models to predict the effect of climate change on species distributions. Glob. Change Biol. 12, 2272–2281 (2006).ADS 
Article 

Google Scholar 
17.Broennimann, O. & Guisan, A. Predicting current and future biological invasions: both native and invaded ranges matter. Biol. Lett. 4, 585–589 (2008).PubMed 
PubMed Central 
Article 

Google Scholar 
18.Godefroid, M., Meurisse, N., Groenen, F., Kerdelhué, C. & Rossi, J.-P. Current and future distribution of the invasive oak processionary moth. Biol. Invasions 22, 523–534 (2020).Article 

Google Scholar 
19.Crall, A. W. et al. Citizen science contributes to our knowledge of invasive plant species distributions. Biol. Invasions 17, 2415–2427 (2015).Article 

Google Scholar 
20.Petrovan, S. O., Vale, C. G. & Sillero, N. Using citizen science in road surveys for large-scale amphibian monitoring: are biased data representative for species distribution?. Biodivers. Conserv. 29, 1767–1781 (2020).Article 

Google Scholar 
21.Hannah, L. J. Climate Change Biology (Academic Press, 2015).
Google Scholar 
22.Buisson, L., Thuiller, W., Casajus, N., Lek, S. & Grenouillet, G. Uncertainty in ensemble forecasting of species distribution. Glob. Change Biol. 16, 1145–1157 (2010).ADS 
Article 

Google Scholar 
23.Porfirio, L. L. et al. Improving the use of species distribution models in conservation planning and management under climate change. PLoS ONE 9, e113749 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
24.Hamilton, G. C., Ahn, J. J., Bu, W., Leskey, T. C., Nielsen, A. L., Park, Y.-L., Rabitsch, W. & Hoelmer, K.A. Halyomorpha halys (Stål). In Invasive stink bugs and related species (Pentatomoidea): biology, higher systematics, semiochemistry, and management (ed McPherson, J. E.) 243–292 (CRC Press, Taylor & Francis, Boca Raton, 2018).25.Bergmann, E. J., Venugopal, P. D., Martinson, H. M., Raupp, M. J. & Shrewsbury, P. M. Host plant use by the invasive Halyomorpha halys (Stål) on woody ornamental trees and shrubs. PLoS ONE 11, e0149975 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
26.Gapon, D. A. First records of the brown marmorated stink bug Halyomorpha halys (Stål, 1855) (Heteroptera, Pentatomidae) in Russia, Abkhazia, and Georgia. Entomol. Rev. 96, 1086–1088 (2016).Article 

Google Scholar 
27.Faúndez, E. I. & Rider, D. A. The brown marmorated stink bug Halyomorpha halys (Stål, 1855) (Heteroptera: Pentatomidae) in Chile. Arq. Entomolóxicos 17, 305–307 (2017).
Google Scholar 
28.McPherson, J. E., ed. Invasive stink bugs and related species (Pentatomoidea): biology, higher systematics, semiochemistry, and management (CRC Press, Taylor & Francis, Boca Raton, 2018).29.Maistrello, L. et al. Halyomorpha halys in Italy: first results of field monitoring in fruit orchards. Integr. Prot. Fruit Crops IOBC-WPRS Bull. 112, 1–5 (2016).
Google Scholar 
30.Bariselli, M., Bugiani, R. & Maistrello, L. Distribution and damage caused by Halyomorpha halys in Italy. EPPO Bull. 46, 332–334 (2016).Article 

Google Scholar 
31.Zhu, G., Bu, W., Gao, Y. & Liu, G. Potential geographic distribution of brown marmorated stink bug invasion (Halyomorpha halys). PLoS ONE 7, e31246 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
32.Kriticos, D. J. et al. The potential global distribution of the brown marmorated stink bug, Halyomorpha halys, a critical threat to plant biosecurity. J. Pest Sci. 90, 1033–1043 (2017).Article 

Google Scholar 
33.Kistner, E. J. Climate change impacts on the potential distribution and abundance of the brown marmorated stink bug (Hemiptera: Pentatomidae) with special reference to North America and Europe. Environ. Entomol. 46, 1212–1224 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
34.R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, Vienna, 2021).35.Vaclavik, T., Kanaskie, A., Hansen, E. M., Ohmann, J. L. & Meentemeyer, R. K. Predicting potential and actual distribution of sudden oak death in Oregon: prioritizing landscape contexts for early detection and eradication of disease outbreaks. For. Ecol. Manag. 260, 1026–1035 (2010).Article 

Google Scholar 
36.Lobo, J. M., Jiménez-Valverde, A. & Hortal, J. The uncertain nature of absences and their importance in species distribution modelling. Ecography 33, 103–114 (2010).Article 

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

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

Google Scholar 
39.Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231–259 (2006).Article 

Google Scholar 
40.Elith, J. et al. Novel methods improve prediction of species’ distributions from occurrence data. Ecography 29, 129–151 (2006).Article 

Google Scholar 
41.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 
42.Wu, T. et al. The Beijing climate center climate system model (BCC-CSM): the main progress from CMIP5 to CMIP6. Geosci. Model Dev. 12, 1573–1600 (2019).ADS 
Article 

Google Scholar 
43.Voldoire, A. et al. Evaluation of CMIP6 DECK experiments With CNRM-CM6-1. J. Adv. Model. Earth Syst. 11, 2177–2213 (2019).ADS 
Article 

Google Scholar 
44.Séférian, R. et al. Evaluation of CNRM earth system model, CNRM-ESM2-1: role of earth system processes in present-day and future climate. J. Adv. Model. Earth Syst. 11, 4182–4227 (2019).ADS 
Article 

Google Scholar 
45.Swart, N. C. et al. The Canadian earth system model version 5 (CanESM5.0.3). Geosci. Model Dev. 12, 4823–4873 (2019).ADS 
CAS 
Article 

Google Scholar 
46.Hajima, T. et al. Development of the MIROC-ES2L Earth system model and the evaluation of biogeochemical processes and feedbacks. Geosci. Model Dev. 13, 2197–2244 (2020).ADS 
Article 

Google Scholar 
47.Tatebe, H. et al. Description and basic evaluation of simulated mean state, internal variability, and climate sensitivity in MIROC6. Geosci. Model Dev. 12, 2727–2765 (2019).ADS 
CAS 
Article 

Google Scholar 
48.Meinshausen, M. et al. The shared socio-economic pathway (SSP) greenhouse gas concentrations and their extensions to 2500. Geosci. Model Dev. 13, 3571–3605 (2020).ADS 
CAS 
Article 

Google Scholar 
49.Guisan, A., Thuiller, W. & Zimmermann, N. E. Habitat Suitability and Distribution Models with Applications in R (Cambridge University Press, 2017).Book 

Google Scholar 
50.Broennimann, O. et al. Evidence of climatic niche shift during biological invasion. Ecol. Lett. 10, 701–709 (2007).CAS 
PubMed 
Article 
PubMed Central 

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

Google Scholar 
52.Varela, S., Anderson, R. P., García-Valdés, R. & Fernández-González, F. Environmental filters reduce the effects of sampling bias and improve predictions of ecological niche models. Ecography https://doi.org/10.1111/j.1600-0587.2013.00441.x (2014).Article 

Google Scholar 
53.VanDerWal, J., Shoo, L. P., Graham, C. & Williams, S. E. Selecting pseudo-absence data for presence-only distribution modeling: How far should you stray from what you know?. Ecol. Model. 220, 589–594 (2009).Article 

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

Google Scholar 
55.Legendre, P. & Legendre, L. Numerical Ecology (Elsevier, 2012).MATH 

Google Scholar 
56.Godefroid, M., Cruaud, A., Streito, J.-C., Rasplus, J.-Y. & Rossi, J.-P. Xylella fastidiosa: climate suitability of European continent. Sci. Rep. 9, 8844 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
57.Vollering, J., Halvorsen, R. & Mazzoni, S. The MIAmaxent R package: variable transformation and model selection for species distribution models. Ecol. Evol. 9, 12051–12068 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
58.Mazzoni, S., Halvorsen, R. & Bakkestuen, V. MIAT: modular R-wrappers for flexible implementation of MaxEnt distribution modelling. Ecol. Inform. 30, 215–221 (2015).Article 

Google Scholar 
59.Elith, J., Kearney, M. & Phillips, S. The art of modelling range-shifting species: the art of modelling range-shifting species. Methods Ecol. Evol. 1, 330–342 (2010).Article 

Google Scholar 
60.Halvorsen, R., Mazzoni, S., Bryn, A. & Bakkestuen, V. Opportunities for improved distribution modelling practice via a strict maximum likelihood interpretation of MaxEnt. Ecography 38, 172–183 (2015).Article 

Google Scholar 
61.Halvorsen, R. A strict maximum likelihood explanation of MaxEnt, and some implications for distribution modelling. Sommerfeltia 36, 1–132 (2013).Article 

Google Scholar 
62.Boyce, M. S., Vernier, P. R., Nielsen, S. E. & Schmiegelow, F. K. A. Evaluating resource selection functions. Ecol. Model. 157, 281–300 (2002).Article 

Google Scholar 
63.Hirzel, A. H., Le Lay, G., Helfer, V., Randin, C. & Guisan, A. Evaluating the ability of habitat suitability models to predict species presences. Ecol. Model. 199, 142–152 (2006).Article 

Google Scholar 
64.Jiménez, L. & Soberón, J. Leaving the area under the receiving operating characteristic curve behind: an evaluation method for species distribution modeling applications based on presence-only data. Methods Ecol. Evol. https://doi.org/10.1111/2041-210X.13479 (2020).Article 

Google Scholar 
65.Chartois, M., Streito, J.-C., Pierre, E., Armand, J.-M., Gaudin, J., Rossi, J.-P. A crowdsourcing approach to track the expansion of the brown marmorated stinkbug Halyomorpha halys (Stål, 1855) in France. Biodivers. Data J. 9, e66335. https://doi.org/10.3897/BDJ.9.e66335 (2021)66.Maurel, J.-P., Blaye G., Valladares L., Roinel, E. & Cochard, P.-O. Halyomorpha halys (Stål, 1855), la punaise diabolique en France, à Toulouse (Heteroptera ; Pentatomidae). Carnets Nat. 3, 21–25 (2016).67.Cherpitel, T. & Casset, L. Halyomorpha halys (Stål, 1855), la Punaise diabolique, atteint la façade atlantique (Heteroptera Pentatomidae). L’Entomologiste 75, 59–60 (2018).
Google Scholar 
68.Pagola-Carte, S. & Zabalegui, I. D. hemípteros asiáticos nuevos para Gipuzkoa, norte de la Península Ibérica (Hemiptera: Pentatomidae, Cicadellidae). Heteropterus Rev. Entomol. 19, 355–360 (2019).
Google Scholar 
69.Streito, J. C., Rossi, J.-P., Haye, T., Hoelmer, K. & Tassus, X. La punaise diabolique à la conquête de la France. Phytoma 677, 26–30 (2014).70.Maistrello, L., Dioli, P., Bariselli, M., Mazzoli, G. L. & Giacalone-Forini, I. Citizen science and early detection of invasive species: phenology of first occurrences of Halyomorpha halys in Southern Europe. Biol. Invasions 18, 3109–3116 (2016).Article 

Google Scholar 
71.Stoeckli, S., Felber, R. & Haye, T. Current distribution and voltinism of the brown marmorated stink bug, Halyomorpha halys, in Switzerland and its response to climate change using a high-resolution CLIMEX model. Int. J. Biometeorol. https://doi.org/10.1007/s00484-020-01992-z (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
72.Leskey, T. C., Lee, D.-H., Glenn, D. M. & Morrison, W. R. Behavioral responses of the invasive Halyomorpha halys (Stål) (Hemiptera: Pentatomidae) to light-based stimuli in the laboratory and field. J. Insect Behav. 28, 674–692 (2015).Article 

Google Scholar 
73.Inkley, D. B. Characteristics of home invasion by the brown marmorated stink bug (Hemiptera: Pentatomidae). J. Entomol. Sci. 47, 125–130 (2012).Article 

Google Scholar 
74.Cambridge, J., Payenski, A. & Hamilton, G. C. The distribution of overwintering brown marmorated stink bugs (Hemiptera: Pentatomidae) in college dormitories. Fla. Entomol. 98, 1257–1259 (2015).Article 

Google Scholar 
75.Hancock, T. J., Lee, D.-H., Bergh, J. C., Morrison, W. R. & Leskey, T. C. Presence of the invasive brown marmorated stink bug Halyomorpha halys (Stål) (Hemiptera: Pentatomidae) on home exteriors during the autumn dispersal period: results generated by citizen scientists: presence of H. halys during the autumn dispersal. Agric. For. Entomol. 21, 99–108 (2019).Article 

Google Scholar 
76.Streito, J.-C., Chartois, M., Pierre, É. & Rossi, J.-P. Beware the brown marmorated stink bug!. IVES Tech Rev. Vine Wine https://doi.org/10.20870/IVES-TR.2020.3304 (2020).Article 

Google Scholar 
77.Haye, T. et al. Range expansion of the invasive brown marmorated stinkbug, Halyomorpha halys: an increasing threat to field, fruit and vegetable crops worldwide. J. Pest Sci. 88, 665–673 (2015).Article 

Google Scholar 
78.Zhu, G., Gariepy, T. D., Haye, T. & Bu, W. Patterns of niche filling and expansion across the invaded ranges of Halyomorpha halys in North America and Europe. J. Pest Sci. 90, 1045–1057 (2017).Article 

Google Scholar 
79.Shabani, F., Kumar, L. & Ahmadi, M. A comparison of absolute performance of different correlative and mechanistic species distribution models in an independent area. Ecol. Evol. 6, 5973–5986 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
80.Leskey, T. C. & Nielsen, A. L. Impact of the invasive brown marmorated stink bug in North America and Europe: history, biology, ecology, and management. Annu. Rev. Entomol. 63, 599–618 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
81.Thuiller, W., Guéguen, M., Renaud, J., Karger, D. N. & Zimmermann, N. E. Uncertainty in ensembles of global biodiversity scenarios. Nat. Commun. 10, 1446 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
82.Pearman, P. B., Guisan, A., Broennimann, O. & Randin, C. F. Niche dynamics in space and time. Trends Ecol. Evol. 23, 149–158 (2008).PubMed 
Article 
PubMed Central 

Google Scholar 
83.Jump, A. S. & Penuelas, J. Running to stand still: adaptation and the response of plants to rapid climate change. Ecol. Lett. 8, 1010–1020 (2005).Article 

Google Scholar 
84.Urvois, T., Auger-Rozenberg, M. A., Roques, A., Rossi, J. P. & Kerdelhue, C. Climate change impact on the potential geographical distribution of two invading Xylosandrus ambrosia beetles. Sci. Rep. 11, 1339 (2021).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
85.Wiens, J. J. & Graham, C. H. Niche conservatism: integrating evolution, ecology, and conservation biology. Annu. Rev. Ecol. Evol. Syst. 36, 519–539 (2005).Article 

Google Scholar  More

1.Murphy, P. G. & Lugo, A. E. Ecology of tropical dry forest. Annu. Rev. Ecol. Syst. 17, 67–88 (1986).Article 

Google Scholar 
2.Kishimoto-Yamada, K. & Itioka, T. How much have we learned about seasonality in tropical insect abundance since Wolda (1988)?. Entomol. Sci. 18, 407–419. https://doi.org/10.1111/ens.12134 (2015).Article 

Google Scholar 
3.dos Santos, J. P. D., Iserhard, C. A., Carreira, J. Y. O. & Freitas, A. V. L. Monitoring fruit-feeding butterfly assemblages in two vertical strata in seasonal Atlantic Forest: Temporal species turnover is lower in the canopy. J. Trop. Ecol. 33, 345–355. https://doi.org/10.1017/s0266467417000323 (2017).Article 

Google Scholar 
4.Bonebrake, T. C., Ponisio, L. C., Boggs, C. L. & Ehrlich, P. R. More than just indicators: A review of tropical butterfly ecology and conservation. Biol. Conser. 143, 1831–1841 (2010).Article 

Google Scholar 
5.Molleman, F. Moving beyond phenology: New directions in the study of temporal dynamics of tropical insect communities. Curr. Sci. 114, 982 (2018).Article 

Google Scholar 
6.Frith, C. B. & Frith, D. W. Seasonality of insect abundance in an Australian upland tropical rainforest. Aust. J. Ecol. 10, 237–248 (1985).Article 

Google Scholar 
7.Braby, M. Seasonal-changes in relative abundance and spatial-distribution of Australian lowland tropical satyrine butterflies. Aust. J. Zool. 43, 209–229 (1995).Article 

Google Scholar 
8.Muniz, D. G., Freitas, A. V. & Oliveira, P. S. Phenological relationships of Eunica bechina (Lepidoptera: Nymphalidae) and its host plant, Caryocar brasiliense (Caryocaraceae), in a Neotropical savanna. Stud. Neotrop. Fauna Environ. 47, 111–118 (2012).Article 

Google Scholar 
9.Wolda, H. Insect seasonality: Why?. Annu. Rev. Ecol. Syst. 19, 1–18 (1988).Article 

Google Scholar 
10.Yonow, T. et al. Modelling the population dynamics of the Queensland fruit fly, Bactrocera (Dacus) tryoni: A cohort-based approach incorporating the effects of weather. Ecol. Model. 173, 9–30. https://doi.org/10.1016/s0304-3800(03)00306-5 (2004).Article 

Google Scholar 
11.Baker, R. et al. Bactrocera dorsalis pest report to support ranking of EU candidate priority pests. EFSA https://doi.org/10.5281/zenodo.2786921 (2019).12.Valtonen, A. et al. Tropical phenology: Bi-annual rhythms and interannual variation in an Afrotropical butterfly assemblage. Ecosphere https://doi.org/10.1890/es12-00338.1 (2013).Article 

Google Scholar 
13.Hernández, C. X. P. & Caballero, S. Z. Temporal variation in the diversity of Cantharidae (Coleoptera), in seven assemblages in tropical dry forest in Mexico. Trop. Conserv. Sci. 9, 439–464 (2016).Article 

Google Scholar 
14.Marchioro, C. A. & Foerster, L. A. Biotic factors are more important than abiotic factors in regulating the abundance of Plutella xylostella L., Southern Brazil. Rev. Bras. Entomol. 60, 328–333 (2016).Article 

Google Scholar 
15.Meats, A. The bioclimatic potential of the Queensland fruit fly, Dacus tryoni, Australia. Proc. Ecol. Soc. Aust. 11, 1–61 (1981).
Google Scholar 
16.Sutherst, R. W. & Yonow, T. The geographical distribution of the Queensland fruit fly, Bactrocera (Dacus) tryoni, in relation to climate. Aust. J. Agric. Res. 49, 935–954 (1998).Article 

Google Scholar 
17.Choudhary, J. S. et al. Potential changes in number of generations of oriental fruit fly, Bactrocera Dorsalis (Diptera: Tephritidae) on mango in India in response to climate change scenarios. J. Agrometeorol. 19, 200–206 (2017).
Google Scholar 
18.Clarke, A. R. Biology and Management of Bactrocera and Related Fruit Flies (CABI, 2019).Book 

Google Scholar 
19.Sakai, S. et al. Plant reproductive phenology over four years including an episode of general flowering in a lowland dipterocarp forest, Sarawak, Malaysia. Am. J. Bot. 86, 1414–1436 (1999).CAS 
PubMed 
Article 

Google Scholar 
20.Land, K. C., Yang, Y. & Zeng, Y. Mathematical demography. Handbook of Population 659–717 (Springer, 2005).21.Carey, J. R. & Roach, D. A. Biodemography: An Introduction to Concepts and Methods (Princeton University Press, 2020).MATH 
Book 

Google Scholar 
22.Carey, J. R. Applied Demography for Biologists: With Special Emphasis on Insects (Oxford University Press, 1993).
Google Scholar 
23.Carey, J. R. Insect biodemography. Annu. Rev. Entomol. 46, 79–110 (2001).CAS 
PubMed 
Article 

Google Scholar 
24.Southwood, T. R. E. Ecological Methods: With Particular Reference to the Study of Insect Populations. xviii + 391 (Methuen, London, 1966).25.Udevitz, M. S. & Ballachey, B. E. Estimating survival rates with age-structure data. J. Wildl. Manag. 62, 779–792 (1998).Article 

Google Scholar 
26.Müller, H. G. et al. Demographic window to aging in the wild: constructing life tables and estimating survival functions from marked individuals of unknown age. Aging Cell 3, 125–131 (2004).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
27.Zajitschek, F., Zajitschek, S. & Bonduriansky, R. Senescence in wild insects: Key questions and challenges. Funct. Ecol. 34, 26–37 (2020).Article 

Google Scholar 
28.Carey, J. R. et al. Age structure changes and extraordinary lifespan in wild medfly populations. Aging Cell 7, 426–437. https://doi.org/10.1111/j.1474-9726.2008.00390.x (2008).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
29.Rao, A. S. S. & Carey, J. R. Generalization of Carey’s equality and a theorem on stationary population. J. Math. Biol. 71, 583–594 (2015).MathSciNet 
MATH 
Article 

Google Scholar 
30.Carey, J. R. Biodemography of the Mediterranean fruit fly: Aging, longevity and adaptation in the wild. Exp. Gerontol. 46, 404–411. https://doi.org/10.1016/j.exger.2010.09.009 (2011).Article 
PubMed 

Google Scholar 
31.Muller, H. G., Wang, J. L., Yu, W., Delaigle, A. & Carey, J. R. Survival and aging in the wild via residual demography. Theor. Popul. Biol. 72, 513–522. https://doi.org/10.1016/j.tpb.2007.07.003 (2007).Article 
PubMed 
PubMed Central 
MATH 

Google Scholar 
32.Vaupel, J. Life lived and left: Carey’s equality. Demogr Res 20, 7–10. https://doi.org/10.4054/DemRes.2009.20.3 (2009).Article 

Google Scholar 
33.Carey, J. R., Papadopoulos, N. T., Papanastasiou, S., Diamantidis, A. & Nakas, C. T. Estimating changes in mean population age using the death distributions of live-captured medflies. Ecol. Entomol. 37, 359–369. https://doi.org/10.1111/j.1365-2311.2012.01372.x (2012).Article 

Google Scholar 
34.Papadopoulos, N. T. et al. Seasonality of post-capture longevity in a medically-important mosquito (Culex pipiens). Front. Ecol. Evol https://doi.org/10.3389/fevo.2016.00063 (2016).Article 

Google Scholar 
35.Behrman, E. L., Watson, S. S., O’Brien, K. R., Heschel, M. S. & Schmidt, P. S. Seasonal variation in life history traits in two Drosophila species. J Evol Biol 28, 1691–1704. https://doi.org/10.1111/jeb.12690 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
36.Drew, R. A. The tropical fruit flies (Diptera: Tephritidae: Dacinae) of the Australasian and Oceanian regions. Mem. Queensland Museum 26, 1 (1989).
Google Scholar 
37.Dominiak, B. C. Components of a systems approach for the management of Queensland fruit fly Bactrocera tryoni (Froggatt) in a post dimethoate fenthion era. Crop prot. 116, 56–67 (2019).Article 

Google Scholar 
38.Boulter, S. L., Kitching, R. L. & Howlett, B. G. Family, visitors and the weather: patterns of flowering in tropical rain forests of northern Australia. J. Ecol. 94, 369–382. https://doi.org/10.1111/j.1365-2745.2005.01084.x (2006).Article 

Google Scholar 
39.Dominiak, B. C. & Mapson, R. Revised distribution of Bactrocera tryoni in eastern Australia and effect on possible incursions of Mediterranean fruit fly: Development of Australia’s eastern trading block. J. Econ. Entomol. 110, 2459–2465 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
40.Bateman, M. Adaptations to temperature in geographic races of the Queensland fruit fly Dacus (Strumenta) tryoni. Aust. J. Zool. 15, 1141–1161 (1967).Article 

Google Scholar 
41.Bateman, M. Determinants of abundance in a population of the Queensland fruit fly. In: Southwood, T.R.E. (ed.) Insect abundance 119–131 (Blackwell Scientific Publications, London, 1968).42.Drew, R., Zalucki, M. & Hooper, G. Ecological studies of eastern Australian fruit flies (Diptera: Tephritidae) in their endemic habitat. Oecologia 64, 267–272 (1984).ADS 
CAS 
PubMed 
Article 

Google Scholar 
43.Muthuthantri, S., Maelzer, D., Zalucki, M. P. & Clarke, A. R. The seasonal phenology of Bactrocera tryoni (Froggatt) (Diptera: Tephritidae) in Queensland. Aust. J. Entomol. 49, 221–233. https://doi.org/10.1111/j.1440-6055.2010.00759.x (2010).Article 

Google Scholar 
44.Lloyd, A. C. et al. Area-wide management of fruit flies (Diptera: Tephritidae) in the Central Burnett district of Queensland. Aust. J. Crop Prot. 29, 462–469. https://doi.org/10.1016/j.cropro.2009.11.003 (2010).CAS 
Article 

Google Scholar 
45.Pritchard, G. The ecology of a natural population of Queensland fruit fly, Dacus tryoni III. The maturation of female flies in relation to temperature. Aust. J. Zool. 18, 77–89 (1970).Article 

Google Scholar 
46.Clarke, A. R., Merkel, K., Hulthen, A. D. & Schwarzmueller, F. Bactrocera tryoni (Froggatt) (Diptera: Tephritidae) overwintering: an overview. Aust. Entomol. 58, 3–8. https://doi.org/10.1111/aen.12369 (2019).Article 

Google Scholar 
47.Merkel, K. et al. Temperature effects on “overwintering” phenology of a polyphagous, tropical fruit fly (Tephritidae) at the subtropical/temperate interface. J. Appl. Entomol. 143, 754–765 (2019).CAS 
Article 

Google Scholar 
48.Raghu, S., Clarke, A. R., Drew, R. A. & Hulsman, K. Impact of habitat modification on the distribution and abundance of fruit flies (Diptera: Tephritidae) in Southeast Queensland. Popul. Ecol. 42, 153–160 (2000).Article 

Google Scholar 
49.Novotny, V., Clarke, A. R., Drew, R. A., Balagawi, S. & Clifford, B. Host specialization and species richness of fruit flies (Diptera: Tephritidae) in a New Guinea rain forest. J. Trop. Ecol. 21, 67–77 (2005).Article 

Google Scholar 
50.Fletcher, B. Temperature-regulated changes in the ovaries of overwintering females of the Queensland Fruit Fly, Dacus tryoni. Aust. J. Zool. 23, 91–102 (1975).Article 

Google Scholar 
51.Meats, A. & Fay, H. The effect of acclimation on mating frequency and mating competitiveness in the Queensland fruit fly, Dacus tryoni, in optimal and cool mating regimes. Physiol. Entomol. 1, 207–212 (1976).Article 

Google Scholar 
52.Balagawi, S. Comparative ecology of Bactrocera Cucumis (French) and Bactrocera Tryoni (Froggatt) (Diptera: Tephritidae)—Understanding the life history consequences of host selection and oviposition behavior. Unpublished Thesis, Griffith University (2006).53.Lee, K. P. et al. Lifespan and reproduction in Drosophila: New insights from nutritional geometry. Proc. Natl. Acad. Sci. 105, 2498–2503 (2008).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
54.Carey, J. R., Liedo, P., Müller, H.-G., Wang, J.-L. & Vaupel, J. W. Dual modes of aging in Mediterranean fruit fly females. Science 281, 996–998 (1998).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
55.Fanson, B. G. & Taylor, P. W. Protein: carbohydrate ratios explain life span patterns found in Queensland fruit fly on diets varying in yeast: Sugar ratios. Age 34, 1361–1368 (2012).CAS 
PubMed 
Article 

Google Scholar 
56.McElderry, R. M. Seasonal life history trade-offs in two leafwing butterflies: Delaying reproductive development increases life expectancy. J. Insect Physiol. 87, 30–34 (2016).CAS 
PubMed 
Article 

Google Scholar 
57.Werfel, J., Ingber, D. E. & Bar-Yam, Y. Theory and associated phenomenology for intrinsic mortality arising from natural selection. PLoS ONE 12, e0173677 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
58.Kozeretska, I. A., Serga, S. V., Koliada, A. K. & Vaiserman, A. M. Epigenetic regulation of longevity in insects. Adv. Insect Physiol. 53, 87–114 (2017).Article 

Google Scholar 
59.Meats, A. Critical periods for developmental acclimation to cold in the Queensland fruit fly. Dacus tryoni. J. Insect Physiol. 29, 943–946 (1983).Article 

Google Scholar 
60.Kumaran, N. et al. Plant-mediated female transcriptomic changes post-mating in a tephritid fruit fly, Bactrocera tryoni. Genome Biol. Evol. 10, 94–107 (2018).CAS 
PubMed 
Article 

Google Scholar 
61.Dominiak, B. C., Sundaralingam, S., Jiang, L., Jessup, A. & Barchia, I. Production levels and life history traits of mass reared Queensland fruit fly Bactrocera tryoni (Froggatt) (Diptera: Tephritidae) during 1999/2002 in Australia. Plant Prot. Q. 23, 131–135 (2008).
Google Scholar 
62.Fanson, B., Sundaralingam, S., Jiang, L., Dominiak, B. & D’arcy, G. A review of 16 years of quality control parameters at a mass-rearing facility producing Queensland fruit fly, Bactrocera tryoni. Entomol. Exp. Appl. 151, 152–159 (2014).Article 

Google Scholar 
63.Papadopoulos, N., Katsoyannos, B., Carey, J. & Kouloussis, N. Seasonal and annual occurrence of the Mediterranean fruit fly (Diptera: Tephritidae) in northern Greece. Ann. Entomol. Soc. Am. 94, 41–50 (2001).Article 

Google Scholar 
64.Brakefield, P. M. & Reitsma, N. Phenotypic plasticity, seasonal climate and the population biology of Bicyclus butterflies (Satyridae) in Malawi. Ecol. Entomol. 16, 291–303 (1991).Article 

Google Scholar 
65.Molleman, F., Zwaan, B., Brakefield, P. & Carey, J. Extraordinary long life spans in fruit-feeding butterflies can provide window on evolution of life span and aging. Exp. Gerontol. 42, 472–482 (2007).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
66.Denlinger, D. L. Dormancy in tropical insects. Ann. Rev. Entomol 31, 239–264 (1986).CAS 
Article 

Google Scholar 
67.Canzano, A. A., Jones, R. E. & Seymour, J. E. Diapause termination in two species of tropical butterfly, Euploea core (Cramer) and Euploea sylvester (Fabricius) (Lepidoptera: Nymphalidae). Aust. J. Entomol 42, 352–356 (2003).Article 

Google Scholar 
68.Lankinen, P. & Forsman, P. Independence of genetic geographical variation between photoperiodic diapause, circadian eclosion rhythm, and Thr-Gly repeat region of the period gene in Drosophila littoralis. J. Biol. Rhythms 21, 3–12 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
69.Kouloussis, N. A. et al. Seasonal trends in Ceratitis capitata reproductive potential derived from live-caught females in Greece. Entomol. Exp. Appl. 140, 181–188. https://doi.org/10.1111/j.1570-7458.2011.01154.x (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
70.Kouloussis, N. A. et al. Life table assay of field-caught Mediterranean fruit flies, Ceratitis capitata, reveals age bias. Entomol. Exp. Appl. 132, 172–181. https://doi.org/10.1111/j.1570-7458.2009.00879.x (2009).Article 
PubMed 
PubMed Central 

Google Scholar 
71.Tasnin, M. S., Silva, R., Merkel, K. & Clarke, A. R. Response of male Queensland fruit fly (Diptera: Tephritidae) to host fruit odors. J. Econ. Entomol. 113, 1888–1893 (2020).PubMed 
Article 

Google Scholar 
72.Clarke, A. R., Powell, K. S., Weldon, C. W. & Taylor, P. W. The ecology of Bactrocera tryoni (Diptera: Tephritidae): What do we know to assist pest management?. Ann. Appl. Biol. 158, 26–54 (2011).Article 

Google Scholar 
73.Chinajariyawong, A., Drew, R., Meats, A., Balagawi, S. & Vijaysegaran, S. Multiple mating by females of two Bactrocera species (Diptera: Tephritidae: Dacinae). Bull. entomol. research 100, 325 (2010).CAS 
Article 

Google Scholar 
74.Pike, N. & Meats, A. Potential for mating between Bactrocera tryoni (Froggatt) and Bactrocera neohumeralis (hardy) (Diptera: Tephritidae). Aust. J. Entomol. 41, 70–74 (2002).Article 

Google Scholar 
75.Tasnin, M. S., Merkel, K. & Clarke, A. R. Effects of advanced age on olfactory response of male and female Queensland fruit fly, Bactrocera tryoni (Froggatt) (Diptera: Tephritidae). J. Insect Physiol. 122, 104024 (2020).CAS 
PubMed 
Article 

Google Scholar 
76.Perez-Staples, D., Prabhu, V. & Taylor, P. W. Post-teneral protein feeding enhances sexual performance of Queensland fruit flies. Physiol. Entomol. 32, 225–232 (2007).Article 

Google Scholar  More

1.Beal, A. P., Kiszka, J. J., Wells, R. S. & Eirin-Lopez, J. M. The Bottlenose dolphin Epigenetic Aging Tool (BEAT): a molecular age estimation tool for small cetaceans. Front. Mar. Sci. https://doi.org/10.3389/fmars.2019.00561 (2019).2.Garde, E., Heide-Jørgensen, M. P., Hansen, S. H., Nachman, G. & Forchhammer, M. C. Age-specific growth and remarkable longevity in narwhals (Monodon monoceros) from West Greenland as estimated by aspartic acid racemization. J. Mammal. 88, 49–58 (2007).Article 

Google Scholar 
3.Matkin, C. O., Ward Testa, J., Ellis, G. M. & Saulitis, E. L. Life history and population dynamics of southern Alaska resident killer whales (Orcinus orca). Mar. Mammal. Sci. 30, 460–479 (2014).Article 

Google Scholar 
4.Olesiuk, P., Bigg, M. & Ellis, G. Life history and population dynamics of resident killer whales (Orcinus orca) in the coastal waters of British Columbia and Washington State. Report of the International Whaling Commission. Special 12, 209–243 (1990).
Google Scholar 
5.Wells, R. S. Primates and Cetaceans: Field Research and Conservation of Complex Mammalian Societies, Primatology Monographs (eds. J. Yamagiwa, & Karczmarski, L.) p. 149–172 (Springer, 2014).6.Robeck, T. R., Willis, K., Scarpuzzi, M. R. & O’Brien, J. K. Survivorship pattern inaccuracies and inappropriate anthropomorphism in scholarly pursuits of killer whale (Orcinus orca) life history: a response to Franks et al.(2016). J. Mammal. 97, 899–905 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
7.Ellis, S. et al. Analyses of ovarian activity reveal repeated evolution of post-reproductive lifespans in toothed whales. Sci. Rep. 8, 1–10 (2018).CAS 
Article 

Google Scholar 
8.Croft, D. P., Brent, L. J., Franks, D. W. & Cant, M. A. The evolution of prolonged life after reproduction. Trends Ecol. Evol. 30, 407–416 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
9.Wursig, B. & Jefferson, T. A. Methods of photo-identification for small cetaceans. Rep. Int. Whal. Comm. 12, 43–52 (1990).
Google Scholar 
10.Perrin, W. F. & Myrick, A. C. Age Determination Of Toothed Whales And Sirenians (International Whaling Commission, 1980).11.Bryden, M. Research on Dolphins (eds. Bryden, M. M. & Harrison, R. J.) p. 211–224 (Clarendon Press Oxford, 1986).12.Myrick, A. C., Yochem, P. K. & Cornell, L. H. Toward calibrating dentinal layers in captive killer whales by use of tetracycline labels. Rit Fiskid. 11, 285–296 (1988).
Google Scholar 
13.Best, P., Meÿer, M. & Lockyer, C. Killer whales in South African waters—a review of their biology. Afr. J. Mar. Sci. 32, 171–186 (2010).Article 

Google Scholar 
14.Foote, A. D., Newton, J., Piertney, S. B., Willerslev, E. & Gilbert, M. T. P. Ecological, morphological and genetic divergence of sympatric North Atlantic killer whale populations. Mol. Ecol. 18, 5207–5217 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
15.Ford, J. K. et al. Shark predation and tooth wear in a population of northeastern Pacific killer whales. Aquat. Biol. 11, 213–224 (2011).Article 

Google Scholar 
16.Hohn, A. A. & Fernandez, S. Biases in dolphin age structure due to age estimation technique. Mar. Mammal. Sci. 15, 1124–1132 (1999).Article 

Google Scholar 
17.Lockyer, C. A report on patterns of deposition of dentine and cement in teeth of pilot whales, genus Globicephala. Rep. Int. Whal. Comm. 14, 137–161 (1993).
Google Scholar 
18.Waugh, D. A., Suydam, R. S., Ortiz, J. D. & Thewissen, J. Validation of Growth Layer Group (GLG) depositional rate using daily incremental growth lines in the dentin of beluga (Delphinapterus leucas (Pallas, 1776)) teeth. PLoS ONE 13, e0190498 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
19.Sergeant, D. E. Age Determination In Odontocete Whales From Dentinal Growth Layers (Norwegian Whaling Gazette, 1959).20.Brodie, P. F. Mandibular layering in Delphinapterus leucas and age determination. Nature 221, 956–958 (1969).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
21.Goren, A. D. et al. Growth layer groups (GLGs) in the teeth of an adult belukha whale (Delphinapterus leucas) of known age: evidence for two annual layers. Mar. Mammal. Sci. 3, 14–21 (1987).Article 

Google Scholar 
22.Brodie, P. & Haulena, M. Dentinal growth layer counts of captive, known-age, mother and daughter belugas (Delphinapterus leucas): confirming two growth layer groups (GLG/2) per year; consequences for recovery and management. J Cetacean. Res Manag. 18, 23–31 (2018).
Google Scholar 
23.Brodie, P., Ramirez, K. & Haulena, M. Growth and maturity of belugas (Delphinapterus leucas) in Cumberland Sound, Canada, and in captivity: evidence for two growth layer groups (GLGs) per year in teeth. J. Cetacean Res. Manag. 13, 1–18 (2013).
Google Scholar 
24.Lockyer, C., Hohn, A. A., Doidge, D. W., Heide-Jørgensen, M. P. & Suydam, R. Age determination in belugas (Delphinapterus leucas in Belugas): a quest for validation of dentinal layering. Aquat. Mamm. 33, 293–304 (2007).Article 

Google Scholar 
25.Stewart, R., Campana, S., Jones, C. & Stewart, B. Bomb radiocarbon dating calibrates beluga (Delphinapterus leucas) age estimates. Can. J. Zool. 84, 1840–1852 (2006).Article 

Google Scholar 
26.Brodie, P. A reconsideration of aspects of growth, reproduction, and behavior of the white whale (Delphinapterus leucas), with reference to the Cumberland Sound, Baffin Island, population. J. Fish. Board Can. 28, 1309–1318 (1971).Article 

Google Scholar 
27.Brodie, P. F., Parsons, J. L. & Sergeant, D. E. Present status of the white whale (Delphinapterus leucas) in Cumberland Sound, Baffin Island.Rep. Int. Whal. Comm. 31, 579–582 (1981).
Google Scholar 
28.Robeck, T. R. et al. Reproduction, growth and development in captive beluga (Delphinapterus leucas). Zoo Biol. 24, 29–49 (2005).Article 

Google Scholar 
29.Bada, J., Brown, S. & Masters, P. Age determination of marine mammals based on aspartic acid racemization in the teeth and lens nucleus. Age Determination of Toothed Whales and Sirenians. p. 113–118 (Report of the International Whaling Commission, Special, 1980).30.George, J. C. et al. Age and growth estimates of bowhead whales (Balaena mysticetus) via aspartic acid racemization. Can. J. Zool. 77, 571–580 (1999).Article 

Google Scholar 
31.Pleskach, K. et al. Use of mass spectrometry to measure aspartic acid racemization for ageing beluga whales. Mar. Mammal. Sci. 32, 1370–1380 (2016).CAS 
Article 

Google Scholar 
32.Garde, E., Peter Heide-Jørgensen, M., Ditlevsen, S. & Hansen, S. H. Aspartic acid racemization rate in narwhal (Monodon monoceros) eye lens nuclei estimated by counting of growth layers in tusks. Polar Res. https://doi.org/10.3402/polar.v31i0.15865 (2012).33.Herman, D. P. et al. Assessing age distributions of killer whale Orcinus orca populations from the composition of endogenous fatty acids in their outer blubber layers. Mar. Ecol. Prog. Ser. 372, 289–302 (2008).CAS 
Article 

Google Scholar 
34.Herman, D. P. et al. Age determination of humpback whales Megaptera novaeangliae through blubber fatty acid compositions of biopsy samples. Mar. Ecol. Prog. Ser. 392, 277–293 (2009).Article 

Google Scholar 
35.Marcoux, M., Lesage, V., Thiemann, G. W., Iverson, S. J. & Ferguson, S. H. Age estimation of belugas, Delphinapterus leucas, using fatty acid composition: a promising method. Mar. Mammal. Sci. 31, 944–962 (2015).Article 

Google Scholar 
36.Olsen, M. T., Berube, M., Robbins, J. & Palsboll, P. J. Empirical evaluation of humpback whale telomere length estimates; quality control and factors causing variability in the singleplex and multiplex qPCR methods. BMC Genet. 13, 77 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Broer, L. et al. Meta-analysis of telomere length in 19 713 subjects reveals high heritability, stronger maternal inheritance and a paternal age effect. Eur. J. Hum. Genet. 21, 1163–1168 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
38.Dunshea, G. et al. Telomeres as age markers in vertebrate molecular ecology. Mol. Ecol. Resour. 11, 225–235 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
39.Polanowski, A. M., Robbins, J., Chandler, D. & Jarman, S. N. Epigenetic estimation of age in humpback whales. Mol. Ecol. Resour. 14, 976–987 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
40.Tanabe, A. et al. Age estimation by DNA methylation in the Antarctic minke whale. Fish. Sci. 86, 35–41 (2020).CAS 
Article 

Google Scholar 
41.Smith, Z. D. & Meissner, A. DNA methylation: roles in mammalian development. Nat. Rev. Genet. 14, 204–220 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
42.Rakyan, V. K. et al. Human aging-associated DNA hypermethylation occurs preferentially at bivalent chromatin domains. Genome Res. 20, 434–439 (2010).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
43.Teschendorff, A. E. et al. Age-dependent DNA methylation of genes that are suppressed in stem cells is a hallmark of cancer. Genome Res. 20, 440–446 (2010).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
44.Horvath, S. & Raj, K. DNA methylation-based biomarkers and the epigenetic clock theory of ageing. Nat. Rev. Genet. 19, 371–384 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
45.Field, A. E. et al. DNA methylation clocks in aging: categories, causes, and consequences. Mol. Cell 71, 882–895 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
46.Horvath, S. DNA methylation age of human tissues and cell types. Genome Biol. 14, R115 (2013).PubMed 
PubMed Central 
Article 

Google Scholar 
47.Bell, C. G. et al. DNA methylation aging clocks: challenges and recommendations. Genome Biol. 20, 1–24 (2019).Article 

Google Scholar 
48.Petkovich, D. A. et al. Using DNA methylation profiling to evaluate biological age and longevity interventions. Cell Metab. 25, 954–960.e956 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Cole, J. J. et al. Diverse interventions that extend mouse lifespan suppress shared age-associated epigenetic changes at critical gene regulatory regions. Genome Biol. 18, 1–16 (2017).Article 
CAS 

Google Scholar 
50.Wang, T. et al. Epigenetic aging signatures in mice livers are slowed by dwarfism, calorie restriction and rapamycin treatment. Genome Biol. 18, 57 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
51.Stubbs, T. M. et al. Multi-tissue DNA methylation age predictor in mouse. Genome Biol. 18, 1–14 (2017).Article 
CAS 

Google Scholar 
52.Thompson, M. J. et al. A multi-tissue full lifespan epigenetic clock for mice. Aging 10, 2832 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Meer, M. V., Podolskiy, D. I., Tyshkovskiy, A. & Gladyshev, V. N. A whole lifespan mouse multi-tissue DNA methylation clock. Elife 7, e40675 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
54.Ito, T., Teo, T. V., Evans, S. A., Neretti, N. & Sedivy, J. Regulation of cellular senescence by polycomb chromatin modifiers through distinct DNA damage- and histone methylation-dependent pathways. Cell Rep. 22, 3480–3492 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
55.St Aubin, D., Deguise, S., Richard, P., Smith, T. & Geraci, J. Hematology and plasma chemistry as indicators of health and ecological status in beluga whales, Delphinapterus leucas. Arctic 54, 317–331 (2001).56.Norman, S. A. et al. Seasonal hematology and serum chemistry of wild beluga whales (Delphinapterus leucas) in Bristol Bay, Alaska, USA. J. Wildl. Dis. 48, 21–32 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
57.Frost, K. J. & Suydam, R. S. Subsistence harvest of beluga or white whales (Delphinapterus leucas) in northern and western Alaska 1987–2006. J. Cetacea. Res. Manag. 11, 293–299 (2010).
Google Scholar 
58.Rosen, A. D. et al. DNA methylation age is accelerated in alcohol dependence. Transl. Psychiatry 8, 182 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
59.Zhang, Q. et al. Improved precision of epigenetic clock estimates across tissues and its implication for biological ageing. Genome Med. 11, 54 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
60.Gronniger, E. et al. Aging and chronic sun exposure cause distinct epigenetic changes in human skin. PLoS Genet. 6, e1000971 (2010).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
61.Doi, A. et al. Differential methylation of tissue-and cancer-specific CpG island shores distinguishes human induced pluripotent stem cells, embryonic stem cells and fibroblasts. Nat. Genet. 41, 1350–1353 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
62.Vandiver, A. R. et al. Age and sun exposure-related widespread genomic blocks of hypomethylation in nonmalignant skin. Genome Biol. 16, 1–15 (2015).CAS 
Article 

Google Scholar 
63.Li, Q. S., Sun, Y. & Wang, T. Epigenome-wide association study of Alzheimer’s disease replicates 22 differentially methylated positions and 30 differentially methylated regions. Clin. Epigenet. 12, 1–14 (2020).Article 
CAS 

Google Scholar 
64.Sun, L., Zhang, X., Wang, T., Chen, M. & Qiao, H. Association of ANK1 variants with new‑onset type 2 diabetes in a Han Chinese population from northeast China. Exp. Ther. Med. 14, 3184–3190 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
65.Luoma, L. M. & Berry, F. B. Molecular analysis of NPAS3 functional domains and variants. BMC Mol. Biol. 19, 1–19 (2018).Article 
CAS 

Google Scholar 
66.Cosgrove, D. et al. Genes influenced by MEF2C contribute to neurodevelopmental disease via gene expression changes that affect multiple types of cortical excitatory neurons. bioRxiv https://doi.org/10.1101/2019.12.16.877837 (2019).67.Decourcelle, A. et al. O-GlcNAcylation links nutrition to the epigenetic downregulation of UNC5A during colon carcinogenesis. Cancers 12, 3168 (2020).CAS 
PubMed Central 
Article 

Google Scholar 
68.Yang, T., Zhang, X.-B., Li, X.-N., Sun, M.-Z. & Gao, P.-Z. Homeobox C4 promotes hepatocellular carcinoma progression by the transactivation of Snail. Neoplasma 68, 23–30 (2020).69.Yeung, B., Law, A. & Wong, C. K. Evolution and roles of stanniocalcin. Mol. Cell. Endocrinol. 349, 272–280 (2012).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
70.Chen, C., Jamaluddin, M. S., Yan, S., Sheikh-Hamad, D. & Yao, Q. Human stanniocalcin-1 blocks TNF-α–induced monolayer permeability in human coronary artery endothelial cells. Arterioscler. Thromb. Vasc. Biol. 28, 906–912 (2008).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
71.Jourdain, E. & Karoliussen, R. Identification catalogue of Norwegian killer whales: 2007–2018. Figshare https://doi.org/10.608/m9.figshare.4205226 (2018).72.Kuningas, S., Similä, T. & Hammond, P. S. Population size, survival and reproductive rates of northern Norwegian killer whales (Orcinus orca) in 1986-2003. J. Mar. Biol. Assoc. UK 94, 1277 (2014).Article 

Google Scholar 
73.Christensen, I. Growth and reproduction of killer whales, Orcinus orca, in Norwegian coastal waters. Rep. Int. Whal. Commn 6, 253–258 (1984).
Google Scholar 
74.Jourdain, E., Vongraven, D., Bisther, A. & Karoliussen, R. First longitudinal study of seal-feeding killer whales (Orcinus orca) in Norwegian coastal waters. PLoS ONE 12, e0180099 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
75.Arneson, A. et al. A mammalian methylation array for profiling methylation levels at conserved sequences. bioRxiv https://doi.org/10.1101/2021.01.07.425637 (2021).76.Zhou, W., Triche, T. J. Jr., Laird, P. W. & Shen, H. SeSAMe: reducing artifactual detection of DNA methylation by Infinium BeadChips in genomic deletions. Nucleic Acids Res. 46, e123 (2018).PubMed 
PubMed Central 

Google Scholar 
77.Friedman, J., Hastie, T. & Tibshirani, R. Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33, 1 (2010).PubMed 
PubMed Central 
Article 

Google Scholar 
78.Shao, J. Linear model selection by cross-validation. J. Am. Stat. Assoc. 88, 486–494 (1993).Article 

Google Scholar 
79.Zhang, P. Model selection via multifold cross validation. Ann. Statist. 21, 299–313 (1993).80.Team, R. C. R.: A language and environment for statistical computing (2020). More

Persian Gulf water and sediment samples along the oil pollution continuumWater and sediment samples were collected along the circulation current of the Persian Gulf from Hormuz Island [HW (SAMN12878178) and HS (SAMN12878113)), Asaluyeh area (AW (SAMN12878179) and AS (SAMN12878114)), and Khark Island (KhW (SAMN12878180) and KhS (SAMN12878115)] (Fig. 1). Physicochemical characteristics and Ionic content of the collected samples are presented in the Supplementary Table S1. The GC-FID analyses showed high TPH and polyaromatic hydrocarbon (PAH) concentrations in the Khark sediment (KhS) (Supplementary Table S2). The GC-SimDis analysis showed that C25–C38 HCs were dominant in the KhS (~ 60%), followed by  > C40 HCs (~ 14%) (Supplementary Fig. S1). Chrysene, fluoranthene, naphthalene, benzo(a)anthracene and phenanthrene were respectively the most abundant PAHs in KhS. This pollution could originate from oil spillage due to Island airstrikes during the imposed war (1980–1988), sub-sea pipeline failures, and discharge of oily wastewater or ballast water of oil tankers (ongoing for ~ 50 years)12. The TPH of other water and sediment samples was below the detection limit of our method ( 2%) representatives were enriched in KhS (Fig. 2B). The co-presence of the orders Methanosarcinals, Alteromanadales, and Thermotogae (Petrotogales) in the KhS hints at potential oil reservoir seepage around the sampling site since these taxa are expected to be present in oil reservoirs38. The main HC pollutants in the Asalouyeh are low molecular weight aromatic compounds that mainly influence the prokaryotic population in the water column and rarely precipitate into sediments hence the similarity of AS to HS microbial composition as they both experience low pollution rates.Apart from oil-degrading Proteobacteria (e.g. Alteromonadales, Rhodobacterales, and Oceanospirillales), a diversity of sulfur/ammonia-oxidizing chemolithoautotrophic Proteobacteria were present in these sediments although at lower abundances e.g., (Acidithiobacillales (KhS 1.8%), Chromatiales (HS 1.5, AS 1.1, KhS 0.85%), Ectothiorhodospirales (HS 3.75, AS 2.3, KhS 1.7%), Halothiobacillales (KhS 2.6%), Thiotrichales (HS 1.5, AS 1.1, KhS 0.3%), Thiohalorhabdales (HS 0.7, AS 1.2, KhS 0.5%), Thiomicrospirales (KhS 1.5%)) (Fig. 2B).Sulfate-reducing bacteria (SRB) in HS comprised up to 16.2% of the community (Desulfobacterales, NB1-j, Myxococcales, Syntrophobacterales, and Thermodesulfovibrionia). Similar groups along with Desulfarculales, comprised the SRB functional guild of the AS (~ 18.9%). In comparison, Desulfuromonadales and Desulfobacterales were the SRB representatives in KhS with a total abundance of only ~ 3.3%. The lower phylogenetic diversity and community contribution of SRBs in KhS hint at the potential susceptibility of some SRBs to oil pollution or that HC degraders might outcompete them (e.g., Deferribacterales). Additionally, KhS was gravel-sized sediment (particles ≥ 4 mm diameter), whereas HS and AS samples were silt and sand-sized sediments39. The higher oxygen penetration in gravel particles of KhS hampers anaerobic metabolism of sulfate/nitrate-reducing bacteria hence their lower relative abundance in this sample (Fig. 2B).Whereas in water samples, sulfur/ammonia-oxidizing chemolithoautotrophs such as Thiomicrospirales and sulfate/nitrate-reducing bacteria such as Desulfobacterales, NB1-j, Deferribacterales, Anaerolineales, Nitrosococcales, Nitrosopumilales, and Pirellulales were present in very small quantities (lower than 0.5% in each sample).Chronic exposure to oil pollution shapes similar prokaryotic communities as oil spill eventsWe analyzed the prokaryotic community composition of 41 oil-polluted marine water metagenomes (different depths in the water column) from Norway (Trondheimsfjord), Deepwater Horizon (Gulf of Mexico), the northern part of the Gulf of Mexico (dead zone) and Coal Oil Point of Santa Barbara; together with 65 oil exposed marine sediment metagenomes (beach sand, surface sediments and deep-sea sediments) originating from DWH Sediment (Barataria Bay), Municipal Pensacola Beach (USA) and a hydrothermal vent in Guaymas Basin (Gulf of California) in comparison with the PG water and sediment samples (in total 112 datasets) (Supplementary Table S3). This extensive analysis allowed us to get a comparative overview of the impact of chronic oil pollution on the prokaryotic community composition.Hydrocarbonoclastic bacteria affiliated to Oceanospirillales, Cellvibrionales (Porticoccaceae family), and Alteromonadales40 comprised a significant proportion of the prokaryotic community in samples with higher aliphatic compounds pollution e.g. DWHW.BD3 (sampled six days after the incubation of unpolluted water with Macondo oil), DWHW.he1, and DWHW.he2 (oil-polluted water samples incubated with hexadecane), DWHW.BM1, DWHW.BM2, DWHW.OV1 and DWHW.OV2 (sampled immediately after the oil spill in the Gulf of Mexico) (Fig. 3). Samples treated with Macondo oil, hexadecane, naphthalene, phenanthrene, and those taken immediately after the oil spill in the Gulf of Mexico had a significantly lower proportion of SAR11 due to the dominance of bloom formers and potential susceptibility of SAR11 to oil pollutants (Fig. 3).Figure 3The abundance of unassembled 16S rDNA reads from unassembled metagenomes of different oil-polluted water samples (41). Row names are microbial taxa at the order level. For taxa with lower frequency, the higher taxonomic level is shown (47 taxa in total). The right-hand dendrogram represents the clustering of rows based on the Pearson correlation. Columns are the name of water samples. Samples are clustered based on Pearson correlation and the color scale on the top left represents the row Z-score. Figure was plotted using “circlize” and “ComplexHeatmap” packages in R.Full size imageFlavobacteriales and Rhodobacterales were present in relatively high abundance in almost all oil-polluted water samples except for those with recent pollution. Samples named NTW5, NTW6, NTW11, NTW12, which were incubated with MC252 oil for 32–64 days, represented similar prokaryotic composition dominating taxa that are reportedly involved in degrading recalcitrant compounds like PAHs in the middle-to-late stages of the oil degradation process (Alteromonadales, Cellvibrionales, Flavobacteriales, and Rhodobacterales). Whereas at the earlier contamination stages, samples represented a different community composition with a higher relative abundance of Oceanospirillales (e.g., NTW8, NTW9, NTW15, NTW16, and NTW17 sampled after 0–8 days incubation) (Fig. 3).The non-metric multidimensional analysis of the prokaryotic community of 106 oil-polluted water and sediment samples, together with the PG samples, is represented in Fig. 4. Water and sediment samples expectedly represented distinct community compositions. The AW sample was placed near samples treated with phenanthrene and naphthalene in the NMDS plot showing the impact of aromatic compounds on its microbial community. The KhW sample was located near NTW13 in the plot, both of which had experienced recent oil pollution.Figure 4Non-metric multidimensional scaling (NMDS) of the Persian Gulf water and sediment metagenomes along with oil-polluted marine water and sediment metagenomes based on Bray–Curtis dissimilarity of the abundance of 16S rDNA reads in unassembled metagenomes at the order level. Samples with different geographical locations are shown in different colors. PG water and sediment samples are shown in red. Water and sediment samples are displayed by triangle and square shapes, respectively. Figure was plotted using “vegan” library in R.Full size imageThe orders Oceanospirillales, Alteromonadales, and Pseudomonadales were present in relatively high abundances in all oil-polluted water samples except for HW (PG input water) and samples collected from the northern Gulf of Mexico dead zone (GOMDZ) (Fig. 3). Persian Gulf was located in the proximity of the developing oxygen minimum zone (OMZ) of the Arabian Sea that is slowly expanding towards the Gulf of Oman41. Potential water exchange with OMZ areas could be the cause of higher similarity to the GOMDZ microbial community42.Our results suggest that water samples with similar contaminants and exposure time to oil pollution enrich for similar phylogenetic diversity in their prokaryotic communities (Fig. 3). Marine prokaryotes represent vertical stratification with discrete community composition across the depth profile. According to our analyses, the prokaryotic communities of the oil-polluted areas are consistently dominated by similar taxa regardless of sampling depth or geographical location. We speculate that the high nutrient input due to crude oil intrusion into the water presumably disturbs this stratification and HC degrading microorganisms are recruited to the polluted sites where their populations flourish.The inherent heterogeneity of the sediment prokaryotic communities is retained even after exposure to oil pollution, reflected in their higher alpha diversity (Supplementary Fig. S3). However, similar taxa dominate the community in response to oil pollution (Fig. 5).Figure 5The abundance of unassembled 16S rDNA reads from unassembled metagenomes of different oil-polluted sediment samples (65). Row names are microbial taxa at the order level. For taxa with lower frequency, the higher taxonomic level is shown (77 taxa in total). The right-hand dendrogram represents the clustering of rows based on the Pearson correlation. Columns are the name of sediment samples. Samples are clustered based on Pearson correlation and the color scale on the top left represents the row Z-score. Figure was plotted using “circlize” and “ComplexHeatmap” packages in R.Full size imageIn sediment samples, Deltaproteobacteria had the highest abundance, followed by Gammaproteobacteria representatives. Ectothiorhodospirales, Rhizobiales, Desulfobacterales, Myxococcales, and Betaproteobacteriales representatives were present in almost all samples at relatively high quantities (Fig. 5). Sulfate/nitrate-reducing bacteria were major HC degraders in sediment, showing substrate specificity for anaerobic HC degradation43. Desulfobacterales and Myxococcales were ubiquitous sulfate-reducers, present in almost all oil-polluted sediment samples44. Sulfate-reducing Deltaproteobacteria play a key role in anaerobic PAH degradation, especially in sediments containing recalcitrant HC types45. Members of Rhizobiales are involved in nitrogen fixation, which accelerates the HC removal process in the sediment samples46, and therefore their abundance increase in response to oil pollution (Fig. 5).Prokaryotes involved in nitrogen/sulfur cycling of sediments are defined by factors such as trace element composition, temperature, pressure, and more importantly, depth and oxygen availability. In oil-polluted sediment samples, the simultaneous reduction of available oxygen with an accumulation of recalcitrant HCs along the depth profile complicates the organic matter removal. However, anaerobic sulfate-reducing HC degrading bacteria will cope with this complexity47. Prokaryotic communities of HS and AS samples represented similar phylogenetic diversity (Figs. 4, 5). Their prokaryotic community involved in the nitrogen and sulfur cycling resembles the community of DWHS samples. The KhS sample had a similar prokaryotic community to deeper sediment samples collected from 30 to 40 cm depth (USFS3, USFS11, and USFS12) which could be due to our sampling method using a grab sampling device.Our results show that the polluted sediments’ sampling depth (surface or subsurface) defines the dominant microbial populations. Hydrocarbon degrading microbes had the ubiquitous distribution in almost all oil-polluted water and sediment samples including Oceanospirillales, Cellvibrionales, Alteromonadales, Flavobacteriales, Pseudomonadales, and Rhodobacterales. Mentioned orders along with Ectothiorhodospirales, Rhizobiales, Desulfobacterales, Myxococcales, and Betaproteobacteriales and also representatives of Deltaproteobacteria phylum dominated in sediment samples. However, their order of frequency varies depending on the type of oil pollution present at the sampling location and the exposure time.Genome-resolved metabolic analysis of the Persian Gulf’s prokaryotic community along the pollution continuumA total of 82 metagenome-assembled genomes (MAGs) were reconstructed from six sequenced metagenomes of the PG (completeness ≥ 40% and contamination ≤ 5%). Amongst them, eight MAGs belonged to domain Archaea and 74 to domain bacteria. According to GTDB-tk assigned taxonomy (release89) (https://data.gtdb.ecogenomic.org/releases/release89/), reconstructed MAGs were affiliated to Gammaproteobacteria (36.6%), Alphaproteobacteria (12.2%), Flavobacteriaceae (9.7%), Thermoplasmatota (5%) together with some representatives of other phyla (MAG stats in Supplementary Table S4).A collection of reported enzymes involved in the degradation of different aromatic and aliphatic HCs under both aerobic and anaerobic conditions was surveyed in the annotated MAGs of this study43,48,49,50. The KEGG orthologous accession numbers (KOs) of genes involved in HC degradation were collected, and the distribution of KEGG orthologues detected at least in one MAG (n = 76 genes) is represented in Fig. 6.Figure 6Hydrocarbon degrading enzymes present in recovered MAGs from the PG water and sediment metagenomes. Row names represent the taxonomy of recovered MAGs and their completeness is provided as a bar plot on the right side. The color indicates the MAG origin. The size of dots indicates the presence or absence of each enzyme in each recovered MAG. Columns indicate the type of hydrocarbon and in the parenthesis is the name of the enzyme hydrolyzing this compound followed by its corresponding KEGG orthologous accession number. Figure was plotted using “reshape2” and “ggplot2” packages in R.Full size imageA combination of different enzymes runs the oil degradation process. Mono- or dioxygenases are the main enzymes triggering the HC degradation process under aerobic conditions. Under anaerobic conditions, degradation is mainly started by the addition of fumarate or in some cases, by carboxylation of the substrate. Therefore, bacteria containing these genes will potentially initiate the degradation process that will be continued by other heterotrophs. Enzymes such as decarboxylase, hydroxylase, dehydrogenase, hydratase, and isomerases act on the products of initiating enzymes mentioned above through a series of oxidation/reduction reactions.Various microorganisms cooperate to cleave HCs into simpler compounds that could enter common metabolic pathways. Mono- or dioxygenases which are involved in the degradation of alkane (alkane 1-monooxygenase, alkB/alkM), cyclododecane (cyclododecanone monooxygenase, cddA), Biphenyl (Biphenyl 2, 3-dioxygenase subunit alpha/beta, bphA1/A2, Biphenyl-2, 3-diol 1, 2-dioxygenase, bphC), phenol (phenol 2-monooxygenase, pheA), toluene (benzene 1, 2-dioxygenase subunit alpha/beta todC1/C2, hydroxylase component of toluene-4-monooxygenase, todE), xylene (toluate/benzoate 1,2-dioxygenase subunit alpha/beta/electron transport component, xylX/Y/Z, hydroxylase component of xylene monooxygenase, xylM) and naphthalene/phenanthrene (catechol 1,2 dioxygenase, catA, a shared enzyme between naphthalene/phenanthrene /phenol degradation) were detected in recovered MAGs of the PG.The key enzymes including Alkylsuccinate synthase (I)/(II) (assA1/A2), benzylsuccinate synthase (BssA)/benzoyl-CoA reductase (BcrA), ethylbenzene dehydrogenase (EbdA), and 6-oxo-cyclohex-1-ene-carbonyl-CoA hydrolase (BamA) that are responsible for initiating the degradation of alkane, toluene, ethylbenzene and benzoate exclusively under anaerobic conditions were not detected in reconstructed MAGs of this study. Consequently, recovered MAGs of this study are not initiating anaerobic degradation via known pathways while they have the necessary genes to continue the degradation process started by other microorganisms.The MAG KhS_63 affiliated to Immundisolibacter contained various types of mono- or dioxygenases and had the potential to degrade a diverse range of HCs such as alkane, cyclododecane, toluene, and xylene (Fig. 6). Members of this genus have been reported to degrade high molecular weight PAHs51.Lutimaribacter representatives have been isolated from seawater and reported to be capable of degrading cyclohexylacetate52. We also detected enzymes responsible for alkane, cycloalkane (even monooxygenase enzymes), and naphthalene degradation under aerobic conditions and alkane, ethylbenzene, toluene, and naphthalene degradation under anaerobic conditions in KhS_39 affiliated to this genus (Fig. 6).The MAGs KhS_15 and KhS_26 affiliated to Roseovarius had the enzymes for degrading alkane (alkane monooxygenase, aldehyde dehydrogenase), cycloalkane, naphthalene, and phenanthrene under aerobic and toluene and naphthalene under anaerobic condition. PAHs degradation has been reported for other representatives of this taxa as well53.The MAGs KhS_11 (a representative of Rhodobacteraceae) and KhS_53 (Marinobacter) had alkB/alkM, KhS_27 (GCA-2701845), KhS_29 (UBA5862) and KhS_40 (from Porticoccaceae family) had cddA, KhS_13 and KhS_21 (UBA5335) and KhS_38 (Oleibacter) had both alkB/alkM and xylM genes. They were among microbes that were initiating the degradation of alkane, cycloalkane and xylene compounds. Other MAGs recovered from Khark sediment were involved in the continuation of the degradation pathway. For example, KhS_1 was affiliated to the genus Halomonas and had different enzymes to degrade intermediate compounds. Halomonas representatives have been frequently isolated from oil-polluted environments54. The phylum Krumholzibacteria has been first introduced in 2019 and reported to contain heterotrophic nitrite reducers55. Two MAGs, KhS_5 and KhS_10, were affiliated to this phylum and contained enzymes involved in the anaerobic degradation of toluene, phenol, and naphthalene (Fig. 6).The MAGs KhS_12 and KhW_31 affiliated to the genus Flexistipes, in Deferribacterales order, were reconstructed from both KhW and KhS samples. Deferribacterales are reported to be present in the medium to high-temperature oil reservoirs with HC degradation activity and also in high-temperature oil-degrading consortia56. The type strain of this species was isolated from environments with a minimum salinity of 3% and a temperature of 45–50 °C57. The presence of this genus in KhS could be due to natural oil seepage from the seabed as PG reservoirs mainly have medium to high temperature and high salinity. Enzymes involved in the degradation of alkane, phenol, toluene and naphthalene under anaerobic conditions were present in MAGs KhS_12 and KhW_31.As mentioned earlier, Flavobacteriales are potent marine indigenous HC degraders that bloom in response to oil pollution58. Flavobacteriales affiliated MAGs (KhW_2, KhW_3, AW_21, and AW_33) were recovered from KhW and AW and mostly contained enzymes that participate in the degradation of aromatic compounds under anaerobic conditions. KhW_2 and KhW_3 also had both alkB/M (alkane monooxygenase) and xylM enzyme, which initiates the alkane and xylene bioremediation in Khark water. Among other recovered MAGs from KhW sample, KhW_18 (UBA724), KhW_24 (clade SAR86), KhW_43 (UBA3478) had alkB/M, and xylM, KhW_24 (clade SAR86) had alkB/M and cddA, and KhW_28 (from Rhodobacteraceae family) had alkB/M and pheA genes in their genome to initiate the degradation process (Fig. 6).Marinobacter (KhW_15) was another MAG reconstructed from KhW sample. This genus is one of the main cultivable genera that play a crucial role in the bioremediation of a wide range of oil derivatives in polluted marine ecosystems54.Marine Group II (MGII) and Poseidonia representatives of Thermoplasmatota that have been reported to be nitrate-reducing Archaea59, were recovered from AW sample (AW_40, AW_45) and contained several enzymes contributing in alkane (alkane monooxygenase, aldehyde dehydrogenase) and naphthalene/phenanthrene/phenol/xylene degradation (decarboxylase) under aerobic conditions. The HC degradation potential of representatives of this phylum has been previously reported60.In the Asalouyeh water sample, MAGs AW_25 (UBA4421) and AW_38 (UBA8337) had cddA, AW_21 (UBA8444) had catA, AW_11 (Poseidonia) and AW_17 (from Rhizobiales order) had both alkB/M and xylM, and AW_4 (UBA8337) had catA and pheA genes and had potential to trigger the breakdown of their corresponding oil derivatives.Other recovered genomes had the potential to metabolize the product of initiating enzymes. For instance, AW_23 contained enzymes involved in the degradation of naphthalene, phenol and cyclododecane and was affiliated to the genus Alteromonas (Fig. 6).Three recovered MAGs of HW affiliated to Pseudomonadales (HW_23), Poseidoniales (HW_24), and Flavobacteriales (HW_30) contained some initiating enzymes to degrade cyclododecane/biphenyl/toluene, alkane/xylene, and alkane/xylene/naphthalene/phenanthrene, respectively. A representative of Heimdallarchaeia that are mainly recovered from sediment samples was reconstructed from the Hormuz water sample (HW_28). It had a completeness of 81% and contained enzymes involved in anaerobic degradation of alkanes. This archaeon could potentially be an input from the neighboring OMZ as this phyla include representatives adopted to microoxic niches61. Containing genes with the potential to initiate the oil derivative degradation in the input water with no oil exposure reiterates the intrinsic ability of marine microbiota for HC degradations and oil bioremediation.While 16S rRNA provides an overview of the community, MAGs provide the possibility to inspect the metabolic capability of the microbiota. We decided to provide both in this manuscript as we believe they are complementary. Having the full picture provided by the combination of these analyses allows for a better understanding of the community structure and their metabolic capabilities. This is even more evident for sediment samples as they are highly diverse, and reconstructing MAGs from sediment metagenomes is still a bottlenecks. In this case, we rely more on the 16S rRNA to provide an overall view of the community composition.This said, we see similar taxonomic distribution in the MAGs and 16S rRNA e.g., the prevalence of Flavobacteriales and Rhodobacterales in KhW and KhS, Synechococcales, and Desulfobacteriales and Flavobacteriales in HW, HS and AW samples, respectively.Additionally, some rare microbiota representatives were recovered among reconstructed MAGs. For example, the Immundisolibacterales showed an abundance of only 0.8% in the KhS sample based on 16S rRNA but the recovered KhS_63 MAG was affiliated to this taxon. Notably, this MAG contained many genes involved in hydrocarbon degradation having the highest potential in hydrocarbon degradation. More

Model-guided procedure guides the exploration of butyrate production landscapesWe aimed to explore the butyrate production landscape as a function of community composition to decipher microbial interactions shaping butyrate production. Exploring the butyrate production functional landscape is a major challenge because the number of sub-communities increases exponentially with the number of species43. To investigate the landscape, we developed a modeling framework to guide the iterative design of informative experiments (Fig. 1a, b). Microbial interactions can impact growth or metabolite production by influencing the availability of ecological niches or facilitating metabolite degradation. To capture these two types of interactions, we implemented a two-stage modeling framework to determine the contributions of microbial interactions to species growth and community assembly or metabolite production. In the first stage, a dynamic ecological model, referred to as the generalized Lotka–Volterra model (gLV), predicts community assembly. The second stage predicts metabolite production as a function of the resulting community composition (Fig. 1b). The gLV model is a set of coupled ordinary differential equations that capture the temporal change in species abundances due to monospecies growth parameters and inter-species growth interactions (see the “Methods” section)16. To estimate parameters for the gLV model, we use Bayesian parameter inference techniques to determine the uncertainty in our parameters based on biological and technical variability in the experimental data44.Fig. 1: Iterative modeling framework to predict microbial community assembly and metabolic function.a Two-stage modeling framework for predicting community assembly and function. The generalized Lotka–Volterra model (gLV) represents community dynamics. A Bayesian Inference approach was used to determine parameter uncertainties due to biological and technical variability. A linear regression model with interactions maps assembled community composition to metabolite concentration. Combining these two models enables prediction of a probability distribution of metabolite concentration from initial species abundances. b Design–Test–Learn cycle for model development. First, we use our model to explore the design space of possible experiments (i.e. different initial conditions of species presence/absence) and design communities that span a desired range of metabolite concentrations. Next, we use high-throughput experiments to measure species abundance and metabolite concentration. Finally, we evaluate the model’s predictive capability and infer an updated set of parameters based on the new experimental measurements. c Phylogenetic tree of the synthetic human gut microbiome composed of 25 highly prevalent and diverse species. Branch color indicates phylum and underlined species denote butyrate producers.Full size imageOur metabolite production model consists of a linear regression model with interaction terms mapping community composition (i.e. abundance of each species) at a specific time point to the concentration of an output metabolite at that time. This model was based on a phenomenological model of metabolite production used in bioprocess engineering expanded to microbial communities (see the “Methods” section). In the regression model, the first-order terms capture the monospecies production per unit biomass and the interaction terms represent the impact of inter-species interactions on metabolite production per unit biomass (i.e. deviations from constant metabolite production per unit biomass19). To estimate parameters for the regression model, we use Lasso regression to identify the most impactful interactions. Altogether, the composite gLV and regression model predicts the probability distribution of the metabolite concentration given an initial condition of species abundances (Fig. 1b, see the “Methods” section).In metabolic and protein engineering, a design–test–learn cycle (DTL) has been used to design biomolecules45 or metabolic pathways46 with properties that satisfy desired performance specifications. We hypothesized that this engineering-inspired approach could be used to explore community design spaces and understand the composition–function mapping for butyrate production. Each cycle consisted of: (1) a design phase wherein we used our model informed by experimental observations to simulate a vast number of potential community compositions to identify sub-communities that satisfied biological objectives (i.e. desired butyrate concentrations), (2) a test phase wherein the selected sub-communities were assembled and species abundance and butyrate concentration were measured, and (3) a learn phase wherein patterns in our experimental data were used to estimate model parameters and to extract information about the key microbial interactions influencing community assembly and butyrate production.Two-stage model enables efficient exploration of low richness community design spaceTo develop a system of microbes representing major metabolic functions in the gut, we selected 25 prevalent bacterial species from all major phyla in the human gut microbiome47 (Fig. 1c, Supplementary Data 1). This community contained five butyrate-producing Firmicutes which have been shown to play important roles in human health and protection from diseases (Fig. 1c, Supplementary Data 1). Due to the lack of a defined medium that universally supports the growth of gut microbes, most in vitro studies use undefined media, making it difficult to interrogate the effects of unknown components on community behaviors48. To maximize our knowledge of the substrates available to the communities, we developed a chemically defined medium to grow the synthetic communities (see the “Methods” section).Based on previous studies using pairwise communities to predict higher richness community behaviors16,18,49, we hypothesized that training our model on single and pairwise community measurements would provide an informative starting point for mapping composition–function relationships determining butyrate production. To do so, we first measured time-resolved growth of single species and observed a wide variety of growth dynamics within each phylum, including disparate growth rates and carrying capacities (Supplementary Fig. 1). We assembled each pairwise community containing at least one butyrate producer (the focal species of our system50) and measured species abundance and the concentrations of organic acid fermentation products (including butyrate, lactate, succinate, and acetate) after 48 h. The pairwise consortia displayed a broad range of butyrate concentrations of 0–50 mM (Fig. 2a).Fig. 2: Exploring the predicted butyrate production of 3–5 member communities with a model trained on 1–2 species communities.a Categorical scatter plot of butyrate production in 1–2 species and 24–25 species communities. Solid datapoints indicate the mean of the biological replicates which are represented by transparent datapoints connected to the mean with transparent lines. The colors indicate which butyrate producer was present in the community with green indicating the presence of multiple butyrate producers. DP− and AC− indicate the 24-species communities lacking Desulfovibrio piger (DP) and Anaerostipes caccae (AC), respectively. b Predicted medians (black line) and 60 percent confidence intervals (gray bars) of butyrate concentration for all 3–5 member communities containing at least one butyrate producer (46,591 community predictions). Colored lines indicate median and 60 percent confidence interval of butyrate production of communities chosen for the experimental design with the color indicating the number of species in the community (156 communities). Subplots separate groups of communities based on the identities of the combination of butyrate producers specified. c Scatter plot of measured butyrate versus predicted butyrate for 3–5 species communities. Colors indicate which butyrate producer was present in the community as in a. Biological replicates (n = 1–5, depending on the community, exact values in source data) are indicated by transparent squares connected to the corresponding mean, which is represented by the large data point. Prediction error bars (x-axis) indicate the 60% confidence interval of the predicted butyrate distribution as in b, with the center being the median prediction. Dashed line indicates the linear regression between the mean measured butyrate and the median predicted butyrate. Indicated statistics are for Pearson correlation (two-sided). Source data are available in the Source Data file.Full size imageSingle-species deletion communities have been used to investigate the contributions of individual species to a community function13,16. Therefore, we characterized the full 25-species community and each single-species deletion sub-community (i.e. 24-member consortia). In stark contrast to the pairwise communities, the 24- and 25-species communities exhibited similar low butyrate production (~2–22 mM Butyrate). The absence of only two species Desulfovibrio piger (DP) (~22 mM Butyrate) and Anaerostipes caccae (AC) (~2 mM Butyrate) resulted in a significant increase or decrease in butyrate concentration compared to the remaining 24-member and 25-member communities (Fig. 2a, Supplementary Fig. 2a). In addition, the concentrations of all measured organic acids spanned a much smaller range in the 24 and 25-member communities than the single and pairwise consortia (Supplementary Fig. 2b). These results suggest that high richness communities may trend towards a similar low butyrate-producing state that is difficult to change by the deletion of most single species and motivates a model-guided design strategy for exploring how community richness shapes butyrate production.To determine whether individual and pairwise communities could predict community composition and butyrate production of low richness communities (i.e. 3–5 species), we estimated the parameters of our model based on experimental measurements. Our initial model was informed only by pairwise communities that contained at least one butyrate producer (Supplementary Data 2, M1) and was thus naïve to all interactions between non-butyrate producers. We assumed that the unobserved growth interactions could be predicted based on trends in measured interactions across phylogenetic relatedness (see the “Methods” section)16. However, the resulting model was unable to predict butyrate production in the 24-and 25-member communities (Supplementary Fig. 3), which we attributed to missing information about non-butyrate producer interactions in our training data. Thus, we used our model to explore a low richness design space of 3–5 species communities based on the assumption that pairwise interactions would be more observable in low than high richness (i.e. >10 species) communities to identify an improved parameter estimate for non-butyrate producer interactions.We used our initial M1 model to predict the probability distributions of butyrate production for all 3–5 species communities containing at least one butyrate producer (46,591 communities). The predicted butyrate production varied substantially based on the combination of butyrate producers present in each community (Fig. 2b). In addition, we observed variations in the shapes of the probability distributions based on how the uncertainty in growth prediction propagated through the regression model. For instance, the butyrate concentration in the AC, Roseburia intestinalis (RI) pairwise community was lower than the AC monoculture, even though RI was low abundance, resulting in a high magnitude negative parameter in the regression model for a production interaction between AC and RI (Supplementary Data 3). Due to the uncertainty in the growth parameters, the model predicted that RI would grow substantially in a subset of the 3–5 member simulations containing both AC and RI. The variability in predicted RI growth combined with the high magnitude negative interaction parameter between AC and RI resulted in distributions where the median butyrate concentration was high (i.e. for simulations where RI did not grow substantially), and the 60 percent confidence interval extended to 0 mM butyrate (i.e. when RI grew substantially) (Fig. 2b). In sum, these results demonstrate that the shape of the predicted probability distributions can provide information about the uncertainty in species growth based on experimental observations.Based on the simulations, we selected 156 communities that spanned a broad range of predicted butyrate concentrations across the butyrate producer groups to evaluate experimentally (Fig. 2b). The model prediction exhibited good agreement with the rank order of butyrate production (Spearman rho = 0.84, p = 9*10−43) (Fig. 2c) and species abundance (Spearman rho = 0.76, p = 3*10−122) (Supplementary Fig. 4a–d), demonstrating that our initial model could predict a wide range of butyrate production in low richness communities.Composition–function landscape predicts contributions of growth and production interactionsEncouraged by our model’s predictive ability, we sought to explore composition–function relationships in higher richness communities (i.e. >10 species) using a model with updated parameters based on measurements of the 3–5 member communities (Supplementary Data 2, M2). Since the human gut microbiome exhibits functional redundancy in butyrate pathways51, we first used model M2 to simulate the assembly of all communities containing all five butyrate producers (5-butyrate producer or 5BP, 1,048,576 total) to map the composition–function landscape for butyrate production (Fig. 3a). In addition, we simulated the assembly of all communities containing the four butyrate producers excluding AC (4-butyrate producer or 4BP, 1,048,576 total) to understand how the composition–function landscape changes in the absence of the most productive butyrate producer (Fig. 3b). The majority of 5BP communities were predicted to have higher butyrate concentration than any of the 4BP communities (Fig. 3a, b), consistent with the substantial decrease in butyrate in the AC deletion community observed previously (Fig. 2a).Fig. 3: Community composition–function landscapes reveal key role of production interactions on A. caccae and negative impact of D. piger on butyrate production.a Scatter plot of predicted total butyrate producer abundance versus predicted butyrate concentration for all possible communities in which all five butyrate producers are present (1,048,576 communities). Histograms indicate the butyrate concentration distribution across the given axis. Communities are colored according to the presence (red) or absence (blue) of D. piger (DP). Blue and red dashed lines indicate the linear regression of communities with (red, y = −2.3x + 25.8, r = −0.34) or without (blue, y = 3.1x + 28.0, r = 0.76) DP. The white star indicates the full 25-member community and black star indicates the community of five butyrate producers alone. Large data points indicate communities chosen for experimental validation. Black triangles indicate leave-one-out communities, black circles indicate designed communities, and gray squares indicate random communities, with open/closed symbols indicate the absence/presence of DP. b Scatter plot of predicted total butyrate producer abundance versus predicted butyrate concentration for all possible communities in which all four butyrate producers excluding AC are present (1,048,576 communities). Histograms indicate the butyrate concentration distribution. Gray dashed line indicates the mean predicted butyrate concentration across all communities. Black dashed line indicates the linear regression of all communities (y = 13.2x−0.1, r = 0.64). The white star indicates the full 24-member community and the black star indicates the four butyrate producers alone. Large data points indicate communities chosen for experimental validation. c, d Scatter plots of experimental measurements of total butyrate producer abundance versus butyrate concentration for communities with (c) and without (d) AC. Data point shapes correspond to the legends in (a) and (b) and represent the mean of biological replicates, which are shown as small datapoints connected to the corresponding mean with lines (n = 2 except for 5 BP, 24 and 25 species communities where n = 5–8). Dashed line in (d) indicates the linear regression (y = 8.9x−0.5). e Comparison of butyrate concentration in communities from a and b with and without DP for both the designed and random experimental sets for the 5BP communities and designed set for the 4BP communities. Data point shapes correspond to the legends in a and b and represent the mean of biological replicates, which are shown as small datapoints connected to the corresponding mean with lines. Box and whisker plots represent the median (center line), quartiles (box), and range (whiskers) of the mean butyrate concentration for each community, excluding outliers (points outside 1.5 times the interquartile range). Indicated p-values are from a Mann–Whitney U test (5BP Designed: n = 28 for DP+ and n = 54 for DP−; 5BP Random: n = 27 for DP+ and n = 55 for DP−; 4BP: n = 42 for DP+ and n = 42 for DP−). f Butyrate concentration per unit biomass as a function of sulfide concentration after 24 h of growth. Butyrate concentration per biomass was normalized to the no sulfide condition. Circles indicate the mean of biological replicates, with individual replicates shown as transparent squares (n = 4). Inset: Butyrate concentration per biomass (mM OD600−1) for AC with and without the addition of 1.6 mM sulfide across time (n = 3). Endpoint sulfide concentrations were higher in the data shown in the inset than in the main figure (Supplementary Fig. 6). Source data are available in the Source Data file.Full size imageThe relationship between butyrate producer abundance and butyrate can provide insight into the contributions of growth and production interactions in the presence and absence of AC (Fig. 3a, b). If butyrate producer abundance correlates with butyrate, then growth interactions drive butyrate production, whereas the contributions of production interactions would reduce the strength of this correlation. The 5BP communities were predicted to have a large contribution of production interactions as evidenced by a weak correlation between butyrate concentration and butyrate producer abundance (Spearman rho = 0.17, p  More

1.Bailey, D. L., Humm, J. L., Todd-Pokropek, A. & van Aswegen, A. Nuclear Medicine Physics: A Handbook for Teachers and Students. International Atomic Energy Agency (International Atomic Energy Agency, 2014).2.McGraw-Hill. McGraw-Hill encyclopedia of science & technology. (McGraw-Hill, 2007).3.Medawar, P. B. An unsolved problem of biology. in The uniqueness of the individual (ed. Medawar, P. B.) 44–70 (Basic Books, Inc., 1952).4.Leike, A. Demonstration of the exponential decay law using beer froth. Eur. J. Phys. 23, 21–26 (2002).Article 

Google Scholar 
5.Pauly, D. On the interrelationships between natural mortality, growth parameters, and mean environmental temperature in 175 fish stocks. ICES J. Mar. Sci. 39, 175–192 (1980).Article 

Google Scholar 
6.Vetter, E. F. Estimation of natural mortality in fish stocks: a review. Fish. Bull. 86, 25–43 (1988).
Google Scholar 
7.Gosselin, J., Zedrosser, A., Swenson, J. E. & Pelletier, F. The relative importance of direct and indirect effects of hunting mortality on the population dynamics of brown bears. Proc. R. Soc. B Biol. Sci. 282, 1–9 (2015).8.Nowak, D. J., Kuroda, M. & Crane, D. E. Tree mortality rates and tree population projections in Baltimore, Maryland, USA. Urban . Urban Green. 2, 139–147 (2004).Article 

Google Scholar 
9.Hoenig, J. M. et al. The logic of comparative life history studies for estimating key parameters, with a focus on natural mortality rate. ICES J. Mar. Sci. 73, 2453–2467 (2016).Article 

Google Scholar 
10.Krebs, C. J. Ecology: The Experimental Analysis of Distribution and Abundance. (Pearson Education Limited, 2014).11.Myers, R. A., Bowen, K. G. & Barrowman, N. J. Maximum reproductive rate of fish at low population sizes. Can. J. Fish. Aquat. Sci. 56, 2404–2419 (1999).
Google Scholar 
12.Simpfendorfer, C. A., Bonfil, R. & Latour, R. J. Mortality estimation. in. FAO Fish. Tech. Pap. 474, 127 (2005).
Google Scholar 
13.Cortés, E. Perspectives on the intrinsic rate of population growth. Methods Ecol. Evol. 7, 1136–1145 (2016).Article 

Google Scholar 
14.IUCN. Guidelines for Using the IUCN Red List Categories and Criteria. Version 14. Prepared by the Standards and Petitions Committee. Geographical 14, 1–113 (2019).15.Myers, R. A. & Worm, B. Extinction, survival or recovery of large predatory fishes. Philos. Trans. R. Soc. B Biol. Sci. 360, 13–20 (2005).Article 

Google Scholar 
16.Conde, D. A. et al. Data gaps and opportunities for comparative and conservation biology. Proc. Natl. Acad. Sci. U. S. A. 116, 9658–9664 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
17.Gavrilov, L. & Gavrilova, N. The biology of life span: a quantitative approach. (Harwood Academic Publishers, 1991).18.Sekharan, K. Estimates of the stocks of oil sardine and mackerel in the present fishing grounds off the West coast of India. Indian J. Fish. 21, 177–182 (1974).
Google Scholar 
19.Alagaraja, K. Simple methods for estimation of parameters for assessing exploited fish stocks. Indian J. Fish. 31, 177–208 (1984).
Google Scholar 
20.Cadima, E. L. Fish stock assessment manual. FAO Fish. Tech. Pap. 393, 161 (2003).
Google Scholar 
21.Hewitt, D. A. & Hoenig, J. M. Comparison of two approaches for estimating natural mortality based on longevity. Fish. Bull. 103, 433–437 (2005).
Google Scholar 
22.Dureuil, M. et al. Unified natural mortality estimation for teleosts and elasmobranchs. Mar. Ecol. Prog. Ser. https://doi.org/10.3354/meps13704 (accepted).23.Litzgus, J. D. Sex differences in longevity in the spotted turtle (Clemmys guttata). Copeia 2, 281–288 (2006).Article 

Google Scholar 
24.Calder, W. A. III Body size, mortality, and longevity. J. Theor. Biol. 102, 135–144 (1983).PubMed 
Article 
PubMed Central 

Google Scholar 
25.Botkin, D. B., Janak, J. F. & Wallis, J. R. Some ecological consequences of a computer model of forest growth. J. Ecol. 60, 849–872 (1972).Article 

Google Scholar 
26.Holt, S. J. A note on the relation between the mortality rate and the duration of life in an exploited fish population. Int. Comm. Northwest Atl. Fish. Res. Bull. 2, 73–75 (1965).
Google Scholar 
27.Hoenig, J. M. Should natural mortality estimators based on maximum age also consider sample size? Trans. Am. Fish. Soc. 146, 136–146 (2017).Article 

Google Scholar 
28.Williams, G. C. Pleiotropy, natural selection, and the evolution of senescence. Evolution 11, 398–411 (1957).Article 

Google Scholar 
29.Hamilton, W. D. The moulding of senescence by natural selection. J. Theor. Biol. 12, 12–45 (1966).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
30.Kirkwood, T. B. L. Evolution of ageing. Nature 270, 301–304 (1977).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
31.Kirkwood, T. B. L. & Rose, M. R. Evolution of senescence: late survival sacrificed for reproduction. Philos. Trans. R. Soc. Lond., B 332, 15–24 (1991).CAS 
Article 

Google Scholar 
32.Froese, R. & Pauly, D. FishBase. World Wide Web Electronic Publication (2019). Available at: www.fishbase.org. (accessed: 6th February 2018)33.I. C. E. S. Herring (Clupea harengus) in Subarea 4 and divisions 3.a and 7.d, autumn spawners (North Sea, Skagerrak and Kattegat, eastern English Channel). in Report of the ICES Advisory Committee, 2019. ICES Advice 2019, her.27.3a47d 11 (2019).34.Caswell, H. & Shyu, E. Senescence, selection gradients and mortality. in The Evolution of Senescence in the Tree of Life (eds. Shefferson, R. P., Jones, O. R. & Salguero-Gómez, R.) 56–82 (Cambridge University Press, 2017). https://doi.org/10.1017/9781139939867.00435.Promislow, D. E. L. Senescence in natural populations of mammals: a comparative study. Evolution 45, 1869–1887 (1991).PubMed 
Article 
PubMed Central 

Google Scholar 
36.Sibly, R. M., Collett, D., Promislow, D. E. L., Peacock, D. J. & Harvey, P. H. Mortality rates of mammals. J. Zool. 243, 1–12 (1997).Article 

Google Scholar 
37.Blumstein, D. T. & Møller, A. P. Is sociality associated with high longevity in North American birds? Biol. Lett. 4, 146–148 (2008).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
38.Nussey, D. H., Froy, H., Lemaitre, J.-F., Gaillard, J.-M. & Austad, S. N. Senescence in natural populations of animals: Widespread evidence and its implications for bio-gerontology. Ageing Res. Rev. 12, 214–225 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
39.Salguero-Gómez, R. & Jones, O. R.. Life history trade-offs modulate the speed of senescence. in The Evolution of Senescence in the Tree of Life (eds. Shefferson, R. P., Jones, O. R. & Salguero-Gómez, R.) 403–421 (Cambridge University Press, 2017). https://doi.org/10.1017/9781139939867.02040.Hoekstra, L. A., Schwartz, T. S., Sparkman, A. M., Miller, D. A. W. & Bronikowski, A. M. The untapped potential of reptile biodiversity for understanding how and why animals age. Funct. Ecol. 34, 38–54 (2020).PubMed 
Article 
PubMed Central 

Google Scholar 
41.Bonduriansky, R. & Brassil, C. E. Rapid and costly ageing in wild male flies. Nature 420, 377 (2002).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
42.Zajitschek, F., Zajitschek, S. & Bonduriansky, R. Senescence in wild insects: Key questions and challenges. Funct. Ecol. 34, 26–37 (2020).Article 

Google Scholar 
43.Roach, D. A. & Smith, E. F. Life-history trade-offs and senescence in plants. Funct. Ecol. 34, 17–25 (2020).Article 

Google Scholar 
44.Jones, O. R. et al. Diversity of ageing across the tree of life. Nature 505, 169–173 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
45.Ruby, J. G., Smith, M. & Buffenstein, R. Naked mole-rat mortality rates defy Gompertzian laws by not increasing with age. Elife 7, 1–18 (2018).Article 

Google Scholar 
46.Keller, L. & Genoud, M. Extraordinary lifespans in ants: a test of evolutionary theories of ageing. Nature 389, 958–960 (1997).CAS 
Article 

Google Scholar 
47.Cooke, G. M., Tonkins, B. M. & Mather, J. A. Care and Enrichment for Captive Cephalopods. in The Welfare of Invertebrate Animals (eds. Carere, C. & Mather, J.). 179–208 (Springer, Cham, 2019). https://doi.org/10.1007/978-3-030-13947-6_848.Baudisch, A. et al. The pace and shape of senescence in angiosperms. J. Ecol. 101, 596–606 (2013).Article 

Google Scholar 
49.Halley, J. M., Van Houtan, K. S. & Mantua, N. How survival curves affect populations’ vulnerability to climate change. PLoS One 13, 1–18 (2018).Article 
CAS 

Google Scholar 
50.Gompertz, B. On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. Philos. Trans. R. Soc. Lond. 115, 513–583 (1825).
Google Scholar 
51.Makeham, W. M. On the law of mortality and the construction of annuity tables. Assur. Mag. J. Inst. Actuar. 8, 301–310 (1860).Article 

Google Scholar 
52.Finch, C. E. & Pike, M. C. Maximum life span predictions from the Gompertz mortality model. J. Gerontol. Biol. Sci. 51A, 183–194 (1996).Article 

Google Scholar 
53.Reznick, D. N., Bryant, M. J., Roff, D., Ghalambor, C. K. & Ghalambor, D. E. Effect of extrinsic mortality on the evolution of senescence in guppies. Nature 431, 1095–1099 (2004).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
54.Kirkwood, T. B. L. Deciphering death: a commentary on Gompertz (1825) ‘On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies’. Philos. Trans. R. Soc. B Biol. Sci. 370, 1–8 (2015).Article 

Google Scholar 
55.Gavrilov, L. A. & Gavrilova, N. S. New trend in old-age mortality: gompertzialization of mortality trajectory. Gerontology 65, 451–457 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
56.Ohsumi, S. Interspecies relationships among some biological parameters in cetaceans and estimation of the natural mortality coefficient of the Southern Hemisphere minke whale. Rep. Int. Whal. Comm. 29, 397–406 (1979).
Google Scholar 
57.Mizroch, S. A. On the relationship between mortality rate and length in baleen whales. Rep. Int. Whal. Comm. 35, 505–510 (1985).
Google Scholar  More