Ecology

1.
van Belkum A, Soriaga LB, LaFave MC, Akella S, Veyrieras JB, Barbu EM, et al. Phylogenetic distribution of CRISPR-Cas systems in antibiotic-resistant Pseudomonas aeruginosa. MBio. 2015;6:e01796–15.
2.
Makarova KS, Wolf YI, Iranzo J, Shmakov SA, Alkhnbashi OS, Brouns SJ, et al. Evolutionary classification of CRISPR-Cas systems: a burst of class 2 and derived variants. Nat Rev Microbiol. 2019;18:1–17.

3.
Portillo MC, Gonzalez JM. CRISPR elements in the Thermococcales: evidence for associated horizontal gene transfer in Pyrococcus furiosus. J Appl Genet. 2009;50:421–30.
CAS  PubMed  Article  Google Scholar 

4.
Watson BN, Staals RH, Fineran, PC. CRISPR-Cas-mediated phage resistance enhances horizontal gene transfer by transduction. MBio. 2018;9:e02406–17.

5.
Varble A, Meaden S, Barrangou R, Westra ER, Marraffini LA. Recombination between phages and CRISPR-cas loci facilitates horizontal gene transfer in staphylococci. Nat Microbiol. 2019;4:956–63.
CAS  PubMed  PubMed Central  Article  Google Scholar 

6.
Levin BR. Nasty viruses, costly plasmids, population dynamics, and the conditions for establishing and maintaining CRISPR-mediated adaptive immunity in bacteria. PLoS Genet. 2010;6:e1001171.
PubMed  PubMed Central  Article  CAS  Google Scholar 

7.
van Houte S, Buckling A, Westra ER. Evolutionary ecology of prokaryotic immune mechanisms. Microbiol Mol Biol Rev. 2016;80:745–63.
PubMed  PubMed Central  Article  Google Scholar 

8.
Rollie C, Chevallereau A, Watson BN, Chyou TY, Fradet O, McLeod I, et al. Targeting of temperate phages drives loss of type I CRISPR-Cas systems. Nature. 2020;578:149–53.
CAS  PubMed  PubMed Central  Article  Google Scholar 

9.
Westra ER, van Houte S, Oyesiku-Blakemore S, Makin B, Broniewski JM, Best A, et al. Parasite exposure drives selective evolution of constitutive versus inducible defense. Curr Biol. 2015;25:1043–9.
CAS  PubMed  Article  Google Scholar 

10.
Vale PF, Lafforgue G, Gatchitch F, Gardan R, Moineau S, Gandon S. Costs of CRISPR-Cas-mediated resistance in Streptococcus thermophilus. Proc R Soc B: Biol Sci. 2015;282:20151270.
Article  CAS  Google Scholar 

11.
van Houte S, Ekroth AK, Broniewski JM, Chabas H, Ashby B, Bondy-Denomy J, et al. The diversity-generating benefits of a prokaryotic adaptive immune system. Nature. 2016;532:385–8.
PubMed  PubMed Central  Article  CAS  Google Scholar 

12.
Stern A, Keren L, Wurtzel O, Amitai G, Sorek R. Self-targeting by CRISPR: gene regulation or autoimmunity? Trends Genet. 2010;26:335–40.
CAS  PubMed  PubMed Central  Article  Google Scholar 

13.
Levy A, Goren MG, Yosef I, Auster O, Manor M, Amitai G, et al. CRISPR adaptation biases explain preference for acquisition of foreign DNA. Nature. 2015;520:505–10.
CAS  PubMed  PubMed Central  Article  Google Scholar 

14.
Vercoe RB, Chang JT, Dy RL, Taylor C, Gristwood T, Clulow JS, et al. Cytotoxic chromosomal targeting by CRISPR/Cas systems can reshape bacterial genomes and expel or remodel pathogenicity islands. PLoS Genet. 2013;9:e1003454.
CAS  PubMed  PubMed Central  Article  Google Scholar 

15.
Staals RH, Jackson SA, Biswas A, Brouns SJ, Brown CM, Fineran PC. Interference-driven spacer acquisition is dominant over naive and primed adaptation in a native CRISPR–Cas system. Nat Commun. 2016;7:1–13.
Article  CAS  Google Scholar 

16.
Weissman JL, Stoltzfus A, Westra ER, Johnson PL. Avoidance of Self during CRISPR immunization. Trends Microbiol. 2020;28:543–53.

17.
Koonin EV. Open questions: CRISPR biology. BMC Biol. 2018;16:1–3.
Article  CAS  Google Scholar 

18.
Patterson AG, Yevstigneyeva MS, Fineran PC. Regulation of CRISPR–Cas adaptive immune systems. Curr Opin Microbiol. 2017;37:1–7.
CAS  PubMed  Article  Google Scholar 

19.
Quax TE, Voet M, Sismeiro O, Dillies MA, Jagla B, Coppee JY, et al. Massive activation of archaeal defense genes during viral infection. J Virol. 2013;87:8419–28.
CAS  PubMed  PubMed Central  Article  Google Scholar 

20.
Young JC, Dill BD, Pan C, Hettich RL, Banfield JF, Shah M, et al. Phage-induced expression of CRISPR-associated proteins is revealed by shotgun proteomics in Streptococcus thermophilus. PloS ONE. 2012;7:e38077.
CAS  PubMed  PubMed Central  Article  Google Scholar 

21.
Agari Y, Sakamoto K, Tamakoshi M, Oshima T, Kuramitsu S, Shinkai A. Transcription profile of Thermus thermophilus CRISPR systems after phage infection. J Mol Biol. 2010;395:270–81.
CAS  PubMed  Article  Google Scholar 

22.
Chabas H, van Houte S, Høyland-Kroghsbo NM, Buckling A, Westra ER. Immigration of susceptible hosts triggers the evolution of alternative parasite defence strategies. Proc R Soc B: Biol Sci. 2016;283:20160721.
Article  Google Scholar 

23.
Pawluk A, Staals RH, Taylor C, Watson BN, Saha S, Fineran PC, et al. Inactivation of CRISPR-Cas systems by anti-CRISPR proteins in diverse bacterial species. Nat Microbiol. 2016;1:1–6.
Article  CAS  Google Scholar 

24.
Bondy-Denomy J, Garcia B, Strum S, Du M, Rollins MF, Hidalgo-Reyes Y, et al. Multiple mechanisms for CRISPR-Cas inhibition by anti-CRISPR proteins. Nature. 2015;526:136–9.
CAS  PubMed  PubMed Central  Article  Google Scholar 

25.
Stanley SY, Borges AL, Chen KH, Swaney DL, Krogan NJ, Bondy-Denomy J, et al. Anti-CRISPR-associated proteins are crucial repressors of anti-CRISPR transcription. Cell. 2019;178:1452–64.
CAS  PubMed  PubMed Central  Article  Google Scholar 

26.
Cady KC, Bondy-Denomy J, Heussler GE, Davidson AR, O’Toole GA. The CRISPR/Cas adaptive immune system of Pseudomonas aeruginosa mediates resistance to naturally occurring and engineered phages. J Bacteriol. 2012;194:5728–38.
CAS  PubMed  PubMed Central  Article  Google Scholar 

27.
Chevallereau A, Meaden S, Fradet O, Landsberger M, Maestri A, Biswas A, et al. Exploitation of the cooperative behaviors of anti-CRISPR phages. Cell Host Microbe. 2020;27:189–98.
CAS  PubMed  PubMed Central  Article  Google Scholar 

28.
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011;17:10–12.
Article  Google Scholar 

29.
Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol. 2012;19:455–77.
CAS  PubMed  PubMed Central  Article  Google Scholar 

30.
Zhang J, Kobert K, Flouri T, Stamatakis A. PEAR: a fast and accurate Illumina Paired-End reAd mergeR. Bioinformatics. 2014;30:614–20.
CAS  PubMed  Article  Google Scholar 

31.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.
CAS  PubMed  Article  Google Scholar 

32.
Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7:335–6.
CAS  PubMed  PubMed Central  Article  Google Scholar 

33.
Mahé F, Rognes T, Quince C, de Vargas C, Dunthorn M. Swarm v2: highly-scalable and high-resolution amplicon clustering. PeerJ. 2015;3:e1420.
PubMed  PubMed Central  Article  Google Scholar 

34.
Chao A. Nonparametric estimation of the number of classes in a population. Scand J Stat. 1984;11:265–70.

35.
Biswas A, Staals RH, Morales SE, Fineran PC, Brown CM. CRISPRDetect: a flexible algorithm to define CRISPR arrays. BMC Genom. 2016;17:1–14.
Article  CAS  Google Scholar 

36.
Biswas A, Gagnon JN, Brouns SJ, Fineran PC, Brown CM. CRISPRTarget: bioinformatic prediction and analysis of crRNA targets. RNA Biol. 2013;10:817–27.
CAS  PubMed  PubMed Central  Article  Google Scholar 

37.
Joshi N. & Fass J. Sickle: A sliding-window, adaptive, quality-based trimming tool for fastQ files (version 1.33) [software]. 2011. https://github.com/najoshi/sickle.

38.
Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009;25:1754–60.
CAS  PubMed  PubMed Central  Article  Google Scholar 

39.
Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31:166–9.
CAS  PubMed  PubMed Central  Article  Google Scholar 

40.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
PubMed  PubMed Central  Article  CAS  Google Scholar 

41.
Galili T. dendextend: an R package for visualizing, adjusting and comparing trees of hierarchical clustering. Bioinformatics. 2015;31:3718–20.
CAS  PubMed  PubMed Central  Article  Google Scholar 

42.
Richter C, Dy RL, McKenzie RE, Watson BN, Taylor C, Chang JT, et al. Priming in the Type IF CRISPR-Cas system triggers strand-independent spacer acquisition, bi-directionally from the primed protospacer. Nucleic Acids Res. 2014;42:8516–26.
CAS  PubMed  PubMed Central  Article  Google Scholar 

43.
Fineran PC, Gerritzen MJ, Suárez-Diez M, Künne T, Boekhorst J, van Hijum SA, et al. Degenerate target sites mediate rapid primed CRISPR adaptation. Proc Natl Acad Sci USA. 2014;111:E1629–38.
CAS  PubMed  Article  Google Scholar 

44.
Xue C, Seetharam AS, Musharova O, Severinov K, Brouns J, Severin SJ, et al. CRISPR interference and priming varies with individual spacer sequences. Nucleic Acids Res. 2015;43:10831–47.
CAS  PubMed  PubMed Central  Article  Google Scholar 

45.
Wortham BW, Oliveira MA, Patel CN. Polyamines in bacteria: pleiotropic effects yet specific mechanisms. In: (eds Perry RD, Fetherston JD). The genus Yersinia. New York, NY: Springer; 2007. p. 106–15.

46.
Bru JL, Rawson B, Trinh C, Whiteson K, Høyland-Kroghsbo NM, Siryaporn A. PQS produced by the Pseudomonas aeruginosa stress response repels swarms away from bacteriophage and antibiotics. J Bacteriol. 2019;201:e00383–19.

47.
Doron S, Fedida A, Hernández-Prieto MA, Sabehi G, Karunker I, Stazic D, et al. Transcriptome dynamics of a broad host-range cyanophage and its hosts. ISME J. 2016;10:1437–55.
CAS  PubMed  Article  Google Scholar 

48.
Landsberger M, Gandon S, Meaden S, Rollie C, Chevallereau A, Chabas H, et al. Anti-CRISPR phages cooperate to overcome CRISPR-Cas immunity. Cell. 2018;174:908–16.
CAS  PubMed  PubMed Central  Article  Google Scholar 

49.
Wei Y, Terns RM, Terns MP. Cas9 function and host genome sampling in Type II-A CRISPR-Cas adaptation. Genes Dev. 2015;29:356–61.
CAS  PubMed  PubMed Central  Article  Google Scholar 

50.
Fokine A, Rossmann MG. Common evolutionary origin of procapsid proteases, phage tail tubes, and tubes of bacterial type VI secretion systems. Structure. 2016;24:1928–35.
CAS  PubMed  PubMed Central  Article  Google Scholar 

51.
Modell JW, Jiang W, Marraffini LA. CRISPR-Cas systems exploit viral DNA injection to establish and maintain adaptive immunity. Nature. 2017;544:101–4.
CAS  PubMed  PubMed Central  Article  Google Scholar 

52.
van Sluijs L, van Houte S, van Der Oost J, Brouns SJ, Buckling A, Westra ER. Addiction systems antagonize bacterial adaptive immunity. FEMS Microbiol Lett. 2019;366:fnz047.
CAS  PubMed  PubMed Central  Google Scholar 

53.
Wagemans J, Blasdel BG, Van den Bossche A, Uytterhoeven B, De Smet J, Paeshuyse J, et al. Functional elucidation of antibacterial phage ORFans targeting Pseudomonas aeruginosa. Cell Microbiol. 2014;16:1822–35.
CAS  PubMed  Article  Google Scholar 

54.
van Houte S, Ros VI, van Oers MM. Walking with insects: molecular mechanisms behind parasitic manipulation of host behaviour. Mol Ecol. 2013;22:3458–75.
PubMed  Article  Google Scholar 

55.
Taylor TB, Buckling A. Bacterial motility confers fitness advantage in the presence of phages. J Evolut Biol. 2013;26:2154–60.
CAS  Article  Google Scholar 

56.
Heussler GE, Cady KC, Koeppen K, Bhuju S, Stanton BA, O’Toole GA. Clustered regularly interspaced short palindromic repeat-dependent, biofilm-specific death of Pseudomonas aeruginosa mediated by increased expression of phage-related genes. MBio. 2015;6:e00129–15.

