Ecology

Allendorf, F. W. & Lundquist, L. L. Introduction: Population biology, evolution, and control of invasive species. Conserv. Biol. 17, 24–30 (2003).Article 

Google Scholar 
Kopf, R. K. et al. Confronting the risks of large-scale invasive species control. Nat. Ecol. Evol. 1, 0172 (2017).Article 

Google Scholar 
Luque, G. M. et al. The 100th of the world’s worst invasive alien species. Biol. Invasions 16, 981–985 (2014).Article 

Google Scholar 
Genovesi, P. Eradications of invasive alien species in Europe: A review. Biol. Invasions 7, 127–133 (2005).Article 

Google Scholar 
Simberloff, D. How much information on population biology is needed to manage introduced species?. Conserv. Biol. 17, 83–92 (2003).Article 

Google Scholar 
Takeshita, F. & Maekawa, T. Paratectonatica tigrina (Gastropoda: Naticidae) adjusts its predation tactics depending on the chosen prey and their shell weight relative to its own. J. Mar. Biol. Assoc. U. K. 100, 921–926 (2020).Article 

Google Scholar 
Wiltse, W. I. Effects of Polinices duplicatus (Gastropoda: Naticidae) on infaunal community structure at Barnstable Harbor, Massachusetts, USA. Mar. Biol. 56, 301–310 (1980).Article 

Google Scholar 
Commito, J. A. Effects of Lunatia heros predation on the population dynamics of Mya arenaria and Macoma balthica in Maine, USA. Mar. Biol. 69, 187–193 (1982).Article 

Google Scholar 
Ansell, A. D. Experimental studies of a benthic predator–prey relationship. I. Feeding, growth, and egg-collar production in long-term cultures of the gastropod drill Polinices alderi (Forbes) feeding on the bivalve Tellina tenuis (da Costa). J. Exp. Mar. Biol. Ecol. 56, 235–255 (1982).Article 

Google Scholar 
Savazzi, E. & Reyment, R. A. Subaerial hunting behaviour in Natica gualteriana (naticid gastropod). Palaeogeogr. Palaeoclimatol. Palaeoecol. 74, 355–364 (1989).Article 

Google Scholar 
Gonor, J. J. Predator–prey reactions between two marine prosobranch gastropods. The Veliger 7, 228–232 (1969).
Google Scholar 
Hughes, R. N. Predatory behaviour of Natica unifasciata feeding intertidally on gastropods. J. Molluscan Stud. 51, 331–335 (1985).Article 

Google Scholar 
Pahari, A. et al. Subaerial naticid gastropod drilling predation by Natica tigrina on the intertidal molluscan community of Chandipur, Eastern Coast of India. Palaeogeogr. Palaeoclimatol. Palaeoecol. 451, 110–123 (2016).Article 

Google Scholar 
Sakai, K. Predation of the moon snail Neverita didyma, on the Manila clam Ruditapes philippinarum, at the culture ground in Mangoku-ura Inlet. Bull. Miyagi Prefect. Fish. Res. Dev. Cent. 16, 109–111 (2000) (in Japanese).
Google Scholar 
Tomiyama, T. et al. Unintentional introduction and the distribution of the nonindigenous moonsnail Euspira fortunei in Matsukawaura Lagoon, Japan. Nippon Suisan Gakkaishi 77, 1020–1026 (2011) (in Japanese with English abstract).Article 

Google Scholar 
Okoshi, K. Alien species introduced with imported clams: The clam-eating moon snail Euspira fortunei and other unintentionally introduced species. Japanese J. Benthol. 59, 74–82 (2004) (in Japanese with English abstract).
Google Scholar 
Okoshi, K. & Sato-Okoshi, W. Euspira fortunei: Biology and fisheries science of an invasive species (Kouseishakouseikaku Press, 2011).
Google Scholar 
Tomiyama, T. Lethal and non-lethal effects of an invasive naticid gastropod on the production of a native clam. Biol. Invasions 20, 2005–2014 (2018).Article 

Google Scholar 
Sato, S., Chiba, T. & Hasegawa, H. Long-term fluctuations in mollusk populations before and after the appearance of the alien predator Euspira fortunei on the Tona coast, Miyagi Prefecture, northern Japan. Fish. Sci. 78, 589–595 (2012).CAS 
Article 

Google Scholar 
Kinoshita, K., Sasaki, N., Seki, A., Matsumasa, M. & Takehara, A. Distribution of the invasive snail Laguncula pulchella and the effect of its predation on the mollusk populations in the Orikasa River Estuary after the 2011 off the Pacific Coast of Tohoku Earthquake. Japanese J. Benthol. 72, 61–70 (2018) (in Japanese with English abstract).
Google Scholar 
Sato, T., Iwasaki, T., Narita, K. & Matsumoto, I. Occurrence of the moonsnail Euspira fortunei after the earthquake in Matsukawaura Lagoon Japan. Bull. Fukushima Prefect. Fish. Exp. Stn. 79–82 (2016). (in Japanese).Chiba, T. & Sato, S. Size-selective predation and drillhole-site selectivity in Euspira fortunei (Gastropoda: Naticidae): Implications for ecological and palaeoecological studies. J. Molluscan Stud. 78, 205–212 (2012).Article 

Google Scholar 
Tanabe, T. Relationship between the shell height of the predatory moon snail Euspira fortunei and drilled hole diameter on the prey shell of manila clam Ruditapes philippinarum. Nippon Suisan Gakkaishi 78, 37–42 (2012) (in Japanese with English abstract).Article 

Google Scholar 
Hasegawa, H. & Sato, S. Predatory behaviour of the naticid Euspira fortunei: Why does it drill the left shell valve of Ruditapes philippinarum?. J. Molluscan Stud. 75, 147–151 (2009).Article 

Google Scholar 
Kingsley-Smith, P. R., Richardson, C. A. & Seed, R. Size-related and seasonal patterns of egg collar production in Polinices pulchellus (Gastropoda: Naticidae) Risso 1826. J. Exp. Mar. Biol. Ecol. 295, 191–206 (2003).Article 

Google Scholar 
Tomiyama, T. Timing and frequency of egg-collar production of the moonsnail Euspira fortunei. Fish. Sci. 79, 905–910 (2013).CAS 
Article 

Google Scholar 
Sakai, K. & Suto, A. Early development and behavior of the moon snail Neverita didyma. Miyagi Prefect. Rep. Fish. Sci. 5, 55–58 (2005) (in Japanese).
Google Scholar 
Cook, N. & Bendell-Young, L. Determining the ecological role of Euspira lewisii: Part I: Feeding ecology. J. Shellfish Res. 29, 223–232 (2010).Article 

Google Scholar 
Peitso, E., Hui, E., Hartwick, B. & Bourne, N. Predation by the naticid gastropod Polinices lewisii (Gould) on littleneck clams Protothaca staminea (Conrad) in British Columbia. Can. J. Zool. 72, 319–325 (1994).Article 

Google Scholar 
Johannesson, K., Saltin, S. H., Duranovic, I., Havenhand, J. N. & Jonsson, P. R. Indiscriminate males: Mating behaviour of a marine snail compromised by a sexual conflict?. PLoS ONE 5, e12005 (2010).ADS 
Article 

Google Scholar 
Vaughn, D., Turnross, O. R. & Carrington, E. Sex-specific temperature dependence of foraging and growth of intertidal snails. Mar. Biol. 161, 75–87 (2014).Article 

Google Scholar 
Watson, P. J., Arnqvist, G. & Stallmann, R. R. Sexual conflict and the energetic costs of mating and mate choice in water striders. Am. Nat. 151, 46–58 (1998).CAS 
Article 

Google Scholar 
Schlacher, T. A. & Wooldridge, T. H. Patterns of selective predation by juvenile, benthivorous fish on estuarine macrofauna. Mar. Biol. 125, 241–247 (1996).Article 

Google Scholar 
Sudo, H. & Azeta, M. Selective predation on mature male Byblis japonicus (Amphipoda: Gammaridea) by the barface cardinalfish, Apogon semilineatus. Mar. Biol. 114, 211–217 (1992).Article 

Google Scholar 
Vadas, R. L., Burrows, M. T. & Hughes, R. N. Foraging strategies of dogwhelks, Nucella lapillus (L.): Interacting effects of age, diet and chemical cues to the threat of predation. Oecologia 100, 439–450 (1994).ADS 
Article 

Google Scholar 
Donelan, S. C. & Trussell, G. C. Sex-specific differences in the response of prey to predation risk. Funct. Ecol. 34, 1235–1243 (2020).Article 

Google Scholar 
Yoshida, K., Sato, T., Narita, K. & Tomiyama, T. Abundance and body size of the moonsnail Laguncula pulchella in the Misuji River estuary, Seto Inland Sea, Japan: Comparison with a population in northern Japan. Plankt. Benthos Res. 12, 53–60 (2017).Article 

Google Scholar 
Sato, T., Ogata, Y., Nemoto, Y. & Shimamura, S. Status review and current concerns of the fishery for the short-neck clam, Ruditapes philippinarum, in Matsukawaura Lagoon, Fukushima Prefecture. Bull. Fukushima Prefect. Fish. Exp. Stn. 14, 57–67 (2007) (in Japanese).
Google Scholar 
Tomiyama, T. & Sato, T. Effects of translocation on the asari clam Ruditapes philippinarum at small spatial scales in Matsukawaura, Japan. Bull. Mar. Sci. 97, 647–664 (2021).Article 

Google Scholar 
Chiba, T. & Sato, S. Invasion of Laguncula pulchella (Gastropoda: Naticidae) and predator–prey interactions with bivalves on the Tona coast, Miyagi prefecture, northern Japan. Biol. Invasions 15, 587–598 (2013).Article 

Google Scholar  More

Doubleday, Z. A. et al. Global proliferation of cephalopods. Curr. Biol. 26, R406–R407 (2016).CAS 
PubMed 
Article 

Google Scholar 
Jereb, P. et al. Cephalopod biology and fisheries in Europe: II. Species Accounts. ICES Cooperative Research Report No vol. 325 (2015).ICES. ICES WGCEPH REPORT 2015 Interim Report of the Working Group on Cephalopod Fisheries and Life History (WGCEPH). 8–11 (2019).Quetglas, A. et al. Long-term spatiotemporal dynamics of cephalopod assemblages in the Mediterranean sea. Sci. Mar. 83, 33–42 (2019).Article 

Google Scholar 
Martins, H. R. Biological studies of the exploited stock of Loligo forbesi (Mollusca: Cephalopoda) in the Azores. J. Mar. Biol. Assoc. United Kingdom 62, 799–808 (1982).Article 

Google Scholar 
Guerra, A. & Rocha, F. The life history of Loligo vulgaris and Loligo forbesi (Cephalopoda: Loliginidae) in Galician waters (NW Spain). Fish. Res. 21, 43–69 (1994).Article 

Google Scholar 
Pierce, G. J. & Boyle, P. R. Empirical modelling of interannual trends in abundance of squid (Loligo forbesi) in Scottish waters. Fish. Res. 59, 305–326 (2003).Article 

Google Scholar 
Lishchenko, F. et al. A review of recent studies on the life history and ecology of European cephalopods with emphasis on species with the greatest commercial fishery and culture potential. Fish. Res. 236, 105847 (2021).Article 

Google Scholar 
Laptikhovsky, V. et al. Identification of benthic egg masses and spawning grounds in commercial squid in the English Channel and Celtic Sea: Loligo vulgaris vs L. forbesii. Fish. Res. 241, 106004 (2021).Article 

Google Scholar 
Souza, H. V. et al. Analysis of the mitochondrial COI gene and its informative potential for evolutionary inferences in the families Coreidae and Pentatomidae (Heteroptera). Genet. Mol. Res. 15, 1–14 (2016).CAS 

Google Scholar 
Brierley, A. S. et al. Genetic variation in the neritic squid Loligo forbesi (Myopsida: Loliginidae) in the northeast Atlantic Ocean. Mar. Biol. 122, 79–86 (1995).Article 

Google Scholar 
Shaw, P. W. et al. Subtle population structuring within a highly vagile marine invertebrate, the veined squid Loligo forbesi, demonstrated with microsatellite DNA markers. Mol. Ecol. 8, 407–417 (1999).CAS 
Article 

