Ecology

1.Sikorski, J. A. & Gruys, K. J. Understanding glyphosate’s molecular mode of action with EPSP synthase: evidence favoring an allosteric inhibitor model. Acc. Chem. Res. 30, 2–8 (1997).CAS 
Article 

Google Scholar 
2.Duke, S. O. & Powles, S. B. Glyphosate: a once‐in‐a‐century herbicide. Pest Manag. Sci. 64, 319–325 (2008).CAS 
PubMed 
Article 

Google Scholar 
3.Siehl, D. L. Inhibitors of EPSP synthase, glutamine synthetase and histidine synthesis. In Herbicide Activity: Toxicology, Biochemistry and Molecular Biology, vol. 1 (eds. Michael Roe, R., Burton, J. D. & Kuhr, R. J.) 37 (IOS Press, 1997).4.Shilo, T., Zygier, L., Rubin, B., Wolf, S. & Eizenberg, H. Mechanism of glyphosate control of Phelipanche aegyptiaca. Planta 244, 1095–1107 (2016).CAS 
PubMed 
Article 

Google Scholar 
5.Tzin, V. & Galili, G. New Insights into the shikimate and aromatic amino acids biosynthesis pathways in plants. Mol. Plant 3, 956–972 (2010).CAS 
PubMed 
Article 

Google Scholar 
6.McFall-Ngai, M. et al. Animals in a bacterial world, a new imperative for the life sciences. Proc. Natl. Acad. Sci. USA 110, 3229–3236 (2013).CAS 
PubMed 
Article 

Google Scholar 
7.Hacker, S. D. & Gaines, S. D. Some implications of direct positive interactions for community species diversity. Ecology 78, 1990–2003 (1997).Article 

Google Scholar 
8.van den Bosch, T. J. M. & Welte, C. U. Detoxifying symbionts in agriculturally important pest insects. Microb. Biotechnol. 10, 531–540 (2017).PubMed 
Article 
CAS 

Google Scholar 
9.Lemoine, M. M., Engl, T. & Kaltenpoth, M. Microbial symbionts expanding or constraining abiotic niche space in insects. Curr. Opin. Insect Sci. 39, 14–20 (2020).PubMed 
Article 

Google Scholar 
10.Feldhaar, H. Bacterial symbionts as mediators of ecologically important traits of insect hosts. Ecol. Entomol. 36, 533–543 (2011).Article 

Google Scholar 
11.Moran, N. A. Symbiosis. Curr. Biol. 16, R866–R871 (2006).CAS 
PubMed 
Article 

Google Scholar 
12.Moran, N. A. & Telang, A. Bacteriocyte-associated symbionts of insects. Bioscience 48, 295–304 (1998).Article 

Google Scholar 
13.Oliver, K. M. & Martinez, A. J. How resident microbes modulate ecologically-important traits of insects. Curr. Opin. Insect Sci. 4, 1–7 (2014).PubMed 
Article 

Google Scholar 
14.Douglas, A. E. The microbial dimension in insect nutritional ecology. Funct. Ecol. 23, 38–47 (2009).Article 

Google Scholar 
15.Douglas, A. E. The B vitamin nutrition of insects: the contributions of diet, microbiome and horizontally acquired genes. Curr. Opin. Insect Sci. 23, 65–69 (2017).PubMed 
Article 

Google Scholar 
16.Vigneron, A. et al. Insects recycle endosymbionts when the benefit is over. Curr. Biol. 24, 2267–2273 (2014).CAS 
PubMed 
Article 

Google Scholar 
17.Andersen, S. O. Cuticular sclerotization and tanning. In Insect Molecular Biology and Biochemistry (ed. Gilbert, L. I.) 167–192 (Elsevier, 2012).18.Anbutsu, H. & Fukatsu, T. Symbiosis for insect cuticle formation. In Cellular Dialogues in the Holobiont (eds. Bosch, T. C. G. & Hadfield, M. G.) 201–216 (CRC Press, 2020).19.Anbutsu, H. et al. Small genome symbiont underlies cuticle hardness in beetles. Proc. Natl. Acad. Sci. USA 114, E8382–E8391 (2017).CAS 
PubMed 
Article 

Google Scholar 
20.Li, A. P. & Long, T. J. An evaluation of the genotoxic potential of glyphosate. Fundam. Appl. Toxicol. 10, 537–546 (1988).CAS 
PubMed 
Article 

Google Scholar 
21.Smith, E. A. & Oehme, F. W. The biological activity of glyphosate to plants and animals: a literature review. Vet. Hum. Toxicol. 34, 531–543 (1992).CAS 
PubMed 

Google Scholar 
22.Smith, D. F. Q. et al. Glyphosate inhibits melanization and increases insect susceptibility to infection. bioRxiv (2020).23.Torretta, V., Katsoyiannis, I., Viotti, P. & Rada, E. Critical review of the effects of glyphosate exposure to the environment and humans through the food supply chain. Sustainability 10, 950 (2018).Article 
CAS 

Google Scholar 
24.Snyder, A. K. & Rio, R. V. M. “Wigglesworthia morsitans” folate (Vitamin B 9) biosynthesis contributes to tsetse host fitness. Appl. Environ. Microbiol. 81, 5375–5386 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
25.Motta, E. V. S., Raymann, K. & Moran, N. A. Glyphosate perturbs the gut microbiota of honey bees. Proc. Natl. Acad. Sci. USA 115, 10305–10310 (2018).CAS 
PubMed 
Article 

Google Scholar 
26.Motta, E. V. S. et al. Oral or topical exposure to glyphosate in herbicide formulation impacts the gut microbiota and survival rates of honey bees. Appl. Environ. Microbiol. 86, e01150–20 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
27.Klein, A. et al. A novel intracellular mutualistic bacterium in the invasive ant Cardiocondyla obscurior. ISME J 10, 376–388 (2016).CAS 
PubMed 
Article 

Google Scholar 
28.Wu, D. et al. Metabolic complementarity and genomics of the dual bacterial symbiosis of sharpshooters. PLoS Biol. 4, e188 (2006).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
29.Dunne, J. A. & Williams, R. J. Cascading extinctions and community collapse in model food webs. Philos. Trans. R. Soc. B Biol. Sci 364, 1711–1723 (2009).Article 

Google Scholar 
30.Dunne, J. A., Williams, R. J. & Martinez, N. D. Network structure and biodiversity loss in food webs: robustness increases with connectance. Ecol. Lett. 5, 558–567 (2002).Article 

Google Scholar 
31.Memmott, J. et al. Biodiversity loss and ecological network structure. In Ecological Networks: Linking Structure to Dynamics in Food Webs (eds Pascual, M. & Dunne, J. A.) 325–347 (Oxford University Press, 2005).32.Liao, C., Upadhyay, A., Liang, J., Han, Q. & Li, J. 3,4-Dihydroxyphenylacetaldehyde synthase and cuticle formation in insects. Dev. Comp. Immunol. 83, 44–50 (2018).CAS 
PubMed 
Article 

Google Scholar 
33.Muthukrishnan, S., Merzendorfer, H., Arakane, Y. & Kramer, K. J. Chitin metabolism in insects. In Insect Molecular Biology and Biochemistry (ed. Gilbert, L. I.) 193–235 (Elsevier, 2012).34.Wirtz, R. A. & Hopkins, T. L. Tyrosine and phenylalanine concentrations in haemolymph and tissues of the American cockroach, Periplaneta americana, during metamorphosis. J. Insect Physiol. 20, 1143–1154 (1974).CAS 
PubMed 
Article 

Google Scholar 
35.Gibbs, A. G. & Rajpurohit, S. Cuticular lipids and water balance. In Insect Hydrocarbons (eds Blomquist, G. J. & Bagneres, A. -G.) 100–120 (Cambridge University Press, 2010).36.Hackman, R. H. Chemistry of the insect cuticle. in The Physiology of Insecta (ed. Rodstein, M.) 215–270 (Academic Press, 1974).37.Mattson, W. J. Herbivory in relation to plant nitrogen content. Annu. Rev. Ecol. Syst. 11, 119–161 (1980).Article 

Google Scholar 
38.Kumar, V. et al. Amino acids distribution in economical important plants: a review. Biotechnol. Res. Innov 3, 197–207 (2019).Article 

Google Scholar 
39.Noh, M. Y., Muthukrishnan, S., Kramer, K. J. & Arakane, Y. Cuticle formation and pigmentation in beetles. Curr. Opin. Insect Sci. 17, 1–9 (2016).PubMed 
Article 

Google Scholar 
40.Sterkel, M. et al. Tyrosine detoxification is an essential trait in the life history of blood-feeding arthropods. Curr. Biol. 26, 2188–2193 (2016).CAS 
PubMed 
Article 

Google Scholar 
41.Herrmann, K. M. & Weaver, L. M. The shikimate pathway. Annu. Rev. Plant Biol. 50, 473–503 (1999).CAS 
Article 

Google Scholar 
42.Engl, T. et al. Ancient symbiosis confers desiccation resistance to stored grain pest beetles. Mol. Ecol. 27, 2095–2108 (2018).CAS 
PubMed 
Article 

Google Scholar 
43.Hirota, B. et al. A novel, extremely elongated, and endocellular bacterial symbiont supports cuticle formation of a grain pest beetle. MBio 8, 1–16 (2017).Article 

Google Scholar 
44.Boyer, S., Zhang, H. & Lempérière, G. A review of control methods and resistance mechanisms in stored-product insects. Bull. Entomol. Res. 102, 213 (2012).CAS 
PubMed 
Article 

Google Scholar 
45.Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P. & Tyson, G. W. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 25, 1043–1055 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
46.Moran, N. A., McCutcheon, J. P. & Nakabachi, A. Genomics and evolution of heritable bacterial symbionts. Annu. Rev. Genet. 42, 165–190 (2008).CAS 
PubMed 
Article 

Google Scholar 
47.McCutcheon, J. P. & Moran, N. A. Extreme genome reduction in symbiotic bacteria. Nat. Rev. Microbiol. 10, 13–26 (2012).CAS 
Article 

Google Scholar 
48.Van Leuven, J. T., Meister, R. C., Simon, C. & McCutcheon, J. P. Sympatric speciation in a bacterial endosymbiont results in two genomes with the functionality of one. Cell 158, 1270–1280 (2014).PubMed 
Article 
CAS 

Google Scholar 
49.Campbell, M. A., Łukasik, P., Simon, C. & McCutcheon, J. P. Idiosyncratic genome degradation in a bacterial endosymbiont of periodical cicadas. Curr. Biol. 27, 3568–3575.e3 (2017).CAS 
PubMed 
Article 

Google Scholar 
50.Campbell, M. A. et al. Genome expansion via lineage splitting and genome reduction in the cicada endosymbiont Hodgkinia. Proc. Natl. Acad. Sci. USA 112, 10192–10199 (2015).CAS 
PubMed 
Article 

Google Scholar 
51.Chen, Y. C., Liu, T., Yu, C. H., Chiang, T. Y. & Hwang, C. C. Effects of GC bias in next-generation-sequencing data on de novo genome assembly. PLoS ONE 8, e62856 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
52.Kozarewa, I. et al. Amplification-free Illumina sequencing-library preparation facilitates improved mapping and assembly of (G+ C)-biased genomes. Nat. Methods 6, 291–295 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Quail, M. A. et al. A tale of three next generation sequencing platforms: comparison of Ion Torrent, Pacific Biosciences and Illumina MiSeq sequencers. BMC Genomics 13, 1–13 (2012).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
54.Treangen, T. J. & Salzberg, S. L. Repetitive DNA and next-generation sequencing: computational challenges and solutions. Nat. Rev. Genet. 13, 36–46 (2012).CAS 
Article 

Google Scholar 
55.Sloan, D. B. et al. Parallel histories of horizontal gene transfer facilitated extreme reduction of endosymbiont genomes in sap-feeding insects. Mol. Biol. Evol. 31, 857–871 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
56.Zucko, J. et al. Global genome analysis of the shikimic acid pathway reveals greater gene loss in host-associated than in free-living bacteria. BMC Genomics 11, 628 (2010).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
57.Tokuda, G. et al. Maintenance of essential amino acid synthesis pathways in the Blattabacterium cuenoti symbiont of a wood-feeding cockroach. Biol. Lett. 9, 20121153 (2013).PubMed 
PubMed Central 
Article 

Google Scholar 
58.Kinjo, Y. et al. Parallel and gradual genome erosion in the Blattabacterium endosymbionts of Mastotermes darwiniensis and Cryptocercus Wood Roaches. Genome Biol. Evol. 10, 1622–1630 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
59.Menzel, R. & Roth, J. Purification of the putA gene product. A bifunctional membrane-bound protein from Salmonella typhimurium responsible for the two-step oxidation of proline to glutamate. J. Biol. Chem. 256, 9755–9761 (1981).CAS 
PubMed 
Article 

Google Scholar 
60.Zhou, Y., Zhu, W., Bellur, P. S., Rewinkel, D. & Becker, D. F. Direct linking of metabolism and gene expression in the proline utilization a protein from Escherichia coli. Amino Acids 35, 711–718 (2008).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
61.Sabree, Z. L., Kambhampati, S. & Moran, N. A. Nitrogen recycling and nutritional provisioning by Blattabacterium, the cockroach endosymbiont. Proc. Natl. Acad. Sci. USA 106, 19521–19526 (2009).CAS 
PubMed 
Article 

Google Scholar 
62.McCutcheon, J. P., McDonald, B. R. & Moran, N. A. Convergent evolution of metabolic roles in bacterial co-symbionts of insects. Proc. Natl. Acad. Sci. USA 106, 15394–15399 (2009).CAS 
PubMed 
Article 

Google Scholar 
63.Sabree, Z. L., Huang, C. Y., Okusu, A., Moran, N. A. & Normark, B. B. The nutrient supplying capabilities of Uzinura, an endosymbiont of armoured scale insects. Environ. Microbiol. 15, 1988–1999 (2013).CAS 
PubMed 
Article 

Google Scholar 
64.Rosas-Pérez, T., Rosenblueth, M., Rincón-Rosales, R., Mora, J. & Martínez-Romero, E. Genome Sequence of “Candidatus Walczuchella monophlebidarum” the Flavobacterial endosymbiont of Llaveia axin axin (Hemiptera: Coccoidea: Monophlebidae). Genome Biol. Evol. 6, 714–726 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
65.Kuriwada, T. et al. Biological role of Nardonella endosymbiont in its weevil host. PLoS ONE 5, e13101 (2010).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
66.Okude, G. et al. Novel bacteriocyte-associated pleomorphic symbiont of the grain pest beetle Rhyzopertha dominica (Coleoptera: Bostrichidae). Zool. Lett. 3, 13 (2017).Article 

Google Scholar 
67.Hirota, B., Meng, X.-Y. & Fukatsu, T. Bacteriome-sssociated rndosymbiotic bacteria of Nosodendron tree sap beetles (Coleoptera: Nosodendridae). Front. Microbiol. 11, 2556 (2020).Article 

Google Scholar 
68.Hopkins, T. L. & Kramer, K. J. Insect cuticle sclerotization. Annu. Rev. Entomol. 37, 273–302 (1992).CAS 
Article 

Google Scholar 
69.Andersen, S. O. Insect cuticular sclerotization: a review. Insect Biochem. Mol. Biol. 40, 166–178 (2010).CAS 
PubMed 
Article 

Google Scholar 
70.Cao, G. et al. A novel 5-enolpyruvylshikimate-3-phosphate synthase shows high glyphosate tolerance in Escherichia coli and tobacco plants. PLoS ONE 7, e38718 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
71.Moran, N. A. & Bennett, G. M. The tiniest tiny genomes. Annu. Rev. Microbiol. 68, 195–215 (2014).CAS 
PubMed 
Article 