57.
Hernandez CA, Koskella B. Phage resistance evolution in vitro is not reflective of in vivo outcome in a plant‐bacteria‐phage system. Evolution. 2019;73:2461–75.
CAS  PubMed  Article  Google Scholar 

58.
Alseth EO, Pursey E, Luján AM, McLeod I, Rollie C, Westra ER. Bacterial biodiversity drives the evolution of CRISPR-based phage resistance. Nature. 2019;574:549–52.
CAS  PubMed  PubMed Central  Article  Google Scholar 

59.
Westra E, Levin BR. How important is CRISPR-Cas for protecting natural populations of bacteria against infections with badass DNAs? BioRxiv. 2020. https://www.biorxiv.org/content/10.1101/2020.02.05.935965v2. More

1.
Bergquist PR. Sponges. London, United Kingdom: Hutchinson and Co. Ltd; 1978.
2.
Brain CKB, Prave AR, Hoffmann KH, Fallick AE, Botha A, Herd DA, et al. The first animals: Ca. 760-million-year-old sponge-like fossils from Namibia. S Afr J Sci. 2012;108:1–8.
Article  Google Scholar 

3.
Hooper JNA, van Soest RWM. Systema Porifera: a guide to the classification of sponges. New York: Kluwer Academic/Plenum Publishers; 2002.
Google Scholar 

4.
Webster NS. Sponge disease: a global threat? Environ Microbiol. 2007;9:1363–75.
CAS  PubMed  Article  PubMed Central  Google Scholar 

5.
Bell JJ. The functional roles of marine sponges. Estuar Coast Shelf Sci. 2008;79:341–53.
Article  Google Scholar 

6.
de Goeij JM, van Oevelen D, Vermeij MJA, Middelburg JJ, Osinga R, de Goeij AFPM, et al. Surviving in a Marine Desert: the sponge loop retains resources within coral reefs. Science. 2013;342:108–10.
PubMed  Article  CAS  PubMed Central  Google Scholar 

7.
Taylor MW, Radax R, Steger D, Wagner M. Sponge-associated microorganisms: evolution, ecology, and biotechnological potential. Microbiol Mol Biol Rev. 2007;71:295–347.
CAS  PubMed  PubMed Central  Article  Google Scholar 

8.
Wehrl M, Steinert M, Hentschel U. Bacterial uptake by the marine sponge Aplysina aerophoba. Microb Ecol. 2007;53:355–65.
PubMed  Article  PubMed Central  Google Scholar 

9.
Thomas T, Rusch D, DeMaere MZ, Yung PY, Lewis M, Halpern A, et al. Functional genomic signatures of sponge bacteria reveal unique and shared features of symbiosis. ISME J. 2010;4:1557–67.
CAS  PubMed  Article  PubMed Central  Google Scholar 

10.
Nguyen MTHD, Liu M, Thomas T. Ankyrin-repeat proteins from sponge symbionts modulate amoebal phagocytosis. Mol Ecol. 2013;23:1635–45.
PubMed  Article  CAS  Google Scholar 

11.
Reynolds D, Thomas T. Evolution and function of eukaryotic-like proteins from sponge symbionts. Mol Ecol. 2016;25:5242–53.
CAS  PubMed  Article  PubMed Central  Google Scholar 

12.
Jahn MT, Arkhipova K, Markert SM, Stigloher C, Lachnit T, Pita L, et al. A phage protein aids bacterial symbionts in eukaryote immune evasion. Cell Host Microbe. 2019;26:542–50.e5.
CAS  PubMed  PubMed Central  Article  Google Scholar 

13.
Degnan SM. The surprisingly complex immune gene repertoire of a simple sponge, exemplified by the NLR genes: a capacity for specificity? Dev Comp Immunol. 2015;48:269–74.
CAS  PubMed  Article  PubMed Central  Google Scholar 

14.
Hentschel U, Usher KM, Taylor MW. Marine sponges as microbial fermenters. FEMS Microbiol Ecol. 2006;55:167–77.
CAS  PubMed  Article  PubMed Central  Google Scholar 

15.
Moitinho-Silva L, Steinert G, Nielsen S, Hardoim CCP, Wu YC, McCormack GP, et al. Predicting the HMA-LMA status in marine sponges by machine learning. Front Microbiol. 2017;8:1–14.
Article  Google Scholar 

16.
Moitinho-Silva L, Nielsen S, Amir A, Gonzalez A, Ackermann GL, Cerrano C, et al. The sponge microbiome project. Gigascience. 2017;6:1–7.
CAS  PubMed  Article  PubMed Central  Google Scholar 

17.
Hentschel U, Hopke J, Horn M, Friedrich AB, Wagner M, Hacker J, et al. Molecular evidence for a uniform microbial community in sponges from different oceans. Appl Environ Microbiol. 2002;68:4431–40.
CAS  PubMed  PubMed Central  Article  Google Scholar 

18.
Simister RL, Deines P, Botté ES, Webster NS, Taylor MW. Sponge-specific clusters revisited: a comprehensive phylogeny of sponge-associated microorganisms. Environ Microbiol. 2012;14:517–24.
CAS  PubMed  Article  PubMed Central  Google Scholar 

19.
Taylor MW, Tsai P, Simister RL, Deines P, Botte E, Ericson G, et al. “Sponge-specific” bacteria are widespread (but rare) in diverse marine environments. ISME J. 2013;7:438–43.
CAS  PubMed  Article  PubMed Central  Google Scholar 

20.
Croué J, West NJ, Escande M-L, Intertaglia L, Lebaron P, Suzuki MT. A single betaproteobacterium dominates the microbial community of the crambescidine-containing sponge Crambe crambe. Sci Rep. 2013;3:2583.
PubMed  PubMed Central  Article  Google Scholar 

21.
Fan L, Reynolds D, Liu M, Stark M, Kjelleberg S, Webster NS, et al. Functional equivalence and evolutionary convergence in complex communities of microbial sponge symbionts. PNAS. 2012;109:E1878–87.
CAS  PubMed  Article  Google Scholar 

22.
Simister RL, Taylor MW, Rogers KM, Schupp PJ, Deines P. Temporal molecular and isotopic analysis of active bacterial communities in two New Zealand sponges. FEMS Microbiol Ecol. 2013;85:195–205.
PubMed  Article  Google Scholar 

23.
Gantt SE, López-Legentil S, Erwin PM. Stable microbial communities in the sponge Crambe crambe from inside and outside a polluted Mediterranean harbor. FEMS Microbiol Lett. 2017;364:1–7.
Article  CAS  Google Scholar 

24.
Bin JeongJ, Kim KH, Park JS. Sponge-specific unknown bacterial groups detected in marine sponges collected from Korea through barcoded pyrosequencing. J Microbiol Biotechnol. 2015;25:1–10.
Article  Google Scholar 

25.
Thiel V, Neulinger SC, Staufenberger T, Schmaljohann R, Imhoff JF. Spatial distribution of sponge-associated bacteria in the Mediterranean sponge Tethya aurantium. FEMS Microbiol Ecol. 2006;59:47–63.
PubMed  Article  CAS  Google Scholar 

26.
Fieth RA, Gauthier M-EA, Bayes J, Green KM, Degnan SM. Ontogenetic changes in the bacterial symbiont community of the tropical demosponge Amphimedon queenslandica: metamorphosis is a new beginning. Front Mar Sci. 2016;3:1–20.
Article  Google Scholar 

27.
Batani G. Fluorescence in situ hybridisation for the localisation and culturing of marine bacteria: co-localisation of symbionts in sponges (Unpublished PhD thesis chapter). UNSW Sydney. 2018.

28.
Waterworth SC, Jiwaji M, Kalinski JCJ, Parker-Nance S, Dorrington RA. A place to call home: an analysis of the bacterial communities in two Tethya rubra Samaai and Gibbons 2005 populations in algoa bay, South Africa. Mar Drugs. 2017;15:95.
PubMed Central  Article  CAS  Google Scholar 

29.
Wu S, Ou H, Liu T, Wang D, Zhao J. Structure and dynamics of microbiomes associated with the marine sponge Tedania sp. during its life cycle. FEMS Microbiol Ecol. 2018;94:1–9.
Article  CAS  Google Scholar 

30.
Matcher GF, Waterworth SC, Walmsley TA, Matsatsa T, Parker-Nance S, Davies-Coleman MT, et al. Keeping it in the family: coevolution of latrunculid sponges and their dominant bacterial symbionts. Microbiologyopen. 2017;6:1–13.
Article  CAS  Google Scholar 

31.
Jackson SA, Flemer B, McCann A, Kennedy J, Morrissey JP, O’Gara F, et al. Archaea appear to dominate the microbiome of Inflatella pellicula deep sea sponges. PLoS ONE. 2013;8:1–8.
Google Scholar 

32.
Steinert G, Taylor MW, Deines P, Simister RL, de Voogd NJ, Hoggard M, et al. In four shallow and mesophotic tropical reef sponges from Guam the microbial community largely depends on host identity. PeerJ. 2016;4:e1936.
PubMed  PubMed Central  Article  CAS  Google Scholar 

33.
Gauthier M-EA, Watson JR, Degnan SM. Draft genomes shed light on the dual bacterial symbiosis that dominates the microbiome of the coral reef sponge Amphimedon queenslandica. Front Mar Sci. 2016;3:1–18.
Article  Google Scholar 

34.
Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Cold Spring Harb Lab Press. 2014;25:1043–55.
Google Scholar 

35.
Mukherjee S, Stamatis D, Bertsch J, Ovchinnikova G, Verezemska O, Isbandi M, et al. Genomes OnLine database (GOLD) v.6: data updates and feature enhancements. Nucleic Acids Res. 2017;45:D446–56.
CAS  PubMed  Article  Google Scholar 

36.
Huntemann M, Ivanova NN, Mavromatis K, James Tripp H, Paez-Espino D, Palaniappan K, et al. The standard operating procedure of the DOE-JGI microbial genome annotation pipeline (MGAP v.4). Stand Genom Sci. 2015;10:4–9.
Article  CAS  Google Scholar 

37.
Eddy S. Profile hidden Markov models. Bioinformatics. 1998;14:755–63.
CAS  PubMed  Article  Google Scholar 

38.
Pruesse E, Peplies J, Glöckner FO. SINA: accurate high-throughput multiple sequence alignment of ribosomal RNA genes. Bioinformatics. 2012;28:1823–9.
CAS  PubMed  PubMed Central  Article  Google Scholar 

39.
Ludwig W, Strunk O, Westram R, Richter L, Meier H, Yadhukumar, et al. ARB: a software environment for sequence data. Nucleic Acids Res. 2004;32:1363–71.
CAS  PubMed  PubMed Central  Article  Google Scholar 

40.
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 2013;41:590–6.
Article  CAS  Google Scholar 

41.
Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30:1312–3.
CAS  PubMed  PubMed Central  Article  Google Scholar 

42.
Letunic I, Bork P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic Acids Res. 2016;44:W242–5.
CAS  PubMed  PubMed Central  Article  Google Scholar 

43.
Chaumeil P-A, Mussig AJ, Hugenholtz P, Parks DH. GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database. Bioinformatics. 2020;36:1925–7.
Google Scholar 

44.
Parks DH, Chuvochina M, Waite DW, Rinke C, Skarshewski A, Chaumeil PA, et al. A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life. Nat Biotechnol. 2018;36:996–1004.
CAS  PubMed  Article  PubMed Central  Google Scholar 

45.
Konstantinidis KT, Tiedje JM. Towards a genome-based taxonomy for prokaryotes. J Bacteriol. 2005;187:6258–64.
CAS  PubMed  PubMed Central  Article  Google Scholar 

46.
Rodriguez-R LM, Konstantinidis KT. The enveomics collection: a toolbox for specialized analyses of microbial genomes and metagenomes. Peer J Prepr. 2016;4:e1900v1.
Google Scholar 

47.
Thomas T, Moitinho-Silva L, Lurgi M, Björk JR, Easson C, Astudillo-García C, et al. Diversity, structure and convergent evolution of the global sponge microbiome. Nat Commun. 2016;7:11870.
CAS  PubMed  PubMed Central  Article  Google Scholar 

48.
Díez-Vives C, Moitinho-Silva L, Nielsen S, Reynolds D, Thomas T. Expression of eukaryotic-like protein in the microbiome of sponges. Mol Ecol. 2017;26:1432–51.
Article  CAS  Google Scholar 

49.
Öztürk B, De Jaeger L, Smidt H, Sipkema D. Culture-dependent and independent approaches for identifying novel halogenases encoded by Crambe crambe (marine sponge) microbiota. Sci Rep. 2013;3:1–9.
Article  Google Scholar 

50.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:1–16.
Article  Google Scholar 

51.
Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.
PubMed  PubMed Central  Article  CAS  Google Scholar 

52.
Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN. RNA-Seq gene expression estimation with read mapping uncertainty. Bioinformatics. 2010;26:493–500.
PubMed  Article  CAS  PubMed Central  Google Scholar 

53.
Wagner GP, Kin K, Lynch VJ. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci. 2012;131:281–5.
CAS  PubMed  Article  PubMed Central  Google Scholar 

54.
Webster NS, Hill RT. The culturable microbial community of the great barrier reef sponge Rhopaloeides odorabile is dominated by an α-Proteobacterium. Mar Biol. 2001;138:843–51.
CAS  Article  Google Scholar 

55.
Engelberts JP, Robbins SJ, de Goeij JM, Aranda M, Bell SC, Webster NS. Characterization of a sponge microbiome using an integrative genome-centric approach. ISME J. 2020;14:1100–10.
CAS  PubMed  Article  PubMed Central  Google Scholar 