Google Scholar 
Ellegren, H. Microsatellites: Simple sequences with complex evolution. Nat. Rev. Genet. 5, 435–445 (2004).CAS 
PubMed 
Article 

Google Scholar 
Begg, G. A. & Waldman, J. R. An holistic approach to fish stock identification. Fish. Res. 43, 35–44 (1999).Article 

Google Scholar 
Shaw, P. W. Polymorphic microsatellite markers in a cephalopod: The veined squid Loligo forbesi. Mol. Ecol. 6, 297–298 (1997).CAS 
Article 

Google Scholar 
Emery, A. M. et al. New microsatellite markers for assessment of paternity in the squid Loligo forbesi (Mollusca: Cephalopoda). Mol. Ecol. 9, 110–112 (2000).CAS 
PubMed 
Article 

Google Scholar 
Butler, J. M. Advanced Topics in Forensic DNA Typing: Interpretation (Elsevier Academic Press, 2015).
Google Scholar 
Park, S. D. E. Trypanotolerance in West African Cattle and the Population Genetics Effects of Selection. Trinity Coll. (2001).Nei, M. Molecular Evolutionary Genetics (Columbia University Press, 1987).
Google Scholar 
Hedrick, P. W. Genetics of Populations (Science Books International, 1983).
Google Scholar 
Weir, B. S. & Cockerham, C. C. Estimating F statistics for Population Structure. Evolution 38, 1358–1370 (1984).CAS 
PubMed 

Google Scholar 
Raymond, M. & Rousset, F. GENEPOP (version 1.2): Population genetics software for exact tests and ecumenicism. J. Hered. 86, 248–249 (1995).Article 

Google Scholar 
Rousset, F. GENEPOP’007: A complete re-implementation of the GENEPOP software for Windows and Linux. Mol. Ecol. Resour. 8, 103–106 (2008).PubMed 
Article 

Google Scholar 
Kalinowski, S. T. HP-RARE 1.0: A computer program for performing rarefaction on measures of allelic richness. Mol. Ecol. Notes 5, 187–189 (2005).CAS 
Article 

Google Scholar 
Excoffier, L. et al. Arlequin (version 3.0): An integrated software package for population genetics data analysis. Evol. Bioinforma. 1, 117693430500100 (2005).Article 

Google Scholar 
Pritchard, J. K. Inference of population structure using multilocus genotype data. Genetics 155, 945–959 (2000).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Gilbert, K. J. et al. Recommendations for utilizing and reporting population genetic analyses: The reproducibility of genetic clustering using the program structure. Mol. Ecol. 21, 4925–4930 (2012).PubMed 
Article 

Google Scholar 
Porras-Hurtado, L. et al. An overview of STRUCTURE: Applications, parameter settings, and supporting software. Front. Genet. 4, 1–13 (2013).CAS 
Article 

Google Scholar 
Evanno, G. et al. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol. Ecol. 14, 2611–2620 (2005).CAS 
PubMed 
Article 

Google Scholar 
Kopelman, N. M. et al. Clumpak: A program for identifying clustering modes and packaging population structure inferences across K. Mol. Ecol. Resour. 15, 1179–1191 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Folmer, O. et al. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol. Mar. Biol. Biotechnol. 3, 294–299 (1994).CAS 
PubMed 

Google Scholar 
Hall, T. A. BioEdit: A user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp. Ser. 41, 95–98 (1999).CAS 

Google Scholar 
Anderson, F. E. Phylogeny and historical biogeography of the loliginid squids (Mollusca: Cephalopoda) based on mitochondrial DNA sequence data. Mol. Phylogenet. Evol. 15, 191–214 (2000).CAS 
PubMed 
Article 

Google Scholar 
Gebhardt, K. & Knebelsberger, T. Identification of cephalopod species from the North and Baltic Seas using morphology, COI and 18S rDNA sequences. Helgol. Mar. Res. 69, 259–271 (2015).ADS 
Article 

Google Scholar 
Lobo, J. et al. Enhanced primers for amplification of DNA barcodes from a broad range of marine metazoans. BMC Ecol. 13, 1–8 (2013).Article 
CAS 

Google Scholar 
de Luna Sales, J. B. et al. New molecular phylogeny of the squids of the family Loliginidae with emphasis on the genus Doryteuthis Naef ,1912: Mitochondrial and nuclear sequences indicate the presence of cryptic species in the southern Atlantic Ocean. Mol. Phylogenet. Evol. 68, 293–299 (2013).Article 

Google Scholar 
Tatulli, G. et al. A rapid colorimetric assay for on-site authentication of cephalopod species. Biosensors 10, 3–10 (2020).Article 
CAS 

Google Scholar 
Velasco, A. et al. A new rapid method for the authentication of common octopus (Octopus vulgaris) in seafood products using recombinase polymerase amplification (rpa) and lateral flow assay (lfa). Foods 10, 1825 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Luz, A. & Keskin, E. Building Reference Library for Marine Fish Species of Azores Archipelago and Bio-monitoring via DNA Metabarcoding. https://www.ncbi.nlm.nih.gov/nuccore/MT491734 (2020).BoldSystems. https://boldsystems.org/index.php/Public_RecordView?processid=AZB030-20 (2018). (Accessed 2 May 2022).Tamura, K. et al. MEGA6: Molecular evolutionary genetics analysis version 6.0. Mol. Biol. Evol. 30, 2725–2729 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ronquist, F. & Huelsenbeck, J. P. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574 (2003).CAS 
PubMed 
Article 

Google Scholar 
Rambaut, A. et al. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67, 901–904 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Bandelt, H.-J. et al. Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48 (2009).Article 

Google Scholar 
Librado, P. & Rozas, J. DnaSP v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452 (2009).CAS 
PubMed 
Article 

Google Scholar 
Schlitzer, R. Ocean Data View. (2013).Shaw, P. W. & Boyle, P. R. Multiple paternity within the brood of single females of Loligo forbesi (Cephalopoda: Loliginidae), demonstrated with microsatellite DNA markers. Mar. Ecol. Prog. Ser. 160, 279–282 (1997).ADS 
Article 

Google Scholar 
Emery, A. M. et al. Assignment of paternity groups without access to parental genotypes: Multiple mating and developmental plasticity in squid. Mol. Ecol. 10, 1265–1278 (2001).CAS 
PubMed 
Article 

Google Scholar 
Catarino, D. et al. The role of the Strait of Gibraltar in shaping the genetic structure of the Mediterranean Grenadier, Coryphaenoides mediterraneus, between the Atlantic and Mediterranean Sea. PLoS ONE 12, 1–24 (2017).
Google Scholar 
Gonzalez, E. G. & Zardoya, R. Relative role of life-history traits and historical factors in shaping genetic population structure of sardines (Sardina pilchardus). BMC Evol. Biol. 7, 1–12 (2007).Article 
CAS 

Google Scholar 
Reichow, D. & Smith, M. J. Microsatellites reveal high levels of gene flow among populations of the California squid Loligo opalescens. Mol. Ecol. 10, 1101–1109 (2001).CAS 
PubMed 
Article 

Google Scholar 
Shaw, P. W. et al. DNA markers indicate that distinct spawning cohorts and aggregations of Patagonian squid, Loligo gahi, do not represent genetically discrete subpopulations. Mar. Biol. 144, 961–970 (2004).CAS 
Article 

Google Scholar 
Göpel, A. Populationsgenetik und Phylogeographie des Nordischen Kalmars Loligo forbesii Steenstrup, 1856 in Europäischen Gewässern. Masterthesis, Univ. Rostock in German, 76pp (2020).Oesterwind, D. et al. Biology and meso-scale distribution patterns of North Sea cephalopods. Fish. Res. 106, 141–150 (2010).Article 

Google Scholar 
Sauer, W. H. H. et al. Tag recapture studies of the chokka squid Loligo vulgaris reynaudii d’Orbigny, 1845 on inshore spawning grounds on the south-east coast of South Africa. Fish. Res. 45, 283–289 (2000).ADS 
Article 

Google Scholar 
Knowlton, N. & Weigt, L. A. New dates and new rates for divergence across the Isthmus of Panama. Proc. R. Soc. B Biol. Sci. 265, 2257–2263 (1998).Article 

Google Scholar 
Pérez-Losada, M. et al. Testing hypotheses of population structuring in the Northeast Atlantic Ocean and Mediterranean Sea using the common cuttlefish Sepia officinalis. Mol. Ecol. 16, 2667–2679 (2007).PubMed 
Article 

Google Scholar 
O’Dor, R. K. Can understanding squid life-history strategies and recruitment improve management?. South African J. Mar. Sci. 7615, 193–206 (1998).Article 

Google Scholar 
Izquierdo, A. et al. Modelling in the Strait of Gibraltar: From operational oceanography to scale interactions. Fundam. i Prikl. Gidrofiz. 9, 15–24 (2016).
Google Scholar 
Clarke, M. & Hart, M. Treatise Online no. 102: Part M, Chapter 11: Statoliths and coleoid evolution. Treatise Online (2018).Hsü, K. J. et al. Late Miocene desiccation of the mediterranean. Nature 242, 240–244 (1973).ADS 
Article 

Google Scholar 
Garcia-Castellanos, D. et al. Catastrophic flood of the Mediterranean after the Messinian salinity crisis. Nature 462, 778–781 (2009).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Thunell, R. C. et al. Atlantic-mediterranean water exchange during the late neocene. Paleoceanography 2(6), 661 (1987).ADS 
Article 

Google Scholar 
Green, C. P. et al. Combining statolith element composition and fourier shape data allows discrimination of spatial and temporal stock structure of arrow squid (Nototodarus gouldi). Can. J. Fish. Aquat. Sci. 72, 1609–1618 (2015).Article 

Google Scholar  More

Paull, S. H. et al. From superspreaders to disease hotspots: Linking transmission across hosts and space. Front. Ecol. Environ. 10, 75–82 (2012).PubMed 
Article 

Google Scholar 
Sorensen, A., Van Beest, F. M. & Brook, R. K. Impacts of wildlife baiting and supplemental feeding on infectious disease transmission risk: A synthesis of knowledge. Prev. Vet. Med. 113, 356–363 (2014).PubMed 
Article 

Google Scholar 
Gortázar, C., Acevedo, P., Ruíz-Fons, F. & Vicente, J. Disease risks and overabundance of game species. Eur. J. Wildl. Res. 52, 81–87 (2006).Article 

Google Scholar 
Brittingham, M. C. & Temple, S. A. Avian disease and winter bird feeding. Passeng. Pigeon 50, (1998).Franz, M., Kramer-Schadt, S., Greenwood, A. D. & Courtiol, A. Sickness-induced lethargy can increase host contact rates and pathogen spread in water-limited landscapes. Funct. Ecol. 32, 2194–2204 (2018).Article 

Google Scholar 
Galbraith, J. A., Stanley, M. C., Jones, D. N. & Beggs, J. R. Experimental feeding regime influences urban bird disease dynamics. J. Avian Biol. 48, 700–713 (2017).Article 

Google Scholar 
Moyers, S. C., Adelman, J. S., Farine, D. R., Thomason, C. A. & Hawley, D. M. Feeder density enhances house finch disease transmission in experimental epidemics. Philos. Trans. R. Soc. B Biol. Sci. 373(1745), 20170090 (2018).Article 
CAS 

Google Scholar 
Keesing, F., Holt, R. D. & Ostfeld, R. S. Effects of species diversity on disease risk. Ecol. Lett. 9, 485–498 (2006).CAS 
PubMed 
Article 

Google Scholar 
Mathiasson, M. E. & Rehan, S. M. Status changes in the wild bees of north-eastern North America over 125 years revealed through museum specimens. Insect Conserv. Divers. 12, 278–288 (2019).
Google Scholar 
Vanbergen, A. J. & Initiative, I. P. Threats to an ecosystem service: pressures on pollinators. Front. Ecol. Env. 11, 251–259 (2013).Article 

Google Scholar 
Buhk, C. et al. Flower strip networks offer promising long term effects on pollinator species richness in intensively cultivated agricultural areas. BMC Ecol. 18(1), 1–13 (2018).Article 