Google Scholar 
72.McCutcheon, J. P., Boyd, B. M. & Dale, C. The life of an insect endosymbiont from the cradle to the grave. Curr. Biol. 29, R485–R495 (2019).CAS 
PubMed 
Article 

Google Scholar 
73.Salem, H. et al. Drastic genome reduction in an herbivore’s pectinolytic symbiont. Cell 171, 1520–1531 (2017).CAS 
PubMed 
Article 

Google Scholar 
74.Reis, F. et al. Bacterial symbionts support larval sap feeding and adult folivory in (semi-) aquatic reed beetles. Nat. Commun. 11, 1–15 (2020).
Google Scholar 
75.Salem, H., Florez, L., Gerardo, N. & Kaltenpoth, M. An out-of-body experience: the extracellular dimension for the transmission of mutualistic bacteria in insects. Proc. R. Soc. B Biol. Sci. 282, 20142957 (2015).Article 

Google Scholar 
76.Salem, H. et al. Symbiont digestive range reflects host plant breadth in herbivorous beetles. Curr. Biol. 30, 2875–2886 (2020).CAS 
PubMed 
Article 

Google Scholar 
77.Hansen, A. K., Pers, D. & Russell, J. A. Symbiotic solutions to nitrogen limitation and amino acid imbalance in insect diets. In Mechanisms Underlying Microbial Symbiosis, vol. 58 (ed. Kerry M. Oliver, J. A. R.) 161–205 (Academic Press, 2020).78.Tanner, J. J. Structural biology of proline catabolism. Amino Acids 35, 719–730 (2008).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
79.Adams, E. & Frank, L. Metabolism of proline and the hydroxyprolines. Annu. Rev. Biochem. 49, 1005–61 (1980).CAS 
PubMed 
Article 

Google Scholar 
80.Bursell, E. The role of proline in energy metabolism.In Energy Metabolism in Insects (ed. Downer R.G.H.) 135–154 (Springer, Boston, 1981).81.Engl, T., Schmidt, T. H. P., Kanyile, S. N. & Klebsch, D. Metabolic cost of a nutritional symbiont manifests in delayed reproduction in a grain pest beetle. Insects 11, 717 (2020).PubMed Central 
Article 
PubMed 

Google Scholar 
82.José de Souza, D., Devers, S. & Lenoir, A. Blochmannia endosymbionts and their host, the ant Camponotus fellah: cuticular hydrocarbons and melanization. C. R. Biol. 334, 737–741 (2011).PubMed 
Article 
CAS 

Google Scholar 
83.Zientz, E., Beyaert, I., Gross, R. & Feldhaar, H. Relevance of the endosymbiosis of Blochmannia floridanus and carpenter ants at different stages of the life cycle of the host. Appl. Environ. Microbiol. 72, 6027–6033 (2006).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
84.Oakeson, K. F. et al. Genome degeneration and adaptation in a nascent stage of symbiosis. Genome Biol. Evol. 6, 76–93 (2013).Article 

Google Scholar 
85.Chong, R. A. & Moran, N. A. Evolutionary loss and replacement of Buchnera, the obligate endosymbiont of aphids. ISME J 12, 898–908 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
86.McCutcheon, J. P. & Moran, N. A. Functional convergence in reduced genomes of bacterial symbionts spanning 200 My of evolution. Genome Biol. Evol. 2, 708–718 (2010).PubMed 
PubMed Central 
Article 

Google Scholar 
87.Gerth, M., Gansauge, M. T., Weigert, A. & Bleidorn, C. Phylogenomic analyses uncover origin and spread of the Wolbachia pandemic. Nat. Commun. 5, 1–7 (2014).Article 
CAS 

Google Scholar 
88.Santos-Garcia, D., Silva, F. J., Morin, S., Dettner, K. & Kuechler, S. M. The all-rounder Sodalis: a new bacteriome-associated endosymbiont of the Lygaeoid bug Henestaris halophilus (Heteroptera: Henestarinae) and a critical examination of its evolution. Genome Biol. Evol. 9, 2893–2910 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
89.Motta, E. V. S. & Moran, N. A. Impact of glyphosate on the honey bee gut microbiota: effects of intensity, duration, and timing of exposure. Msystems 5, e00268–20 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
90.Helander, M., Pauna, A., Saikkonen, K. & Saloniemi, I. Glyphosate residues in soil affect crop plant germination and growth. Sci. Rep. 9, 19653 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
91.Kiers, E. T., Rousseau, R. A., West, S. A. & Denlson, R. F. Host sanctions and the legume-rhizobium mutualism. Nature 425, 78–81 (2003).CAS 
PubMed 
Article 

Google Scholar 
92.Whiteside, M. D., Digman, M. A., Gratton, E. & Treseder, K. K. Organic nitrogen uptake by arbuscular mycorrhizal fungi in a boreal forest. Soil Biol. Biochem. 55, 7–13 (2012).CAS 
Article 

Google Scholar 
93.Faita, M. R., Cardozo, M. M., Amandio, D. T. T., Orth, A. I. & Nodari, R. O. Glyphosate-based herbicides and Nosema sp. microsporidia reduce honey bee (Apis mellifera L.) survivability under laboratory conditions. J. Apic. Res. 59, 1–11 (2020).Article 

Google Scholar 
94.Wilson, A. C. C. et al. Genomic insight into the amino acid relations of the pea aphid, Acyrthosiphon pisum, with its symbiotic bacterium Buchnera aphidicola. Insect Mol. Biol. 19, 249–258 (2010).CAS 
PubMed 
Article 

Google Scholar 
95.Sánchez-Bayo, F. & Wyckhuys, K. A. G. Worldwide decline of the entomofauna: a review of its drivers. Biol. Conserv. 232, 8–27 (2019).Article 

Google Scholar 
96.Wagner, D. L. Insect declines in the anthropocene. Annu. Rev. Entomol. 65, 457–480 (2020).CAS 
Article 

Google Scholar 
97.Desneux, N., Decourtye, A. & Delpuech, J.-M. The sublethal effects of pesticides on beneficial arthropods. Annu. Rev. Entomol. 52, 81–106 (2007).CAS 
Article 

Google Scholar 
98.Goulson, D. The insect apocalypse, and why it matters. Curr. Biol. 29, R967–R971 (2019).CAS 
PubMed 
Article 

Google Scholar 
99.Hallmann, C. A. et al. More than 75 percent decline over 27 years in total flying insect biomass in protected areas. PLoS ONE 12, e0185809 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
100.Hayes, T. B. & Hansen, M. From silent spring to silent night: agrochemicals and the anthropocene. Elem. Sci. Anthropol. 5, (2017).101.Bowler, D. E., Heldbjerg, H., Fox, A. D., Jong, M. & Böhning‐Gaese, K. Long‐term declines of European insectivorous bird populations and potential causes. Conserv. Biol. 33, 1120–1130 (2019).PubMed 
Article 

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

Google Scholar 
103.Wood, D. E. & Salzberg, S. L. Kraken: ultrafast metagenomic sequence classification using exact alignments. Genome Biol. 15, 1–12 (2014).Article 

Google Scholar 
104.Wood, D. E., Lu, J. & Langmead, B. Improved metagenomic analysis with Kraken 2. Genome Biol. 20, 257 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
105.Bankevich, A. et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J. Comput. Biol. 19, 455–477 (2015).Article 
CAS 

Google Scholar 
106.Laczny, C. C. et al. BusyBee Web: metagenomic data analysis by bootstrapped supervised binning and annotation. Nucleic Acids Res. 45, W171–W179 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
107.Aziz, R. K. et al. The RAST Server: rapid annotations using subsystems technology. BMC Genomics 9, 1–15 (2008).Article 
CAS 

Google Scholar 
108.Arkin, A. P. et al. KBase: The United States department of energy systems biology knowledgebase. Nat. Biotechnol. 36, 566 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
109.Krzywinski, M. et al. Circos: an information aesthetic for comparative genomics. Genome Res. 19, 1639–1645 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
110.Price, M. N., Dehal, P. S. & Arkin, A. P. FastTree 2–approximately maximum-likelihood trees for large alignments. PLoS ONE 5, e9490 (2010).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
111.Tatusov, R. L., Galperin, M. Y., Natale, D. A. & Koonin, E. V. The COG database: a tool for genome-scale analysis of protein functions and evolution. Nucleic Acids Res. 28, 33–36 (2000).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
112.Edgar, R. C. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797 (2004).CAS 
PubMed 
PubMed Central 

Google Scholar 
113.Guindon, S. et al. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst. Biol. 59, 307–321 (2010).CAS 
PubMed 
Article 

Google Scholar 
114.Brettin, T. et al. RASTtk: a modular and extensible implementation of the RAST algorithm for building custom annotation pipelines and annotating batches of genomes. Sci. Rep. 5, 8365 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
115.Li, L., Stoeckert, C. J. & Roos, D. S. OrthoMCL: identification of ortholog groups for eukaryotic genomes. Genome Res. 13, 2178–2189 (2003).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
116.Kanehisa, M., Sato, Y. & Morishima, K. BlastKOALA and GhostKOALA: KEGG tools for functional characterization of genome and metagenome sequences. J. Mol. Biol. 428, 726–731 (2016).CAS 
Article 
PubMed 

Google Scholar 
117.Overbeek, R. et al. The SEED and the Rapid Annotation of microbial genomes using Subsystems Technology (RAST). Nucleic Acids Res. 42, D206–14 (2014).CAS 
PubMed 
Article 

Google Scholar 
118.Weiss, B. & Kaltenpoth, M. Bacteriome-localized intracellular symbionts in pollen-feeding beetles of the genus Dasytes (Coleoptera, Dasytidae). Front. Microbiol. 7, 1486 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
119.Dunn, O. J. Multiple comparisons using rank sums. Technometrics 6, 241–252 (1964).Article 

Google Scholar 
120.Tanahashi, M. Natsumushi: Image measuring software for entomological studies. Entomol. Sci. 21, 347–360 (2018).Article 

Google Scholar 
121.Pérez-Palacios, T., Barroso, M. A., Ruiz, J. & Antequera, T. A rapid and accurate extraction procedure for analysing free amino acids in meat samples by GC–MS. Int. J. Anal. Chem. 2015, 209214 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
122.Miller, R. G. Simultaneous Statistical Inference (Springer, 1981).123.Engl, T., Kiefer, J.S.T. Data from: Inhibition of a nutritional endosymbiont by glyphosate abolishes mutualistic benefit on cuticle synthesis in Oryzaephilus surinamenis. Max Planck Soc. https://doi.org/10.17617/3.5l (2021). More