56.
Moitinho-Silva L, Díez-Vives C, Batani G, Esteves AIS, Jahn MT, Thomas T. Integrated metabolism in sponge-microbe symbiosis revealed by genome-centered metatranscriptomics. ISME J. 2017;11:1651–66.
CAS  PubMed  PubMed Central  Article  Google Scholar 

57.
Lobb B, Tremblay BJM, Moreno-Hagelsieb G, Doxey AC. An assessment of genome annotation coverage across the bacterial tree of life. Microb Genom. 2020;6:e000341.
PubMed Central  Google Scholar 

58.
Lopez JV, Ranzer LK, Ledger A, Schoch B, Duckworth A, Mccarthy PJ, et al. Comparison of bacterial diversity within the coral reef sponge, axinella corrugata, the encrusting coral erythropodium caribaeorum. Proc. 11 Int. Coral Reef Symp. 2008;2:1362–6.
Google Scholar 

59.
Yarza P, Yilmaz P, Pruesse E, Glöckner FO, Ludwig W, Schleifer KH, et al. Uniting the classification of cultured and uncultured bacteria and archaea using 16S rRNA gene sequences. Nat Rev Microbiol. 2014;12:635–45.
CAS  PubMed  Article  PubMed Central  Google Scholar 

60.
Konstantinidis KT, Rosselló-Móra R, Amann R. Uncultivated microbes in need of their own taxonomy. ISME J. 2017;11:2399–406.
PubMed  PubMed Central  Article  Google Scholar 

61.
Bowers RM, Kyrpides NC, Stepanauskas R, Harmon-smith M, Doud D, Jarett J, et al. Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea. Nat Biotechnol. 2017;35:725–31.
CAS  PubMed  PubMed Central  Article  Google Scholar 

62.
Chuvochina M, Rinke C, Parks DH, Rappé MS, Tyson GW, Yilmaz P, et al. The importance of designating type material for uncultured taxa. Syst Appl Microbiol. 2019;42:15–21.
PubMed  Article  PubMed Central  Google Scholar 

63.
Tripp E. Crowell’s handbook of classical mythology. New York: Thomas Y. Crowell Company; 1970.

64.
Athanassakis AN. Hesiod: Theogony, Works and Days, Shield. 2nd ed. Baltimore and London: The John Hopkins University Press; 2005.

65.
Webster NS, Negri AP, Munro MMHG, Battershill CN. Diverse microbial communities inhabit Antarctic sponges. Environ Microbiol. 2004;6:288–300.
PubMed  Article  PubMed Central  Google Scholar 

66.
Coelho FJRC, Cleary DFR, Gomes NCM, Pólonia ARM, Huang YM, Liu LL, et al. Sponge prokaryote communities in Taiwanese coral reef and shallow hydrothermal vent ecosystems. Microb Ecol. 2018;75:239–54.
CAS  PubMed  Article  PubMed Central  Google Scholar 

67.
Sipkema D, de Caralt S, Morillo JA, Al-Soud WAB, Sørensen SJ, Smidt H, et al. Similar sponge-associated bacteria can be acquired via both vertical and horizontal transmission. Environ Microbiol. 2015;17:3807–21.
CAS  PubMed  Article  PubMed Central  Google Scholar 

68.
Bjork JR, Diez-Vives C, Astudillo-Garcia C, Archie E, Montoya JM. Vertical transmission of sponge microbiota is inconsistent and unfaithful. Nat Ecol Evol. 2019;3:1172–83.
PubMed  PubMed Central  Article  Google Scholar 

69.
Gonzalez-Zapata FL, Bongaerts P, Ramírez-Portilla C, Adu-Oppong B, Walljasper G, Reyes A, et al. Holobiont diversity in a reef-building coral over its entire depth range in the mesophotic zone. Front Mar Sci. 2018;5:1–13.
Article  Google Scholar 

70.
Yang S, Sun W, Zhang F, Li Z. Phylogenetically diverse denitrifying and ammonia-oxidizing bacteria in corals Alcyonium gracillimum and Tubastraea coccinea. Mar Biotechnol. 2013;15:540–51.
CAS  PubMed  Article  Google Scholar 

71.
Høj L, Levy N, Baillie BK, Clode PL, Strohmaier RC, Siboni N, et al. Crown-of-thorns sea star Acanthaster cf. solaris has tissue—characteristic microbiomes with potential roles in health and reproduction. Appl Environ Microbiol. 2018;84:1–18.
Article  Google Scholar 

72.
Schmitt S, Deines P, Behnam F, Wagner M, Taylor MW. Chloroflexi bacteria are more diverse, abundant, and similar in high than in low microbial abundance sponges. FEMS Microbiol Ecol. 2011;78:497–510.
CAS  PubMed  Article  PubMed Central  Google Scholar 

73.
Bergquist PR, Kelly-Borges M. An evaluation of the genus Tethya (Porifera: Demospongiae: Hadromerjda) with descriptions of new species from the southwest Pacific. Beagle Rec Mus Art Galleries North Territ. 1991;8:37–72.
Google Scholar 

74.
Hoshino T, Yilmaz LS, Noguera DR, Daims H, Wagner M. Quantification of target molecules needed to detect microorganisms by fluorescence in situ hybridization (FISH) and catalyzed reporter deposition-FISH. Appl Environ Microbiol. 2008;74:5068–77.
CAS  PubMed  PubMed Central  Article  Google Scholar 

75.
Moran NA, Wernegreen JJ. Lifestyle evolution in symbiotic bacteria: insights from genomics. Trends Ecol Evol. 2000;15:321–6.
CAS  PubMed  Article  PubMed Central  Google Scholar 

76.
Vacelet J, Donadey C. Electron microscope study of the association between some sponges and bacteria. J Exp Mar Bio Ecol. 1977;30:301–14.
Article  Google Scholar 

77.
Maldonado M. Intergenerational transmission of symbiotic bacteria in oviparous and viviparous demosponges, with emphasis on intracytoplasmically-compartmented bacterial types. J Mar Biol Assoc UK. 2007;87:1701–13.
Article  Google Scholar 

78.
Ilan M, Abelson A. The life of a sponge in a sandy lagoon. Biol Bull. 1995;189:363–9.
CAS  PubMed  Article  PubMed Central  Google Scholar 

79.
Borisov VB, Gennis RB, Hemp J, Verkhovsky MI. The cytochrome bd respiratory oxygen reductases. Biochim Biophys Acta Bioenerg. 2011;1807:1398–413.
CAS  Article  Google Scholar 

80.
Lavy A, Keren R, Yahel G, Ilan M. Intermittent hypoxia and prolonged suboxia measured in situ in a marine sponge. Front Mar Sci. 2016;3:1–11.
Article  Google Scholar 

81.
Hoffmann F, Larsen O, Thiel V, Rapp HT, Pape T, Michaelis W, et al. An anaerobic world in sponges. Geomicrobiol J. 2005;22:1–10.
Article  Google Scholar 

82.
Moreno-Vivian C, Cabello C, Martinez-Luque M, Blasco R, Castillo F. Prokaryotic nitrate reduction: molecular properties and functional distinction among bacterial nitrate reductases. J Bacteriol. 1999;181:6573–84.
CAS  PubMed  PubMed Central  Article  Google Scholar 

83.
Toffanin A, Wu Q, Maskus M, Casella S, Abruña HD, Shapleigh JP. Characterization of the gene encoding nitrite reductase and the physiological consequences of its expression in the nondenitrifying Rhizobium “hedysari” strain HCNT1. Appl Environ Microbiol. 1996;62:4019–25.
CAS  PubMed  PubMed Central  Article  Google Scholar 

84.
Jones PM, George AM. The ABC transporter structure and mechanism: perspectives on recent research. Cell Mol Life Sci. 2004;61:682–99.
CAS  PubMed  Article  PubMed Central  Google Scholar 

85.
Mulligan C, Fischer M, Thomas GH. Tripartite ATP-independent periplasmic (TRAP) transporters in bacteria and archaea. FEMS Microbiol Rev. 2011;35:68–86.
CAS  PubMed  Article  PubMed Central  Google Scholar 

86.
Vastermark A, Wollwage S, Houle ME, Rio R, Saier MH. Expansion of the APC superfamily of secondary carriers. Proteins Struct Funct Bioinform. 2014;82:2797–811.
CAS  Article  Google Scholar 

87.
Pao SS, Paulsen IT, Saier MH. Major facilitator superfamily. Microbiol Mol Biol Rev. 1998;62:1–34.
CAS  PubMed  PubMed Central  Article  Google Scholar 

88.
Ziegler C, Bremer E, Krämer R. The BCCT family of carriers: from physiology to crystal structure. Mol Microbiol. 2010;78:13–34.
CAS  PubMed  PubMed Central  Google Scholar 

89.
Andrade SLA, Einsle O. The Amt/Mep/Rh family of ammonium transport proteins (Review). Mol Membr Biol. 2007;24:357–65.
CAS  PubMed  Article  PubMed Central  Google Scholar 

90.
Kamke J, Sczyrba A, Ivanova N, Schwientek P, Rinke C, Mavromatis K, et al. Single-cell genomics reveals complex carbohydrate degradation patterns in poribacterial symbionts of marine sponges. ISME J. 2013;7:2287–300.
CAS  PubMed  PubMed Central  Article  Google Scholar 

91.
Clifford EL, Hansell DA, Varela MM, Nieto-Cid M, Herndl GJ, Sintes E. Crustacean zooplankton release copious amounts of dissolved organic matter as taurine in the ocean. Limnol Oceanogr. 2017;62:2745–58.
CAS  PubMed  PubMed Central  Article  Google Scholar 

92.
Karimi E, Keller-Costa T, Slaby BM, Cox CJ, da Rocha UN, Hentschel U, et al. Genomic blueprints of sponge-prokaryote symbiosis are shared by low abundant and cultivatable Alphaproteobacteria. Sci Rep. 2019;9:1–15.
Article  CAS  Google Scholar 

93.
Dahl C, Franz B, Hensen D, Kesselheim A, Zigann R. Sulfite oxidation in the purple sulfur bacterium Allochromatium vinosum: Identification of SoeABC as a major player and relevance of SoxYZ in the process. Microbiology. 2013;159:2626–38.
CAS  PubMed  Article  PubMed Central  Google Scholar 

94.
Bardischewsky F, Quentmeier A, Friedrich CG. The flavoprotein SoxF functions in chemotrophic thiosulfate oxidation of Paracoccus pantotrophus in vivo and in vitro. FEMS Microbiol Lett. 2006;258:121–6.
CAS  PubMed  Article  PubMed Central  Google Scholar 

95.
Gregersen LH, Bryant DA, Frigaard NU. Mechanisms and evolution of oxidative sulfur metabolism in green sulfur bacteria. Front Microbiol. 2011;2:1–14.
Article  CAS  Google Scholar 

96.
Lavy A, Keren R, Yu K, Thomas BC, Alvarez-Cohen L, Banfield JF, et al. A novel Chromatiales bacterium is a potential sulfide oxidizer in multiple orders of marine sponges. Environ Microbiol. 2018;20:800–14.
CAS  PubMed  Article  PubMed Central  Google Scholar 

97.
Tian R-M, Wang Y, Bougouffa S, Gao ZM, Cai L, Bajic V, et al. Genomic analysis reveals versatile heterotrophic capacity of a potentially symbiotic sulfur-oxidizing bacterium in sponge. Environ Microbiol. 2014;16:3548–61.
CAS  PubMed  Article  PubMed Central  Google Scholar 

98.
Keller MD, Kiene RP, Matrai PA, Bellows WK. Production of glycine betaine and dimethylsulfoniopropionate in marine phytoplankton. II. N-limited chemostat cultures. Mar Biol. 1999;135:249–57.
CAS  Article  Google Scholar 

99.
Diaz MR, Visscher PT, Taylor BF. Metabolism of dimethylsulfoniopropionate and glycine betaine by a marine bacterium. FEMS Microbiol Lett. 1992;96:61–5.
CAS  Article  Google Scholar 

100.
Sun J, Steindler L, Thrash JC, Halsey KH, Smith DP, Carter AE, et al. One carbon metabolism in SAR11 pelagic marine bacteria. PLoS ONE. 2011;6:e23973.
CAS  PubMed  PubMed Central  Article  Google Scholar 

101.
Equar MY, Tani Y, Mihara H. Purification and properties of glycine oxidase from pseudomonas putida KT2440. J Nutr Sci Vitaminol. 2015;61:506–10.
CAS  PubMed  Article  Google Scholar 

102.
Karimi E, Slaby BM, Soares AR, Blom J, Hentschel U, Costa R. Metagenomic binning reveals versatile nutrient cycling and distinct adaptive features in alphaproteobacterial symbionts of marine sponges. FEMS Microbiol Ecol. 2018;94:1–18.
Article  CAS  Google Scholar 

103.
Moreno-Pino M, Cristi A, Gillooly JF, Trefault N. Characterizing the microbiomes of Antarctic sponges: a functional metagenomic approach. Sci Rep. 2020;10:1–12.
Article  CAS  Google Scholar 

104.
Mayer J, Huhn T, Habeck M, Denger K, Hollemeyer K, Cook AM. 2,3-Dihydroxypropane-1-sulfonate degraded by Cupriavidus pinatubonensis JMP134: purification of dihydroxypropanesulfonate 3-dehydrogenase. Microbiology. 2010;156:1556–64.
CAS  PubMed  Article  Google Scholar 

105.
Michael AJ. Polyamines in eukaryotes, bacteria, and archaea. J Biol Chem. 2016;291:14896–903.
CAS  PubMed  PubMed Central  Article  Google Scholar 