Google Scholar 
Morandin, L. A. & Kremen, C. Hedgerow restoration promotes pollinator populations and exports native bees to adjacent fields. Ecol. Appl. 23, 829–839 (2013).PubMed 
Article 

Google Scholar 
Williams, N. M. et al. Native wildflower plantings support wild bee abundance and diversity in agricultural landscapes across the United States. Ecol. Appl. 25, 2119–2131 (2015).PubMed 
Article 

Google Scholar 
Graystock, P. et al. Dominant bee species and floral abundance drive parasite temporal dynamics in plant-pollinator communities. Nat. Ecol. Evol. 4, 1358–1367 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Adler, L. S. et al. Disease where you dine: Plant species and floral traits associated with pathogen transmission in bumble bees. Ecology 99, 2535–2545 (2018).PubMed 
Article 

Google Scholar 
Alger, S. A., Burnham, P. A. & Brody, A. K. Flowers as viral hot spots: Honey bees (Apis mellifera) unevenly deposit viruses across plant species. PLoS ONE 14(9), e0221800 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
McNeil, D. J. et al. Bumble bees in landscapes with abundant floral resources have lower pathogen loads. Sci. Rep. 10, 1–12 (2020).Article 
CAS 

Google Scholar 
Daughenbaugh, K. F. et al. Metatranscriptome analysis of sympatric bee species identifies bee virus variants and a new virus, andrena-associated bee virus-1. Viruses 13, 291 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Alger, S. A., Alexander Burnham, P., Boncristiani, H. F. & Brody, A. K. RNA virus spillover from managed honeybees (Apis mellifera) to wild bumblebees (Bombus spp.). PLoS ONE 14, e0217822 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ravoet, J. et al. Widespread occurrence of honey bee pathogens in solitary bees. J. Invertebr. Pathol. 122, 55–58 (2014).PubMed 
Article 

Google Scholar 
Hayes, S. E., Tuiwawa, M., Stevens, M. I. & Schwarz, M. P. A recipe for weed disaster in islands: A super-generalist native pollinator aided by a ‘Parlourmaid’ plant welcome new arrivals in Fiji. Biol. Invasions 21, 1643–1655 (2019).Article 

Google Scholar 
Levenson, H. & Tarpy, D. R. Pollinator community response to planted pollinator habitat in agroecosystems over time. Authorea https://doi.org/10.22541/au.164191433.37143936/v1 (2022).Article 

Google Scholar 
Graystock, P., Yates, K., Darvill, B., Goulson, D. & Hughes, W. O. H. Emerging dangers: Deadly effects of an emergent parasite in a new pollinator host. J. Invertebr. Pathol. 114, 114–119 (2013).PubMed 
Article 

Google Scholar 
Genersch, E., Yue, C., Fries, I. & De Miranda, J. R. Detection of Deformed wing virus, a honey bee viral pathogen, in bumble bees (Bombus terrestris and Bombus pascuorum) with wing deformities. J. Invertebr. Pathol. 91, 61–63 (2006).PubMed 
Article 

Google Scholar 
Müller, U., McMahon, D. P. & Rolff, J. Exposure of the wild bee Osmia bicornis to the honey bee pathogen Nosema ceranae. Agric. For. Entomol. 21, 363–371 (2019).Article 

Google Scholar 
Strobl, V., Yañez, O., Straub, L., Albrecht, M. & Neumann, P. Trypanosomatid parasites infecting managed honeybees and wild solitary bees. Int. J. Parasitol. 49, 605–613 (2019).PubMed 
Article 

Google Scholar 
Gisder, S. et al. Rapid gastrointestinal passage may protect Bombus terrestris from becoming a true host for Nosema ceranae. Appl. Environ. Microbiol. 86(12), e00629-20 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Tehel, A., Streicher, T., Tragust, S. & Paxton, R. J. Experimental infection of bumblebees with honeybee-associated viruses: No direct fitness costs but potential future threats to novel wild bee hosts. R. Soc. Open Sci. 7(7), 200480 (2020).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Reynaldi, F. J., Sguazza, G. H., Albicoro, F. J., Pecoraro, M. R. & Galosi, C. M. First molecular detection of co-infection of honey bee viruses in asymptomatic Bombus atratus in South America. Braz. J. Biol. 73, 797–800 (2013).CAS 
PubMed 
Article 

Google Scholar 
Schoonvaere, K. et al. Unbiased RNA shotgun metagenomics in social and solitary wild bees detects associations with eukaryote parasites and new viruses. PLoS ONE 11(12), e0168456 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Melathopoulos, A. et al. Viruses of managed alfalfa leafcutting bees (Megachille rotundata Fabricus) and honey bees (Apis mellifera L.) in Western Canada: Incidence, impacts, and prospects of cross-species viral transmission. J. Invertebr. Pathol. 146, 24–30 (2017).PubMed 
Article 

Google Scholar 
Schoonvaere, K., Smagghe, G., Francis, F. & de Graaf, D. C. Study of the metatranscriptome of eight social and solitary wild bee species reveals novel viruses and bee parasites. Front. Microbiol. 9, 177 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Payne, A. N., Shepherd, T. F. & Rangel, J. The detection of honey bee (Apis mellifera)-associated viruses in ants. Sci. Rep. 10(1), 1–8 (2020).Article 
CAS 

Google Scholar 
Dalmon, A. et al. Possible spillover of pathogens between bee communities foraging on the same floral resource. Insects 12(2), 122 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Kojima, Y. et al. Infestation of Japanese native honey bees by tracheal mite and virus from non-native European honey Bees in Japan. Microb. Ecol. 62, 895–906 (2011).PubMed 
Article 

Google Scholar 
Graystock, P. et al. The Trojan hives: Pollinator pathogens, imported and distributed in bumblebee colonies. J. Appl. Ecol. 50, 1207–1215 (2013).Article 

Google Scholar 
Plischuk, S. et al. South American native bumblebees (Hymenoptera: Apidae) infected by Nosema ceranae (Microsporidia), an emerging pathogen of honeybees (Apis mellifera). Environ. Microbiol. Rep. 1, 131–135 (2009).PubMed 
Article 

Google Scholar 
Evison, S. E. et al. Pervasiveness of parasites in pollinators. PLoS ONE 7(1), e30641 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Graystock, P., Goulson, D. & Hughes, W. O. H. The relationship between managed bees and the prevalence of parasites in bumblebees. PeerJ 2, e522 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Graystock, P., Goulson, D. & Hughes, W. O. Parasites in bloom: Flowers aid dispersal and transmission of pollinator parasites within and between bee species. Proc. R. Soc. B Biol. Sci. 282(1813), 20151371 (2015).Article 

Google Scholar 
Tripodi, A. D., Szalanski, A. L. & Strange, J. P. Novel multiplex PCR reveals multiple trypanosomatid species infecting North American bumble bees (Hymenoptera: Apidae: Bombus). J. Invertebr. Pathol. 153, 147–155 (2018).PubMed 
Article 

Google Scholar 
Singh, R. et al. RNA viruses in hymenopteran pollinators: evidence of inter-taxa virus transmission via pollen and potential impact on non-Apis hymenopteran species. PLoS ONE 5(12), e14357 (2010).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Peng, W. et al. Host range expansion of honey bee Black Queen Cell Virus in the bumble bee, Bombus huntii. Apidologie 42, 650–658 (2011).Article 

Google Scholar 
Levitt, A. L. et al. Cross-species transmission of honey bee viruses in associated arthropods. Virus Res. 176, 232–240 (2013).CAS 
PubMed 
Article 

Google Scholar 
Fürst, M. A., McMahon, D. P., Osborne, J. L., Paxton, R. J. & Brown, M. J. F. Disease associations between honeybees and bumblebees as a threat to wild pollinators. Nature 506, 364–366 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Gamboa, V. et al. Bee pathogens found in Bombus atratus from Colombia: A case study. J. Invertebr. Pathol. 129, 36–39 (2015).PubMed 
Article 

Google Scholar 
Radzevičiūtė, R. et al. Replication of honey bee-associated RNA viruses across multiple bee species in apple orchards of Georgia, Germany and Kyrgyzstan. J. Invertebr. Pathol. 146, 14–23 (2017).PubMed 
Article 
CAS 

Google Scholar 
Murray, E. A. et al. Viral transmission in honey bees and native bees, supported by a global black queen cell virus phylogeny. Environ. Microbiol. 21, 972–983 (2019).PubMed 
Article 

Google Scholar 
Dobelmann, J., Felden, A. & Lester, P. J. Genetic strain diversity of multi-host RNA viruses that infect a wide range of pollinators and associates is shaped by geographic origins. Viruses 12, 13–15 (2020).Article 
CAS 

Google Scholar 
Olgun, T., Everhart, S. E., Anderson, T. & Wu-Smart, J. Comparative analysis of viruses in four bee species collected from agricultural, urban, and natural landscapes. PLoS ONE 15(6), e0234431 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Fearon, M. L. & Tibbetts, E. A. Pollinator community species richness dilutes prevalence of multiple viruses within multiple host species. Ecology 102(5), e03305 (2021).PubMed 
Article 

Google Scholar 
Sokół, R., Michalczyk, M. & Michołap, P. Preliminary studies on the occurrence of honeybee pathogens in the national bumblebee population. Ann. Parasitol. 64, 385–390 (2018).PubMed 

Google Scholar 
Bravi, M. E. et al. Wild bumble bees (Hymenoptera: Apidae: Bombini) as a potential reservoir for bee pathogens in northeastern Argentina. J. Apic. Res. 58, 710–713 (2019).Article 

Google Scholar 
Mazzei, M. et al. Detection of replicative Kashmir Bee Virus and Black Queen Cell Virus in Asian hornet Vespa velutina (Lepelieter 1836) in Italy. Sci. Rep. 9, 1–9 (2019).CAS 
Article 

Google Scholar 
Li, J. et al. Cross-species infection of deformed wing virus poses a new threat to pollinator conservation. J. Econ. Entomol. 104, 732–739 (2011).PubMed 
Article 

Google Scholar 
Sachman-Ruiz, B., Narváez-Padilla, V. & Reynaud, E. Commercial Bombus impatiens as reservoirs of emerging infectious diseases in central México. Biol. Invasions 17, 2043–2053 (2015).Article 

Google Scholar 
Jones, L. J., Ford, R. P., Schilder, R. J. & López-Uribe, M. M. Honey bee viruses are highly prevalent but at low intensities in wild pollinators of cucurbit agroecosystems. J. Invertebr. Pathol. 185, 107667 (2021).CAS 
PubMed 
Article 

Google Scholar 
Dolezal, A. G. et al. Honey bee viruses in wild bees: Viral prevalence, loads, and experimental inoculation. PLoS ONE 11(11), e0166190 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Mazzei, M. et al. First detection of replicative deformed wing virus (DWV) in Vespa velutina nigrithorax. Bull. Insectology 71, 211–216 (2018).
Google Scholar 
Plischuk, S. et al. Parasites and pathogens associated with native bumble bees (Hymenoptera: Apidae: Bombus spp.) from highlands in Bolivia and Peru. Stud. Neotrop. Fauna Environ. Stud. https://doi.org/10.1080/01650521.2020.1743551 (2020).Article 

Google Scholar 
McMahon, D. P. et al. A sting in the spit: Widespread cross-infection of multiple RNA viruses across wild and managed bees. J. Anim. Ecol. 84, 615–624 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
Bailes, E. J. et al. First detection of bee viruses in hoverfly (syrphid) pollinators. Biol. Lett. 14(2), 20180001 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Pritchard, Z. A. et al. Do viruses from managed honey bees (Hymenoptera: Apidae) endanger wild bees in native prairies?. Environ. Entomol. 50, 455–466 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Danforth, B. N., Mitchell, P. L. & Packer, L. Mitochondrial DNA differentiation between two cryptic Halictus (Hymenoptera: Halictidae) species. Ann. Entomol. Soc. Am. 91, 387–391 (1998).CAS 
Article 

Google Scholar 
Grozinger, C. M. & Flenniken, M. L. Bee viruses: Ecology, pathogenicity, and impacts. Annu. Rev. Entomol. 64, 205–226 (2019).CAS 
PubMed 
Article 