Isolation and characterisation of viruses infecting the picoeukaryote Bathycoccus Clade BIIBathycoccus BII isolates RCC716 and RCC715 used in our experiments were originally cultured from a nutrient-limited region in the Indian Ocean. Clade BII as a whole has been reported extensively in warm oligotrophic ocean gyres based on metagenome analyses [22,23,24]. Peak abundances occurr when well-developed deep chlorophyll maxima are present, or throughout the photic zone during mixing periods at Station ALOHA of the Hawaii Ocean Time-series [12]. We targeted BATS for viral isolation in springtime because Bathycoccus has been observed at relatively high abundance during this period using qPCR [74]. Here, three viruses were isolated against RCC716 [12] using seawater flown from BATS/Bermuda to the laboratory, obviating bringing this finicky strain into the field for use as a viral host. We then purified the viruses by serial dilutions and sequenced the partial PolB gene to determine whether they were evolutionarily different from other cultured viruses. BLASTn and preliminary phylogenetic analysis using GenBank nr sequences indicated they were distinct from described viruses with deposited sequences, with best BLASTn hits to Bathycoccus prasinos viruses (62–74% nucleotide identity). Transmission electron microscopy (TEM) revealed that all three viruses have similar morphology to other characterised prasinoviruses [75], with icosahedral capsids diameter ranging between 120 and 140 nm (Fig. 1A).Fig. 1: Morphology and evolutionary relationships of newly discovered Bathycoccus viruses.A Transmission electron micrographs of BII-V1, BII-V2 and BII-V3 (scale bar, 50 nm). The capsid diameters (n = 6 virions) measured 138 ± 2 nm (BII-V1), 150 ± 5 nm (BII-V2) and 152 ± 11 nm (BII-V3). B Maximum Likelihood (ML) phylogenetic reconstruction of green algal viruses inferred from a concatenated alignment of 22 core proteins shared among the viruses (7,001 positions) under the LG + G + F model. Node support was calculated from 1000 bootstrap (BS) replicates, with all branches acquired support values of 100% (white dots). Viruses infecting Chlorella were used as an outgroup and the branch connecting the prasinoviruses to the outgroup was truncated for display purpose. The new Bathycoccus viruses isolated against Bathycoccus Clade II (sensu [12]) isolate RCC716 (named as species Bathycoccus calidus herein, see below) are in bold. Colours reflect different host species within each genus. Letters alongside vertical lines (a and b) correspond to Bathycoccus viral clades. C Venn diagram of the shared and unique protein-encoding genes in the genome sequences of the new Bathycoccus viruses.Full size imageGenomic sequencing and multi-gene evolutionary analysesAssembly of DNA sequences from the viral isolates after deep sequencing by Illumina rendered one complete dsDNA genome sequence (BII-V3), and two others may still be partial (Table 1). The BII-V2 genome, which was in one contig, was similar in size (~208 kb) to that of BII-V3 (~212 kb). The BII-V1 genome assembly was ~174 kb and comprised of four linear dsDNA scaffolds. The viral concentrate was deeply sequenced ( >50x coverage) and minor fragmentation of the genome was partially related to repeats that were not resolved during assembly. The total number of putative open reading frames (ORFs) varied from 220 in BII-V1 to 235 in BII-V2 (Table 1). Gene synteny was globally well-conserved across the BII-Vs and the BpV1 and BpV2 viruses of B. prasinos (Fig. S1), with limited genomic rearrangements. Other genome characteristics such as the coding proportion (~90%) and G + C % (~36%) were similar to other described prasinoviruses infecting Mamiellophyceae [64, 75], for which the reported number of proteins range from 203 to 268 and G + C % from 37 to 45%. However, the full-length PolB gene from the genome assemblies differed for BII-V3 from the other two, in having a 329 amino acid intein (Supplementary information table S3). Likewise, inteins have been reported at the same PolB position in uncultivated prasinoviruses from the subtropical Pacific Ocean [76], where Bathycoccus BII is abundant [12].Table 1 Genomic characteristics of the three Bathycoccus viruses (BII-Vs) isolated against Clade BII isolate RCC716.Full size tableTo reconstruct a robust phylogeny for the new viruses, we employed 22 proteins previously identified as being shared across all available green algal virus genomes, including both prasinoviruses and chloroviruses [65]. We found all 22 in the predicted coding sequences of BII-V1; however, DNA helicase (SNF2) was not found in BII-V2 or -V3, FAD-dependent thymidylate synthase (thy1) and the topoisomerase IV were not found in BII-V2, nor was the prolyl 4-hydroxylase in the BII-V3 genome. Additional searches with tBLASTn did not recover these genes or fragments of them, suggesting they have been lost. Phylogenomic reconstruction grouped the three BII-Vs with the two BpVs [32], in a fully supported clade that branched adjacent to a large group of viruses that infect various species of Ostreococcus and Micromonas (Fig. 1B). The clade of Bathycoccus viruses was segregated in two subclades with BII-V2 and BII-V3 clustering together adjacent to BII-V1 and BpVs (Fig. 1B). While better resolution of the position of BII-V1 awaits greater taxonomic sampling, our results demonstrated that the three new viruses branch adjacent or basally to BpVs.Variation in prasinovirus gene content and functions encodedThe three Bathycoccus Clade BII viruses had 72–77% of their proteins held in common, and ~30 unique proteins as well as a few proteins shared by just two of the three viruses (Fig. 1C). The 170 shared proteins had higher amino acid identities between BII-V2 and BII-V3 (73% aa identity) than to BII-V1 (69% and 68%, respectively). Generally, only 19–21% of Bathycoccus viral genes could be assigned a functional category, based on EggNOG classification. Similar functional category distributions were observed across both prasinoviruses and chloroviruses, including lipid metabolism, RNA processing and modification, and nucleotide metabolism and transport (Fig. 2A). Other functional categories were more variable, such as cell wall/membrane/envelope biogenesis genes prevalent in chloroviruses (potentially related to their enveloped nature), as well as genes involved in modification of the capsid with compounds such as with chitin and hyaluronan [77, 78] that are absent from prasinoviruses sequenced to date (Fig. 2A). Within prasinoviruses, most of the unique proteins in the Bathycoccus viruses lack defined functional categories. Among those with functional assignments, all five Bathycoccus viruses encoded a P2X receptor in the intracellular trafficking and secretion category, and both BII-V2 and -V3 encode two proteins putatively involved in degrading the aromatic compound 4-hydroxy-2-oxopentanoate to acetyl-CoA (secondary metabolite category), that otherwise are only encoded by one other prasinovirus, MpV1 [32]. Similar to the phylogenetic relationships, the functional category distributions of BII-V1 were closer to those of BpVs than to BII-Vs. The primary difference was in carbohydrate metabolism, where BII-V2 and -V3 each encodes ribulose-phosphate 3-epimerase (involved in the pentose phosphate pathway and carbon fixation; not found in any other available virus genomes, but encoded by B. prasinos) and TDP-glucose 4,6-dehydratase (involved in biosynthesis of rhamnose and encoded by most other chloroviruses and prasinoviruses [79]). Notably, the putative high-affinity phosphate transporter (PHO4, also termed HAPT) was only present in BII-V1 and BpV1, as well as OtV2 (isolated against the Ostreococcus Clade OII ecotype), and most sequenced viruses of O. lucimarinus (Supplementary information table S3). This gene is hypothesised to enhance phosphate uptake during infection under phosphorus‐limited host growth [25], as observed for the PstS phosphate transport system expressed by cyanophages [80], mitigating limitation of this key component of viral genomes. However, most isolated prasinovirus genomes come from waters that are not considered phosphate-limited, hence presence of this gene may connect to poising the host for responding to sudden availability of other nutrients, such as nitrogen, which is often limiting in the ecosystems from which these viruses were isolated. Studies of virus-cell responses under various limiting nutrients are required to understand the retention of this host-derived HGT.Fig. 2: Distribution of functions and orthologous protein families across genome-sequenced prasinoviruses.A Functional category distributions across 21 genome-sequenced prasinoviruses and chloroviruses based on EggNOG categorisation. Viruses are clustered by similarity in their distribution of the functional categories on the y-axis and the frequency of each category across the viral genomes determines clustering along the x-axis ordering. Genes with homology to proteins in the EggNOG database but could not be assigned a function are in the “function unknown” category. B Orthogroups presence/absence patterns ordered along the x-axis by ranking according to the total number of genes in the orthogroup. For inclusion, the orthogroup was required to include protein sequences from at least two different viral genomes. Viruses are ordered along the vertical by their presence/absence pattern reconstructed by hierarchical clustering (topology on the left). Top histogram: frequency of each orthogroup in sequenced prasinoviruses. C Genes in each virus (number) not assigned to any orthogroup, with viruses in the same vertical order as B.Full size imageHierarchical clustering of orthologous proteins revealed patterns across prasinoviruses that generally corresponded with phylogenetic relationships. The BII- and Bp-viruses shared 130 orthologous proteins and hierarchical clustering (Fig. 2B) followed the clade structure of the phylogenomic reconstruction (Fig. 1B) with the exception of BII-V1 that grouped with BII-Vs, as well as OtV6, which grouped with Micromonas viruses. These orthologous proteins had on average 72% amino acid identity between BII-V2 and BII-V3, and 88% between the two B. prasinos viruses, but between 65 to 67% when comparing members of these two groups (Table 2). BII-V1 orthologs also had 67% and 66% amino acid identity to BII-V2 and BII-V3, respectively, while they had 83% and 80% identity to BpV1 and BpV2, respectively. Collectively, these results indicate that BII-V2 and -V3 diverged from BpVs prior to the divergence of BII-V1.Table 2 Average percent amino acid identity of the orthologous proteins between the five Bathycoccus viruses.Full size tableOf the 130 orthologous Bathycoccus virus proteins, 37% were assigned putative functions revealing core components of this viral group (Supplementary information table S3). These included genes involved in DNA replication and transcription, including PolB (type II), a DNA topoisomerase, a transcription factor S-II, mRNA capping enzymes, ribonucleases, a ribonucleotide reductase, and a dUTPase. Several others are necessary for viral particle synthesis, such as genes encoding structural elements for assembling the virion, including capsid proteins (5–6 copies per genome), as well as transcriptional regulators connected to the replication cycle. The BII viruses showed a number of differences among orthologous protein families. In addition to each having “unique” protein sets, there was a set of BII-V specific orthogroups, as well as some shared with BpVs, and/or other prasinoviruses (Fig. 1C and Supplementary information table S3). First, six predicted proteins showed orthologs across the three BII-Vs, but were not present in other prasinoviruses sequenced to date. Only one of these six was assigned putative function, belonging to the XRE family of transcriptional regulators. Additionally, all BII viruses harboured a tandem duplication of the FstH gene, while other sequenced prasinoviruses (including the two Clade BI viruses) have one copy (Supplementary information table S3). This ATP-dependent metalloprotease has been shown to be involved in photosystem II repair in cyanobacteria [81], and is present in genomes of photosynthetic eukaryotes, including all Mamiellophyceae [15, 16]. In Arabidopsis and Chlamydomonas it has been shown to be involved in protein quality control in the thylakoid membranes [82]. A gene of unknown function was also duplicated in the BII-virus genomes, that is a single copy in BpVs and absent from other sequenced prasinoviruses. Genes putatively encoding a glucose-1-phosphate adenylyltransferase, a glycosyltransferase and a thiamine pyrophosphate-requiring enzyme involved in amino acid biosynthesis were sporadically found in BII-viruses.Considering the two Bathycoccus virus subclades (Fig. 1B), there is one predicted protein of unknown function exclusive to BpV1, BpV2 and BII-V1 and six predicted proteins shared only by BII-V2 and BII-V3. Among the latter, one belonged to the Ribulose-5-Phosphate-3-Epimerase (RPE) family, which catalyses the interconversion of D-ribulose 5-phosphate (Ru5P) into d-xylulose 5-phosphate, as part of the Calvin cycle (although no transit peptide was detected using TargetP) and in the oxidative pentose phosphate pathway. The ortholog analyses further showed that among prasinoviruses, 9, 17 and 18 genes were unique to BII-V1, BII-V2 and BII-V3, respectively (Fig. 2B). Apart from one nucleotidyltransferase and one glycosyltransferase (group 1) in BII-V1, none of these unique genes had known functions.To study the evolutionary aspects of the shared prasinovirus proteins, we constructed and examined 130 phylogenies of orthogroups shared between Bathycoccus viruses. Nine showed a topology where all three BII-Vs grouped together with full support (100% bootstrap support), separate from the BpV orthologs, and in contrast to the multi-gene phylogeny where BII-V1 grouped with BpVs (Fig. 1B). The average amino acid similarities within each of these nine protein ortholog groups ranged from 85 to 88% between BII-Vs proteins, while they were 77 to 81% between BII-Vs and BpVs, different from overall amino acid similarity averages (Table 2). Interestingly, proteins from three of these nine ortholog groups, all lacking known functions, were adjacent to each other in the genome, or separated by only one gene. This synteny and co-location likely reflects the acquisition of these genes before co-infecting viruses diverged via viral recombination [83].Infection dynamics of Bathycoccus virusesGeneral host specificity of BII-viruses was assessed using two B. prasinos isolates (CCMP1898 and RCC4222, Clade BI), the two available Clade BII isolates (RCC715 and RCC716), four Ostreococcus species and one Micromonas species (Table 3). None were able to infect the B. prasinos, Ostreococcus or Micromonas isolates tested, suggesting BII-V specificity for Bathycoccus Clade BII. Similar host specificity has been observed in O. lucimarinus viruses, none of which infect O. tauri [64], and other viruses of eukaryotic and prokaryotic algae [84, 85]. Some other prasinoviruses appear to have broader host ranges [85,86,87], or their host species are less divergent than the two known Bathycoccus clades. For example, generalist viruses isolated against Micromonas commoda can infect M. bravo [85]. Further investigations are necessary to determine the extent to which the six shared proteins in BII-Vs (absent from BpVs), are responsible for the differences in host and virus specificity of interactions, versus variations in the shared Bathycoccus virus proteins (65–83% similarity). Importantly, host specificity tests for the new viruses described herein were limited by weak sampling of Bathycoccus diversity (in culture; all that we could acquire were tested).Table 3 Results of cross infectivity tests of BII-V1, BII-V2 and BII-V3 against isolates representing various picoprasinophyte species within the Class Mamiellophyceae.Full size tableAlthough specific for the BII clade, the three BII-Vs exhibited variations in infectivity of the two cultured BII strains available, despite their isolation from the same sample and having identical ITS1 and ITS2 sequences. BII-V1 lysed and cleared RCC715 and RCC716 cultures after four days (Table 3). The same was true for BII-V2 and BII-V3, when incubated with RCC716. Different from results for BII-V1, we found that while BII-V2 and -V3 initially lysed RCC715 cultures, resistant populations became evident at day 7 of infectivity tests, and measureable lysis of RCC715 could not be achieved thereafter. These results underscored the need to further examine host-virus interactions for the three new viruses.Infection dynamics over time course experiments further illuminated differences in BII-V impacts on hosts. In these experiments, growth rates of the uninfected (control) RCC715 and RCC716 cultures were 0.45 ± 0.04 day−1 and 0.49 ± 0.06 per day, respectively, similar to rates during the pre-experiment acclimation period (T-test, p  > 0.05). Host and virus dynamics were similar for RCC715 and RCC716 infected with BII-V1 (Fig. S2 and Fig. 3), with cell numbers starting to diverge from control abundances 10 h after inoculation (T-test, p  More

1.Francis, R. C., Hare, S. R., Hollowed, A. B. & Wooster, W. S. Effects of interdecadal climate variability on the oceanic ecosystems of the NE Pacific. Fish. Oceanogr. 7, 1–21. https://doi.org/10.1046/j.1365-2419.1998.00052.x (1998).Article 

Google Scholar 
2.Hollowed, A. B. & Wooster, W. S. Variability of winter ocean conditions and strong year classes of Northeast Pacific groundfish. ICES Mar. Sci. Symp. 195, 433–444 (1992).
Google Scholar 
3.Mantua, N. J. & Hare, S. R. The Pacific decadal oscillation. J. Oceanogr. 58, 35–44. https://doi.org/10.1023/A:1015820616384 (2002).Article 

Google Scholar 
4.Di Lorenzo, E. et al. North Pacific gyre oscillation links ocean climate and ecosystem change. Geophys. Res. Lett.https://doi.org/10.1029/2007GL032838 (2008).Article 

Google Scholar 
5.Di Lorenzo, E. et al. Synthesis of pacific ocean climate and ecosystem dynamics. Oceanography 26, 68–81. https://doi.org/10.5670/oceanog.2013.76 (2013).Article 

Google Scholar 
6.Anderson, P. J. & Piatt, J. F. Community reorganization in the Gulf of Alaska following ocean climate regime shift. Mar. Ecol. Prog. Ser. 189, 117–123 (1999).ADS 
Article 

Google Scholar 
7.Polovina, J. J., Mitchum, G. T. & Evans, G. T. Decadal and basin-scale variation in mixed layer depth and the impact on biological production in the Central and North Pacific, 1960–88. Deep Sea Res. Part I Oceanogr. Res. Pap. 42, 1701–1716 (1995).8.Litzow, M. A. & Mueter, F. J. Assessing the ecological importance of climate regime shifts: an approach from the North Pacific Ocean. Prog. Oceanogr. 120, 110–119. https://doi.org/10.1016/j.pocean.2013.08.003 (2014).ADS 
Article 

Google Scholar 
9.Möllmann, C., Folke, C., Edwards, M. & Conversi, A. Marine regime shifts around the globe: theory, drivers and impacts. Philos. Trans. R. Soc. B Biol. Sci. 370, 20130260. https://doi.org/10.1098/rstb.2013.0260 (2015).10.Litzow, M. A., Mueter, F. J. & Hobday, A. J. Reassessing regime shifts in the North Pacific: incremental climate change and commercial fishing are necessary for explaining decadal-scale biological variability. Glob. Change Biol. 20, 38–50. https://doi.org/10.1111/gcb.12373 (2014).ADS 
Article 

Google Scholar 
11.Goen, J., & Erikson, L. Fishery Statistics. Technical Report. IPHC-2018-AM094-05, International Pacific Halibut Commission (2017).12.Johnson, K. F. et al. Status of the U.S. Sablefish Resource in 2015. Technical Report. Pacific Fishery Management Council (2016).13.Pacific Fishery Management Council. Pacific Coast Groundfish Fishery Management Plan. Technical Report, NOAA (2016).14.NPFMC. Fishery Management Plan for Groundfish of the Gulf of Alaska. Technical Report, North Pacific Fishery Management Council (2017).15.Pennoyer, S. & Balsiger, J. Groundfish Total Allowable Catch Specifications and Prohibited Species Catch Limits Under the Authority of the Fishery Management Plans for the Groundfish Fishery of the Bering Sea and Aleutian Islands Area and Groundfish of the Gulf of Alaska: Final Supplemental Environmental Impact Statement. Technical Report, United States National Marine Fisheries Service Alaska Regional Office, Juneau, Alaska (1998).16.Rodgveller, C. J., Lunsford, C. R. & Fujioka, J. T. Evidence of hook competition in longline surveys. Fish. Bull. 106, 364–374 (2008).
Google Scholar 
17.Hutchings, J. A. Collapse and recovery of marine fishes. Nature 406, 882–885. https://doi.org/10.1038/35022565 (2000).ADS 
CAS 
Article 
PubMed 

Google Scholar 
18.Moore, J. A. Deep-sea finfish fisheries: lessons from history. Fisheries 24, 16–21 (1999).Article 

Google Scholar 
19.Moore, J. & Mace, P. Challenges and prospects for deep-sea finfish fisheries. Fisheries 24, 22–23 (1999).Article 

Google Scholar 
20.Rijnsdorp, A. D., Peck, M. A., Engelhard, G. H., Möllmann, C. & Pinnegar, J. K. Resolving the effect of climate change on fish populations. ICES J. Mar. Sci. J. Conseil 66, 1570–1583. https://doi.org/10.1093/icesjms/fsp056 (2009).Article 

Google Scholar 
21.Xia, Y. & Li, W. K. On single-index coefficient regression models. J. Am. Stat. Assoc. 94, 1275–1285. https://doi.org/10.1080/01621459.1999.10473880 (1999).MathSciNet 
Article 
MATH 

Google Scholar 
22.Kammann, E. E. & Wand, M. P. Geoadditive models. J. R. Stat. Soc. Ser. C (Appl. Stat.) 52, 1–18. https://doi.org/10.1111/1467-9876.00385 (2003).MathSciNet 
Article 
MATH 

Google Scholar 
23.Lu, Z., Steinskog, D. J., Tjøstheim, D. & Yao, Q. Adaptively varying-coefficient spatiotemporal models. J. R. Stat. Soc. Ser. B (Stat. Methodol.) 71, 859–880. https://doi.org/10.1111/j.1467-9868.2009.00710.x (2009).24.Ruppert, D., Wand, M. P. & Carroll, R. J. Semiparametric Regression. Cambridge Series in Statistical and Probabilistic Mathematics (Cambridge University Press, 2003).25.Scheipl, F., Staicu, A.-M. & Greven, S. Functional additive mixed models. J. Comput. Graph. Stat. 24, 477–501. https://doi.org/10.1080/10618600.2014.901914 (2015).MathSciNet 
Article 
PubMed 
MATH 

Google Scholar 
26.Wood, S. N. Generalized Additive Models. Texts in Statistical Science Series (Chapman & Hall/CRC, 2006). An Introduction with (R).27.Wood, S. N., Scheipl, F. & Faraway, J. J. Straightforward intermediate rank tensor product smoothing in mixed models. Stat. Comput. 23, 341–360. https://doi.org/10.1007/s11222-012-9314-z (2013).MathSciNet 
Article 
MATH 

Google Scholar 
28.Conn, P. B., Johnson, D. S. & Boveng, P. L. On extrapolating past the range of observed data when making statistical predictions in ecology. PLoS ONE 10, 1–16. https://doi.org/10.1371/journal.pone.0141416 (2015).CAS 
Article 

Google Scholar 
29.Hodges, J. S. & Reich, B. J. Adding spatially-correlated errors can mess up the fixed effect you love. Am. Stat. 64, 325–334. https://doi.org/10.1198/tast.2010.10052 (2010).MathSciNet 
Article 
MATH 

Google Scholar 
30.Kim, M. & Wang, L. Generalized spatially varying coefficient models. J. Comput. Graph. Stat.https://doi.org/10.1080/10618600.2020.1754225 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
31.Mu, J., Wang, G. & Wang, L. Estimation and inference in spatially varying coefficient models. Environmetrics 29, e2485. https://doi.org/10.1002/env.2485 (2018).MathSciNet 
Article 

Google Scholar 
32.Brumback, B. A. & Rice, J. A. Smoothing spline models for the analysis of nested and crossed samples of curves. J. Am. Stat. Assoc. 93, 961–976. https://doi.org/10.1080/01621459.1998.10473755 (1998).MathSciNet 
Article 
MATH 

Google Scholar 
33.Augustin, N. H., Trenkel, V. M., Wood, S. N. & Lorance, P. Space-time modelling of blue ling for fisheries stock management. Environmetrics 24, 109–119. https://doi.org/10.1002/env.2196 (2013).MathSciNet 
Article 

Google Scholar 
34.Finley, A. O. Comparing spatially-varying coefficients models for analysis of ecological data with non-stationary and anisotropic residual dependence. Methods Ecol. Evol. 2, 143–154. https://doi.org/10.1111/j.2041-210X.2010.00060.x (2011).Article 

Google Scholar 
35.Al-Sulami, D., Jiang, Z., Lu, Z. & Zhu, J. Estimation for semiparametric nonlinear regression of irregularly located spatial time-series data. Econom. Stat. 2, 22–35. https://doi.org/10.1016/j.ecosta.2017.01.002 (2017).MathSciNet 
Article 

Google Scholar 
36.Gelfand, A. E., Kim, H.-J., Sirmans, C. F. & Banerjee, S. Spatial modeling with spatially varying coefficient processes. J. Am. Stat. Assoc. 98, 387–396. https://doi.org/10.1198/016214503000170 (2003).MathSciNet 
Article 
MATH 

Google Scholar 
37.Feng, S. & Xue, L. Model detection and estimation for single-index varying coefficient model. J. Multivariate Anal. 139, 227–244. https://doi.org/10.1016/j.jmva.2015.03.008 (2015).MathSciNet 
Article 
MATH 

Google Scholar 
38.Zhao, P. & Xue, L. Variable selection for semiparametric varying coefficient partially linear models. Stat. Probab. Lett. 79, 2148–2157 (2009).MathSciNet 
Article 
MATH 

Google Scholar 
39.Guisan, A. Jr. & Hastie, T. Generalized linear and generalized additive models in studies of species distributions: setting the scene. Ecol. Model. 157, 89–100. https://doi.org/10.1016/S0304-3800(02)00204-1 (2002).Article 

Google Scholar 
40.Zhang, L. & Gove, J. H. Spatial assessment of model errors from four regression techniques. For. Sci. 51, 334–346. https://doi.org/10.1093/forestscience/51.4.334 (2005).Article 

Google Scholar 
41.Cai, A., Tsay, R. S. & Chen, R. Variable selection in linear regression with many predictors. J. Comput. Graph. Stat. 18, 573–591 (2009).MathSciNet 
Article 

Google Scholar 
42.Tibshirani, R. Regression shrinkage and selection via the lasso. J. R. Stat. Soc. Ser. B (Stat. Methodol.) 58, 267–288. https://doi.org/10.1111/j.2517-6161.1996.tb02080.x (1996).43.Ledolter, J. Penalty-Based Variable Selection in Regression Models with Many Parameters (LASSO), chap. 6, 71–82 (Wiley, 2013). https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781118596289.ch6.44.Yuan, M. & Lin, Y. Model selection and estimation in regression with grouped variables. J. R. Stat. Soc. Ser. B (Stat. Methodol.) 68, 49–67. https://doi.org/10.1111/j.1467-9868.2005.00532.x (2006).45.Fan, J., Yao, Q. & Cai, Z. Adaptive varying-coefficient linear models. J. R. Stat. Soc. Ser. B (Stat. Methodol.) 65, 57–80 (2003).46.Matsui, H. & Misumi, T. Variable selection for varying-coefficient models with the sparse regularization. Comput. Stat. 30, 43–55 (2015).MathSciNet 
Article 
MATH 

Google Scholar 
47.Wang, H. & Xia, Y. Shrinkage estimation of the varying coefficient model. J. Am. Stat. Assoc. 104, 747–757. https://doi.org/10.1198/jasa.2009.0138 (2009).MathSciNet 
CAS 
Article 
MATH 

Google Scholar 
48.Xue, L. & Qu, A. Variable selection in high-dimensional varying-coefficient models with global optimality. J. Mach. Learn. Res. 13, 1973–1998 (2012).MathSciNet 
MATH 

Google Scholar 
49.Feng, S. & Xue, L. Variable selection for single-index varying-coefficient model. Front. Math. China 8, 541–565. https://doi.org/10.1007/s11464-013-0284-z (2013).MathSciNet 
Article 
MATH 

Google Scholar 
50.Song, Y., Jian, L. & Lin, L. Robust exponential squared loss-based variable selection for high-dimensional single-index varying-coefficient model. J. Comput. Appl. Math. 308, 330–345. https://doi.org/10.1016/j.cam.2016.05.030 (2016).MathSciNet 
Article 
MATH 

Google Scholar 
51.Yang, J. & Yang, H. Robust modal estimation and variable selection for single-index varying-coefficient models. Commun. Stat. Simul. Comput 46, 2976–2997. https://doi.org/10.1080/03610918.2015.1069346 (2017).MathSciNet 
Article 
MATH 

Google Scholar 
52.Wei, F., Huang, J. & Li, H. Variable selection and estimation in high-dimensional varying-coefficient models. Stat. Sin. 21, 1515–1540. https://doi.org/10.5705/ss.2009.316 (2011).MathSciNet 
Article 
PubMed 
PubMed Central 
MATH 

Google Scholar 
53.Sun, W., Bindele, H. F., Abebe, A. & Correia, H. E. Robust functional coefficient selection for the single-index varying coefficients regression model. J. Stat. Comput. Simul. 1, 17. https://doi.org/10.1080/00949655.2020.1867548 (2021).Article 

Google Scholar 
54.Alaska Fisheries Science Center. AFSC/ABL: Longline Sablefish Survey. https://noaa-fisheries-afsc.data.socrata.com/dataset/AFSC-ABL-Longline-Sablefish-Survey/itxd-qjvg/data (2019). Accessed 14 Apr 2014.55.Sigler, M. F. & Lunsford, C. R. Survey Protocol for the Alaska Sablefish Longline Survey. Technical Report, Alaska Fisheries Science Center (2009).56.National Data Buoy Center. Meteorological and oceanographic data collected from the National Data Buoy Center Coastal-Marine Automated Network (C-MAN) and moored (weather) buoys. https://accession.nodc.noaa.gov/NDBC-CMANWx (2018).57.Alaska Fisheries Science Center. AFSC/RACE/GAP: RACEBASE Database. Online: http://www.afsc.noaa.gov/RACE/groundfish/survey_data/default.htm (2019).58.O’Brien, T. D. COPEPOD: The Global Plankton Database. A Review of the 2007 Database Contents and New Quality Control Methodology. Technical Report. NOAA Tech. Memo. NMFS-F/ST-34, U.S. Dep. Commerce (2007).59.Boyer, T. P. et al. World Ocean Database 2013. Technical Report. National Oceanographic Data Center, Ocean Climate Laboratory, NOAA (2013). https://doi.org/10.7289/V5NZ85MT.60.Neter, J., Kutner, M. H., Nachtsheim, C. J. & Wasserman, W. Applied Linear Statistical Models, vol. 4 (Irwin Chicago, 1996).61.O’Brien, R. M. A caution regarding rules of thumb for variance inflation factors. Qual. Quant. 41, 673–690 (2007).Article 

Google Scholar 
62.Li, L., Losser, T., Yorke, C. & Piltner, R. Fast inverse distance weighting-based spatiotemporal interpolation: a web-based application of interpolating daily fine particulate matter pm2.5 in the contiguous U.S. using parallel programming and k–d tree. Int. J. Environ. Res. Public Health 11, 9101–9141. https://doi.org/10.3390/ijerph110909101 (2014).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
63.Melo, C. & Melo, O. geosptdb: Spatio-Temporal Inverse Distance Weighting and Radial Basis Functions with Distance-Based Regression (2015). R package version 0.5-0.64.Whitney, F. A. Nutrient variability in the mixed layer of the subarctic Pacific Ocean, 1987–2010. J. Oceanogr. 67, 481–492. https://doi.org/10.1007/s10872-011-0051-2 (2011).CAS 
Article 

Google Scholar 
65.Whitney, F. A., Bograd, S. J. & Ono, T. Nutrient enrichment of the subarctic Pacific Ocean pycnocline. Geophys. Res. Lett. 40, 2200–2205. https://doi.org/10.1002/grl.50439 (2013).ADS 
CAS 
Article 

Google Scholar 
66.Brodeur, R. D. & Ware, D. M. Long-term variability in zooplankton biomass in the subarctic Pacific Ocean. Fish. Oceanogr. 1, 32–38. https://doi.org/10.1111/j.1365-2419.1992.tb00023.x (1992).Article 

Google Scholar 
67.Chiba, S., Tadokoro, K., Sugisaki, H. & Saino, T. Effects of decadal climate change on zooplankton over the last 50 years in the western subarctic North Pacific. Glob. Change Biol. 12, 907–920. https://doi.org/10.1111/j.1365-2486.2006.01136.x (2006).ADS 
Article 

Google Scholar 
68.Childers, A. R., Whitledge, T. E. & Stockwell, D. A. Seasonal and interannual variability in the distribution of nutrients and chlorophyll a across the Gulf of Alaska shelf: 1998–2000. Deep Sea Res. Part II Top. Stud. Oceanogr. 52, 193–216. https://doi.org/10.1016/j.dsr2.2004.09.018 (2005). U.S. GLOBEC Biological and Physical Studies of Plankton, Fish and Higher Trophic Level Production, Distribution, and Variability in the Northeast Pacific.69.Sackmann, B., Mack, L., Logsdon, M. & Perry, M. J. Seasonal and inter-annual variability of SeaWiFS-derived chlorophyll a concentrations in waters off the Washington and Vancouver Island coasts, 1998–2002. Deep Sea Res. Part II Top. Stud. Oceanogr. 51, 945–965 (2004).ADS 
CAS 
Article 

Google Scholar 
70.Wong, C. et al. Seasonal cycles of nutrients and dissolved inorganic carbon at high and mid latitudes in the North Pacific Ocean during the Skaugran cruises: determination of new production and nutrient uptake ratios. Deep Sea Res. Part II Top. Stud. Oceanogr. 49, 5317–5338 (2002).ADS 
CAS 
Article 

Google Scholar 
71.Houde, E. D. chap. Recruitment variability. In Fish Reproductive Biology: Implications for Assessment and Management (Eds. Jakobsen, T., Fogarty, M. J., Megrey, B. A. & Moksness, E.) (Wiley, 2016).72.Mason, J. C., Beamish, R. J. & McFarlane, G. A. Sexual maturity, fecundity, spawning, and early life history of sablefish (Anoplopoma fimbria) off the Pacific Coast of Canada. Can. J. Fish. Aquat. Sci. 40, 2126–2134. https://doi.org/10.1139/f83-247 (1983).Article 

Google Scholar 
73.Stark, J. W. Geographic and seasonal variations in maturation and growth of female Pacific cod (Gadus macrocephalus) in the Gulf of Alaska and Bering Sea. Fish. Bull. 105, 396–407 (2007).
Google Scholar 
74.Clark, W. G., Hare, S. R., Parma, A. M., Sullivan, P. J. & Trumble, R. J. Decadal changes in growth and recruitment of Pacific halibut (Hippoglossus stenolepis). Can. J. Fish. Aquat. Sci. 56, 242–252. https://doi.org/10.1139/f98-163 (1999).Article 

Google Scholar 
75.Echeverria, T. W. Thirty-four species of California rockfishes: maturity and seasonality of reproduction. Fish. Bull. 85, 229–250 (1987).
Google Scholar 
76.Di Lorenzo, E. North Pacific Gyre Oscillation (2018). NPGO index.77.NOAA ESRL Physical Sciences Division. Multivariate ENSO Index Version 2 (MEI.v2) (2019). ENSO index.78.Mantua, N. J. & JISAU, University of Washington. The Pacific Decadal Oscillation. http://research.jisao.washington.edu/pdo/ (2016).79.Di Lorenzo, E. et al. Central Pacific El Niño and decadal climate change in the North Pacific Ocean. Nat. Geosc. 3, 762–765. https://doi.org/10.1038/ngeo984 (2010).ADS 
CAS 
Article 

Google Scholar 
80.Ladd, C. & Stabeno, P. J. Stratification on the Eastern Bering Sea shelf revisited. Deep Sea Res. Part II Top. Stud. Oceanogr. 65-70, 72–83. https://doi.org/10.1016/j.dsr2.2012.02.009 (2012). Understanding Ecosystem Processes in the Eastern Bering Sea.81.Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B (Methodol.) 57, 289–300 (1995).82.Frainer, A. et al. Climate-driven changes in functional biogeography of Arctic marine fish communities. Proc. Natl. Acad. Sci. 114, 12202–12207, https://doi.org/10.1073/pnas.1706080114 (2017).83.Hewitt, J. E., Ellis, J. I. & Thrush, S. F. Multiple stressors, nonlinear effects and the implications of climate change impacts on marine coastal ecosystems. Glob. Change Biol. 22, 2665–2675. https://doi.org/10.1111/gcb.13176 (2016).ADS 
Article 

Google Scholar 
84.Liu, H. et al. Nonlinear dynamic features and co-predictability of the Georges Bank fish community. Mar. Ecol. Prog. Ser. 464, 195–207 (2012).ADS 
Article 

Google Scholar 
85.Echave, K., Rodgveller, C. & Shotwell, S. K. Calculation of the Geographic Area Sizes Used to Create Population Indices for the Alaska Fisheries Science Center Longline Survey. Technical Report. NOAA Tech. Memo. NMFS-AFSC-253, U.S. Department of Commerce (2013).86.Webster, R. A., Soderlund, E., Dykstra, C. L. & Stewart, I. J. Monitoring change in a dynamic environment: spatio-temporal modelling of calibrated data from different types of fisheries surveys of Pacific halibut. Can. J. Fish. Aquat. Sci.https://doi.org/10.1139/cjfas-2019-0240 (2020).Article 

Google Scholar 
87.Spencer, P. D., Hollowed, A. B., Sigler, M. F., Hermann, A. J. & Nelson, M. W. Trait-based climate vulnerability assessments in data-rich systems: an application to eastern Bering Sea fish and invertebrate stocks. Glob. Change Biol. 25, 3954–3971. https://doi.org/10.1111/gcb.14763 (2019).ADS 
Article 

Google Scholar 
88.Sogard, S. M. & Olla, B. L. Growth and behavioral responses to elevated temperatures by juvenile sablefish Anoplopoma fimbria and the interactive role of food availability. Mar. Ecol. Prog. Ser. 217, 121–134 (2001).ADS 
Article 

Google Scholar 
89.Stoner, A. W. & Sturm, E. A. Temperature and hunger mediate sablefish (Anoplopoma fimbria) feeding motivation: implications for stock assessment. Can. J. Fish. Aquat. Sci. 61, 238–246. https://doi.org/10.1139/f03-170 (2004).Article 

Google Scholar 
90.Sogard, S. Interannual variability in growth rates of early juvenile sablefish and the role of environmental factors. Bull. Mar. Sci.https://doi.org/10.5343/bms.2010.1045 (2011).Article 

Google Scholar 
91.Shotwell, S. K., Hanselman, D. H. & Belkin, I. M. Toward biophysical synergy: investigating advection along the Polar Front to identify factors influencing Alaska sablefish recruitment. Deep Sea Res. Part II Top. Stud. Oceanogr. 107, 40–53. https://doi.org/10.1016/j.dsr2.2012.08.024 (2014). Fronts, Fish and Top Predators.92.Tolimieri, N., Haltuch, M. A., Lee, Q., Jacox, M. G. & Bograd, S. J. Oceanographic drivers of sablefish recruitment in the California current. Fish. Oceanogr. 27, 458–474. https://doi.org/10.1111/fog.12266 (2018).Article 

Google Scholar 
93.Harrison, P. J., Whitney, F. A., Tsuda, A., Saito, H. & Tadokoro, K. Nutrient and plankton dynamics in the NE and NW gyres of the subarctic Pacific Ocean. J. Oceanogr. 60, 93–117 (2004).CAS 
Article 

Google Scholar 
94.Coffin, B. & Mueter, F. Environmental covariates of sablefish (Anoplopoma fimbria) and Pacific ocean perch (Sebastes alutus) recruitment in the Gulf of Alaska. Deep Sea Res. Part II Top. Stud. Oceanogr. 132, 194–209. https://doi.org/10.1016/j.dsr2.2015.02.016 (2016). Understanding Ecosystem Processes in the Gulf of Alaska: Volume 1.95.Hagen, P. T. & Quinn, T. J. Long-term growth dynamics of young Pacific halibut: evidence of temperature-induced variation. Fish. Res. 11, 283–306. https://doi.org/10.1016/0165-7836(91)90006-2 (1991). Fish Population Dynamics: Solving Fishery Management Problems.96.Hurst, T. P., Spencer, M. L., Sogard, S. M. & Stoner, A. W. Compensatory growth, energy storage and behavior of juvenile Pacific halibut Hippoglossus stenolepis following thermally induced growth reduction. Mar. Ecol. Prog. Ser. 293, 233–240 (2005).ADS 
Article 

Google Scholar 
97.Holsman, K. K., Aydin, K., Sullivan, J., Hurst, T. & Kruse, G. H. Climate effects and bottom-up controls on growth and size-at-age of Pacific halibut (Hippoglossus stenolepis) in Alaska (USA). Fish. Oceanogr. 28, 345–358. https://doi.org/10.1111/fog.12416 (2019).Article 

Google Scholar 
98.Lynam, C. P. et al. Interaction between top-down and bottom-up control in marine food webs. Proc. Natl. Acad. Sci. 114, 1952–1957. https://doi.org/10.1073/pnas.1621037114 (2017).ADS 
CAS 
Article 
PubMed 

Google Scholar 
99.Desmit, X., Ruddick, K. & Lacroix, G. Salinity predicts the distribution of chlorophyll a spring peak in the southern North Sea continental waters. J. Sea Res. 103, 59–74. https://doi.org/10.1016/j.seares.2015.02.007 (2015).ADS 
Article 

Google Scholar 
100.Benson, A. J. & Trites, A. W. Ecological effects of regime shifts in the Bering Sea and eastern North Pacific Ocean. Fish Fish. 3, 95–113. https://doi.org/10.1046/j.1467-2979.2002.00078.x (2002).Article 

Google Scholar 
101.Feng, J. et al. Contrasting correlation patterns between environmental factors and chlorophyll levels in the global ocean. Glob. Biogeochem. Cycles 29, 2095–2107. https://doi.org/10.1002/2015GB005216 (2015).ADS 
CAS 
Article 

Google Scholar 
102.Kahru, M. et al. Global correlations between winds and ocean chlorophyll. J. Geophys. Res. Oceanshttps://doi.org/10.1029/2010JC006500 (2010).Article 

Google Scholar 
103.Sadorus, L. L., Mantua, N. J., Essington, T., Hickey, B. & Hare, S. Distribution patterns of Pacific halibut (Hippoglossus stenolepis) in relation to environmental variables along the continental shelf waters of the US West Coast and southern British Columbia. Fish. Oceanogr. 23, 225–241. https://doi.org/10.1111/fog.12057 (2014).Article 

Google Scholar 
104.Barbeaux, S. et al. Gulf of Alaska Stock Assessments. Technical Report, North Pacific Fishery Management Council, Anchorage, AK (2018).105.Barbeaux, S. J. & Hollowed, A. B. Ontogeny matters: climate variability and effects on fish distribution in the eastern Bering Sea. Fish. Oceanogr. 27, 1–15. https://doi.org/10.1111/fog.12229 (2018).Article 

Google Scholar 
106.Yang, Q. et al. How “The Blob” affected groundfish distributions in the Gulf of Alaska. Fish. Oceanogr. 28, 434–453. https://doi.org/10.1111/fog.12422 (2019).107.Hagens, M. & Middelburg, J. J. Attributing seasonal pH variability in surface ocean waters to governing factors. Geophys. Res. Lett. 43, 12528–12537. https://doi.org/10.1002/2016GL071719 (2016).ADS 
CAS 
Article 

Google Scholar 
108.Fry, C. H., Tyrrell, T., Hain, M. P., Bates, N. R. & Achterberg, E. P. Analysis of global surface ocean alkalinity to determine controlling processes. Mar. Chem. 174, 46–57 (2015).CAS 
Article 

Google Scholar 
109.Bromhead, D. et al. The potential impact of ocean acidification upon eggs and larvae of yellowfin tuna (Thunnus albacares). Deep Sea Res. Part II Top. Stud. Oceanogr. 113, 268–279, https://doi.org/10.1016/j.dsr2.2014.03.019 (2015). Impacts of climate on marine top predators.110.Doney, S. C. et al. Impact of anthropogenic atmospheric nitrogen and sulfur deposition on ocean acidification and the inorganic carbon system. Proc. Natl. Acad. Sci. 104, 14580–14585. https://doi.org/10.1073/pnas.0702218104 (2007).111.Napp, J. M. & Hunt, G. L. Anomalous conditions in the south-eastern Bering Sea 1997: linkages among climate, weather, ocean, and biology. Fish. Oceanogr. 10, 61–68. https://doi.org/10.1046/j.1365-2419.2001.00155.x (2001).Article 

Google Scholar 
112.Noakes, D. J. & Beamish, R. J. Synchrony of marine fish catches and climate and ocean regime shifts in the North Pacific Ocean. Mar. Coast. Fish. 1, 155–168. https://doi.org/10.1577/C08-001.1 (2009).Article 

Google Scholar  More

Data setsRadioactivity measurement data for several species of wild game mammals and birds in Fukushima Prefecture from May 2011 to March 2018 were released to the public by the Fukushima Prefecture Government (https://emdb.jaea.go.jp/emdb/en/portals/1040501000/). We extracted the data for wild boar (Sus scrofa), 1404 samples, and Asian black bear (Ursus thibetanus), 422 samples. The resulting boar and bear data sets contained total radioactive cesium activity (134Cs + 137Cs isotopes) values (in Bq/kg) from animals captured at different times and locations within Fukushima Prefecture. The data were imported for analysis into R 4.0.3 software21.We ln-transformed the cesium activity values to bring their distribution closer to normal, creating the variable LnCsTot. To facilitate regression analyses (described below), we removed instances of missing data and cesium levels below detection: 20 samples (1.4%) for boar and 15 samples (3.3%) for bears. The time when each sample was taken (labeled “Day of collection” in the Fukushima Prefecture Government data set) was converted to years since the Fukushima accident (since March 11, 2011), assuming that 1 year = 365.25 days. This time of sample collection in years was called variable T.Since for each sample some time passed between sample collection and radioactivity measurement (labeled “Result found Date”, called Tr in our notation), we needed to correct the reported LnCsTot values for physical decay over this time, which was different for different samples. The procedure used to perform this correction is described in Supplementary methods. The data with corrected total cesium values (LnCsc) are provided in Supplementary data (Supplementary_Dataset_File_Full).Mathematical modelTo describe the data on ln-transformed total radioactive cesium levels (LnScc) in each species as function of time after the accident (T), we developed the following simple mathematical model (Eqs. 1A, 1B):$${LnCs}_{c}=X+Q-mu times {T}^{nu }+Atimes mathrm{sin}left[2times pi times left(T+Pright)right], $$
(1A)
$$X=mathrm{ln}left[mathrm{exp}left(LnCs{134}_{t{0}_{r}}right)times {2}^{-frac{T}{{Th}_{Cs134}}}+mathrm{exp}left(LnCs{137}_{t{0}_{r}}right)times {2}^{-frac{T}{{Th}_{Cs137}}}right]$$
(1B)
Here the term X represents the estimated average radioactive cesium level in the studied area, based on the intercepts (LnCs134t0r for 134Cs and LnCs137t0r for 137Cs, respectively) from robust regression discussed in Supplementary methods, and taking into account only physical decay for each isotope (with half-lives of ThCs134 for 134Cs and ThCs134 for 137Cs, respectively). The terms Q, µ, ν, A and P represent adjustable model parameters. Parameter Q represents the fitted relationship between radioactive cesium levels in the animal (Bq/kg), relative to the external environment (Bq/m2). Parameter µ represents the net rate of radioactive cesium reduction in animal tissues over time due to all processes except physical decay (e.g. decrease in bioavailability due to migration of cesium into deeper soil layers, human-mediated cleanup efforts, etc.). Parameter ν is a potential power dependence for these processes. By default, ν was set to ν = 1, but exploratory calculations using ν = 2 or treating ν as a freely adjustable parameter (≥ 0.1) were performed as well. Parameters A and P in the sine function represent a sinusoidal approximation for seasonal changes in radioactive cesium levels in animal tissues (e.g. due to seasonal variations in diet and life style), where A is the amplitude of the oscillations, P is the phase shift, and the period is set to 1 year. For simplicity, these parameters were assumed to be the same for both studied cesium isotopes. The descriptions of each parameter are also presented in Table 1.Table 1 The meanings of all parameters used in our mathematical model (Eq. 1A, 1B) for radioactive cesium levels in wild boar (Sus scrofa) and Asian black bear (Ursus thibetanus).Full size tableModel fitting approachesInitially, we used nonlinear ordinary least squares (OLS) regression (nls R function) to fit the model (Eq. 1A, 1B) to the data. To find the global optimum fit, we repeated the fitting procedure 2000 times with slightly different random initial parameter values and recorded the solution with the smallest root mean squared error (RMSE). Diagnostics on this regression included checking of convergence criteria and analyses of residuals (by scatter plot and histogram, regressing residuals as function of T, visualizing the QQ plot, autocorrelation and partial autocorrelation functions with 95% confidence intervals, performing the Shapiro–Wilk normality test, and calculating skewness and kurtosis). For boar data, diagnostics revealed problems with convergence (both X-convergence and relative convergence) and non-normality of residuals: e.g. Shapiro–Wilk p-value = 1.476 × 10–7, skewness = − 0.37, kurtosis = 3.50. For black bear data similar problems occurred with convergence, but residuals were closer to the normal distribution (perhaps due to smaller sample size): e.g. Shapiro–Wilk p-value = 0.0526, skewness = − 0.058, kurtosis = 2.45.Due to these issues, we used robust nonlinear regression (nlrob R package) to reduce the effects of “outlier” data points. To find the global optimum, we repeated the fitting procedure 2000 times with slightly different random initial parameter values and selected the solution with the smallest absolute value of median residuals. The best-fit parameters for OLS and robust regressions were somewhat different for both boar and bear data. For boar data, the minimum robustness weight was 0.339 and the median was 0.762, and the corresponding values for black bear data were 0.557 and 0.821, respectively.For each species, we compared the performances of model variants with different assumptions about parameter ν: (1) The default case with ν = 1, which represents an exponential rate of radioactive cesium decrease due to processes other than physical decay. (2) The case with ν = 2, which represents quadratic decay. (3) The case with ν being freely adjustable (≥ 0.1). The comparisons were based on Akaike information criterion (AIC)22,23. The purpose of these calculations was to better assess the shape of the time course for non-physical factors involved in radioactive cesium level decline in animal tissues over time after the accident.In addition to analyzing the full data set for each species, we also performed separate analyses on subsets of data from specific locations: from those districts of Fukushima Prefecture where the mean radioactive cesium levels in animal samples were the highest, and where a sufficiently large number of samples was present. For wild boar, the two selected districts for this subset analysis were Soso and Kenpoku (819 samples), and for black bear they were Kenpoku and Kenchu (163 samples).To further assess the sensitivity of model results to geographical and temporal factors, we also constructed a separate subset of data for each species. This subset excluded data from the Aizu and Minamiaizu districts, which are separated by mountains from the Fukushima Daiichi Nuclear Power Plant, and excluded data collected ≤ 6 months after the accident. These restrictions were intended to determine model performance on data collected in a more geographically contiguous area after the initial abrupt changes in contamination levels were completed and the system entered the phase of more stable kinetics. The purpose of all these analyses was to assess whether the time course of radioactive cesium levels in the bodies of each species differed between locations with high contamination vs. those with lower contamination, and as function of time after the accident.We were interested in quantifying not only the center of the distribution of radioactive cesium values in each species over time, but also in assessing the lower and upper tails of this distribution. For this purpose, we fitted the model (Eq. 1A, 1B) for each species using quantile regression (nlrq function in quantreg R package) for the median (50th percentile), and also for the 25th and 75th percentiles. Initial parameter estimates for the quantile regressions were taken from best-fit parameters from robust regression described above. The 25th and 75th percentiles were selected instead of more extreme values (e.g. 5th and 95th) because the latter resulted in poor quality fits due to limited amounts of data at the fringes of the distribution.To assess the variability of model parameters by location in more detail, we used mixed effects modeling (nlme R package) on the data from each species. Since original OLS fits suggested substantial deviations of residuals from the normality assumption, we performed mixed effects modeling on data with some outlier data points removed. The OutlierDetection package in R removed 43 boar samples and 5 bear samples. These outliers are listed in the Supplementary_outlier_data_points file. The remaining samples were used for mixed effects model fitting, but model performance metrics like coefficient of determination (R2) and RMSE were assessed on the full data set (with outliers included) for each species.Since the Fligner-Killeen test of homogeneity of variances by district generated low p-values for both species (4.6 × 10–14 for boar and 0.018 for black bear), we allowed modelled variances to differ by district (using the weights option in nlme). We investigated several random effects structures for some or all model parameters, with randomness by district only, or by district and municipality within district. Model diagnostics were the same as for fixed effects OLS modeling described above, and also included boxplots of model residuals by district. The mixed effects model variants with different random effects structures were compared using the anova function in R, and also by assessing convergence criteria, normality of residuals, skewness, and kurtosis. Consequently, preferred mixed effects model variants were selected for the full data as well as for the subset of two districts with high radioactive cesium levels, separately for each species.Model extrapolation from training to testing dataTo investigate how the robust and quantile regression fits of our model could extrapolate beyond the time range that was used for model fitting, we split the data for each species into “training” (early times) and “testing” (later times) parts. The split was done based on time since the accident (T variable), so that approximately ½ of the samples were assigned to the training and testing sets, respectively. For wild boar data, the training set included times between 0.20 and 3.45 years after the accident, and the testing set included times between 3.45 and 7.03 years. For black bear data, the training set included times between 0.42 and 3.46 years after the accident, and the testing set included times between 3.46 and 6.87 years.We also evaluated an alternative approach to splitting the data, where the split was done randomly instead of by time. In other words, any data point regardless of time had an equal probability of being assigned to either the training or the testing data set. Both the training and testing data subsets generated by this random split included the complete time range. This approach was implemented in context of the sensitivity analysis described above.For each species, robust and quantile regressions were fitted to training data, and their predictions were calculated for testing data. For robust regression, RMSE was calculated on testing data for two scenarios: (1) for the model fitted to training data only, and (2) for the model fitted over the entire data range (training + testing). These RMSE values for conditions 1 and 2 were compared to assess the quality of model extrapolation. Extrapolation performance for robust and quantile regressions was also assessed visually by plotting the model predictions and data.Application of the model to wild boar data from the Chernobyl accident areaTo compare the results of our analysis of wild boar contamination with radioactive cesium in the area affected by the Fukushima accident with data from another location, we also analyzed wild boar data from the Chernobyl accident area. These data were published by Gulakov14 and contain summaries of 137Cs contamination levels in the muscles of 188 boar collected between 1991 and 2008 (i.e. from 5 to 22 years after the 1986 accident). Sampling was carried out in three zones with different land contamination levels with 137Cs. This data set provides important information on radioactive cesium contamination in wild boar in the Chernobyl area. Unfortunately, 137Cs measurements in each sampled boar were not provided by Gulakov14, and only summary statistics are available for each zone and year after the accident (Tables 1–3 in reference14): number of animals, mean, minimum and maximum 137Cs levels.We could not apply the full model (Eq. 1A, 1B) to these summary data which lacked seasonality information and 134Cs data. However, we were able to perform a weighted linear regression to quantify the ecological half-life of 137Cs in Chernobyl boar and the relationship between 137Cs levels in the animals (Bq/kg), relative to the external environment (Bq/m2). The data used for this analysis, derived from Gulakov14, are provided in Supplementary data (Supplementary_Dataset_File_Full). They contain the following variables. Zone = location of sample collection (Alienation, Permanent control or Periodic control). Time = time in years after the Chernobyl accident. LnMeanCs = ln-transformed mean 137Cs level in boar muscle (Bq/kg). LnMeanCs_c = LnMeanCs − X, where X is ln-transformed 137Cs land contamination (Bq/m2) in the given zone, corrected for physical decay of 137Cs. Weight = weighting of each data point used for regression. Weight = number of animals/(ln[maximum 137Cs level] − ln[minimum 137Cs level])2. These approximately inverse-variance weights were normalized by the overall mean across all data points, so that the mean weight across all data points was set to 1.These data were analyzed by weighted linear regression in R, where LnMeanCs_c was allowed to depend on Time and Zone variables. Model variants containing all possible combinations and pairwise interactions between these predictor variables were fitted and their performances were compared using the Akaike information criterion with correction for small sample size (AICc). These calculations were performed using the glmulti R package. Multimodel inference (MMI) was performed on this collection of fitted model variants. It resulted in the calculation of model-averaged parameter estimates, 95% CIs and importance scores, corrected for model selection uncertainty. Of main interest here were the intercept parameter, which is analogous to parameter Q in the full model (Eq. 1A, 1B), and the Time parameter, which is analogous to parameter µ in the full model. The ecological half-life for 137Cs was calculated based on the Time parameter. More

1.Dong Y, Silbermann M, Speiser A, Forieri I, Linster E, Poschet G, et al. Sulfur availability regulates plant growth via glucose-TOR signaling. Nat Commun. 2017;8:1174.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
2.Freney JR, Melville GE, Williams CH. Soil organic matter fractions as sources of plant-available sulphur. Soil Biol Biochem. 1975;7:217–21.CAS 
Article 

Google Scholar 
3.Kopittke PM, Dalal RC, Finn D, Menzies NW. Global changes in soil stocks of carbon, nitrogen, phosphorus, and sulphur as influenced by long‐term agricultural production. Global Change Biol. 2017;23:2509–19.Article 

Google Scholar 
4.Ciaffi M, Paolacci AR, Celletti S, Catarcione G, Kopriva S, Astolfi S. Transcriptional and physiological changes in the S assimilation pathway due to single or combined S and Fe deprivation in durum wheat (Triticum durum L.) seedlings. J Exp Bot. 2013;64:1663.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Ma Q, Luo Y, Wen Y, Hill PW, Chadwick DR, Wu L, et al. Carbon and sulphur tracing from soil organic sulphur in plants and soil microorganisms. Soil Biol Biochem. 2020;150:107971.CAS 
Article 

Google Scholar 
6.Piotrowska-Długosz A, Siwik-Ziomek A, Długosz J, Gozdowski D. Spatio-temporal variability of soil sulfur content and arylsulfatase activity at a conventionally managed arable field. Geoderma. 2017;295:107–18.Article 
CAS 

Google Scholar 
7.Fitzgerald JW, Watwood ME. Amino-acid metabolism in forest soil—Isolation and turnover of organic matter covalently labelled with 35S-methionine. Soil Biol Biochem. 1988;20:833–8.CAS 
Article 

Google Scholar 
8.Ma Q, Wen Y, Pan W, Macdonald A, Hill PW, Chadwick DR, et al. Soil carbon, nitrogen, and sulphur status affects the metabolism of organic S but not its uptake by microorganisms. Soil Biol Biochem. 2020;149:107943.CAS 
Article 

Google Scholar 
9.Jan MT, Roberts P, Tonheim SK, Jones DL. Protein breakdown represents a major bottleneck in nitrogen cycling in grassland soils. Soil Biol Biochem. 2009;41:2272–82.CAS 
Article 

Google Scholar 
10.Farrell M, Macdonald LM, Hill PW, Wanniarachchi SD, Farrar J, Bardgett RD, et al. Amino acid dynamics across a grassland altitudinal gradient. Soil Biol Biochem. 2014;76:179–82.CAS 
Article 

Google Scholar 
11.Wilkinson A, Hill PW, Farrar JF, Jones DL, Bardgett RD. Rapid microbial uptake and mineralization of amino acids and peptides along a grassland productivity gradient. Soil Biol Biochem. 2014;72:75–83.CAS 
Article 

Google Scholar 
12.Hill PW, Jones DL. Plant–microbe competition: does injection of isotopes of C and N into the rhizosphere effectively characterise plant use of soil N? N Phytol. 2018;221:796–806.Article 
CAS 

Google Scholar 
13.Godwin CM, Cotner JB. Aquatic heterotrophic bacteria have highly flexible phosphorus content and biomass stoichiometry. ISME J. 2015;9:2324–27.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
14.Hartman WH, Ye R, Horwath WR, Tringe SG. A genomic perspective on stoichiometric regulation of soil carbon cycling. ISME J. 2017;11:2652–65.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Manzoni S, Jackson RB, Trofymow JA, Porporato A. The global stoichiometry of litter nitrogen mineralization. Science. 2008;321:684–6.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
16.Cui J, Zhu Z, Xu X, Liu S, Jones DL, Kuzyakov Y, et al. Carbon and nitrogen recycling from microbial necromass to cope with C:N stoichiometric imbalance by priming. Soil Biol Biochem. 2020;142:107720.CAS 
Article 

Google Scholar 
17.Wei X, Zhu Z, Liu Y, Luo Y, Deng Y, Xu X, et al. C:N:P stoichiometry regulates soil organic carbon mineralization and concomitant shifts in microbial community composition in paddy soil. Biol Fert Soils. 2020;56:1093–107.Article 
CAS 

Google Scholar 
18.Maria M, Wolfgang W, Sophie ZB, Andreas R. Stoichiometric imbalances between terrestrial decomposer communities and their resources: mechanisms and implications of microbial adaptations to their resources. Front Microbiol. 2014;5:22.
Google Scholar 
19.Qiao N, Xu XL, Hu YH, et al. Carbon and nitrogen additions induce distinct priming effects along an organic-matter decay continuum. Sci Rep-UK. 2016;6:19865.CAS 
Article 

Google Scholar 
20.Kirkby CA, Kirkegaard JA, Richardson AE, Wade LJ, Blanchard C, Batten G. Stable soil organic matter: a comparison of C:N:P:S ratios in Australian and other world soils. Geoderma. 2011;163:197–208.CAS 
Article 

Google Scholar 
21.Manzoni S, Čapek P, Mooshammer M, Lindahl BD, Richter A, Šantrůčková H. Optimal metabolic regulation along resource stoichiometry gradients. Ecol Lett. 2017;20:1182–91.PubMed 
Article 
PubMed Central 

Google Scholar 
22.Mooshammer M, Wanek W, Hämmerle I, Fuchslueger L, Hofhansl F, Knoltsch A, et al. Adjustment of microbial nitrogen use efficiency to carbon:nitrogen imbalances regulates soil nitrogen cycling. Nat Commun. 2014;5:3694.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
23.Delgado-Baquerizo M, Reich PB, Khachane AN, Campbell CD, Thomas N, Freitag TE, et al. It is elemental: soil nutrient stoichiometry drives bacterial diversity. Environ Microbiol. 2017;19:1176–88.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
24.Mayor JR, Mack MC, Schuur EAG. Decoupled stoichiometric, isotopic, and fungal responses of an ectomycorrhizal black spruce forest to nitrogen and phosphorus additions. Soil Biol Biochem. 2015;88:247–56.CAS 
Article 

Google Scholar 
25.Mooshammer M, Hofhansl F, Frank AH, Wanek W, Hämmerle I, Leitner S, et al. Decoupling of microbial carbon, nitrogen, and phosphorus cycling in response to extreme temperature events. Sci Adv. 2017;3:e1602781.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
26.Chen J, Seven J, Zilla T, Dippold MA, Blagodatskaya E, Kuzyakov Y. Microbial C:N:P stoichiometry and turnover depend on nutrients availability in soil: A 14C, 15N and 33P triple labelling study. Soil Biol Biochem. 2019;131:206–16.CAS 
Article 

Google Scholar 
27.Ma Q, Wu L, Wang J, Ma J, Zheng N, Hill PW, et al. Fertilizer regime changes the competitive uptake of organic nitrogen by wheat and soil microorganisms: An in-situ uptake test using 13C, 15N labelling, and 13C-PLFA analysis. Soil Biol Biochem. 2018;125:319–27.CAS 
Article 

Google Scholar 
28.Jiang G, Zhang W, Xu M, Kuzyakov Y, Zhang X, Wang J, et al. Manure and Mineral Fertilizer Effects on Crop Yield and Soil Carbon Sequestration: a Meta‐Analysis and Modeling Across China. Glob Biogeochem Cy. 2018;32:1659–72.CAS 
Article 

Google Scholar 
29.Liu S, Wang J, Pu S, Blagodatskaya E, Kuzyakov Y, Razavi BS. Impact of manure on soil biochemical properties: a global synthesis. Sci Total Environ. 2020;745:141003.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
30.Cheng L, Zhang N, Yuan M, Xiao J, Qin Y, Deng Y, et al. Warming enhances old organic carbon decomposition through altering functional microbial communities. ISME J. 2017;11:1825–35.PubMed 
PubMed Central 
Article 

Google Scholar 
31.Ma Q, Wen Y, Wang D, Sun X, Hill PW, Macdonald A, et al. Farmyard manure applications stimulate soil carbon and nitrogen cycling by boosting microbial biomass rather than changing its community composition. Soil Biol Biochem. 2020;144:107760.CAS 
Article 

Google Scholar 
32.Ma Q, Wen Y, Ma J, Macdonald A, Hill PW, Chadwick DR, et al. Long-term farmyard manure application affects soil organic phosphorus cycling: a combined metagenomic and 33P/14C labelling study. Soil Biol Biochem. 2020;149:107959.CAS 
Article 

Google Scholar 
33.Jones DL, Hill PW, Smith AR, Farrell M, Ge T, Banning NC, et al. Role of substrate supply on microbial carbon use efficiency and its role in interpreting soil microbial community-level physiological profiles (CLPP). Soil Biol Biochem. 2018;123:1–6.CAS 
Article 

Google Scholar 
34.Jones DL, Magthab EA, Gleeson DB, Hill PW, Sánchez-Rodríguez AR, Roberts P, et al. Microbial competition for nitrogen and carbon is as intense in the subsoil as in the topsoil. Soil Biol Biochem. 2018;117:72–82.CAS 
Article 

Google Scholar 
35.Mariano E, Jones DL, Hill PW, Trivelin PCO. Mineralisation and sorption of dissolved organic nitrogen compounds in litter and soil from sugarcane fields. Soil Biol Biochem. 2016;103:522–32.CAS 
Article 

Google Scholar 
36.Corre M, Brumme R, Veldkamp EF. Changes in nitrogen cycling and retention processes in soils under spruce forests along a nitrogen enrichment gradient in Germany. Glob Change Biol. 2010;13:1509–27.Article 

Google Scholar 
37.Spohn M, Kuzyakov Y. Phosphorus mineralization can be driven by microbial need for carbon. Soil Biol Biochem. 2013;61:69–75.CAS 
Article 

Google Scholar 
38.Wu J, Joergensen RG, Pommerening B, Chaussod R, Brookes PC. Measurement of soil microbial biomass C by fumigation-extraction—an automated procedure. Soil Biol Biochem. 1990;22:1167–9.CAS 
Article 

Google Scholar 
39.Glanville H, Hill PW, Schnepf A, Oburger E, Jones DL. Combined use of empirical data and mathematical modelling to better estimate the microbial turnover of isotopically labelled carbon substrates in soil. Soil Biol Biochem. 2015;94:154–68.Article 
CAS 

Google Scholar 
40.Greenfield LM, Hill PW, Paterson E, Baggs EM, Jones DL. Methodological bias associated with soluble protein recovery from soil. Sci Rep-UK. 2018;8:11186.Article 
CAS 

Google Scholar 
41.Näsholm T, Kielland K, Ganeteg U. Uptake of organic nitrogen by plants. N Phytol. 2010;182:31–48.Article 
CAS 

Google Scholar 
42.Farrell M, Prendergast-Miller M, Jones DL, Hill PW, Condron LM. Soil microbial organic nitrogen uptake is regulated by carbon availability. Soil Biol Biochem. 2014;77:261–7.CAS 
Article 

Google Scholar 
43.Niknahad-Gharmakher H, Piutti S, Machet JM, Benizri E, Recous S. Mineralization-immobilization of sulphur in a soil during decomposition of plant residues of varied chemical composition and S content. Plant Soil. 2012;360:391–404.CAS 
Article 

Google Scholar 
44.Vong PC, Piutti S, Slezackdeschaumes S, Benizri E, Guckert A. Effects of low-molecular-weight organic compounds on sulphur immobilization and re-mineralization and extraction of immobilized sulphur by hot-water and acid hydrolysis. Eur J Soil Sci. 2010;61:287–97.CAS 
Article 

Google Scholar 
45.Takriti M, Wild B, Schnecker J, Mooshammer M, Knoltsch A, Lashchinskiy N, et al. Soil organic matter quality exerts a stronger control than stoichiometry on microbial substrate use efficiency along a latitudinal transect. Soil Biol Biochem. 2018;121:212–20.CAS 
Article 

Google Scholar 
46.Zhu Z, Ge T, Luo Y, Liu S, Xu X, Tong C, et al. Microbial stoichiometric flexibility regulates rice straw mineralization and its priming effect in paddy soil. Soil Biol Biochem. 2018;121:67–76.CAS 
Article 

Google Scholar 
47.Xu X, Hui D, King AW, Song X, Thornton PE, Zhang L. Convergence of microbial assimilations of soil carbon, nitrogen, phosphorus and sulfur in terrestrial ecosystems. Sci Rep-UK. 2015;5:17445.CAS 
Article 

Google Scholar 
48.Fanin N, Fromin N, Buatois B, Ttenschwiler SH. An experimental test of the hypothesis of non-homeostatic consumer stoichiometry in a plant litter-microbe system. Ecol Lett. 2013;16:764–72.PubMed 
Article 
PubMed Central 

Google Scholar 
49.Peters MK, Hemp A, Appelhans T, Becker JN, Behler C, Classen A, et al. Climate–land-use interactions shape tropical mountain biodiversity and ecosystem functions. Nature. 2019;568:88–92.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
50.Kaiser C, Franklin O, Dieckmann U, Richter A. Microbial community dynamics alleviate stoichiometric constraints during litter decay. Ecol Lett. 2014;17:680–90.PubMed 
PubMed Central 
Article 

Google Scholar 
51.Xu X, Schimel JP, Thornton PE, Song X, Yuan F, Goswami S. Substrate and environmental controls on microbial assimilation of soil organic carbon: a framework for Earth system models. Ecol Lett. 2014;17:547–55.PubMed 
Article 
PubMed Central 

Google Scholar 
52.Apostel C, Dippold M, Kuzyakov Y. Biochemistry of hexose and pentose transformations in soil analyzed by position-specific labeling and 13C-PLFA. Soil Biol Biochem. 2015;80:199–208.CAS 
Article 

Google Scholar 
53.Manzoni S, Taylor P, Richter A, Porporato A, Ågren GI. Environmental and stoichiometric controls on microbial carbon‐use efficiency in soils. N Phytol. 2012;196:79–91.CAS 
Article 

Google Scholar 
54.Fitzgerald JW, Hale DD, Swank WT. Sulphur-containing amino acid metabolism in surface horizons of a hardwood forest. Soil Biol Biochem. 1988;20:825–31.CAS 
Article 

Google Scholar 
55.Akashi H, Gojobori T. Metabolic efficiency and amino acid composition in the proteomes of Escherichia coli and Bacillus subtilis. P Natl Acad Sci USA. 2002;99:3695–700.CAS 
Article 

Google Scholar 
56.Nozaki T, Ali V, Tokoro M. Sulfur-Containing Amino Acid Metabolism in Parasitic Protozoa. Adv Parasit. 2005;60:1–99.Article 

Google Scholar 
57.Takagi H, Ohtsu I. l-Cysteine Metabolism and Fermentation in Microorganisms. Adv Biochem Eng Biotechnol. 2016;159:129–51.
Google Scholar 
58.Bustos I, Miguel AM, Fouad A, Carmen P, Teresa R, CM M. Volatile sulphur compounds-forming abilities of lactic acid bacteria: C-S lyase activities. Int J Food Microbiol. 2011;148:121–7.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
59.Assefa MK, Tucher SV, Schmidhalter U. Soil sulfur availability due to mineralization: soil amended with biogas residues. J Soil Sci Enviro Manag. 2014;5:13–9.CAS 
Article 

Google Scholar  More

Characteristics of the single agrometeorological disaster scenariosSAD-f occurred in 49 out of 57 years at different spatial scales, with a maximum IOC of 0.519 in 2013; SAD-d occurred in 33 years with a maximum IOC value of 0.808 in 1995; SAD-s occurred in 5 years, with a maximum IOC value of 0.115 in 1977 (Fig. 4). SAD-d showed a declining trend over the past 57 years, but the SAD-d frequency was higher than SAD-f and SAD-d. Since the mid-1980s, the frequency of SAD-f has increased, while the frequency and scale of SAD-s were relatively small.Figure 4IOC change curve for single agrometeorological disasters (SAD) in different scenarios.Full size imageA large-scale grade SAD-f event occurred in 2013 in Liaoning Province, and regional SAD-f occurred in 17 years. Five years showed a large-scale SAD-d and 11 years had regional SAD-d. There were no large-scale and regional years for SAD-s (Table 5).Table 5 Occurrence years of large-scale and regional single agrometeorological disaster (SAD) in Liaoning Province.Full size tableThe occurrence of SAD-f was recorded at 15 sites with a frequency greater than 20%, 13 sites with frequency in the 10%–20% range, and 24 sites with a low frequency (P ≤ 10%) in Liaoning Province from 1961 to 2017 (Fig. 5). SAD-d occurred at 12 sites with a frequency higher than 20%, 22 sites with frequency between 10% and 20%, and 18 sites with a low frequency (P ≤ 10%). In the three scenarios, the occurrence frequency and distribution of SAD-f was the highest and SAD-s was the lowest.Figure 5Frequency of single agrometeorological disasters (SAD) in different scenarios (a SAD-f, b SAD-d, c SAD-s). Maps generated in ArcGIS 9.3.Full size imageComparison of the characteristics of single agrometeorological disasters and combined agrometeorological disastersThe maximum IOC of SAD was 0.808 in 1995 and the mean value was 0.294 for the 57 years of the study; the maximum IOC of CAD was 0.654 in 1987, and the mean value was 0.180 over the past 57 years; SAD and CAD occurred in all 57 years (Fig. 6). Both SAD and CAD showed declining trends from 1961 to 2017. The IOC was lower for CAD than for SAD for 42 years and higher than SAD for 14 years.Figure 6Change in the IOC for agrometeorological disasters in rice crops.Full size imageThis paper analysed the mean IOC of SAD and CAD over six decades and found that the interdecadal mean value of the IOC in CAD was lower than that of SAD over five of the periods, but the IOC of SAD was lower than that of CAD in 1971–1980 (Fig. 7). The IOC of SAD showed a decreasing trend from the 1970s to the 2010s but showed an increasing trend after 2011. The IOC of CAD showed a decreasing trend from the 1970s to the 2000s, but showed an increasing trend after 2001 (Fig. 7).Figure 7Interdecadal mean value of IOC for agrometeorological disasters in rice crops.Full size imageThere was one site (Fushun) with a SAD frequency of more than 50% in Liaoning Province from 1961 to 2017, 42 sites with a frequency between 20% and 50%, and nine sites with a frequency lower than 20%. There were four sites (Xinfeng, Jianping, Xinbin and Caohekou) with a CAD frequency higher than 50%, 13 sites with a frequency in the range 20%–50%, and 35 sites with a frequency lower than 20% (Fig. 8). The frequency and range of CAD were less than those of SAD.Figure 8Frequency of agrometeorological disaster in rice crops (a combined agrometeorological disaster (CAD), b single agrometeorological disaster (SAD)). Maps generated in ArcGIS 9.3.Full size imageThere has been little research into the temporal or spatial distribution of CAD for rice and its occurrence characteristics: most research has been on SAD. For example, studies have examined the characteristics of SCD, DCD, FD for rice in northeast China26, 27, 29, and the risk of multiple disasters for rice in northeast China30, 31. Han et al.31 analysed the risk of disaster using the reduction rate of rice yield in Liaoning Province from 1980 to 2011, and found that the high-risk areas were distributed in the west and northeast of Liaoning Province; higher rates of yield reduction in lean years were mainly found in western Liaoning and its surrounding areas. In this study, a higher frequency of CAD was mainly distributed in the northwest of Liaoning Province, while that of SAD occurred in the northeast of Liaoning Province. The median frequency of CAD occurred in the northwest and northeast of Liaoning Province, while that of SAD covered most areas in Liaoning Province. The range of medium and higher frequency occurrence in CAD was consistent with the distribution of high-risk and high yield reduction areas in the study of Han et al.31. Therefore, it can be speculated that the CAD scenarios might magnify the effect of each single disaster, and, therefore, CAD would more easily lead to a higher reduction in the rice yield.Comparison of the occurrence of single agrometeorological disasters and combined agrometeorological disastersDuring the rice growing season in Liaoning Province, there were three scenarios of SAD and six scenarios of CAD. Compared with SAD, CAD had more scenarios and more complex processes, and its effect on rice was more difficult to evaluate. In SAD, the occurrence frequency and distribution of SAD-f and SAD-d were both high, when FD and DCD occurred alone in only one rice growth stage. In CAD, the occurrence frequency and distribution of TD-1, when FD and DCD occurred simultaneously, was the highest in the six scenarios. A single or combined occurrence of FD and DCD was most common disaster for rice in Liaoning Province. The occurrence frequency and distribution of OD-1 were both smaller than that of SAD-f, indicating that the occurrence was lower when FD happened at both the seedling and milk stages. SAD-s and OD-2 had the lowest frequency and range in all scenarios, indicating that DSD rarely appeared in SAD and CAD. The occurrence of SCD was not major disaster in the growth and development of rice in Liaoning Province, but the occurrence of DCD or FD, or both, was.In this study, the occurrence frequency and range of SAD and CAD for rice showed declining trends in most sites over the past 57 years, which was consistent with the results of other studies. Studies on rice DCD and SCD concluded that cold damage events of rice in most areas of northeast China showed decreasing trends26, 27. Because of events such as climate warming, earlier warming in spring, delaying first frost dates and fewer low temperature days in summer, the trend of disasters was lower in rice planting areas30. However, although rice disasters showed a decreasing trend, local disasters may increase because of the frequent occurrence of climate anomalies. SAD-f and OD-1 scenarios in this study showed no significant decreasing trend, and even a partial increasing trend. Jiang et al.29 believed that the possibility of frequent SCD in north-east China was still high. According to Xi et al.32, cold periods would still occur in the growing season of rice in northeast China. Hu et al.33 concluded that the increase of SCD in northeast China was mainly because of the increase of climate variability, and most of the sites with increases were located in areas with decreasing temperature or no obvious trend of temperature increase.Rice is a higher temperature-loving crop, which is mainly restricted by temperature conditions during its growing season. Liaoning Province is in the south of the rice planting area of the colder regions in China. Because of the relatively low latitude, heat conditions during the rice growing season were better than those in Jilin and Heilongjiang to the north of Liaoning Province. The climatic risk of cold damage in the rice growing season was lower than other regions in northeast China34. The occurrence of CAD was generally caused by low temperatures, which were the dominant factor. When two or more disasters occur together, there is a coupling or amplifying effect on rice growth compared with a single disaster.A comparison of the rice yield reduction rates in the years when CAD or SAD occurred in more than 50% of stations in Liaoning Province revealed that the former happened in 5 years, 1969, 1974, 1976, 1980 and 1987, whereas the latter happened in 7 years, 1972, 1973, 1985, 1986, 1990, 1995 and 2013. When CAD was the major occurrence, the average yield reduction rate in the five years was 10.6%. The yield reduction rate in 1969 was 34.6%, which was the highest in the past 57 years. When SAD was the major occurrence, the average yield reduction rate in the seven years was 9.8%. The average yield reduction rate in the years when CAD dominated was greater than in the years when SAD dominated. Therefore, it can be speculated that CAD has a greater effect on rice growing than any single disaster within CAD. However, it is difficult to quantify the effect on rice yield of CAD, and further controlled field experiments should be conducted to verify these. It is difficult to control field experiments that are limited by conditions and facilities.Comparison of the occurrence of agrometeorological disasters in years having rice yield reductionsOn the basis of the rice yield reduction rate in calculations Liaoning Province from 1961 to 2017, a total of 10 years (Table 6) were screened. Six years had large-scale disasters (including SAD and CAD) and four years had regional disasters. In 1969, which showed the highest yield reduction rate (34.6%), 30 sites had TD-1 disasters and the other 22 sites had SAD-f disasters. In 1972, the second highest reduction year (29.1%), 11 sites had MD-1 disasters, i.e. three kinds of disasters occurred, seven sites had TD-1 disasters, one site had a TD-2 disaster, 31 sites had SAD-d disasters, one site had a SAD-f disaster, and only one station had no disaster. The TD-1 disaster, i.e. delayed cold damage and frost injury, was the most frequent CAD over the years, and SAD-d, i.e., delayed cold damage, was the most frequent SAD. The occurrence of single and combined agrometeorological disasters in different regions strongly affected the rice yield. Generally, the larger the disaster range, the higher the yield reduction. However, some years were not completely consistent with this conclusion. The yield reduction rate was also related to the type, severity, occurrence period and geographical location of the disasters.Table 6 Comparison of agrometeorological disasters in years having greater than 10% rice yield reduction rates in Liaoning Province.Full size tableIn every year from 1961 to 2017, CAD or SAD occurred in Liaoning Province, and the rice yields declined in 23 of the 57 years owing to meteorological disasters (Fig. 9). Although meteorological disasters occurred in the other 34 years, there was no reduction in rice production, which may be related to the gradient of the disaster or the spatial distribution of the rice planting areas. The rice yield reduction rates in 1969 and 1976 were 34.6% and 15.6%, respectively. In these two years, CAD occurred at 30 stations and SAD occurred at 22 stations, and TD-1 was the main type of CAD, whereas SAD-d was the main type of SAD. Using statistical data, on the rice planting area of each city in Liaoning Province, the provincial area can be divided into four regions. The first region was Shenyang City, which has the largest rice planting area, accounting for 20%–25% of the total rice planting area; the second region was Panjin City, accounting for 15%–20% of the total rice area; the third region encompassed Tieling and other six cities, accounting for nearly 50% of the total rice area, with each city representing 5%–10%; and the fourth region encompassed Jinzhou and five other cities, accounting for 10%–15% of the total, with each city representing 0–5%. As shown in Fig. 10a,b, TD-1 occurred in the first region in both 1969 and 1976 and in the second region in 1969. SAD-d occurred in the second region in 1976. In the third region, TD-1 occurred at more stations of 1969 than in 1976. The rice area in the first three regions accounted for nearly 80% of the total rice area, and CAD occurred more often than SAD in these regions. Thus, there was a greater yield reduction rate in 1969 than in 1976.Figure 9The IOC change curve of all agrometeorological disasters and the rice yield reduction rate from 1961 to 2017.Full size imageFigure 10Distributions of the types of agrometeorological disasters and the percentages of rice planting areas in different regions of Liaoning Province in 1969 and 1976 (a: 1969; b: 1976). Maps generated in ArcGIS 9.3.Full size imageThe occurrence characteristics of single disasters or the risk of yield reduction were analysed in previous studies, but the quantitative effect on rice production was rarely evaluated. Ji et al.26 reported that the delayed cold damage in 1961, 1962, 1969, 1972, 1976, 1989 and 1995 was so severe that there was a large reduction in rice production. In our paper, we examined the occurrence of not just one disaster, i.e. delayed cold damage, over time, but also other types of disasters including SAD and CAD. For example, in 1972 and 1976, the disaster scenario affecting the largest number of stations was TD-1, i.e., both delayed cold damage and frost damage occurred in the growing season of rice. In 1961, the most widespread damage came from a single disaster—frost damage. According to the records35, Liaoning Province experienced frost damage in 1961, 1962, 1969, 1972, 1976 and 1995, and the rice yield was seriously reduced. Most regions of Liaoning Province experienced both delayed cold damage and frost damage in 1976 and 1995. There was a low temperature during the critical period of rice growth (mid-July to mid-August) in 1995. In 1985, the growing season in most areas was characterized by unusually persistent low temperature and little sunshine. These statistics were basically consistent with the conclusion of this study. In the process of rice production, a variety of disasters occurred caused by low temperature, such as delayed cold damage, frost damage and sterile cold damage. More

Fish and experimental protocolThree-hundred-fifty market-size (755–3865 g) common carp Cyprinus carpio were obtained from six sources at various times of year to for effects of variation in rearing conditions. The weight of collected carp corresponded to the weight of carp normally delivered to the market. Fish were obtained from the Faculty of Fisheries and Protection of Waters of the University of South Bohemia in Ceske Budejovice (FFPW USB), Vodnany and the fisheries Chlumec nad Cidlinou, Blatna, Hodonin, Klatovy, Lnare, and Tabor. Ten fish were collected from each fishery at the spring (March/April), summer (June/July), and autumn harvests (October/November) in 2018 and 2019. Carp were transported live to the laboratory of the FFPW, killed by a blow to the head, weighed, measured, and filleted. Two fillets, one with skin removed, from each fish were individually vacuum packed, immediately frozen, and stored at − 32 °C until chemical analysis.Ethics approvalAll the methods used in the present study followed relevant guidelines and regulations. Also, the competent authority (Ethical Committee for the Protection of Animals in Research of the University of South Bohemia, FFPW Vodnany) approved the fish sampling and protocols of the present study and reporting herein follows the recommendations in the ARRIVE guidelines.Chemical analysisSeven-hundred carp fillets were analysed for basic nutritional composition, dry matter, protein, fat, and ash. All samples were homogenized by grinding before analysis.The determination of dry matter followed ISO 1442:1997 Meat and meat products—Determination of moisture content (Reference method)26. The homogenized samples were dried with sand to constant weight at 103 ± 2 °C in a laboratory oven (Memmert UE 500, Memmert GmbH + Co. KG, Germany).The determination of ash was based on the standard ISO 936:1998 Meat and meat products—Determination of total ash27. The homogenized samples were burned in a muffle furnace (Nabertherm A11/HR, Nabertherm GmbH, Germany) at 550 ± 25 °C to a grey-white colour.The determination of total fat was based on the standard ISO 1443:1973 Meat and meat products—Determination of total fat content28. The homogenized samples were hydrolysed by hydrochloric acid, and fat was extracted by light petroleum in SOXTEC 2050 (FOSS Headquarters, Denmark).The determination of nitrogen used the Kjeldahl method based on the standard method ISO 937:1978 Meat and meat products—Determination of nitrogen content (Reference method)29. The homogenized samples were digested by sulphuric acid and a catalyser in a KjelROC Digestor 20 (OPSIS AB, Sweden) digestion unit at 420 ± 10 °C. Organically bound nitrogen was measured on the KJELTEC 8400 with KJELTEC sampler 8420 (FOSS Headquarters, Denmark). Calculation of protein content from nitrogen used the conversion factor for meat of 6.25.All analysis of dry matter, ash, and total fat were performed in duplicate and analysis of nitrogen (protein) was performed in triplicate for each sample.Calculation of fat-free nitrogen (Nff) in g/100 g used the formula24:$$ N_{ff} = frac{{100 times N { }}}{{100 – F { }}}. $$This formula was applied to nitrogen (N) and fat (F) content for all samples, providing a fat-free nitrogen value for each sample.Fish meat content calculated based on nitrogen factor Nf (total fillet) in g/100 g used the formula9:$$ Fish ;content_{Nf} = frac{N times 100}{{N_{f} }}. $$Fish meat content calculated based on fat-free nitrogen factor (Nff) and DCC (defatted carp content) in g/100 g used formulas11:$$ Fishc; content_{Nff} = DCC + F, $$$$ DCC = frac{N times 100}{{N_{ff} }}. $$Statistical analysisKolmogorov–Smirnov and Bartlett’s tests were applied to assess normal distribution data and the homoscedasticity of variance, respectively. A two-way ANOVA and Tukey’s test was conducted to analyse effects of season, weight, fishery, and difference between fillets with and without skin. The significance level was set at p  More

Uganda has been exceptionally successful in scaling-up coverage of LLINs. Following the mass distribution campaigns to deliver free LLINs in 2013–14 and 2017–18, 90 and 83% of households, respectively reported ownership of at least one LLIN7,14. However, despite this success, the burden of malaria remains high in much of the country. Uganda had the third highest number of malaria cases reported in 2019, with reported case incidence increasing since 20142. If Uganda is to achieve the goals established by the World Health Organization’s Global Technical Strategy for malaria including reducing malaria case incidence by at least 90% by 2030 as compared with 201515, additional malaria control measures will be needed. This report highlights the critical role of IRS in substantially reducing the burden of malaria in areas where transmission remains high despite deployment of LLINs. Withdrawing IRS after 5 years of sustained use in three districts in northern Uganda was associated with a more than fivefold increase in malaria cases within 10 months. Restarting IRS with a single round in nine districts in Northern Uganda ~3 years after IRS had been stopped was associated with a transient but important (more than a fivefold) decrease in malaria cases within 8–12 months, returning to pre-IRS levels after 23 months. Initiating and sustaining IRS in five districts in Eastern Uganda was associated with a gradual reduction in malaria cases reaching almost a sevenfold reduction after 4–5 years.Robust evidence supports the widespread use of LLINs for malaria control. In a systematic review of clinical trials conducted between 1987 and 2001, insecticide treated nets reduced all cause child mortality by 17% and the incidence of uncomplicated P. falciparum malaria by almost half16. However, there is concern that the effectiveness of LLINs may be diminishing due to widespread resistance to pyrethroids which until recently were the only class of insecticides approved for LLINs. Similar to many other African countries, high-level resistance to pyrethroids among the principle Anopheles vectors has been reported recently throughout Uganda17,18,19. In addition, behavioral changes in vector biting activity following the introduction of LLINs have been reported which could present new challenges for malaria control20,21,22. Finally, the effectiveness of LLINs may be further compromised by poor adherence and waning coverage in the setting of free distribution campaigns done intermittently. In Uganda, less than 18% of households reported adequate coverage (defined as at least one LLIN per two residents) 3 years after the 2013–14 distribution campaign23 and adequate coverage decreased from 71% to 51% between 6 and 18 months following the 2017–18 distribution campaign24. Although the World Health Organization recommends mass distribution campaigns every 3 years, mounting evidence suggests that LLINs should be distributed more frequently to sustain high coverage25,26,27,28,29,30,31.Given concerns about the current effectiveness of pyrethroid-based LLINs and the persistently high burden of malaria despite aggressive scale-up of LLINs in countries like Uganda, additional malaria control measures are needed. IRS is an attractive option. Historically, IRS programs were used to dramatically reduce and even eliminate malaria in many parts of the world. Thus, while there is some evidence for the impact of IRS in the absence of LLINs32, it is surprising that the evidence base from contemporary controlled trials on the impact of adding IRS to LLINs for vector control is limited. A recent systematic review of cluster randomized controlled trials conducted in sub-Saharan Africa since 2008, reported that adding IRS using a “pyrethroid-like” insecticide to LLINs did not provide any benefits, while adding IRS with a “non-pyrethroid-like” insecticide produced mixed results5. Among the four trials comparing IRS plus LLINs with LLINs alone, three evaluated IRS with a carbamate (bendiocarb) and one evaluated a long-lasting organophosphate, pirimiphos-methyl (Actellic 300CS®)33,34,35,36. Only two trials (both using bendiocarb) assessed malaria incidence; one from Sudan found a 35% reduction when adding IRS to LLINs34, while another from Benin found no benefit of adding IRS33. All four trials assessed parasite prevalence, with an overall non-significant trend towards a lower prevalence when adding IRS to LLINs (RR = 0.67, 95% CI 0.35–1.28)5. However, when the analyses were restricted to include only the two studies with LLIN usage over 50%, adding IRS reduced parasite prevalence by over 50% (RR = 0.47, 95% CI 0.33–0.67)5. Of note, none of the trials that evaluated the impact of adding IRS with a “non-pyrethroid-like” insecticide assessed outcomes beyond 2 years. More recently, a number of observational studies have reported benefits of using IRS with pirimiphos-methyl (Actellic 300CS®). In the Mopti Region of Mali, delivery of a single round of IRS with Actellic 300CS® was associated with a 42% decrease in the peak incidence of laboratory-confirmed malaria cases reported at public health facilities37. In the Koulikoro Region of Mali, villages that received a single round of IRS with Actellic 300CS® combined with LLINs observed a greater than 50% decrease in the incidence of malaria compared to villages that only received LLINs38. In the Northern Region of Ghana, districts that received IRS with Actellic 300CS® reported 26–58% fewer cases of laboratory-confirmed malaria cases reported at public health facilities over a 2-year period, compared to districts that did not receive IRS39. In Northern Zambia, implementation of IRS with Actellic 300CS® targeting only high burden areas over a 3 year period was associated with a 25% decline in parasite prevalence during the rainy season, but no decline during the dry season40. In Western Kenya, the introduction of a single round of IRS with Actellic 300CS® was associated with a 44–65% decrease in district level malaria case counts over a 10 month period compared to pre-IRS levels41. In addition, several recent reports have documented dramatic resurgences of malaria following the withdrawal of IRS with bendiocarb in Benin42, and the withdrawal of IRS with Actellic 300CS® in Mali and Ghana37,39.The results from this study provides additional support for the critical role IRS can play in reducing the burden of malaria in African countries with high LLINs coverage. A strength of the study was its use of a large, rigorously collected dataset. Data were collected over nearly 7 years through an enhanced health facility-based surveillance system covering 14 districts in Uganda where IRS was being withdrawn, re-started, and initiated. This enhanced surveillance system facilitated laboratory testing and provided prospectively collected, individual-level data, allowing for analyses of quantitative changes in laboratory-confirmed cases of malaria over time, controlling for temporal changes in rainfall, seasonal effects, diagnostic practices, and health seeking behavior. Previous work by our group documented a marked decrease in malaria TPRs after 4 years of sustained IRS with bendiocarb in one district of Northern Uganda followed by a rapid resurgence over an 18-month period after IRS was withdrawn11. In this study we expand on these findings by including data from three districts and covering a 31-month period following the withdrawal of IRS. We were able to quantify more than a fivefold increase in malaria cases which was sustained over the 10–31 months following the withdrawal of IRS. This marked resurgence occurred despite the fact the first universal LLIN distribution campaign was timed to occur right after IRS was withdrawn. Given the dramatic nature of the resurgence, the Ugandan government was able to procure funding for a single round of IRS with Actellic 300CS® ~3 years after IRS was withdrawn in 10 districts of Northern Uganda. In this study, we assessed the impact of this single round in nine of these districts. This single round was associated with over a fivefold decrease in malaria cases after 8–12 months, with malaria cases returning to pre-IRS levels after almost 2 years. These data suggest that IRS with longer-acting formulations such as Actellic 300CS® administered every 2 years could be considered as a strategy for mitigating the risk of resurgence following sustained IRS and/or enabling countries to expand coverage when resources are limited, but formal assessment and a cost-effectiveness analyses are needed. This study also evaluated the impact of 5 years of sustained IRS in five districts of Eastern Uganda, starting first with bendiocarb and then switching to Actellic 300CS® after 18 months. Rounds of IRS were initially associated with marked decreases in malaria cases followed by peaks before subsequent rounds until the fourth and fifth years after IRS was initiated when there was a sustained decrease of almost sevenfold compared to pre-IRS level. Given the before-and-after nature of our study design, it is not clear whether the maximum sustained benefits of IRS seen after 4–5 years were due to the cumulative effect of multiple rounds of IRS, the switch from bendiocarb to Actellic 300CS®, improvements in implementation (although campaigns occurred regularly and coverage was universally high across rounds, see Supplementary Table 4), the second universal LLIN distribution campaign which occurred in this area in 2017, and/or other factors.This study had several limitations. First, we used an observational study design, with measures of impact based on comparisons made before-and-after key changes in IRS policy. Although cluster randomized controlled trials are the gold standard study design for estimating the impact of IRS, it could be argued that withholding IRS would be unethical, given what is known about its impact in Uganda. Second, our estimates of impact could have been confounded by secular trends in factors not accounted for in our analyses. However, we feel that our overall conclusions are robust given the large amount of data available from multiple sites over an extended period with multiple complementary objectives providing consistent findings. Third, we could not assess the impact of IRS independent of LLIN use and did not have access to measures of IRS or LLIN coverage from our study populations. It is possible that some of the impacts we observed were from LLIN distributions in combination with IRS campaigns. However, we were able to provide a “real world” assessment of IRS in a setting where LLIN use is strongly supported by repeated universal distribution campaigns that are becoming increasingly common in sub-Saharan Africa. Similarly, we cannot draw conclusions on the impact of different IRS compounds given all sites received the same formulations consecutively. The results from Objective 3 indicate that malaria incidence dropped substantially in the years that districts stopped receiving bendiocarb and began receiving Actellic 300CS®. However, we cannot conclude whether this reduction was a result of this change or rather the cumulative impact of sustained IRS campaigns, as it has been suggested that in very high transmission settings, several years of IRS may be needed to maximize impact on measures of morbidity.43,44 Finally, our study outcome was limited to case counts of laboratory-confirmed malaria captured at health facilities. Thus, we were unable to measure the impact of IRS on other important indicators such as measures of vector distribution, parasite prevalence, or mortality.There is a growing body of evidence that combining LLINs with IRS using “non-pyrethroid-like” insecticides, especially the long acting organophosphate Actellic 300CS®, is highly effective at reducing the burden of malaria in Uganda, and elsewhere in Africa. Despite these encouraging findings, IRS coverage in Africa has been moving in the wrong direction. The proportion of those at risk protected by IRS in Africa peaked at just over 10% in 2010. However, the spread of pyrethroid resistance has led many control programs to switch to more expensive formulations resulting in a 53% decrease in the number of houses sprayed between years of peak coverage and 2015 across 18 countries supported by the US President’s Malaria Initiative45 and an overall reduction in the proportion protected by IRS in Africa to less than 2% in 20192. Given the lack of recent progress in reducing the global burden of malaria coupled with challenges in funding, renewed commitments are needed to address the “high burden to high impact” approach now being advocated by the World Health Organization2. IRS is a widely available tool that could be scaled up, however demands currently exceed the availability of resources. Additional work is needed to optimize the use of IRS, prevent further spread of insecticide resistance, and better evaluate the cost-effectiveness of IRS in the context of other control interventions. More