106.
Tsukamoto S, Kato H, Hirota H, Fusetani N. Pseudoceratidine: a new antifouling spermidine derivative from the marine sponge Pseudoceratina purpurea. Tetrahedron Lett. 1996;37:1439–40.
CAS  Article  Google Scholar 

107.
Igarashi K, Kashiwagi K. Polyamine transport in bacteria and yeast. Biochem J. 1999;344:633–42.
CAS  PubMed  PubMed Central  Article  Google Scholar 

108.
Shah P, Swiatlo E. A multifaceted role for polyamines in bacterial pathogens. Mol Microbiol. 2008;68:4–16.
CAS  PubMed  Article  Google Scholar 

109.
Schneider BL, Reitzer L. Pathway and enzyme redundancy in putrescine catabolism in Escherichia coli. J Bacteriol. 2012;194:4080–8.
CAS  PubMed  PubMed Central  Article  Google Scholar 

110.
Mou X, Sun S, Rayapati P, Moran MA. Genes for transport and metabolism of spermidine in Ruegeria pomeroyi DSS-3 and other marine bacteria. Aquat Microb Ecol. 2010;58:311–21.
Article  Google Scholar 

111.
Dasu VV, Nakada Y, Ohnishi-Kameyama M, Kimura K, Itoh Y. Characterization and a role of Pseudomonas aeruginosa spermidine dehydrogenase in polyamine catabolism. Microbiology. 2006;152:2265–72.
CAS  PubMed  Article  PubMed Central  Google Scholar 

112.
Tofalo R, Cocchi S, Suzzi G. Polyamines and gut microbiota. Front Nutr. 2019;6:1–5.
Article  CAS  Google Scholar 

113.
Sharfstein ST, Keasling JD. Polyphosphate metabolism in Escherichia coli. Ann N Y Acad Sci. 1994;745:77–91.
CAS  PubMed  Article  PubMed Central  Google Scholar 

114.
Zhang F, Blasiak LC, Karolin JO, Powell RJ, Geddes CD, Hill RT. Phosphorus sequestration in the form of polyphosphate by microbial symbionts in marine sponges. Proc Natl Acad Sci USA. 2015;112:4381–6.
CAS  PubMed  Article  PubMed Central  Google Scholar 

115.
Wang L, Yan J, Wise MJ, Liu Q, Asenso J, Huang Y, et al. Distribution patterns of polyphosphate metabolism pathway and its relationships with bacterial durability and virulence. Front Microbiol. 2018;9:1–10.
Article  Google Scholar 

116.
Romano S, Schulz-Vogt HN, González JM, Bondarev V. Phosphate limitation induces drastic physiological changes, virulence-related gene expression, and secondary metabolite production in Pseudovibrio sp. strain FO-BEG1. Appl Environ Microbiol. 2015;81:3518–28.
CAS  PubMed  PubMed Central  Article  Google Scholar 

117.
Bergquist PR. In: eLS (ed) Porifera (Sponges). John Wiley & Sons, Ltd. 2001.

118.
Kovacs-Simon A, Titball RW, Michell SL. Lipoproteins of bacterial pathogens. Infect Immun. 2011;79:548–61.
CAS  PubMed  Article  PubMed Central  Google Scholar 

119.
Aoki S, Ito M, Iwasaki W. From β- To α-proteobacteria: The origin and evolution of rhizobial nodulation genes nodij. Mol Biol Evol. 2013;30:2494–508.
CAS  PubMed  Article  PubMed Central  Google Scholar 

120.
Darby AC, Cho NH, Fuxelius HH, Westberg J, Andersson SGE. Intracellular pathogens go extreme: genome evolution in the Rickettsiales. Trends Genet. 2007;23:511–20.
CAS  PubMed  Article  PubMed Central  Google Scholar 

121.
Kuo CH, Ochman H. Inferring clocks when lacking rocks: the variable rates of molecular evolution in bacteria. Biol Direct. 2009;4:35.
PubMed  PubMed Central  Article  CAS  Google Scholar  More

1.
Li, Y. H., Hallerman, E. M., Liu, Q. S., Wu, K. M. & Peng, Y. F. The development and status of Bt rice in China. Plant Biotechnol. J. 14, 839–848 (2016).
Article  Google Scholar 
2.
Li, X. et al. Comparison of nutritional quality between chinese indica rice with sck and crylAc genes and its nontransgenic counterpart. J. Food Sci. 72, 420–424 (2007).
Article  Google Scholar 

3.
Liu, Q. S., Hallerman, E., Peng, Y. F. & Li, Y. H. Development of Bt rice and Bt maize in china and their efficacy in target pest control. Int. J. Mol. Sci. 17, 1–15 (2016).
ADS  Google Scholar 

4.
Chen, M., Shelton, A. & Ye, G. Y. Insect-resistant genetically modified rice in china: from research to commercialization. Annu. Rev. Entomol. 56, 81–101 (2011).
CAS  Article  Google Scholar 

5.
Mannakkara, A., Niu, L., Ma, W. H. & Lei, C. L. Zero effect of Bt rice on expression of genes coding for digestion, detoxification and immune responses and developmental performances of brown planthopper Nilaparvata lugens (Stål). J. Insect Physiol. 59, 985–993 (2013).
CAS  Article  Google Scholar 

6.
Tian, J. C. et al. Assessing the effects of Cry1C rice and Cry2A rice to Pseudogonatopus flavifemur, a parasitoid of rice planthoppers. Sci. Rep. 7, 7838. https://doi.org/10.1038/s41598-017-08173-w (2017).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

7.
Tian, J. C. et al. The rice planthopper parasitoid Anagrus nilaparvatae is not at risk when feeding on honeydew derived from Bacillus thuringiensis (Bt) rice. Pest Manag. Sci. 74, 1854–1860 (2018).
CAS  Article  Google Scholar 

8.
Zhang, L., Guo, R., Fang, Z. & Liu, B. Genetically modified rice Bt-Shanyou63 expressing Cry1Ab/c protein does not harm Daphnia magna. Ecotox. Environ. Safe. 132, 196–201 (2016).
CAS  Article  Google Scholar 

9.
Chen, Y. et al. Bt rice expressing Cry1Ab does not stimulate an outbreak of its non-target herbivore, Nilaparvata lugens. Transgen. Res. 21, 279–291 (2012).
CAS  Article  Google Scholar 

10.
Akhtar, Z. R. et al. Impact of six transgenic Bacillus thuringiensis rice lines on four nontarget thrips species attacking rice panicles in the paddy field. Environ. Entomol. 42, 173–180 (2013).
CAS  Article  Google Scholar 

11.
Zhou, X., Cheng, J. A., Yang, H. U. & Lou, Y. G. Effects of transgenic Bt rice on the population development of Nephotettix cincticeps. Chin. J. Rice Sci. 19, 74–78 (2005).
Google Scholar 

12.
Lu, Z. B. et al. Impacts of Bt rice expressing Cry1C or Cry2A protein on the performance of nontarget leafhopper, Nephotettix cincticeps (Hemiptera: Cicadellidae), under laboratory and field conditions. Environ. Entomol. 43, 209–217 (2014).
CAS  Article  Google Scholar 

13.
Hagenbucher, S. et al. Pest trade-offs in technology: reduced damage by caterpillars in Bt cotton benefits aphids. Proc. Biol. Sci. 280, 20130042. https://doi.org/10.1098/rspb.2013.0042 (2013).
CAS  Article  PubMed  PubMed Central  Google Scholar 

14.
Hagenbucher, S., Wackers, F. L. & Romeis, J. Indirect multi-trophic interactions mediated by induced plant resistance: Impact of caterpillar feeding on aphid parasitoids. Biol. Lett. 10, 20130795. https://doi.org/10.1098/rsbl.2013.0795 (2014).
Article  PubMed  PubMed Central  Google Scholar 

15.
Wang, X. Y. et al. Bt rice could provide ecological resistance against nontarget planthoppers. Plant Biotechnol. J. 16, 1748–1755 (2018).
CAS  Article  Google Scholar 

16.
Chen, G. et al. Odor, not performance, dictates Bemisia tabaci’s selection between healthy and virus infected plants. Front. Physiol. 8, 146 (2017).
PubMed  PubMed Central  Google Scholar 

17.
Xu, H. X., He, X. C., Zheng, X. S., Yang, Y. J. & Lu, Z. X. Influence of rice black streaked dwarf virus on the ecological fitness of non-vector planthopper Nilaparvata lugens (Hemiptera: Delphacidae). Insect Sci. 21, 507–514 (2013).
Article  Google Scholar 

18.
Xu, H. X. et al. Effects of transgenic rice infected with SRBSDV on Bt expression and the eological fitness of non-vector brown planthopper Nilaparvata lugens. Sci. Rep. 7, 6328. https://doi.org/10.1038/s41598-017-02218-w (2017).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

19.
Wang, Q. J. et al. Combined influence of Bt rice and rice dwarf virus on biological parameters of a non-target herbivore, Nephotettix cincticeps (Uhler) (Hemiptera: Cicadellidae). PLoS ONE 12, e0181258 (2017).
Article  Google Scholar 

20.
Chen, H. Y. et al. Sequential infection of rice dwarf virus in the internal organs of its insect vector after ingestion of virus. Virus Res. 160, 389–394 (2011).
CAS  Article  Google Scholar 

21.
Honda, K. et al. Retention of rice dwarf virus by descendants of pairs of viruliferous vector insects after rearing for 6 years. Phytopathology 97, 712–716 (2007).
Article  Google Scholar 

22.
Wei, T. Y. & Li, Y. Rice reoviruses in insect vectors. Annu. Rev. Phytopathol. 54, 99–120 (2016).
CAS  Article  Google Scholar 

23.
Zhao, S. S. et al. A viral protein promotes host SAMS1 activity and ethylene production for the benefit of virus infection. elife. 6, e27529. https://doi.org/10.7554/eLife.27529 (2017).
Article  PubMed  PubMed Central  Google Scholar 

24.
Ling, K. C., Tiongco, E. R. & Aguiero, V. M. Rice ragged stunt a new virus disease. Plant. Dis. Rep. 62, 701–705 (1978).
Google Scholar 

25.
Hibino, H. et al. Rice grassy stunt virus: A planthopper-borne circular filament. Phytopathology 75, 894–899 (1985).
Article  Google Scholar 

26.
Shimizu, T. et al. Strong resistance against rice grassy stunt virus is induced in transgenic rice plants expressing double-stranded RNA of the viral genes for nucleocapsid or movement proteins as targets for RNA interference. Phytopathology 103, 513–519 (2013).
Article  Google Scholar 

27.
Dang, C. et al. Does Bt rice pose risks to non-target arthropods? Results of a meta-analysis in China. Plant Biotechnol. J. 15, 1047–1053 (2017).
CAS  Article  Google Scholar 

28.
Ge, L. Q., Wu, J. C., Sun, Y. C., Ouyang, F. & Ge, F. Effects of triazophos on biochemical substances of transgenic Bt rice and its nontarget pest Nilaparvata lugens Stål under elevated CO2. Pest. Biochem. Physiol. 107, 188–199 (2013).
CAS  Article  Google Scholar 

29.
Ge, L. Q., Sun, Y. C., Ouyang, F., Wu, J. C. & Ge, F. The effects of triazophos applied to transgenic Bt rice on the nutritional indexes, Nlvg expression, and population growth of Nilaparvata lugens Stal under elevated CO2. Pest. Biochem. Physiol. 118, 50–57 (2015).
CAS  Article  Google Scholar 

30.
Che, Q. Y., Liu, X. H., Liang, Y. Y., Yang, B. & Ge, F. Effects of transgenic Bt-rice and insecticides on the communitystructure of soil nematodes. J. Plant Protect. 42(5), 724–733 (2015).
Google Scholar 

31.
Yang, Y. J. et al. Impacts of nitrogen fertilizer on major insect pests and their predators in transgenic Bt rice lines T2A–1 and T1C–19. Entomol. Exp. Appl. 160, 281–291 (2016).
CAS  Article  Google Scholar 

32.
Hill, D. S. Agricultural insect pests of temperate regions and their control. Q. Rev. Biol. 63, 344 (1987).
Google Scholar 

33.
Wang, Q. J. et al. Rice dwarf virus infection alters green rice leafhopper host preference and feeding behavior. PLoS ONE 13, e0203364 (2018).
Article  Google Scholar 

34.
Sun, X., Yan, M. J., Zhang, A. J. & Wang, M. Q. Transgenic cry1C gene rough rice line T1C–19 does not change the host preferences of the non-target stored product pest, Rhyzopertha dominica (Fabricius) (Coleoptera: Bostrichidae), and its parasitoid wasp, Anisopteromalus calandrae (Howard) (Hymenoptera: Pteromalidae). Ecotox. Environ. Safe. 120, 449–456 (2015).
CAS  Article  Google Scholar 

35.
Lu, Z. B. et al. Transgenic cry1C or cry2A rice has no adverse impacts on the life-table parameters and population dynamics of the brown planthopper, Nilaparvata lugens (Hemiptera: Delphacidae). Pest Manag. Sci. 71, 937–945 (2015).
ADS  CAS  Article  Google Scholar 

36.
He, X. C. et al. Ecological fitness of non-vector planthopper Sogatella furcifera on rice plants infected with rice black streaked dwarf virus. Rice Sci. 19, 335–338 (2012).
Article  Google Scholar 

37.
Xu, H. X. et al. Effects of SRBSDV-infected rice plants on the fitness of vector and non-vector rice planthoppers. J. Asia-Pac. Entomol. 19, 707–710 (2016).
Article  Google Scholar 