Google Scholar 
Antúnez, K. et al. Immune suppression in the honey bee (Apis mellifera) following infection by Nosema ceranae (Microsporidia). Environ. Microbiol. 11, 2284–2290 (2009).PubMed 
Article 
CAS 

Google Scholar 
Cameron, S. A. et al. Patterns of widespread decline in North American bumble bees. PNAS 108, 662–667 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Leite, G. M., Magan, N. & Medina, A. Comparison of different bead-beating RNA extraction strategies: An optimized method for filamentous fungi. J. Microbiol. Methods 88, 413–418 (2012).CAS 
PubMed 
Article 

Google Scholar 
Simms, D., Cizdziel, P. & Chomczynski, P. TRIzol: A new reagent for optimal single-step isolation of RNA. Focus (Madison) 15, 99–102 (1993).
Google Scholar 
Vandesompele, J. et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 3(7), 1–12 (2002).Article 

Google Scholar 
Mwalili, S. M., Lesaffre, E. & Declerck, D. The zero-inflated negative binomial regression model with correction for misclassification: An example in caries research. Stat. Methods Med. Res. 17, 123–139 (2008).MathSciNet 
PubMed 
MATH 
Article 

Google Scholar 
R Core Team. R: Language and Environment for Statistical Computing. R Foundation for Statistical Computer (2018). Available at: https://www.r-project.org/.Jackman, S. et al. Package ‘pscl’. (2020).Canty, A. & Ripley, B. Package ‘boot’. (2021).Figueroa, L. L. et al. Landscape simplification shapes pathogen prevalence in plant-pollinator networks. Ecol. Lett. 23, 1212–1222 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
National Heritage Program. Species/Community Search. National Heritage Program: Natural and Cultural Resources (2021). Available at: https://ncnhp.org/data/speciescommunity-search.Hatfield, R. et al. IUCN Assessments for North American Bombus spp. (2014).Sersic, A. N., Masco, M. & Noy-Meir, I. Natural hybridization between species of Calceolaria with different pollination syndromes in southern Patagonia, Argentina. Plant Syst. Evol. 230, 111–124 (2001).Article 

Google Scholar 
Otti, O. & Schmid-Hempel, P. Nosema bombi: A pollinator parasite with detrimental fitness effects. J. Invertebr. Pathol. 96, 118–124 (2007).PubMed 
Article 

Google Scholar 
Crabbe, J. C., Wahlsten, D. & Dudek, B. C. Genetics of mouse behavior: Interactions with laboratory environment. Science 284(5420), 1670–1672 (1999).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Wahlsten, D. et al. Different data from different labs: Lessons from studies of gene-environment interaction. J. Nuerobiol. 54, 283–311 (2003).Article 

Google Scholar 
Brownie, J. et al. The elimination of primer-dimer accumulation in PCR. Nucleic Acids Res. 25(16), 3235–3241 (1997).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Boncristiani, H. F. et al. In vitro infection of pupae with israeli acute paralysis virus suggests disturbance of transcriptional homeostasis in honey bees (Apis mellifera). PLoS ONE 8(9), e73429 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar  More

Gerber, B. D., Karpanty, S. M. & Randrianantenaina, J. Activity patterns of carnivores in the rain forests of Madagascar: Implications for species coexistence. J. Mammal. 93, 667–676 (2012).Article 

Google Scholar 
Azevedo, F., Lemos, F., Freitas-Junior, M., Rocha, D. & Azevedo, F. Puma activity patterns and temporal overlap with prey in a human-modified landscape at Southeastern Brazil. J. Zool. 305, 246–255 (2018).Article 

Google Scholar 
Wellborn, G. A., Skelly, D. K. & Werner, E. E. Mechanisms creating community structure across a freshwater habitat gradient. Annu. Rev. Ecol. Syst. 27, 337–363 (1996).Article 

Google Scholar 
Sitvarin, M. I., Rypstra, A. L. & Harwood, J. D. Linking the green and brown worlds through nonconsumptive predator effects. Oikos 125, 1057–1068 (2016).Article 

Google Scholar 
Damien, M. & Tougeron, K. Prey–predator phenological mismatch under climate change. Curr. Opin. Insect Sci. 35, 60–68 (2019).PubMed 
Article 

Google Scholar 
Relyea, R. A. & Werner, E. E. Quantifying the relation between predator-induced behavior and growth performance in larval anurans. Ecology 80, 2117–2124 (1999).Article 

Google Scholar 
Relyea, R. A. Morphological and behavioral plasticity of larval anurans in response to different predators. Ecology 82, 523–540 (2001).Article 

Google Scholar 
Osman, R. W. & Whitlatch, R. B. The control of the development of a marine benthic community by predation on recruits. J. Exp. Mar. Biol. Ecol. 311, 117–145 (2004).Article 

Google Scholar 
Schmidt, B. R., Băncilă, R. I., Hartel, T., Grossenbacher, K. & Schaub, M. Shifts in amphibian population dynamics in response to a change in the predator community. Ecosphere 12, e03528 (2021).Article 

Google Scholar 
Falaschi, M., Melotto, A., Manenti, R. & Ficetola, G. F. Invasive species and amphibian conservation. Herpetologica 76, 216–227 (2020).Article 

Google Scholar 
Gamradt, S. C. & Kats, L. B. Effect of introduced crayfish and mosquito fish on California newts. Conserv. Biol. 10, 1155–1162 (1996).Article 

Google Scholar 
Matthews, K. R., Knapp, R. A. & Pope, K. L. Garter snake distributions in high-elevation aquatic ecosystems: Is there a link with declining amphibian populations and nonnative trout introductions?. J. Herpetol. 36, 16–22 (2002).Article 

Google Scholar 
Dodds, W. K. & Whiles, M. R. Freshwater Ecology: Concepts and Environmental Applications 3rd edn. (Elsevier, 2002).
Google Scholar 
Preisser, E. L., Bolnick, D. I. & Benard, M. F. Scared to death? The effects of intimidation and consumption in predator–prey interactions. Ecology 86, 501–509 (2005).Article 

Google Scholar 
Le Roux, E., Kerley, G. I. & Cromsigt, J. P. Megaherbivores modify trophic cascades triggered by fear of predation in an African savanna ecosystem. Curr. Biol. 28, 2493–2499 (2018).PubMed 
Article 
CAS 

Google Scholar 
Daversa, D. et al. Broadening the ecology of fear: Non-lethal effects arise from diverse responses to predation and parasitism. Proc. R. Soc. B 288, 20202966 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Benard, M. F. Predator-induced phenotypic plasticity in organisms with complex life histories. Annu. Rev. Ecol. Evol. Syst. 35, 651–673 (2004).Article 

Google Scholar 
Van Buskirk, J. & Schmidt, B. R. Predator-induced phenotypic plasticity in larval newts: Trade-offs, selection, and variation in nature. Ecology 81, 3009–3028 (2000).Article 

Google Scholar 
McCollum, S. A. & Van Buskirk, J. Costs and benefits of a predator-induced polyphenism in the gray treefrog Hyla chrysoscelis. Evolution 50, 583–593 (1996).PubMed 

Google Scholar 
Skelly, D. K. Tadpole communities: pond permanence and predation are powerful forces shaping the structure of tadpole communities. Am. Sci. 85, 36–45 (1997).ADS 

Google Scholar 
Sih, A. & Moore, R. D. Delayed hatching of salamander eggs in response to enhanced larval predation risk. Am. Nat. 142, 947–960 (1993).CAS 
PubMed 
Article 

Google Scholar 
Warkentin, K. M. Adaptive plasticity in hatching age: A response to predation risk trade-offs. Proc. Natl. Acad. Sci. 92, 3507–3510 (1995).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Johnson, J. B., Saenz, D., Adams, C. K. & Conner, R. N. The influence of predator threat on the timing of a life-history switch point: Predator-induced hatching in the southern leopard frog (Rana sphenocephala). Can. J. Zool. 81, 1608–1613 (2003).Article 

Google Scholar 
Wilbur, H. M. & Fauth, J. E. Experimental aquatic food webs: Interactions between two predators and two prey. Am. Nat. 135, 176–204 (1990).Article 

Google Scholar 
Laurila, A. Behavioural responses to predator chemical cues and local variation in antipredator performance in Rana temporaria tadpoles. Oikos 88, 159–168 (2000).CAS 
Article 

Google Scholar 
Gomez-Mestre, I. et al. The shape of things to come: Linking developmental plasticity to post-metamorphic morphology in anurans. J. Evol. Biol. 23, 1364–1373 (2010).CAS 
PubMed 
Article 

Google Scholar 
Vieira, E. A., Duarte, L. F. L. & Dias, G. M. How the timing of predation affects composition and diversity of species in a marine sessile community?. J. Exp. Mar. Biol. Ecol. 412, 126–133 (2012).Article 

Google Scholar 
Andrade, M. R., Albeny-Simões, D., Breaux, J. A., Juliano, S. A. & Lima, E. Are behavioural responses to predation cues linked across life cycle stages?. Ecol. Entomol. 42, 77–85 (2017).Article 

Google Scholar 
Knapp, R. A. Effects of nonnative fish and habitat characteristics on lentic herpetofauna in Yosemite National Park, USA. Biol. Conserv. 121, 265–279 (2005).Article 

Google Scholar 
Kiesecker, J. M. & Blaustein, A. R. Population differences in responses of red-legged frogs (Rana aurora) to introduced bullfrogs. Ecology 78, 1752–1760 (1997).Article 

Google Scholar 
Nunes, A. L., Orizaola, G., Laurila, A. & Rebelo, R. Rapid evolution of constitutive and inducible defenses against an invasive predator. Ecology 95, 1520–1530 (2014).PubMed 
Article 

Google Scholar 
Polo-Cavia, N., Gonzalo, A., López, P. & Martín, J. Predator recognition of native but not invasive turtle predators by naïve anuran tadpoles. Anim. Behav. 80, 461–466 (2010).Article 

Google Scholar 
Zhang, F., Zhao, J., Zhang, Y., Messenger, K. & Wang, Y. Antipredator behavioral responses of native and exotic tadpoles to novel predator. Asian Herpetol. Res. 6, 51–58 (2015).
Google Scholar 
Lowe, S., Browne, M., Boudjelas, S. & De Poorter, M. 100 of the World’s Worst Invasive Alien Species: A Selection from the Global Invasive Species Database (Invasive Species Specialist Group, 2000).
Google Scholar 
TTWG. Conservation biology of freshwater turtles and tortoises: A compilation project of the IUCN/SSC tortoise and freshwater turtle specialist group. in Chelonian Research Monographs 7 Turtle of the World: Annotated Checklist and Atlas of Taxonomy, Synonymy, Distribution, and Conversation Status. 8th edn. (eds. Rhodin, A.G.J., Iverson, J.B., van Dijk, P.P., Saumure, R.A., Buhlmann, K.A., Pritchard, P.C.H., Mittermeier, R.A.). 1–292. (Chelonian Research Foundation and Turtle Conservancy, 2017).GISD. Global Invasive Species Database. http://www.issg.org/database (2021).Berec, M., Klapka, V. & Zemek, R. Effect of an alien turtle predator on movement activity of European brown frog tadpoles. Ital. J. Zool. 83, 68–76 (2016).Article 

Google Scholar 
Vodrážková, M., Šetlíková, I. & Berec, M. Chemical cues of an invasive turtle reduce development time and size at metamorphosis in the common frog. Sci. Rep. 10, 1–6 (2020).Article 
CAS 

Google Scholar 
Gibbons, J., Greene, J. & Congdon, J. Life history and ecology of the slider turtle. in Temporal and Spatial Movement Patterns of Sliders and Other Turtles (ed. Gibbons, J.). 201–215. (Smithsonian Institution Press, 1990).Formanowicz, D. R. Anuran tadpole/aquatic insect predator-prey interactions: tadpole size and predator capture success. Herpetologica 42, 367–373 (1986).
Google Scholar 
Semlitsch, R. D. & Gibbons, J. W. Fish predation in size-structured populations of treefrog tadpoles. Oecologia 75, 321–326 (1988).ADS 
PubMed 
Article 