38.
Gómez-Torres, M. L., Nava, D. E. & Parra, J. R. P. Life table of Tamarixia radiata (Hymenoptera: Eulophidae) on Diaphorina citri (Hemiptera: Psyllidae) at different temperatures. J. Econ. Entomol. 105, 338–343 (2012).
Article  Google Scholar 

39.
Tang, W. et al. Development of insect-resistant transgenic indica rice with a synthetic cry1C* gene. Mol. Breed. 18, 1–10 (2006).
CAS  Article  Google Scholar 

40.
Hulting, F. L., Orr, D. B. & Obrycki, J. J. A computer-program for calculation and statistical comparison of intrinsic rates of increase and associated life table parameters. Fla. Entomol. 73, 601–612 (1990).
Article  Google Scholar 

41.
Maia, A. D. N., Luiz, A. J. B. & Campanhola, C. Statistical inference on associated fertility life table parameters using Jackknife technique: computational aspects. J. Econ. Entomol. 93, 511–518 (2000).
Article  Google Scholar  More

Sampling site and regimes
Field work was carried out in the agricultural catchment of Grytelandsbekken (a creek of approx. 2.5 km in length) also known as Skuterud catchment, located in the municipality of Ås (20,000 people), 30 km southeast of Oslo (Fig. 1). The catchment area is of approx. 4.5 km2 and largely consists of farmlands (60%) and forest/marshlands (31%). Grytelandsbekken has previously been studied for spatial variations in the microbial communities of various lotic freshwater ecosystems in different regions of Norway21. Based on that study, it was characterised as a rural creek with the highest diversity and abundance of microbial communities. Thus, in this follow-up study, Grytelandsbekken, specifically, was selected for assessing the seasonal dynamics of lotic freshwater bacterial communities.
All field measurements and samplings were carried out over a 2-year period at the same site of the creek, i.e. at about 0.25 km. In total, there were 16 sampling events, split equally between cold and warm seasonal regimes. The cold season refers to December–March, while the warm season includes June–September. In detail, there were four events during the cold season in the first year, Cold_1 (Cold_1a/Dec.2014, Cold_1b/Jan.2015, Cold_1c/Feb.2015, and Cold_1d/Mar.2015) and four events during the cold season in the second year, Cold_2 (Cold_2a/Dec.2016, Cold_2b/Jan.2017, Cold_2c/Feb.2017, and Cold_2d/Mar.2017). A similar sample assembly was applied to the warm seasonal regimes, i.e. four in the first year, Warm_1 (Warm_1a/Jun.2015, Warm_1b/Jul.2015, Warm_1c/Aug.2015, and Warm_1d/Sep.2015) and four in the second year, Warm_2 (Warm_2a/Jun.2016, Warm_2b/Jul.2016, Warm_2c/Aug.2016, and Warm_2d/Sep.2016). The entire study duration and sampling sets were conceived based on similar settings reported in aquatic research worldwide, profiling microbial communities through seasonal and spatial variations20,35,36.
Environmental measures
Primary physico-chemical parameters that are routinely measured in standard catchment water quality control37,38 were selected for abiotic characteristics of the lotic environment. These included organic matter content (expressed as CODCr), TSS, nutrients (Ptot and Ntot), EC, pH, and Temp. The latter was measured in situ at one of the weather stations, administrated by the Norwegian Institute of Bioeconomy Research (NIBIO) and located at the field measuring/sampling site of Grytelandsbekken. The station is equipped with a number of automatic sensors, providing hourly measurements and registering various climatic data online, which are all open and available at the AgroMetBase hosted by NIBIO (https://lmt.nibio.no/agrometbase/getweatherdata.php). The former parameters were analysed post-sampling in an accredited laboratory of the ALS Laboratory Group Norway AS. These analyses were performed in accordance with ISO and national standards for the respective parameters: CODCr (ISO 15705), TSS (CSN EN 872, NS 4733), Ptot (ISO 6878, ISO 15681-1), Ntot (EN 12260), EC (EN 27 888, SM 2520B, EN 16192), and pH (ISO 10523, EPA 150.1, EN 16192).
Genomic DNA purification: 16S rRNA amplicon library preparation and MiSeq sequencing
The seasonal water samples underwent DNA extraction using the QIAGEN DNeasy PowerWater Kit (QIAGEN GmbH, Hilden, Germany). In practice, 0.5 L of water was subjected to ultrafiltration to obtain a solid mass on a membrane filter (0.45 µm). DNA was then extracted from the collected filters, following the manufacturer’s instruction. The concentration and quality of the purified DNA was analysed on the NanoDrop Spectrophotometer (Thermo Fisher Scientific, Wilmington, USA). Five nanograms of purified DNA was used in PCR to prepare an amplicon library using the NEXTflex 16S V4 Amplicon-Seq Kit 2.0 (Bioo Scientific Corporation, Austin, TX, USA), following the provided protocol, which has previously been described in detail21. The applied specific 16S V4 forward primer (5′-GACGCTCTTCCGATCTTATGGTAATTGTGTGCCAGCMGCCGCGGTAA-3′) and reverse primer (5′-TGTGCTCTTCCGATCTAGTCAGTCAGCCGGACTACHVGGGTWTCTAAT-3′) were provided by the manufacturer and included in the sequencing kit. The library was prepared in triplicate for each sample. The final concentration of each library was measured on a Qubit Fluorometer (Life Technologies, Eugene, OR, USA) using the Quant-IT dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA, USA). Pooling of all prepared libraries was achieved after concentration normalization using the SequalPrep Normalization Plate Kit (Thermo Fisher Scientific, Wilmington, USA). The pooled library was analysed on the Agilent 2100 Bioanalyzer system using the Agilent High Sensitivity DNA Kit (Agilent, CA, USA). It indicated a single band/unique product at 451 bp. The multiplexed library was sequenced on the Illumina MiSeq system using the MiSeq Reagent Kit V3, 600 Cycles (Illumina Inc., San Diego, CA, USA), following the default standard procedures.
Sequence data processing
The output sequence datasets were analysed using Microbial Genomics Module 2.0, added onto the CLC Genomic Workbench 10.1.1 (CLC Bio, QIAGEN Company, Aarhus, Denmark, https://www.qiagenbioinformatics.com/products/clc-genomics-workbench). The processing workflow consisted of four key components: quality filtration, OTU clustering, and alpha and beta diversity measures. Adapter and primer sequences were trimmed. Unqualified reads were trashed when the quality score was less than 20 or a higher number of ambiguous nucleotides (more than two) were detected. The average length after trimming was between 220–230 bp. Chimeric sequences and singletons were detected and discarded. The remaining qualified reads were used to characterize OTUs based on a reference database (Greengenes v_13_5)39 at a 97% identity level. The bacterial alpha diversity of each sample was estimated in rarefaction analysis with a depth cutoff at 50,000 reads. The bacterial beta diversity applied the Euclidean distance criterion (EDC) to estimate the community similarities between the examined samples. All sequence data are available at NCBI Sequence Read Archive, under accession number SRR10835654-669, as part of BioProject PRJNA599104.
Statistical analyses
Alpha diversity differences were tested using a two-tailed Student’s t-test at 0.05 significance. This was performed in the XLSTAT-ECOLOGY statistical software package version 2019.1.1 (Addinsoft 2020, Boston, USA, https://www.xlstat.com). A hierarchical clustering heat map was created to elucidate the bacterial community similarities/relatedness (pairwise) among all examined samples. It was conducted on a subset of top 700 OTUs, based on the EDC using the trimmed mean of M-values and Z-score (standard deviation numbers from the population mean) normalizations. This was further supported by the PERMANOVA analysis, performed to ascertain statistical significance of the clusters (p-value = 0.00001). These tests were executed using a package of functional features included in the Microbial Genomics Module (CLC Bio, QIAGEN Company, Aarhus, Denmark, https://www.qiagenbioinformatics.com/products/clc-genomics-workbench). Furthermore, the LEfSe tool40 was applied to identify the responsible bacterial members, accounting for community discrepancy between cold and warm seasons. The formatted abundance table of bacterial classes was uploaded to the Galaxy/Hutlab application web-based platform (Biostatistics Department, Harvard School of Public Health, Boston, MA, USA, https://huttenhower.sph.harvard.edu/galaxy) for pairwise comparison. Statistical significance of the comparison was determined by the Wilcoxon rank-sum test at the alpha value of 0.05. The identified features characterising the microbial differences among samples were processed using LDA, with a threshold score set at 2.0. Beyond that, RDA was carried out to determine the key abiotic environmental variables driving seasonal changes of the bacterial community based on the Pearson correlation test, with a statistical significance level higher than 95% (p  More

Phenological and climatic data
We used data from the Pan European Phenology Project (PEP725)29, an open-access database with long-term plant phenological observations across 25 European countries (http://www.pep725.eu/). The regional/national network partners of PEP725 are following a consistent guideline for phenological observations30 and prepare the data for submission to the PEP725 database curators29. We selected 30 species for which sufficient observational data were available: 21 deciduous broadleaved trees or shrubs, 6 herbaceous perennials, 2 evergreen coniferous trees, and 1 deciduous coniferous tree (Supplementary Table 3). Particularly, our data set included one fruit tree (Prunus avium) and one nut tree (Corylus avellana) since some of the chilling models are specifically developed for fruit and nut trees. A total of 2,493,644 individual records from 15,533 phenological stations were used. The stations were mainly distributed in moderate climates in Central Europe (Supplementary Fig. 3). Four spring events based on the BBCH code were investigated: BBCH 10, 11, 60, and 69, representing first leaves separated, first leaves unfolded, first flowers open, and end of flowering, respectively31.
We used the E-OBS v19.0eHOM data set32 with a spatial resolution of 0.1 × 0.1° for 1950–2018 for calculating CA and HR of the in situ phenological records. This data set is provided by the European Climate Assessment & Data set project and includes homogenized series of daily mean, minimum, and maximum temperatures. We also use the daily maximum and minimum temperature data from the GHCN data set33 to assess the scale effect. The GHCN data set contains station-based measurements from over 90,000 land-based stations worldwide, but only parts of PEP725 stations match with the GHCN stations.
For future climatic data (2019–2099), we used daily minimum and maximum temperatures simulated by the HADGEM2-ES model (with a spatial resolution of 0.5 × 0.5°) under two climatic scenarios (RCP 4.5 and RCP 8.5). These data have been bias-corrected by applying the method used in the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP)34, which were available on the ISIMIP server (https://esg.pik-potsdam.de/projects/isimip2b/).
Chilling models
We used 12 chilling models to measure the amount of chilling. One type of chilling model is based on several specific temperature thresholds. The most commonly used model, developed in the 1930s and 1940s for peach35, calculates the number of hours or days with temperatures 0 °C for calculating CA11,38. Models C1–C6 were developed based on various combinations of the upper and lower temperature limits. The rate of chilling was 1 for daily temperatures 5,{mathrm{or}},{T} < - 10} hfill end{array}} right.,$$ (2) $${mathop{rm{CU}}nolimits} _3 = left{ {begin{array}{*{20}{l}} 1 hfill & {0 le {mathop{T}nolimits} le 5} hfill \ 0 hfill & {{T} > 5,{mathrm{or}},{T} < 0} hfill end{array}} right.,$$ (3) $${mathrm{CU}}_{mathrm{4}} = left{ {begin{array}{*{20}{c}} 1 & {{mathop{T}nolimits} le 7} \ 0 & {{T} > 7} end{array}} right.,$$
(4)

$${mathop{rm{CU}}nolimits} _5 = left{ {begin{array}{*{20}{l}} 1 hfill & { – 10 le {mathop{T}nolimits} le 7} hfill \ 0 hfill & {{T} > 7,{mathrm{or}},{T} < - 10} hfill end{array}} right.,$$ (5) $${mathop{rm{CU}}nolimits} _6 = left{ {begin{array}{*{20}{l}} 1 hfill & {0 le {mathop{T}nolimits} le 7} hfill \ 0 hfill & {{T} > 7,{mathrm{or}},{T} < 0} hfill end{array}} right.,$$ (6) where CUi is the rate of chilling for Model Ci, and T is the daily mean temperature (°C). Model C7 is also known as the Utah Model39, which assigned different weights to different ranges of temperatures and was first used to measure the chilling requirements of peach (Eq. (7)). The Utah Model was modified to produce Model C8 (Eq. (8)), which removed the negative contributions of warm temperatures to accumulated chilling20. $${mathop{rm{CU}}nolimits} _7 = left{ {begin{array}{*{20}{l}} 0 hfill & {{mathop{T}nolimits} le 1.4} hfill \ {0.5} hfill & {1.4 < {mathop{T}nolimits} le 2.4} hfill \ 1 hfill & {2.4 < {mathop{T}nolimits} le 9.1} hfill \ {0.5} hfill & {9.1 < {mathop{T}nolimits} le 12.4} hfill \ 0 hfill & {12.4 < {mathop{T}nolimits} le 15.9} hfill \ { - 0.5} hfill & {15.9 < {mathop{T}nolimits} le 18} hfill \ { - 1} hfill & {{mathop{T}nolimits} > 18} hfill end{array}} right.,$$
(7)

$${mathop{rm{CU}}nolimits} _8 = left{ {begin{array}{*{20}{l}} 0 hfill & {{mathop{T}nolimits} le 1.4} hfill \ {0.5} hfill & {1.4 < {mathop{T}nolimits} le 2.4} hfill \ 1 hfill & {2.4 < {mathop{T}nolimits} le 9.1} hfill \ {0.5} hfill & {9.1 < {mathop{T}nolimits} le 12.4} hfill \ 0 hfill & {{mathop{T}nolimits} > 12.4} hfill end{array}} right.,$$
(8)

where CUi is the rate of chilling for Model Ci, and T is the daily mean temperature (°C).
Model C9 is a dynamic model developed for peach in Israel and South Africa40,41 and now adopted for apricot cultivars42. The most important characteristic of Model C9 was that a previous intermediate product affected the rate of chilling in the current hour or day. We did not provide equations for Model C9 for simplicity (see the equation in Luedeling et al.20).
Harrington et al.43 summarized published results for chilling units and constructed a chilling function based on a three-parameter Weibull distribution, coded as Model C10 (Eq. (9)). Model C11 has a triangular form, which was fitted by Hänninen44 using previous experimental results for Finnish birch seedlings (Eq. (10)). Zhang et al.45 recently fitted observational data to the triangular model for 24 plant species and found that a mean optimal chilling temperature of 0.2 °C and an upper limit of the chilling temperature of 6.9 °C were most effective. Model C12, therefore, uses the triangular form with parameters of 0.2 and 6.9 °C (Eq. (11)).

$${mathop{rm{CU}}nolimits} _{10} = left{ {begin{array}{*{20}{l}} 1 hfill & {2.5 < {mathop{T}nolimits} < 7.4} hfill \ 0 hfill & {{mathop{T}nolimits} < - 4.7,{mathrm{or}},{mathop{T}nolimits} > 16} hfill \ {3.13left( {frac{{{mathop{T}nolimits} + 4.66}}{{10.93}}} right)^{2.10}{mathop{rm{e}}nolimits} ^{ – left( {frac{{{mathop{T}nolimits} + 4.66}}{{10.93}}} right)^{3.10}}} hfill & {{mathrm{else}}} hfill end{array}} right.,$$
(9)

$${mathrm{CU}}_{11} = left{ {begin{array}{*{20}{l}} 0 hfill & {{mathop{T}nolimits} le – 3.4,{mathop{rm{or}}nolimits} ,{mathop{T}nolimits} ge 10.4} hfill \ {frac{{{mathop{T}nolimits} + 3.4}}{{5 + 3.4}}} hfill & { – 3.4 < {mathop{rm{T}}nolimits} le 5} hfill \ {frac{{{mathop{T}nolimits} - 10.4}}{{5 - 10.4}}} hfill & {5 < {mathop{T}nolimits} < 10.4} hfill end{array}} right.,$$ (10) $${mathop{rm{CU}}nolimits} _{12} = left{ {begin{array}{*{20}{l}} 0 hfill & {{mathop{T}nolimits} le - 6.5,{mathop{rm{or}}nolimits} ,{mathop{T}nolimits} ge 6.9} hfill \ {frac{{{T + }6.5}}{{6.9 - 0.2}}} hfill & { - 6.5 < {mathop{T}nolimits} le 0.2} hfill \ {frac{{6.9 - {T}}}{{6.9 - 0.2}}} hfill & {0.2 < {T} < 6.9} hfill end{array}} right.,$$ (11) where CUi is the rate of chilling for Model Ci, and T is the daily mean temperature in °C. Forcing models Forcing models were used to measure HR for the spring events of plants. The GDD model is the most commonly used forcing model, which assumes that the rate of forcing is linearly correlated with temperature if the temperature is above a particular threshold. We mainly used Model F1 (Eq. (12)), which adopts a temperature threshold of 0 °C46,47,48, for examining the relationship between CA and HR: $${mathrm{FU}}_{mathrm{1}} = {mathrm{max}}({T},0),$$ (12) where FU1 is the rate of forcing for Model F1, and T is the daily mean temperature (°C). We also validated the chilling models by correlating them with HR based on seven other forcing models to assess the impact of the choice of forcing model on the results. Model F2 (Eq. (13)) is also a GDD model but has a temperature threshold of 5 °C12,18: $${mathrm{FU}}_{mathrm{2}} = {mathrm{max}}(T - 5,0),$$ (13) where FU2 is the rate of forcing for Model F2, and T is the daily mean temperature (°C). Piao et al.48 found that leaf onset in the Northern Hemisphere was triggered more by daytime than nighttime temperature. They thus proposed a GDD model using maximum instead of mean temperature. Models F3 (Eq. (14)) and F4 (Eq. (15)) are based on maximum temperature with thresholds of 0 and 5 °C, respectively. $${mathop{rm{FU}}nolimits} _3 = max ({mathop{T}nolimits} _{rm{max}},0),$$ (14) $${mathop{rm{FU}}nolimits} _4 = max ({mathop{T}nolimits} _{rm{max}} - 5,0),$$ (15) where FUi is the rate of forcing for Model Fi, and Tmax is the daily maximum temperature (°C). A recent experiment demonstrated that the impact of daytime temperature on leaf unfolding for temperate trees was approximately threefold higher than the impact of nighttime temperature49. Model F5 (Eq. 16) thus uses two parameters (0.75 and 0.25) to weigh the impact of daytime and nighttime temperatures on HR. $${mathop{rm{FU}}nolimits} _5 = 0.75 times max ({mathop{T}nolimits} _{rm{max}} - 5,0){mathrm{ + }}0.25 times max ({mathop{T}nolimits} _{rm{min}} - 5,0),$$ (16) where FU5 is the rate of forcing for Model F5. Tmax and Tmin are the daily maximum and minimum temperatures (°C), respectively. Many studies have suggested that the rate of forcing followed a logistic function of temperature44,50. Model F6 (Eq. (17)) uses a logistic function proposed by Hänninen44, and Model F7 (Eq. (18)) uses another logistic function proposed by Harrington et al.43. $${mathop{rm{FU}}nolimits} _6 = left{ {begin{array}{*{20}{l}} {frac{{28.4}}{{1 + {mathop{rm{e}}nolimits} ^{ - 0.185({mathop{T}nolimits} - 18.5)}}}} hfill & {{mathop{T}nolimits} > 0} hfill \ 0 hfill & {{mathop{rm{else}}nolimits} } hfill end{array}} right.,$$
(17)

$${mathop{rm{FU}}nolimits} _7 = frac{1}{{1 + {mathop{rm{e}}nolimits} ^{ – 0.47{mathop{T}nolimits} {mathrm{ + 6}}{mathrm{.49}}}}},$$
(18)

where FUi is the rate of forcing for Model Fi, and T is the daily mean temperature (°C).
Model F8 is a growing degree hour (GDH) model, where species have an optimum temperature for growth and where temperatures above or below that optimum have a smaller impact51. Model F8 (Eq. (19)) was first designed for calculating HR at hourly intervals, but we applied it at a daily interval. The stress factor in the original GDH model was ignored, because we assumed that the plants were not under other stresses.

$${mathrm{FU}}_{mathrm{8}}{mathrm{ = }}left{ {begin{array}{*{20}{l}} {mathrm{0}} hfill & {{T < }T_{mathrm{L}},{mathrm{or}},{T > }T_{mathrm{c}}} hfill \ {frac{{T_{mathrm{u}} – T_{mathrm{L}}}}{{mathrm{2}}}left( {{mathrm{1 + cos}}left( {pi {mathrm{ + }}pi frac{{{T} – T_{mathrm{L}}}}{{T_{mathrm{u}} – T_{mathrm{L}}}}} right)} right)} hfill & {T_{mathrm{L}} ge {T} ge T_{mathrm{u}}} hfill \ {left( {T_{mathrm{u}} – T_{mathrm{L}}} right)left( {{mathrm{1 + cos}}left( {frac{{uppi }}{{mathrm{2}}}{mathrm{ + }}frac{{uppi }}{{mathrm{2}}}frac{{{T} – T_{mathrm{u}}}}{{T_{mathrm{c}} – T_{mathrm{u}}}}} right)} right)} hfill & {T_{mathrm{u}}{ < T} le T_{mathrm{c}}} hfill end{array}} right.,$$ (19) where FU8 is the rate of forcing for Model F8, T is the daily mean temperature (°C), TL = 4, Tu = 25, and Tc = 36. Analysis We assessed the ability of each chilling model to represent long-term trends in the chilling conditions by calculating CA using each chilling model for each station for 1951–2018. CA was calculated as the sum of CUi from 1 November in the previous year to 30 April. The trend of CA at each station was visualized as the slope of the linear regression of CA against year. We also calculated Pearson’s r between each pair of chilling models for each station to determine if the chilling models were interrelated. We calculated HR and CA of spring events for each species, station, and year to determine if the chilling models are consistent with the physiological assumption that the reduction in chilling would increase HR (Fig. 1). HR was calculated as the sum of FUi from 1 January to the date of onset of spring events using Model F1, and the performances of the other forcing models (F2–F8) were also tested. We also compared 1 January with the other two starting dates of temperature accumulation (15 January and 1 February) to test any potential difference causing by the date when temperature accumulation begins. CA was calculated as the sum of CUi from 1 November in the previous year to the date of onset of spring events. We chose 1 November as the start date for CA because the endodormancy of temperate trees began around 1 November52. We only tested the linear relationship because the data were better fitted by the linear regression than the exponential model (Fig. 3), even though CA was linearly or nonlinearly negatively correlated with HR17. Pearson’s r between CA and HR for all records was calculated for each species, with a significantly negative Pearson’s r (p  More

Acoustic data and analysis
Audio data were collected from online sound archives (Xeno-Canto—https://www.xeno-canto.org—and Macaulay libraries—https://www.macaulaylibrary.org), creating a pool of over 2000 audio tracks. We assessed the sound quality of these audio tracks by listening and through visual inspection of sound spectrograms. To capture intra-specific variation, we limited audio extraction to one drum per audio track (which also avoided pseudoreplication) and only included species for which at least 3 high-quality drums could be extracted. We retained 736 high-quality drums suitable for further analyses. These drums were distributed among 92 species (out of the 217 recognized species of woodpeckers50 and 22 genera, providing a representative sampling of the phylogenetic diversity found in this family (Fig.1b)). Background noise and other artifacts were reduced by wavelet continuous reconstruction (R ‘WaveletComp’ package72), following the methods and description outlined in previous work73. The full script is available on demand. Finally, 22 acoustic variables were extracted from these filtered sound samples using the R ‘Seewave’ package74. Given the pulsed-like nature of drumming, the chosen variables emphasized the temporal and amplitude-related (all normalized to the maximal amplitude within a given drumming signal) features of the sounds (Supplementary Table 1). These 22 variables were z-scored and then used in all subsequent analyses. Since these variables were partly correlated to varying degrees, we performed a principal component analysis (PCA) to reduce the number of descriptive variables quantifying drumming acoustic structure. This dimensionality reduction was useful for visualization and necessary for regularization (i.e. to prevent overfitting) in many of our analyses. This resulted in six principal components (PCs) with eigenvalues >1 which together explained 75% of the variance (Supplementary Table 11).
We used these variables to evaluate the similarity between species-specific drums, by performing a hierarchical cluster analysis (HCA)75,76 based on Euclidean distances and the ‘Ward.D2’ method (‘NbClust’ R package77). NbClust provides a clustering output resulting from the use of multiple indices (in the case of our analysis, 26 indices were used). The best number of clusters is chosen according to the majority rule, i.e. it is the one supported by the highest number of indices used. This entails creating a vector of acoustic features (22 raw acoustic measures or 22 PCs) for each of the 92 species in our dataset and calculating the Euclidean distance between these vectors to evaluate how close acoustically species were. Note that one can still use Euclidian distances in a non-orthonormal space to calculate the ‘distance’ between signals. The result is a distance metric that might give more weights to measures that co-vary. This could theoretically affect the clustering results. However, when we performed the same analysis using the 22 PCs, we obtained the same grouping (6 clusters) with very minor differences in species grouping and distances between clusters as shown in the relative length of branches (Supplementary Fig. 13). The output of this HCA established an optimal classification of woodpeckers’ drums into 6 main drumming types (Fig. 2a), described as follows:
– Acceleration (AC): Beak strikes decrease in amplitude as they are produced within successively shorter time intervals.
– Regular sequence (RS): Beak strikes are produced in bouts, each comprising a relatively fixed (stereotyped) number of strikes.
– Irregular sequence (IS): Beak strikes are produced in bouts, each comprising a variable number of strikes (as opposed to RS).
– Steady fast (SF): Beak strikes are produced with constant time intervals and at a similar amplitude, with a high pulse rate (on average >20 strikes/s).
– Steady slow (SS): Beak strikes are produced with constant time intervals and at a similar amplitude, with a low pulse rate (on average 2; see models Supplementary Tables 3–6), a likelihood ratio test (LRT) was conducted to test for the specific effect of predictor variables (Supplementary Table 4). Because both models (null and fitted) differ in their fixed effects, model comparison was performed on models fit by maximum likelihood (ML) with the phylogenetic correlation structure (Pagel’s λ) fixed to the estimates obtained from initial fit by REML. The statistics reported for model comparison are likelihood ratios.
PGLS models testing for a relationship between life-history variables and drumming structure included either of PC1 to PC6 as the dependent variable to investigate whether differences exist between these proxies for acoustic structure. Similarly, LDs were used to verify our results with these different loading combinations of drumming acoustic variables. No significant correlations were found between life-history traits and acoustic structure using LDs instead of PCs (Supplementary Table 6), indicating that the combination of structural variation captured by the LDs differed from that of the PCs, while not leading to fundamentally different conclusions. Similarly, no significant correlations were found between life-history traits and information content (no decrease in AICc >2; Supplementary Table 3), overall emphasizing that none of the variables investigated here (and which could have potentially affected drumming structure) seemed to have influenced species-specific information, or at least not directly.
Ancestral states reconstructions
We carried out two types of ancestral state reconstructions: discrete reconstruction of drumming status in Fig. 1b and of drumming types in Fig. 3a (using ‘ace’ from the R ‘ape’ package81), or continuous character reconstruction of drumming acoustic structure based on Brownian motion models (using ‘fastanc’ from the R ‘phytools’ package82 in Supplementary Figs. 4 and 5) and using relaxed Brownian motion model (using ‘rjmcmc.bm’ from the R ‘geiger’83 package in Fig. 4a and Supplementary Figs. 8 and 9).
While evaluating the likelihood that drumming was already present at an early stage of woodpecker’s phylogeny, we tried to represent the most complete tree of the family, based on very recent molecular data50. Note that strictly speaking, we evaluate the state at the root but at the next internal node, i.e. at the node including Picumninae and Picinae (the largest pie-chart in our Fig. 1b), as Wrynecks do no drum, and neither do honeyguides or barbets). To include species with an unknown drumming status in this discrete reconstruction, we attributed equal probability distribution between the 3 states (i.e. when the ‘drummer state’ of a species is unknown, the species is given, prior to ancestral state reconstruction, a 1/3 probability of belonging to each of the three categories ‘drummer’, ‘occasional drummer’ and ‘non-drummer’). Stochastic mapping was performed under an MCMC model, sampling the rate matrix from its posterior distribution for Q (‘Q = mcmc’ in make.simmap function from the R ‘phytools’ package), with an equiprobable default prior at the root, and 200 simulations. Under a symmetrical model for the probability to change among the three states, scaled likelihood on woodpeckers’ ancestral node indicated 56.4%, 38.3% and 5.3% probabilities of being a drummer, an occasional drummer and a non-drummer, respectively. This is in line with the fact that morphological adaptations for drilling (including reinforced rhamphotheca, frontal overhang and processus dorsalis pterygoidei) evolved in the ancestral lineage of Picumninae and Picinae64.
To prevent overfitting, the discrete reconstructions for drumming types were estimated for six different rate models: equal rate model (ER), symmetric rate model (SYM), all rates difference model (ARD) and three sequential transition models based on the normalized MIL as measures of complexity as shown in Supplementary Fig. 7 ((SF leftrightarrow SS leftrightarrow DK leftrightarrow AC leftrightarrow RS leftrightarrow IS)). These three models assumed (1) sequential and equal, (2) sequential and incremental and (3) sequential and reversed transition rates, respectively. The number of parameters for these 6 rate models were 25, 1, 15, 1, 2 and 10. The final regularized likelihoods of each ancestral states were then obtained by model averaging using Akaike weights.
Calculation of information at different evolutionary steps was carried out as an extension of the drumming types reconstruction described above. From the discrete ancestral reconstruction procedure, probability distributions of drumming types were obtained for each node of the phylogenetic tree. We then obtained probability distributions at 20 fixed time intervals (dt = 1 myr) by linear interpolation. Using these probability distributions, we sampled drumming types proportionally from extant species descending the node closest to the time interval to estimate ancestral information values. This bootstrap procedure was repeated 30 times in order to obtain reliable estimates of mean and standard error. In this manner, we obtained information-through-time plots. These plots quantify a putative diversity of drumming signals in the clade at a particular point in time. They are similar in spirit to the disparity-through-time plots that have been used to measure specific morphological diversity in a clade through time using phylogenetic trees based on molecular data in combination with morphological measures in extant species84.
Continuous ancestral character trait reconstruction of drumming acoustic structure was carried out using either the six PCs that explain variation among the 22 drumming acoustic variables, or the six LDs that explain the variation in discriminating potential among the same variables (see above, ‘Acoustic data and analysis’ and ‘Calculation of information’ sections; Supplementary Figs. 4 and 5). The results and conclusions were similar for all PC’s and since the PC1 component has strong loading of multiple acoustic variables and the highest acoustic structure variance explained (Supplementary Table 11) it serves well as an illustrative example. The measure of phylogenetic signal on continuous traits (i.e. the historical contingency between species-specific drums that renders a trait non-randomly distributed along the phylogenetic tree) was made using Pagel’s lamba (Supplementary Table 2).
Reconstructing information content from raw MIL values would not have been biologically relevant since information calculation is based on the number of species involved, a factor that changes as branches merge going backward along the phylogenetic tree. We thus reconstructed MIL based on the normalized MIL values to avoid this pitfall. We used a Bayesian model implemented in the R package ‘Geiger’83 (model ‘rbm’ in the function ‘rjmcmc.bm’) to estimate branch-specific rates of trait evolution (i.e. changes in rates through time and across lineages). In this method, a reversible jump Markov Chain Monte Carlo (MCMC) sampling algorithm is used to detect shifts in rates of continuous traits evolution under a relaxed Brownian motion model85. The results of the model fit were summarized by the branch-specific average rate, estimated from the posterior samples. To obtain relative variations in posterior average rates, drumming structure (PC1-PC6) and MIL were standardized, i.e. these traits were divided by their standard deviation prior to running the ‘rbm’ models.
Analytical simulations of selection for information
In Fig. 3c, we compared that reconstructed evolution of information to what might be expected in different scenarios to further support those conclusions. More specifically, we estimated the ancestral MI for two simulated scenarios using an analytical model that describes species-specific information based on the probability of correct detection and the number of species (see ‘Calculation of information’ section). In the ‘No Diversifying Selection’ scenario (dark brown), the probability of correct detection for the initial pair of species, p2, is first estimated from the data using the approach described in the main text. It is then assumed that that additional species are randomly just as different/similar than these original species pair, yielding a probability of correct detection through time given by (p_c(t) = p_2^{n_s(t) – 1}), where ns(t) is the number of species at a given time. In the ‘Strong Diversifying Selection’ scenario (light brown), the probability of correct detection estimated at the first time point in our reconstruction (−20 M years ago, 3 species) is kept constant, (p_c(t) = p_2). In other words, the only species that survive would be species that can discriminate themselves from all other species equally well than the currently existing species. The reconstructed (actual) scenario is found between these two extreme values, showing that the drumming types are clearly not random but were also not under high evolutionary pressure to increase species-specific information. New drumming types evolved and species within types used signals that were distinct enough to result in the maintenance of normalized MI.
In Supplementary Fig. 6, we showed that the non-normalized reconstructed MI increased more rapidly when new drumming types appeared but that the normalized MI was relatively constant, reflecting the fact that the appearance of novel drumming types could co-occur with rapid radiation and increase in species numbers.
Playback experiments
Initial preparation involved identifying and mapping the areas prone to high densities of great-spotted woodpeckers Dendrocopos major, the study species of this experimental phase, using GIS maps provided by the LPO (French Bird Protection Organization). D. major is commonly found in European forests, ranging from open coniferous to mature deciduous forests. Playback experiments were carried out on wild individuals around Saint-Etienne, France, during this species’ breeding season (February–April 2017). All experiments were performed in accordance with relevant guidelines and regulations including French national guidelines, permits and regulations regarding animal care and experimental use (approval no. D42-218-0901, ENES lab agreement, Direction Départementale de la Protection des Populations, Préfecture du Rhône).
Two sets of experiments were conducted over the course of the breeding season, although we implemented the same general design which consisted in simulating a territorial intrusion. Playback stimuli tracks consisted of eight drums spread unevenly over about 60 s, aiming at representing the variation encountered in natural sequences (ref. 44 and personal observations). The first experiment (Exp. 1) aimed at investigating D. major’s response to conspecific vs. heterospecific drums. The other experiment (Exp. 2) aimed at investigating D. major’s response to drums from conspecifics vs. drums modified through acoustic manipulation (i.e. signal re-synthesis). D. major typically drums with an ‘acceleration’ pattern, which is mainly characterized by a shortening of the inter-strike time interval, a progressive decrease in strikes’ amplitude, and a gradual change in spectral properties as strikes get faster and weaker.
In Exp. 1, we used a paired and randomized order design, presenting each focal individual with one D. major drum and one drum from one out of 4 different species: 2 of which have very different drumming patterns (Picus canus and Dryobates minor, both producing ‘steady fast’ drums), and 2 others which have similar (accelerating) drumming patterns (Dendrocopos syriacus and Dendrocopos hyperythrus). A potentially confounding factor (which is nevertheless in line with our phylogenetic analyses) lies in that the allopatric species producing a drum similar to that of D. major also happened to be closely related to our model species. We carried out 48 playback experiments (testing 24 individuals with one of 4 categories of paired signals).
In Exp. 2, we altered one of the 3 acoustic features described above or all of them together (thus having 4 categories of modified signals), using Praat sound analysis software86. The design was paired so that each focal individual was exposed to one conspecific drum and one modified drum, following a randomized presentation order. This led to 48 playback experiments (24 individuals, each tested with one of 4 categories of paired signals).
Within each of Exp. 1 and Exp. 2, tested individuals were all separated by at least 500 m, ensuring different identities since their territory sizes vary between 200 and 400 m87,88. Upon visual or aural detection of (an) active individual(s), the experimenter set up an Anchor Megavox loudspeaker at about 1–1.5 m from ground level. The speaker was connected to an Edirol R-09 recorder (stimuli tracks were created and stored as WAV files, 44.1 kHz sampling frequency). Playbacks started at about 50 m from where the experimenter last saw or heard the focal individual. Following the work from Schuppe and colleagues89, playback intensity was calibrated and kept at about 80 dB measured 1 m away from the speaker. Behavioural data collection started when the first drum of the stimuli track was broadcasted and lasted 10 min from that moment on. To document focal individuals’ responses, notes were taken manually and continuously, while audio was recorded with a Sennheiser ME67 microphone mounted on a tripod and connected to a digital recorder (Zoom H4N, 44.1 kHz, 16 bit). If a response was elicited from multiple individuals in the area, only the one from a particular individual (ideally the one seen or heard before setting up the experiment) was monitored and used in further analyses. Six behavioural variables were reported, namely the number of screams, the number of drums, the approach (which was divided into three categories: ‘within 25 m’, ‘25–50 m’ and ‘further than 50 m’) as well as the latencies to first scream and drum and the latency to closest approach. When no occurrence was observed for the first three behaviours, latencies were set by default to the maximum value, i.e. the duration of the full experiment (10 min = 600 s). To characterize D. major’s behaviour, a PCA was then performed on scaled/centred data, where we retained the first principal component (‘Playback-PC1’) as an indicator of the behavioural response’s strength. A higher Playback-PC1 score indicates a stronger territorial response, i.e. more screams, a closer approach to the speaker and shorter latencies to these 2 behaviours. A second significant component resulted from this PCA, which represented the drumming’s response (inversely related: a higher Playback-PC2 score indicates fewer drums and a longer latency to drum; see Supplementary Table 13). None of the pairwise comparisons were statistically significant for PC2, besides a stronger drumming response to drums resynthesized without temporal variation than to D. major drums (Supplementary Fig. 15). This can be explained by the fact that birds were tested during their breeding season. At this time, drumming behaviour is likely to occur more consistently and commonly across experiments, independently from the stimulus played back, while screams and approach do not occur unless threat of intrusion is clear. Therefore, we used Playback-PC1 to represent birds’ behavioural response in our analysis (as it is in addition explaining much more variance in the behavioural data than Playback-PC2). Note that, as two playback sets were involved in this study, while we considered them independently in our statistical analysis, for standardization of the behavioural scale, we used the same polynomial equation. More specifically, the linear equation obtained from the loading scores of Exp. 1 was applied to the behavioural data of Exp. 2 for computation of Playback-PC2 scores.
Finally, distances were approximated during continuous note-taking and confirmed post-experimentally using a National Geographic 4*21 rangefinder (measurement accuracy: ±1 m up to 200 m). Sex was not documented as sometimes birds were not seen (but just heard drumming or calling back at our playback), which we nevertheless believe to be negligible since both sexes drum and are territorial in this monogamous species44,90.
Statistical analyses tested for differential responses of focal birds to drums of their own species versus either another species or a modified resynthesized condition. A paired comparison design was used by means of LMMs and contrasts using R software (‘lme4’ and ‘lsmeans’ packages)91,92. LMMs included study day and time, order of presentation and focal bird identity as random factors, and tested for a fixed effect of the interaction between treatment and group of paired condition. Contrasts were then computed between treatments (i.e. conspecific versus non-conspecific drums) for each group (i.e. each paired testing condition, such as D. major versus D. minor for which n = 6 birds were exposed to paired playback presentations—see Fig. 5a, b). Before contrasts and using the ‘lsmeans’ function, a Tukey adjustment for multiple testing was used; two-sided statistics are reported.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

Publication collection and data compilation
The detailed methods of publication search and data collection were described in our related work7. In brief, 10 databases in Web of Science (WoS; 1 January 1900 to 13 December 2016) including BIOSIS Previews, Chinese Science Citation Database, Data Citation Index, Derwent Innovations Index, Inspec, KCI-Korean Journal Database, MEDLINE, Russian Science Citation Index, SciELO Citation Index, and WoS Core Collection were used for searching peer-reviewed publications that reported GCMEs. The 18 keywords for WoS title search were: global change, climate change, free-air carbon dioxide enrichment, free-air CO2 enrichment, elevated carbon dioxide, elevated CO2, elevated atmospheric CO2, CO2 enrichment, eCO2, [CO2], warming, elevated temperature, changing precipitation, increased precipitation, decreased precipitation, nitrogen deposition, nitrogen addition, and nitrogen application. Through these search, 310,177 publication records that might be relevant to our topic were found.
First, we identified all the 310,177 records via reading each title. Second, we read the abstracts of all the records collected in the first step to further screen publications. During the two steps, we excluded 291,436 records because these studies were reviews/meta-analyses or conducted in non-terrestrial ecosystems such as oceans. Third, we read the methods of the remaining 18,741 publications to identify which of them met the following three inclusion criteria:
1.
Publications reported results of outdoor GCMEs which had at least three control and global change treatment plots ( > = 1 m2).

2.
The GCMEs were conducted in terrestrial ecosystems except for croplands and lab incubation studies.

3.
The GCMEs aimed to examine effects of simulated global change drivers on carbon, nitrogen, and water-cycle variables as well as plant and microbial parameters.

During the screening in the third step, 1,290 publications met these defined criteria.
We subsequently cross-checked the list of the 1,290 publications with references cited by the previous review/meta-analysis articles in global change research as well as the 1,290 publications, and collected 756 publications. In addition, 184 studies were collected by searching the websites of ecology laboratories and experiment networks and checking the references of the papers downloaded from these websites. In total, 2,230 publications were collected in the original version of the database7. Moreover, another 12 publications were found when we checked and reorganized all the data extracted from the 2,230 publications8,9,10,11,12,13,14,15,16,17,18,19. This database compiled 11 plant production and ecosystem carbon exchange variables including net primary productivity (NPP), above- and below-ground NPP (ANPP and BNPP), total biomass, aboveground biomass (AGB), root biomass, litter mass, gross and net ecosystem productivity (GEP and NEP), and ecosystem and soil respiration (ER and SR). Data of mean values, standard deviations or standard errors, and sample sizes (number of plot replications) of these variables in the control and treatment (e.g., elevated CO2, nitrogen addition, warming, increased/decreased precipitation, or their combinations) groups were extracted from each publication when possible. The figures were digitized using SigmaScan Pro 5.0 (SPSS, Inc.) and the numerical values were extracted when a publication presented experimental data graphically. Data of the experiments that were conducted over less than one year/growing season were excluded in this database. However, we included short-term data from tundra studies because most of measurements in this ecosystem were performed during July-August. Overall, 5,213 pairs (the control versus global change treatment) of plant production and ecosystem carbon exchange samples were collected in this database, having 2,247, 2,120, 81, and 765 pairs from single-, two-, three-, and four-factor manipulative experiments, respectively (Fig. 1).
Fig. 1

Number of samplings. Number of sample pairs of ecosystem carbon-cycling variables including net primary productivity (NPP), above- and below-ground NPP (ANPP and BNPP), total biomass, aboveground biomass (AGB), root biomass, litter mass, gross and net ecosystem productivity (GEP and NEP), and ecosystem and soil respiration (ER and SR) extracted from publications reporting single-, two-, three-, and four-factor global change manipulative experiments.

Full size image

Environmental metadata: Climate and vegetation
Information on the locations and altitudes of each experimental site, site climate including mean annual temperature (MAT) and precipitation (MAP) as well as wetness index ((frac{{rm{MAP}}}{{rm{MAT+10}}})), ref. 20, and vegetation types were extracted from each of the 2,242 publications. If a study did not report climate characteristics for its experimental site, data of MAT and MAP were downloaded from Climate Model Intercomparison Project phase 5 (CMIP5; https://esgf-node.llnl.gov/projects/cmip5/) based on the site coordinate. The dataset selection in CMIP5 was “historical (simulation of recent past 1850–2005)” and the climate data averaged from 20 (i.e. BCC_CSM1_1, BCC_CSM1_1_M, CANESM2, CCSM4, CMCC_CM, CMCC-CMS, CNRM-CM5, CSIRO_MK3_6_0, GFDL_CM3, GISS_E2_H, HADGEM2_AO, HADGEM2_ES, INMCM4, MIROC_ESM, MIROC_ESM_CHEM, MIROC5, MPI_ESM_LR, MPI_ESM_MR, MRI_ESM1, and NORESM1_M)21, that contained historical climate data, out of the 35 global climate models available in CMIP5 were used in this study. In addition, we downloaded data of climate means at global 1 × 1° land grid cells from Princeton University (http://hydrology.princeton.edu/data/pgf/v3/) to construct global climate space. Moreover, we classified ecosystems subjected to ecosystem manipulative experiments into five typical types: forests (mature forests and tree seedlings), grasslands (grasslands, meadows, short- and tall-grass prairies, temperate/semi-arid steppes, shrublands, savannas, pastures, and old-fields), tundra, wetlands (peatlands, bogs, marshes, and fens), and deserts.
Metadata of experimental facilities and performance
Information on CO2 enrichment and warming facilities were also extracted from the related publications reporting CO2 and warming effects on plant production and ecosystem carbon exchange. Facilities used in elevated CO2 experiments included greenhouse, open-top chamber, free-air CO2 enrichment, and tunnels. Warming experiments primarily used greenhouse, open-top chamber, soil heating cables, infrared radiator, and infrared reflector to elevate vegetation canopy and soil temperature. In addition, the manipulation magnitudes of global change drivers imposed by manipulative experiments, such as the increases in CO2 concentrations (ppm) and temperature (°C), the changes in precipitation amount (mm), and the rates of nitrogen input (g N m−2 yr−1), were also collected and added into this updated database. More

1.
Letchumanan, V., Chan, K. G. & Lee, L. H. Vibrio parahaemolyticus: a review on the pathogenesis, prevalence, and advance molecular identification techniques. Front. Microbiol. 5, 705, https://doi.org/10.3389/fmicb.2014.00705 (2014).
Article  PubMed  PubMed Central  Google Scholar 
2.
DePaola, A., Hopkins, L. H., Peeler, J. T., Wentz, B. & McPhearson, R. M. Incidence of Vibrio parahaemolyticus in U.S. coastal waters and oysters. Appl. Environ. Microbiol. 56, 2299–2302 (1990).
CAS  Article  Google Scholar 

3.
Daniels, N. A. et al. Emergence of a new Vibrio parahaemolyticus serotype in raw oysters: a prevention quandary. JAMA. 284, 1541–1545, https://doi.org/10.1001/jama.284.12.1541 (2000).
CAS  Article  PubMed  Google Scholar 

4.
McLaughlin, J. B. et al. Outbreak of Vibrio parahaemolyticus gastroenteritis associated with Alaskan oysters. N. Engl. J. Med. 353, 1463–1470, https://doi.org/10.1056/NEJMoa051594 (2005).
CAS  Article  PubMed  Google Scholar 

5.
Ralph, A. & Currie, B. J. Vibrio vulnificus and V. parahaemolyticus necrotising fasciitis in fishermen visiting an estuarine tropical northern Australian location. J. Infect. 54, e111–114, https://doi.org/10.1016/j.jinf.2006.06.015 (2007).
Article  PubMed  Google Scholar 

6.
Akther, F. et al. Major tdh(+)Vibrio parahaemolyticus serotype changes temporally in the Bay of Bengal estuary of Bangladesh. Infect. Genet. Evol. 41, 153–159, https://doi.org/10.1016/j.meegid.2016.04.003 (2016).
Article  PubMed  Google Scholar 

7.
Caburlotto, G. et al. Occurrence and molecular characterisation of Vibrio parahaemolyticus in crustaceans commercialised in Venice area, Italy. Int. J. Food Microbiol. 220, 39–49, https://doi.org/10.1016/j.ijfoodmicro.2015.12.007 (2016).
CAS  Article  PubMed  Google Scholar 

8.
Arakawa, E. et al. Emergence and prevalence of a novel Vibrio parahaemolyticus O3:K6 clone in Japan. Jpn. J. Infect. Dis. 52, 246–247 (1999).
CAS  PubMed  Google Scholar 

9.
Leal, N. C. et al. Vibrio parahaemolyticus serovar O3:K6 gastroenteritis in northeast Brazil. J. Appl. Microbiol. 105, 691–697, https://doi.org/10.1111/j.1365-2672.2008.03782.x (2008).
CAS  Article  PubMed  Google Scholar 

10.
Shaw, K. S., Sapkota, A. R., Jacobs, J. M., He, X. & Crump, B. C. Recreational swimmers’ exposure to Vibrio vulnificus and Vibrio parahaemolyticus in the Chesapeake Bay, Maryland, USA. Environ. Int. 74, 99–105, https://doi.org/10.1016/j.envint.2014.09.016 (2015).
Article  PubMed  Google Scholar 

11.
Liu, J. et al. Trends of foodborne diseases in China: lessons from laboratory-based surveillance since 2011. Front. Med. 12, 48–57, https://doi.org/10.1007/s11684-017-0608-6 (2018).
Article  PubMed  Google Scholar 

12.
Paudyal, N. et al. A meta-analysis of major foodborne pathogens in Chinese food commodities between 2006 and 2016. Foodborne Pathog. Dis. 15, 187–197, https://doi.org/10.1089/fpd.2017.2417 (2018).
CAS  Article  PubMed  Google Scholar 

13.
Deng, C., Deng, Y. & Yi, J. Analysis of microbial food poisoning from 2010 to 2016 in Sanya city. Hainan Med. J. 28, 2723–2725 (2017).
Google Scholar 

14.
Xu, X. et al. Prevalence, pathogenicity, and serotypes of Vibrio parahaemolyticus in shrimp from Chinese retail markets. Food Control. 46, 81–85, https://doi.org/10.1016/j.foodcont.2014.04.042 (2014).
Article  Google Scholar 

15.
Xie, T., Wu, Q., Xu, X., Zhang, J. & Guo, W. Prevalence and population analysis of Vibrio parahaemolyticus in aquatic products from South China markets. FEMS Microbiol. Lett. 362, https://doi.org/10.1093/femsle/fnv178 (2015).

16.
Xu, X., Cheng, J., Wu, Q., Zhang, J. & Xie, T. Prevalence, characterization, and antibiotic susceptibility of Vibrio parahaemolyticus isolated from retail aquatic products in North China. BMC Microbiol. 16, 32, https://doi.org/10.1186/s12866-016-0650-6 (2016).
CAS  Article  PubMed  PubMed Central  Google Scholar 

17.
Xie, T., Xu, X., Wu, Q., Zhang, J. & Cheng, J. Prevalence, molecular characterization, and antibiotic susceptibility of Vibrio parahaemolyticus from Ready-to-Eat foods in China. Front. Microbiol. 7, 549, https://doi.org/10.3389/fmicb.2016.00549 (2016).
Article  PubMed  PubMed Central  Google Scholar 

18.
Pang, R. et al. Comparative genomic analysis reveals the potential risk of Vibrio parahaemolyticus isolated from Ready-To-Eat foods in China. Front. Microbiol. 10, 186, https://doi.org/10.3389/fmicb.2019.00186 (2019).
Article  PubMed  PubMed Central  Google Scholar 

19.
Bankevich, A. et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J. Comput. Biol. 19, 455–477, https://doi.org/10.1089/cmb.2012.0021 (2012).
MathSciNet  CAS  Article  PubMed  PubMed Central  Google Scholar 

20.
Seemann, T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 30, 2068–2069, https://doi.org/10.1093/bioinformatics/btu153 (2014).
CAS  Article  PubMed  Google Scholar 

21.
Page, A. J. et al. Roary: rapid large-scale prokaryote pan genome analysis. Bioinformatics. 31, 3691–3693, https://doi.org/10.1093/bioinformatics/btv421 (2015).
CAS  Article  PubMed  PubMed Central  Google Scholar 

22.
Nabil-Fareed, A., Zhemin, Z., Sergeant, M. J. & Achtman, M. A genomic overview of the population structure of Salmonella. PLoS Genet. 14, e1007261, https://doi.org/10.1371/journal.pgen.1007261 (2018).
CAS  Article  Google Scholar 

23.
Nakamura, T., Yamada, K. D., Tomii, K. & Katoh, K. Parallelization of MAFFT for large-scale multiple sequence alignments. Bioinformatics. 34, 2490–2492, https://doi.org/10.1093/bioinformatics/bty121 (2018).
CAS  Article  PubMed  PubMed Central  Google Scholar 

24.
Price, M. N., Dehal, P. S. & Arkin, A. P. FastTree 2–approximately maximum-likelihood trees for large alignments. PLoS One. 5, e9490, https://doi.org/10.1371/journal.pone.0009490 (2010).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

25.
Pang, R. & Wu, Q. A database for risk assessment and comparative genomic analysis of foodborne Vibrio parahaemolyticus in China. figshare https://doi.org/10.6084/m9.figshare.12210287 (2020).

26.
Pang, R. & Wu, Q. Genome assemblies and annotations of food-borne Vibrio parahaemolyticus strains. figshare https://doi.org/10.6084/m9.figshare.12004416 (2020).

27.
Pang, R. et al. Comparative genomic analysis of foodborne Vibrio parahaemolyticus in China. NCBI Sequence Read Archive https://identifiers.org/ncbi/insdc.sra:SRP253458 (2020).

28.
Chen, Y. et al. Foodborne disease outbreaks in 2006 report of the National Foodborne Disease Surveillance Network, China. Wei Sheng Yan Jiu. 39, 331–334 (2010).
PubMed  Google Scholar 

29.
Li, L. et al. Comparative genomic analysis of clinical and environmental strains provides insight into the pathogenicity and evolution of Vibrio parahaemolyticus. BMC Genomics. 15, 1135, https://doi.org/10.1186/1471-2164-15-1135 (2014).
CAS  Article  PubMed  PubMed Central  Google Scholar 

30.
McInerney, J. O., McNally, A. & O’Connell, M. J. Why prokaryotes have pangenomes. Nat. Microbiol. 2, 17040, https://doi.org/10.1038/nmicrobiol.2017.40 (2017).
CAS  Article  PubMed  Google Scholar 

31.
Gonzalez-Escalona, N., Jolley, K. A., Reed, E. & Martinez-Urtaza, J. Defining a core genome multilocus sequence typing scheme for the global epidemiology of Vibrio parahaemolyticus. J. Clin. Microbiol. 55, 1682–1697, https://doi.org/10.1128/JCM.00227-17 (2017).
CAS  Article  PubMed  PubMed Central  Google Scholar  More