Google Scholar 
Teplitsky, C., Piha, H., Laurila, A. & Merilä, J. Common pesticide increases costs of antipredator defenses in Rana temporaria tadpoles. Environ. Sci. Technol. 39, 6079–6085 (2005).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Travis, J. Anuran size at metamorphosis: experimental test of a model based on intraspecific competition. Ecology 65, 1155–1160 (1984).Article 

Google Scholar 
Wilbur, H. M. & Collins, J. P. Ecological aspects of amphibian metamorphosis: Nonnormal distributions of competitive ability reflect selection for facultative metamorphosis. Science 182, 1305–1314 (1973).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Gosner, K. L. A simplified table for staging anuran embryos and larvae with notes on identification. Herpetologica 16, 183–190 (1960).
Google Scholar 
Woodward, G. & Hildrew, A. G. Body-size determinants of niche overlap and intraguild predation within a complex food web. J. Anim. Ecol. 71, 1063–1074 (2002).Article 

Google Scholar 
Relyea, R. A. Getting out alive: How predators affect the decision to metamorphose. Oecologia 152, 389–400 (2007).ADS 
PubMed 
Article 

Google Scholar 
Pujol-Buxó, E., San Sebastián, O., Garriga, N. & Llorente, G. A. How does the invasive/native nature of species influence tadpoles’ plastic responses to predators?. Oikos 122, 19–29 (2013).Article 

Google Scholar 
Phuge, S., Shetye, K. & Pandit, R. Effect of water level on insect-tadpole predator-prey interactions. Acta Oecol. 108, 103649 (2020).Article 

Google Scholar 
Loman, J. Early metamorphosis in common frog Rana temporaria tadpoles at risk of drying: An experimental demonstration. Amphibia-Reptilia 20, 421–430 (1999).Article 

Google Scholar 
Stav, G., Kotler, B. P. & Blaustein, L. Direct and indirect effects of dragonfly (Anax imperator) nymphs on green toad (Bufo viridis) tadpoles. Hydrobiologia 579, 85–93 (2007).Article 

Google Scholar 
Goldberg, T., Nevo, E. & Degani, G. Phenotypic plasticity in larval development of six amphibian species in stressful natural environments. Zool. Stud. 51, 345–361 (2012).
Google Scholar 
Kishida, O., Costa, Z., Tezuka, A. & Michimae, H. Inducible offences affect predator–prey interactions and life-history plasticity in both predators and prey. J. Anim. Ecol. 83, 899–906 (2014).PubMed 
Article 

Google Scholar 
Leips, J. & Travis, J. Metamorphic responses to changing food levels in two species of hylid frogs. Ecology 75, 1345–1356 (1994).Article 

Google Scholar 
Alford, R. A. & Harris, R. N. Effects of larval growth history on anuran metamorphosis. Am. Nat. 131, 91–106 (1988).Article 

Google Scholar 
Loman, J. Temperature, genetic and hydroperiod effects on metamorphosis of brown frogs Rana arvalis and R. temporaria in the field. J. Zool. 258, 115–129 (2002).Article 

Google Scholar 
Laugen, A. T. et al. Quantitative genetics of larval life-history traits in Rana temporaria in different environmental conditions. Genet. Res. 86, 161–170 (2005).CAS 
PubMed 
Article 

Google Scholar 
Brodie, E. D. & Formanowicz, D. R. Prey size preference of predators: Differential vulnerability of larval anurans. Herpetologica 39, 67–75 (1983).
Google Scholar 
Eklöv, P. & Werner, E. E. Multiple predator effects on size-dependent behavior and mortality of two species of anuran larvae. Oikos 88, 250–258 (2000).Article 

Google Scholar 
Urban, M. C. Predator size and phenology shape prey survival in temporary ponds. Oecologia 154, 571–580 (2007).ADS 
PubMed 
Article 

Google Scholar 
Jara, F. G. & Perotti, M. G. Risk of predation and behavioural response in three anuran species: influence of tadpole size and predator type. Hydrobiologia 644, 313–324 (2010).Article 

Google Scholar 
Wassersug, R. J. & Sperry, D. G. The relationships of locomotion to differential predation on Pseudacris triseriata (Anura: Hylidae). Ecology 58, 830–839 (1977).Article 

Google Scholar 
Huey, R. B. Sprint velocity of tadpoles (Bufo boreas) through metamorphosis. Copeia 1980, 537–540 (1980).Article 

Google Scholar 
Laurila, A. & Kujasalo, J. Habitat duration, predation risk and phenotypic plasticity in common frog (Rana temporaria) tadpoles. J. Anim. Ecol. 68, 1123–1132 (1999).Article 

Google Scholar 
Metcalfe, N. B. & Monaghan, P. Compensation for a bad start: grow now, pay later?. Trends Ecol. Evol. 16, 254–260 (2001).PubMed 
Article 

Google Scholar 
Downie, J. & Weir, A. Developmental arrest in Leptodactylus fuscus tadpoles (Anura: Leptodactylidae) III effect of length of arrest period on growth potential. Herpetol. J. 7, 85–92 (1997).
Google Scholar 
Smith, D. C. Adult recruitment in chorus frogs: Effects of size and date at metamorphosis. Ecology 68, 344–350 (1987).Article 

Google Scholar 
Altwegg, R. & Reyer, H. U. Patterns of natural selection on size at metamorphosis in water frogs. Evolution 57, 872–882 (2003).PubMed 
Article 

Google Scholar 
Brunelli, E. et al. Environmentally relevant concentrations of endosulfan impair development, metamorphosis and behaviour in Bufo bufo tadpoles. Aquat. Toxicol. 91, 135–142 (2009).CAS 
PubMed 
Article 

Google Scholar 
Boone, M. D. Juvenile frogs compensate for small metamorph size with terrestrial growth: Overcoming the effects of larval density and insecticide exposure. J. Herpetol. 39, 416–423 (2005).Article 

Google Scholar 
Schmidt, B. R., Hödl, W. & Schaub, M. From metamorphosis to maturity in complex life cycles: Equal performance of different juvenile life history pathways. Ecology 93, 657–667 (2012).PubMed 
Article 

Google Scholar  More

Sampling sitesLongitudinal observation campaigns for I. ricinus nymph activity were carried out at 11 sampling sites in forest areas from seven different tick observatories across France. Tick observatories are located at the following French municipal areas, where the coordinates of the centre of each municipal area and the climatic types29 are also provided as: (1) La Tour de Salvagny (45° 48′ 50.6″ N 4° 42′ 53.2″ E; Mixed climates); (2) Saint-Genès-Champanelle (45° 43′ 23.8″ N 3°01′ 08.0″ E; Mountain climate); (3) Etiolles (48° 37′ 59.9″ N 2° 28′ 00.1″ E; Degraded oceanic climate); (4) Carquefou (47° 17′ 58.5″ N 1° 29′ 26.0″ W; Oceanic climate); (5) Gardouch (43°23′ 25.7″ N 1° 41′ 02.1″ E; South-West Basin climate); (6) Velaine-en-Haye (48° 42′ 13.4″ N 6° 01′ 16.1″ E; Semi-continental climate); (7) Les Bordes (47° 48′ 47.3″ N 2° 24′ 01.3″ E; Degraded-oceanic climate) (Fig. 1). The observation campaigns were carried out from April/June 2014 to May/June 2021 in most observatories, except for Les Bordes, which began in April 2018.Figure 1The map was created using QGIS version 3.8, Zanzibar (https://www.qgis.org). The climatic region types were previously classified by Joly et al.29.The distribution of tick observatories according to the climatic region types of continental France: (1) Etiolles (degraded oceanic); (2) Velaine-en-Haye (semi-continental); (3) Les Bordes (degraded oceanic); (4) Carquefou (oceanic); (5) La Tour de Salvagny (mixed); (6) Saint-Genès-Champanelle (mountain); (7) Gardouch (south-west basin). Phenological patterns observed at each observatory were also indicated.Full size imageEach tick observatory corresponds to one sampling site except La Tour de Salvagny, Gardouch, and Les Bordes (Table S1). In La Tour de Salvagny, we had to withdraw the observations at the original site (La Tour de Salvagny A) in September 2016 because the site became no longer accessible. In April 2017, we continued our observations at a nearby site, approximately 2 km apart (La Tour de Salvagny B). In Gardouch, the activity of questing nymphs was observed both inside and outside the enclosed area of an experimental station on roe deer (Capreolus capreolus), referred to as Gardouch Inside and Gardouch Outside, respectively. The estimated population density of roe deer in Gardouch Inside (50 individuals per 100 ha) was higher than Gardouch Outside (less than 20 individuals per 100 ha) (H. Verheiden, personal communication, 15th October 2021). Furthermore, three sampling sites in Les Bordes, approximately 1.2 km apart, were referred to as Les Bordes A, B, and C, respectively. Additional sampling sites of these observatories were considered and reported as distinct sampling sites in further analyses, resulting in a total to 11 sampling sites from 7 observatories. Furthermore, due to their geographical proximity, meteorological/climatic factors of different sampling sites from the same observatories were considered identical in subsequent statistical analyses, whereas land cover and topography factors could be varied.Field observation campaigns were planned and carried out by local investigators who had been trained on the sampling protocol. The locations of forests, sampling sites, and passages were chosen where their biotopes are known to be suitable for I. ricinus tick populations around each observatory at the time the field observation campaigns started30. The observations were never carried out during the daytime when the weather was highly unfavourable to questing ticks, e.g., heavy rain, snow, or snow cover.Sampling protocol for questing Ixodes ricinus nymphsActivity of questing I. ricinus nymphs was observed by a cloth-dragging sampling technique31. Within a 1-km radius, a 1 m × 1 m white cloth was dragged over 10 observation units of 10 m short-grass vegetative forest floors, called transects. For each transect, a repeated removal sampling design was used27. The cloth-dragging sampling process was successively repeated three times per sampling. All nymphs found on white cloth in each campaign were removed and collected in a vial for subsequent morphological identification32 by the same acarologists at the corresponding laboratories. As a result, the questing nymph activity of each sampling site was monitored as a total number of confirmed I. ricinus nymphs collected from three repeated sampling on 10 transects, equivalent to a surface area of 100 m2. This measure was considered as an indicator for tick abundance on the day of sampling. The same transects were repeatedly sampled throughout the study period at approximately 1-month intervals.Environmental dataWe tested 28 environmental variables to explain the observed I. ricinus nymph activity (Table 1). These variables could be categorized as: (1) Daytime duration and meteorological variables (time-dependent, 9 variables); (2) Land cover, topography, and bioclimatic variables (time-independent, 19 variables).Table 1 Environmental variables (meteorological, land cover, topography, and bioclimatic variables) used to explain I. ricinus nymph counts per 100 m2 in regression analysis.Full size tableDaytime duration and meteorological variablesDaytime duration ((daytime)) from January 2013 to June 2021 at each sampling site was obtained from the corresponding latitude using geosphere package33. Hourly meteorological data (2-m temperature and relative humidity) were recorded locally at each forest. Subsequently, daily mean, minimum, and maximum values of temperature (({T}_{M}), ({T}_{N}), and ({T}_{X}); in °C) and relative humidity (({U}_{M}), ({U}_{N}), and ({U}_{X}); in %) were derived from these hourly records. The meteorological seasons of the temperate area in northern hemisphere are defined as: (1) Spring, 1st March to 31st May; (2) Summer, 1st June to 31st August; (3) Autumn, 1st September to 30th November; (4) Winter, 1st December to 28th or 29th February.Missing values found on these local daily-level variables were imputed by the random forest algorithm in mice package34. External daily meteorological data, i.e., daily average temperature and relative humidity, derived from neighbouring weather stations (Météo-France or INRAE), as well as month and year information, were used as auxiliary variables (Table S2). As a result, the imputation process creates a total of 500 iterated values for each variable. The median values of 500 imputations were used to replace the missing values.The imputed daily meteorological data were subsequently used to calculate the averaged values in different lagged time intervals for further analysis, called interval-average variables15. The interval-average variables were generated to reduce the uncertainty that might arise during the imputation process and to capture the cumulative effects of the meteorological variables, which were mean temperature ({T}_{M}) and minimum relative humidity ({U}_{N}). The interval-average variables were defined as the average values of a meteorological variable (Min) {({T}_{M}), ({U}_{N})} during a period between ({t}_{1}) to ({t}_{2}) month(s) before the sampling, denoted as ({M}^{{t}_{1}:{t}_{2}}), where 1 month consists of 28 days. As temperature conditions affect several ecological processes of tick populations, particularly developmental and questing rates3, the mean temperature ({T}_{M}) was selected for further analysis to reflect the overall temperature effects. While the minimum relative humidity ({U}_{N}) was chosen for the following reasons: (1) the survival of I. ricinus is highly sensitive to desiccation conditions6,7,8. As a result, when compared to mean or maximum relative humidity, minimum relative humidity is a relatively strong indicator of the effects of desiccation stress; (2) the variation of minimum relative humidity among all sites was higher than that of the mean and maximum relative humidity. This high variation allowed us to better describe meteorological characteristics of each sampling site.Here, we hypothesized that interval-average meteorological conditions influence the dynamics of observed nymph activity at different time lags in different manners. Short-term lags may have an impact on immediate responses, such as the probability of questing. At the same time, long-term lags may influence the dynamics of nymph abundance, which is associated with development and survival rates. Therefore, we explored the impact of each meteorological variable at following time lags on the observed nymphs activity in subsequent regression analysis: (1) 1-month moving average condition, ({M}^{0:1}); (2) previous 3-to-6-month moving average condition, ({M}^{3:6}); (3) 6-month moving average condition, ({M}^{0:6}); (4) 12-month moving average condition, ({M}^{0:12}). For instance, ({T}_{M}^{0:1}) denotes 1-month moving average temperature, representing an average of temperature between 0 and 1 months (0–28 days) before the day of sampling.In addition to the interval-average variables, monthly and seasonal average values of mean temperature and minimum relative humidity during the observation period were also calculated to describe the characteristics of meteorological conditions of each sampling site.Land cover, topography, and bioclimatic variablesWe obtained land cover, topography, and bioclimatic data from a 1-km radius buffer area around the center of each sampling site to capture habitat characteristics across all 10 transects. All the variables were handled and obtained by using QGIS version 3.8.035. The digital elevation model (DEM) data derived from the Shuttle Radar Topography Mission (SRTM) database36 was used to describe the topographic features of sampling sites, which included the mean (({mean}_{elv})) and standard deviation (({sd}_{elv})) of the elevation (in m above sea level), the proportion of flat area (({p}_{flat}); defined by the slope ≤ 2.5%37), the proportion of area facing north (({p}_{north})), east (({p}_{east})), west (({p}_{west})), and south (({p}_{south})), and the catchment area ((catchment)) as a proxy variable for moisture. Bioclimatic variables for each site (historical average conditions during 1970–2000) were derived from the WorldClim database38, including the annual mean temperature (({BIO1}_{Temp}); in °C), the mean diurnal range (({BIO2}_{Diur}); in °C), the maximum temperature of the warmest month (({BIO5}_{maxTemp}); in °C), and the annual precipitation (({BIO12}_{Prec}); in mm). The land cover features of each sampling site were described using the CORINE Land Cover (CLC) 201839, while the characteristics of forests were explained by the BD forêt version 2 data40. The forest fragmentation was characterized by the percentage of forest-covering area (({p}_{Forest})), the forest edge density (({ED}_{Forest}); in m/km2), and the number of forest patches (({n}_{Forest})). While the diversities of the land cover types (level-1 and level-2 CLC) and the forest types were calculated by using the Shannon’s diversity index41 ((H)) as (H=sum_{i=1}^{S}{p}_{i}mathrm{ln}{p}_{i}), where (S) is the total number of land cover/forest types and ({p}_{i}) is the proportion of land cover/forest type (i) within the 1-km radius buffer area. The Shannon’s diversity index for level-1 CLC, level-2 CLC, and forest types were denoted as ({H}_{CLC1}), ({H}_{CLC2}), and ({H}_{Forest}), respectively. Finally, the soil pH data (({pH}_{soil})) was retrieved from the European Soil Data Centre (ESDC) database42.Statistical analysisAll the statistical analyses were carried out using the programming language R version 3.6.043. The variations of questing nymph population of each site were described by using (1) baseline annual nymph counts (spatial variation); (2) phenological patterns (seasonal variation). A baseline annual nymph count of site (i) (({{N}_{base}}_{i})) was defined as a summation of monthly median nymph counts ({varvec{tilde{N}}}_{i}={{tilde{N }}_{i,t}}) across all 12 months (tin left{mathrm{1,2},dots ,12right}) and expressed as: ({{N}_{base}}_{i}=sum_{t=1}^{12}{tilde{N }}_{i,t}). Subsequently, the monthly median nymph counts of each site ({varvec{tilde{N}}}_{i}) were transformed into normalized monthly median nymph counts ({varvec{tilde{N}}}_{i}^{*}={{tilde{N }}_{i,t}^{*}}) following Eq. (1) to have a range value of 0 to 1, which allows us to compare phenological patterns among all sites that have different annual baseline nymph counts.$${tilde{N }}_{i,t}^{*}=frac{{tilde{N }}_{i,t}}{mathrm{max}({stackrel{sim }{{varvec{N}}}}_{i})}$$
(1)
The term (mathrm{max}({stackrel{sim }{{varvec{N}}}}_{i})) denoted the maximum monthly median nymph counts. The normalized median nymph count ({tilde{N }}_{i,t}^{*}) of 1 indicates the maximum nymph activity (peak), while the value ({tilde{N }}_{i,t}^{*}) of 0 designates the absence of nymph activity. Afterwards, the phenological patterns were descriptively classified using the following criteria: (1) the season which the peaks of activity arrive; (2) evidence of reduced activity during winter (November–January); (3) the number of activity waves in a year, whether the pattern is unimodal or bimodal. After assigning phenological patterns to each site, the overall trends of different patterns were derived from medians of the normalized monthly median nymph count ({tilde{N }}_{i,t}^{*}) from all sites that belonged to each pattern. Furthermore, the directional changes in the maximum nymph counts were tested using a Spearman’s rank correlation coefficient, a p-value More

Our results show that the effects of the forest disturbance drivers on biodiversity are likely to be different from those simply expected from the baseline proportions of the forest disturbance drivers if we take into account the threatened species’ distributions. The amount of forest habitat is a primary factor for species diversity of many taxa, including mammals, amphibians, reptiles, birds, insects, and plants18. Indeed, our results revealed that threatened forest species have been exposed to a disproportional decrease in their habitat amount globally (i.e., lower proportions of forest with no or minor loss in all regions when species ranges were considered). Although this finding may be intuitive as population size and/or species range are part of the criteria in the IUCN assessment19, the detected pattern supports the validity of our approach of combining a forest disturbance map and species ranges for evaluating the impact of forest disturbances on threatened species. Moreover, we found that the dominant drivers differ among regions: the proportion of forestry, for example, increased in northern regions such as North America and Europe, whereas that of shifting agriculture increased in tropical regions when threatened species’ distributions were considered. These facts indicate although several influential international schemes for conservation have been implemented for regulating forestry20,21, different mechanisms aiming to directly tackle the over land use for local agriculture may increase their importance when we consider conservation in tropical regions. Our findings suggest that the social and economic drivers underlying the forest disturbance that impacts biodiversity differ among regions or nations, and it is important to establish specific conservation strategies in order to be effective.Based on the findings, we further emphasize that the combinations of multiple interacting drivers are likely to vary among regions. For example, the frequency and extent of stand-replacing natural disturbances such as wildfires have clearly been magnified by climate change, particularly in the Northern Hemisphere (e.g.,22). After such natural disturbances, societal demand for timber and/or pest reduction compels forest managers to ‘salvage’ timber by logging before it deteriorates, a common practice even in locations otherwise exempt from conventional green-tree harvesting, such as national parks or wilderness areas23. Thus, salvage logging clearly mediates the interaction between disturbances by forestry and wildfires and is likely to further affect biodiversity under climate change. Especially in regions where infrastructure (e.g., irrigation systems) has not been well developed, unpredictable changes in precipitation due to climate change was reported to increase forest disturbance by unregulated increases of agricultural land use24. Such regions largely overlapped with regions where shifting agriculture was identified as a dominant disturbance driver for threatened species in this study. Moreover, species themselves shift their ranges in response to climate change25, which would also shift major disturbance drivers and influential interactions of drivers to which the species are exposed, given the region-specific driver patterns. These examples clearly suggest the necessity to understand both the region-specific interrelations among multiple drivers and species’ responses for better prediction of land-use change and thus its effects on biodiversity.Shifting agriculture was the most dominant driver in all tropical regions corresponding to the recent estimates suggesting that the cover of regenerating secondary forest is increasing worldwide26. We demonstrated that this tendency is more drastic especially within the range of threatened species. The effect of shifting agriculture per unit area might be more limited than that of commodity-driven deforestation, which permanently alters forests into other land uses, since habitat structure might recover as the forest vegetation regenerates to a secondary state following the abandonment of the small clearings. However, ample evidence shows that many types of agricultural activities significantly degrade the conservation value of primary forest, especially in the tropics27, which often recovers very slowly if ever28 with the loss of irreplaceable conservation values. Therefore, given the wide areas of dominance of shifting agriculture across all tropical regions, its effect is likely to be pervasive. Consistently, our results show that species extinction risk (i.e., IUCN Red List status) is positively related to the proportional coverage of shifting agriculture (Fig. 2). In addition, as expected, a larger current proportion of shifting agriculture within a species range worsens the change rate in IUCN Red List status of the species (Fig. 4b). Furthermore, the effect is anticipated to be magnified for forest specialists because they are exposed to larger proportions of shifting agriculture than are forest generalist (Fig. 2), and they are also reported to recover more slowly than do forest habitat generalists27,28.A guideline for forest restoration suggested that appropriately sized landscapes should contain ≥40% forest cover (higher percentages are likely needed in the tropics), with about 10% in a very large forest patch and the remaining 30% in many evenly dispersed smaller patches and semi-natural wooded elements (e.g., vegetation corridors)29. Importantly, the guideline also suggests that the patches should be embedded in a high-quality matrix. Although younger secondary forest cannot be a substitute for pristine forest until 50 years or more after a disturbance, it can help to improve the quality of matrix in agricultural landscapes30. Indeed, we show that the negative impacts of shifting agriculture and forestry on IUCN status change have improved over time (Fig. 4b, c), presumably corresponding to the forest regenerating and recovery process. In contrast, the pattern of commodity-driven deforestation, a land use accompanied with permanent forest loss, showed a prolonged negative impact on IUCN status change (Fig. 4a). Notably, whether regenerating forests can move towards a highly diverse and structurally complex state or towards a state of low to intermediate levels of biodiversity and structural complexity depends on the amount of remaining intact mature forest in the landscape29. Therefore, a promising direction for future research would be to develop our analysis further to include spatiotemporal relationships among mature forest remnants, secondary forests, disturbance drivers, and threatened species populations.For conserving the core patches of mature forests, the establishment of protected areas (PAs) is one of the most effective legal measures that has been widely used to regulate land use for biodiversity31. On the other hand, for improving matrix quality, balancing conservation and use of the ecosystem would be critically important; shifting agriculture, for example, causes forest degradation, but it also contributes to food supply chains sourced from smallholder farmers and to food security of local communities8. In fact, establishing mechanisms for managing biodiversity-friendly landscapes has been intensively discussed recently, given the large potential influence of these landscapes on conservation32. These mechanisms include setting an international target on OECMs15. Our finding of a disproportional decrease in forest proportions with minor or no loss within species ranges supports the urgency of the discussion. At the same time, our results highlight an opportunity because large portions of the disturbed forests for threatened species are dominated by shifting agriculture at the global scale, especially in the tropics. As suggested above, if manged properly, such landscapes can still retain or improve functions as essential habitats and/or matrix for a variety of forest-dwelling species. Our analytical method provides a tool set to identify and prioritize areas where such attempts are urgently needed.Global demands for natural resources and ecosystem services drive land use in forests33 and thus affect biodiversity. Therefore, connecting the supply chains to the five major drivers of forest disturbance and their spatial overlaps with biodiversity is essential to inform how we should regulate and design material flows from forest ecosystems to keep them sustainable by minimizing the effects on biodiversity. Existing studies examining the impacts of resource consumption on biodiversity through supply chains of various sectors have often been assessed at the country scale (e.g.,12), partly because the availability of statistics needed to estimate material flows in supply chains is usually limited at finer (i.e., subnational) scales (but see34). We believe that our study provides the first basis for filling the resolution gap between trade statistics and local biodiversity effects by identifying patterns of the local co-occurrence of biodiversity and the forest disturbance drivers that can be directly linked to resource production at the national scale. Note, however, that downscaling a remotely sensed global data set into finer scales inevitably propagates errors and biases which include both those in the original maps and those in the processed data produced by analyses. Thus, preparation of more high-resolution data sets is essential, especially for disturbance drivers and threatened species’ distributions in our case, to keep the errors and biases at a reasonable level at focal spatial scales.The effectiveness of area-based conservation measures to regulate land use for conservation including PAs and OECMs also depends strongly on social and ecosystem conditions. For example, a few studies show that the effectiveness of PAs in halting or slowing forest disturbances depends on PA characteristics such as size and history, as well as on the management entities such as subnational governments or indigenous peoples35,36,37. Moreover, there has been no attempt to elucidate whether PAs and OECMs are effective at regulating supply chains as a supply-side measure by balancing resource production, ecosystem services for local communities, and biodiversity conservation; to tackle this issue, it will be necessary to conduct extensive analyses integrating spatial and temporal patterns of biodiversity, forest loss, its drivers, and material flows in global food supply chains. Though it is challenging and beyond the scope of this paper, solving this issue is urgent and raises a promising opportunity for future research. More

Growth in harvests of crops meant for exports, processing and industrial use, together with their higher yields and faster yield gains, stands out globally; at a more granular level, this was driven by specific global regions that are getting increasingly specialized in harvesting crops for these usages.Changes in global-level harvested areasAt the global scale, we find that crops harvested for direct food utilization have the highest area and have been relatively stable over the study period (Fig. 1a). However, as the total harvested hectares have increased globally (Supplementary Table 1), this has translated into decreasing fractions of crops harvested for direct food utilization, from ~51% in the 1960s (average over 1964 to 1968) to ~37% in the 2010s (average over 2009 to 2013), with a similar reduction in feed crop harvests (Table 1). Conversely, there has been a substantial increase in crops for processing, exports and industrial use (Fig. 1a, Table 1 and Supplementary Table 1). The increase in industrial crop harvests occurred after year 2000. Around the same time, harvested hectares for exported crops ramped up and by the 2010s had surpassed those of crops harvested for feed use (Fig. 1a). Crops harvested for seed usage and losses are relatively minor, and we will not discuss them further. If the global trends observed in the past 20 years continue (Fig. 1a), by 2030, crops harvested for exports, processing and industrial use will account for ~ 23%, 17% and 8% of overall harvested hectares, whereas those for food will decrease to ~29% (Table 1).Fig. 1: Sector-based global crop-utilization trends.a–d, Observed total harvested ha (a), average yield in kcal ha−1 per year (b), average yield in protein ha−1 per year (c) and average yield in fat ha−1 per year (d) in the seven sectors of food, feed, processing, export, other uses (non-food/industrial), seed and losses from 1964 to 2013, annually, and projections to 2030 based on the past 20 years. The shading shows the 90% confidence interval for the significant linear model projections.Full size imageTable 1 Sector-based global crop-utilization changesFull size tableChanges in global-level crop yieldsWe find that crops harvested for direct food usage generally have had lower yields than all other sectors at the global scale over the time period of the study (Fig. 1b–d). This is not a new phenomenon, as crops harvested for direct food utilization have always had lower yields relative to other sectors (Supplementary Table 1). What has changed, however, is the ramping up (steeper positive slopes) of industrial, export and processing crop yields (Fig. 1b–d and Table 1). At these rates, caloric yields of industrial-use crops could increase by 28% from the 2010s to 2030 compared with 24% and 21% yield increases of crops harvested for directly consumed food and for feed use (Fig. 1b). Given that caloric yields of industrial-use crops are already substantially higher than food and feed crops (2× and 1.4×, respectively, in the 2010s), the faster caloric yield increases for industrial-use crops will widen this gap (2.1× and 1.5×, respectively). Yield measurements in other units of protein and fat show similar results (Table 1, Fig. 1c,d and Supplementary Table 1).Changes in the spatial patterns of harvested areas and productionWithin country-level information on harvested areas and productivity based on utilization categories is required for developing more locally effective agricultural policies. Over the course of the study time period 1964 to 2013 (Fig. 2a,b and Supplementary Video 1), we find changes in all continents when spatially analysed at the grid-cell level, except for most parts of Africa. Even in Africa, there are locations with fractional reductions in food crop harvests over the study period, such as parts of Angola, Ghana, Nigeria and South Africa. Within these and other countries, the exact location, magnitude and direction of the change varies from one region to the next (that is, compare Fig. 2a with Fig. 2b).Fig. 2: Sector-based spatial changes in crop harvests.a,b, The fraction of a grid cell in one of seven categories—food, feed, processing, export, other (non-food/industrial use), seed and losses—in each period, 1964–1968 (a) and 2009–2013 (b).Full size imageCrops harvested for direct food utilization have been prevalent in Asia, though much has changed since the 1960s (Fig. 2a,b and Supplementary Video 1). In China, there appears to be an imaginary belt, north and west of which harvests of crops used as directly consumed food decreased between the 1960s (Fig. 2a) and 2010s (Fig. 2b), while those for other uses increased. This belt appears to roughly extend from the northern half of Jiangsu (a province on the Yellow Sea in the east), curving westwards and southwards through northern Anhui, southern Henan, central Hubei and the northern tip of Hunan, and then turning sharply south and splitting Guangdong (a province on the South China Sea) through the middle. The sector gaining from the 10–20% fractional food harvest reduction varies. The increase in crops for feed, processing and industrial usage increases as one moves northward, especially north of Jiangsu and Anhui (Fig. 2a,b and Supplementary Video 1).Similarly, in India, there is a north–south zone encompassing eastern Haryana in the north, moving southwards through eastern Rajasthan, western Madhya Pradesh to eastern Maharashtra in the south, where there was a drastic reduction in crops harvested for direct food utilization over the study period (Fig. 2 and Supplementary Video 1); crops harvested for processing primarily increased. Changes in South and Southeast Asia over the study period are primarily away from once-dominant harvests of directly consumed food crops to feed crops, followed by processing crops, export crops and industrial-use crops, as in Myanmar and Thailand. In Malaysia, the growth was in export and industrial-usage crops, whereas in Indonesia, it was export crops and smaller increases in industrial-utilization crops. Central Asian states, especially Kazakhstan and some parts of Russia, witnessed a large reduction in crops harvested for direct food use over the study period, replaced by the crops destined for exports between the two periods (Fig. 2 and Supplementary Video 1).In Australia in the 1960s, food crops were harvested everywhere, accounting for ~10% of the total, which declined to ~5% by the 2010s. This was accompanied by small reductions in crops harvested for feed and export and balanced mainly by increases in crops for processing and industrial utilization (Fig. 2 and Supplementary Video 1).In Europe in the 1960s, crops were dominantly harvested for food and feed, but by the 2010s, this changed to include crops harvested for processing (Fig. 2 and Supplementary Video 1). In France, major reductions in feed crops have been balanced by growth in processing, export and industrial-use crops. In Spain, the primary change is from crops harvested for direct food to those of feed. In Germany, crops harvested for export have replaced those for direct food utilization.Latin America used to dominantly harvest food crops (as in Mexico) or food and feed crops (as in Brazil and Argentina) (Fig. 2 and Supplementary Video 1). Midwestern Brazil used to harvest only food crops, and feed and processing crop harvests were restricted to the Atlantic states (the 1960s; Fig. 2a), but by the 2010s (Fig. 2b), harvests of food crops had become a negligible fraction in Midwestern Brazil (as in Mato Grosso), and crops harvested for processing and exports are dominant now. In the Atlantic states of Brazil, one of the major changes is the increased proportion of harvests for industrial crops. In Argentina, over the study period, the proportion of crops harvested for food and feed has decreased, and this utilization has been mainly replaced by crops harvested for processing; crops harvested for exports changed, but the direction of change was spatially heterogeneous across Argentina (Fig. 2 and Supplementary Video 1). In Mexico, the primary change is the reduction in the fraction of crops harvested for direct food consumption and the increased harvests of crops for feed.Crops harvested for food and feed are also on the decline proportionally in North America. The United States has experienced a change from the dominance of food and feed crops in the 1960s to processing and industrial-usage crops in the 2010s. Detailed changes in the United States and Canada vary from one location to the next (Fig. 2), though the major change is the lower fraction of crops harvested for direct food consumption.Results are similar when viewed through the lens of calories, protein and fat with local-level differences as yields vary based on the measurement units (Supplementary Fig. 1). Further dramatic changes can be expected if observed linear trends from 1994 to 2013 at each grid cell continued until 2030 (Supplementary Fig. 2).Calories harvested in 2030 and achieving UN SDG 2We compare the extra food calories that will potentially be harvested in 2030 (Fig. 3a and Supplementary Data 2) to those required for both the projected extra population and feeding the projected undernourished population in each country (Fig. 3b and Supplementary Data 2). As an extreme case, we also compared whether total calories (all seven utilization sectors) would be sufficient (Fig. 3c and Supplementary Data 2). Altogether, we evaluated 156 countries, of which 86 had reported undernourished populations (Supplementary Data 2). On the basis of the minimum dietary energy requirement (MDER), we find that countries with reported undernourished populations will have a shortfall of ~675.4 trillion kcal per year to nourish the increased population and the expected undernourished from their extra harvested food calories. However, compared with the more realistic average dietary energy requirement (ADER), this shortfall will be ~993.9 trillion kcal per year (or ~70% from requirements) in 2030 (15 additional scenarios of undernourished populations in 2030 (provided in Supplementary Data 3) show global calorie shortfalls may similarly range from ~587.2 trillion kcal per year to ~1,269.3 trillion kcal per year based on the MDER level of nutrition requirement, and ~880.7 trillion kcal per year to ~1,755.6 trillion kcal per year based on the more realistic ADER level of nutrition requirement in 2030).Fig. 3: Meeting UN SDG goal 2 in 2030.a, Same as Fig. 2 but for the projected kcal ha−1 per year in 2030 per utilization sector and then mapping the fraction of total kcal ha−1 per year projected as harvested. b, Shortfall or gap from kcal per year harvested in 2030 as crops for direct food use and those to plug the gap from population growth and/or undernourished population. Computed based on the 2018 to 2020 ADER number for the country. c, Same as b but the kcal per year harvested used for computation is the total across all the seven sectors and shortfall is from whether the total calories harvested were used for direct food consumption (little to no processing).Full size imageCountries reporting undernourishment can, however, meet their requirement of extra calories in 2030 for both population change and those for the undernourished if calories from other utilization sectors are diverted and consumed directly as food calories (Fig. 3c and Supplementary Data 2 and 3). Though at the global scale, it appears that countries with high levels of undernourishment in 2030 can divert just a portion of their total harvested calories and meet some of the requirements of UN’s SDG 2 (ref. 4). In reality, many of the individual countries concentrated in sub-Saharan Africa have limited scope of diversion of calories from other sectors such as feed, processing or exports as crops for direct food use, as they already harvest most crops for direct food consumption (Fig. 2 and Supplementary Figs. 1 and 2). As such, many countries in this region may see deepening reliance on food imports. Note that the UN’s second SDG goal is broader in scope, including efforts to end malnutrition and increase agricultural productivity, among other goals4. Reconfiguration planning19 can use our spatially detailed information (Figs. 2 and 3, Supplementary Figs. 1 and 2 and Supplementary Data 2 and 3) in conjunction with policies that incentivize increased food crop harvests globally and ensure their equitable distribution to undernourished regions when local production is not sufficient20,21. This will require supply chain management22,23 and detailed analysis of optimization scenarios24 with our maps and tables as an important step linking specific production regions with the initial use of that production. More

Virus propagation and enumerationEchovirus-11 (E11, Gregory strain, ATCC VR737) and Coxsackievirus-A9 (CVA9, environmental strain from sewage, kindly provided by the Finnish National Institute for Health and Welfare) stocks were produced by infecting sub-confluent monolayers of BGMK cells as described previously [7]. Viruses were released from infected cells by freezing and thawing the culture flasks three times. To eliminate cell debris, the suspensions were centrifuged at 3000 × g for 5 min. Each stock solution was stored at −20 °C until use. Infectious virus concentrations were enumerated by a most probable number (MPN) infectivity assay as described in the Supplementary Information. The assay limit of detection (LoD), defined as the concentration corresponding to one positive cytopathic effect in the lowest dilution of the MPN assay under the experimental conditions used, corresponding to 2 MPN/mL.Inactivation of enteroviruses by bacterial consortia from lake waterTo study the inactivation of CVA9 and E11 by a bacterial consortium from lake water, four surface water samples were collected from Lake Geneva (Ecublens, Switzerland) during the summer 2021. Each sampling event was conducted on warm and sunny days, to minimize biological variation. Immediately after sampling, large particles of the sample were removed by filtering 500 mL of water on a 8 μm nitrocellulose filter membrane (Merck Millipore, Cork, Ireland). The sample was then filtered through a 0.8 μm nitrocellulose filter membrane (Merck Millipore) to remove large microorganisms such as protists. The resulting water sample corresponds to the bacterial fraction used to study virus inactivation.For inactivation experiments, each virus was spiked into individual 1 mL aliquots of fractionated lake water to a final concentration of 106 MPN/mL, and samples were incubated for 48 h at 30 °C without shaking. Duplicate experiments were conducted for each virus and each lake water sample. Experiments to control for thermal inactivation were conducted using the same procedure but by replacing the fractionated lake water with sterile milliQ water. Viral infectivity at times 0 h and 48 h was determined by MPN as described above. Virus decay was calculated as log10 (C/C0), where C is the residual titer after 48 h of incubation, and C0 is the initial titer. The experimental LoD was approximately 5-log10.These same experiments were conducted for three new water samples in the presence of four protease inhibitors with the following final concentrations: E64—10 μM (E3132, Sigma–Aldrich, Saint-Louis, MO, USA), GM6001—4 μM (CC1010, Sigma–Aldrich), Chymostatin—100 μM (C7268, Sigma–Aldrich), and PMSF—100 μM (P7626, Sigma–Aldrich). Each inhibitor was added to 1 mL of fractionated lake water, vortexed for 30 seconds, and incubated at room temperature for 15 min, before adding the two viral strains under the same conditions as described above.Bacterial isolation, cultivation, and storageBacteria were isolated from two water samples from Lake Geneva’s Ecublens beach, taken in November 2019 (Fall, 89 isolates) and May 2020 (Spring, 47 isolates). Bacteria recovery was performed on R2A agar plate (BD Difco, Franklin Lakes, NJ, USA) as described previously [15]. Briefly, successive dilutions from 10−1 to 10−5 were carried out in sterile water for each sample. For each dilution, a volume of 1 mL was deposited on three separate R2A plates, before being incubated at 22, 30, and 37 °C. After 5 days of incubation, each colony was picked and enriched on a new R2A plate. To ensure purity, each isolate was successively plated five times on R2A plate and incubated at the same temperature as the initial isolation. Each purified isolate was cryopreserved in R2A / 20% glycerol at −80 °C. The isolates were named based on the water body (Lake (L)), isolation temperature, and the isolation order (L-T°C-number).Bacterial identificationThe identification of each isolate was performed by 16 S rRNA gene sequencing using the pair of primers 27 F (5’- AGA GTT TGA TCM TGG CTC AG- 3’, Microsynth AG, Balgach, Switzerland) / 786 R (5’- CTA CCA GGG TAT CTA ATC – 3’, Microsynth AG), following a methodology previously described [15]. The thermocycling conditions and the purification of PCR products are described in the Supplementary Information. The complete list of isolated bacteria and associated accession numbers is given in Supplementary Table 1.Phylogenetic inference and metadata visualizationThe consensus from 16 S rRNA gene sequences of the 136 isolates was aligned using the MUSCLE algorithm [16]. The phylogenetic analysis of 566 bp aligned sequences from the V2-V4 16 S rRNA gene regions (Positions: 152–717) was performed using Molecular Evolutionary Genetics Analysis X software [17]. Phylogeny was inferred by maximum likelihood, with 1000 bootstrap iterations to test the robustness of the nodes. The resulting tree was uploaded and formatted using iTOL [18].Virus incubation with bacterial isolatesFor the preparation of the bacteria before co-incubation, each one was first cultured on R2A agar for 48 h at their initial isolation temperature. Overnight suspensions of each bacterial isolate were grown in R2A broth at room temperature under constant agitation (180 rpm). For co-incubation experiments, 200 μL of each bacterial suspension were mixed with 100 μL of a 105 MPN/mL stock of E11 or CVA9. Then, each condition was supplemented with 600 μL of R2A broth. Incubation was carried out for 96 h at room temperature, without shaking. At the end of the co-incubation, each tube was centrifuged for 15 min at 9000 × g (4 °C) to eliminate bacteria, and the residual infectious viral titer was enumerated by MPN assay as described above [7]. Each co-incubation experiment was carried out in triplicate. Control experiments were performed under the same conditions but using sterile R2A. Virus decay was quantified as log10 (Cexp/Cctrl), where Cexp is the residual titer after a co-incubation for 96 h, and Cctrl is the titer after incubation of the virus in sterile R2A for 96 h. The experimental LoD was 3-log10.Protease activity measurement using casein and gelatin agar platesCasein agar was prepared as follows: 20 g of skim milk (BD Difco), supplemented with 1 g glucose were reconstituted with 200 mL of distilled water. Likewise, a 10% bacteriological agar solution was prepared in a final volume of 200 mL. Finally, a solution consisting of 0.8% NaCl, 0.02% KCl, 0.144% Na2HPO4, and 0.024% KH2PO4 was reconstituted in 600 mL of water. All solutions were autoclaved for 15 min at 110 °C. The solutions were mixed, and 25 mL were poured into each Petri dish. Gelatin agar was composed of 0.4% peptone, 0.1% yeast extract, 1.5% gelatin and 1.5% bacteriological agar. The mixture was autoclaved 15 min at 120 °C, and 25 mL of medium was poured into each Petri dish.For each isolate, an overnight suspension was performed in R2A broth at room temperature, before spotting 15 μL of each suspension at the center of both gelatin and casein agar plates. Each plate was incubated at 22, 30, or 37 °C for 72 h, depending on the initial isolation temperature of the bacteria. Casein-degrading activity (cas), which is exerted by many different protease classes, and gelatin-degrading activity (gel), which is mostly caused by MMPs, were revealed by a hydrolysis halo around the producing bacteria. Hydrolysis diameters were measured in millimeters (mm) to report the extent of the proteolytic effect of each strain on both substrates.Protease activity quantification in cell-free supernatantUsing the same bacterial suspensions as for bacterial/virus co-incubation, 200 μL of each suspension was inoculated into 600 μL of R2A broth and incubated without shaking for 96 h at room temperature. Each culture was centrifuged for 15 min at 9000 × g at 4 °C. The resulting cell-free supernatants (CFS) were stored at −20 °C until use. For each CFS, protease activity was measured using the Protease Activity Assay Kit (ab112152, Abcam, Cambridge, UK), which measures general protease activity (pgen) except MMPs, and the MMP Activity Assay Kit (ab112146, Abcam), which selectively measures MMP activity (mmp). Briefly, for the Protease Activity Assay kit, 50 μL of the substrate was added into each well of a dark-bottom plate containing 50 μL of each CFS. Standard trypsin provided by the kit was used as a positive control. For the MMP Activity Assay kit, 50 μL of each CFS was incubated with 50 μL of 2 mM APMA for 3 h at 37 °C, prior to the activity test. Collagenase I (C0130, Sigma–Aldrich) was used as a positive control. R2A broth was used as a negative control for each assay. Protease activity was measured at time 0 and after 60 min, using a Synergy MX fluorescence reader (BioTek). The excitation and emission wavelengths were set to 485 and 530 nm, respectively. The emitted fluorescence, generated by proteolytic cleavage of the substrate of each kit, was calculated as follows: ∆RFU = RFU (60 min) − RFU (0 min). Proteolytic activity was calculated in mmol/min/μL based on the emitted fluorescence measured for trypsin and collagenase I at known proteolytic activities.Data analysisStatistical analyses to compare inactivation data were performed by one-way t-test or one-way ANOVA with Dunnett’s post-hoc test in GraphPad Prism v.9. An alpha value of 0.05 was used as a threshold for statistical significance. For each dataset we confirmed that data were normally distributed.To analyze a potential correlation between protease activity and viral decay, the decay values for each virus strain was related to the four protease activity tests of this study using a scatterplot combined with a Kernel density estimation. The analyses were performed with R v.3.6.1 using the SmoothScatter function of the R Base package.A Left-Censored Tobit model (CTM) with mixed effects was chosen to investigate interactions between protease activity and the decay measured for each virus strain. Briefly, the CTM with mixed effect was chosen for three reasons: (1) The protocol used to measure viral decay had a limit of quantification of −3-log10, and 152 measurement points reached the detection limit, requiring the use of this value as the left-censored value of the model; (2) The two virus strains used in the study showed distinct responses after exposure to environmental bacteria, preventing the use of a multiple linear regression model; (3) Among biological replicates of co-incubation experiments, inactivation variability was observed, suggesting the concomitant action of random biological effects (e.g., production of other compounds than proteases by bacteria, or differences in protease production rate between replicates for each bacterial isolate). The resulting statistical model was then formulated as follows:$$log left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right) = ; beta _0 + beta _1;{rm I}_{{{{{{{{mathrm{virus}}}}}}}}_i = 2} + beta _2sqrt {left[ {pgen} right]_i} + beta _3sqrt {left[ {mmp} right]_i} + beta _4sqrt {left[ {cas} right]_i} \ + beta _5sqrt {left[ {gel} right]_i} + beta _6I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {pgen} right]_i} + beta _7I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {mmp} right]_i} \ + beta _8I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {cas} right]_i} + beta _9I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {gel} right]_i} + alpha _{{{{{{{{mathrm{id}}}}}}}}_i} + varepsilon _i$$$${{{mbox{where}}}}; log left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right) = left{ {begin{array}{*{20}{c}} { – 3} & {{{{{{{{mathrm{if}}}}}}}};{{{{{{{mathrm{log}}}}}}}}left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right) le – 3} \ {{{{{{{{mathrm{log}}}}}}}}left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right)} & {{{{{{{{mathrm{otherwise}}}}}}}}} end{array}} right.$$$$alpha _{{{{{{{{mathrm{id}}}}}}}}_i}sim {{{{{{{mathrm{i}}}}}}}}.{{{{{{{mathrm{i}}}}}}}}.;{{{{{{{mathrm{d}}}}}}}}.;{rm N}left( {0,;sigma _{{{{{{{{mathrm{id}}}}}}}}}^2} right)$$$${{{{{{{mathrm{for}}}}}}}};i in left{ {1,2, ldots } right}$$for which β0 defines the model intercept, (beta _1{rm I}_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}) corresponds to the main effect of the virus factor on the viral decay, (beta _2,;beta _3,;beta _4,;{{{{{{{mathrm{and}}}}}}}};beta _5) corresponds to the main effects of the different protease activity measurements on viral decay, (beta _6I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2},;beta _7I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2},;beta _8I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2},{{{{{{{mathrm{and}}}}}}}};beta _9I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}) corresponds to the interaction effects between each of these variables and the viral decay, (alpha _{{{{{{{{mathrm{id}}}}}}}}_i}) corresponds to the mixed effect of the model and (varepsilon _i) corresponds to the error term of the model. The selection of the model is further detailed in the Supplementary Information (Supplementary Material and Figs. S1 and S2).The full dataset included in the correlation analysis and the CTM is provided in Supplementary Table 2. A description of the variables used is given in the Supplementary Information. The dataset was analyzed using the censReg package in R [19]. The R code is given in the Supplementary Information. More