Ecology

Mace, G. M., Norris, K. & Fitter, A. H. Biodiversity and ecosystem services: a multilayered relationship. Trends Ecol. Evol. 27, 19–26 (2012).Article 

Google Scholar 
Global Assessment Report on Biodiversity and Ecosystem Services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES, 2019).Williams, B. A. et al. A robust goal is needed for species in the Post-2020 Global Biodiversity Framework. Conserv. Lett. 14, e12778 (2021).Article 

Google Scholar 
Rodrigues, A. S. L. et al. Effectiveness of the global protected area network in representing species diversity. Nature 428, 640–643 (2004).CAS 
Article 

Google Scholar 
Watson, J. E. M., Dudley, N., Segan, D. B. & Hockings, M. The performance and potential of protected areas. Nature 515, 67–73 (2014).CAS 
Article 

Google Scholar 
McCarthy, D. P. et al. Financial costs of meeting global biodiversity conservation targets: current spending and unmet needs. Science 338, 946–949 (2012).CAS 
Article 

Google Scholar 
Pe’er, G. et al. Action needed for the EU Common Agricultural Policy to address sustainability challenges. People Nat. 2, 305–316 (2020).Article 

Google Scholar 
McDonald, R. I. et al. Research gaps in knowledge of the impact of urban growth on biodiversity. Nat. Sustain. 3, 16–24 (2020).Article 

Google Scholar 
Rosenzweig, M. L. Reconciliation ecology and the future of species diversity. Oryx 37, 194–205 (2003).Article 

Google Scholar 
Dunn, R. R., Gavin, M. C., Sanchez, M. C. & Solomon, J. N. The pigeon paradox: dependence of global conservation on urban nature. Conserv. Biol. 20, 1814–1816 (2006).Article 

Google Scholar 
Callaghan, C. T. et al. How to build a biodiverse city: environmental determinants of bird diversity within and among 1581 cities. Biodivers. Conserv. 30, 217–234 (2021).Article 

Google Scholar 
Ives, C. D. et al. Cities are hotspots for threatened species. Glob. Ecol. Biogeogr. 25, 117–126 (2016).Article 

Google Scholar 
Soanes, K. & Lentini, P. E. When cities are the last chance for saving species. Front. Ecol. Environ. 17, 225–231 (2019).Article 

Google Scholar 
Luck, G. W., Davidson, P., Boxall, D. & Smallbone, L. Relations between urban bird and plant communities and human well-being and connection to nature. Conserv. Biol. 25, 816–826 (2011).Article 

Google Scholar 
Maller, C., Mumaw, L. & Cooke, B. in Rewilding (eds Pettorelli, N. et al.) Ch. 9 (Cambridge, Univ. Press, 2019).Jiang, L. & O’Neill, B. C. Global urbanization projections for the Shared Socioeconomic Pathways. Glob. Environ. Change 42, 193–199 (2017).Article 

Google Scholar 
Prévot, A.-C., Cheval, H., Raymond, R. & Cosquer, A. Routine experiences of nature in cities can increase personal commitment toward biodiversity conservation. Biol. Conserv. 226, 1–8 (2018).Article 

Google Scholar 
Berthon, K., Thomas, F. & Bekessy, S. The role of ‘nativeness’ in urban greening to support animal biodiversity. Landsc. Urban Plan. 205, 103959 (2021).Article 

Google Scholar 
van Heezik, Y., Freeman, C., Davidson, K. & Lewis, B. Uptake and engagement of activities to promote native species in private gardens. Environ. Manag. 66, 42–55 (2020).Article 

Google Scholar 
Jorgensen, A. & Keenan, R. Urban Wildscapes (Routledge, 2012).Majewska, A. A. & Altizer, S. Planting gardens to support insect pollinators. Conserv. Biol. 34, 15–25 (2020).Article 

Google Scholar 
Tallamy, D. W. Bringing Nature Home: How You Can Sustain Wildlife with Native Plants (Timber Press, 2007).Burghardt, K. T., Tallamy, D. W. & Gregory Shriver, W. Impact of native plants on bird and butterfly biodiversity in suburban landscapes. Conserv. Biol. 23, 219–224 (2009).Article 

Google Scholar 
Kraljevic, A. & Mitlacher, G. Barometer on CBD’s Target for International Resource Mobilization (WWF, 2020).Eichenberg, D. et al. Widespread decline in Central European plant diversity across six decades. Glob. Change Biol. 27, 1097–1110 (2020).Article 

Google Scholar 
Lughadha, E. N. et al. Extinction risk and threats to plants and fungi. Plants People Planet 2, 389–408 (2020).Article 

Google Scholar 
Metzing, D., Hofbauer, N., Ludwig, G. & Matzke-Hajek, G. Rote Liste Gefährdeter Tiere, Pflanzen und Pilze Deutschlands: Pflanzen/Redaktion: Detlev Metzing, Natalie Hofbauer, Gerhard Ludwig und Günter Matzke-Hajek (Bundesamt für Naturschutz, 2018).Kalusová, V. et al. Naturalization of European plants on other continents: the role of donor habitats. Proc. Natl Acad. Sci. USA 114, 13756–13761 (2017).Article 
CAS 

Google Scholar 
Staude, I. R. et al. Directional turnover towards larger-ranged plants over time and across habitats. Ecol. Lett. 25, 466–482 (2021).Article 

Google Scholar 
Galloway, J. N. et al. The nitrogen cascade. Bioscience 53, 341–356 (2003).Article 

Google Scholar 
Lundholm, J. T. & Richardson, P. J. Mini-Review: Habitat analogues for reconciliation ecology in urban and industrial environments. J. Appl. Ecol. 47, 966–975 (2010).Article 

Google Scholar 
Ellenberg, H. Gefahrdung wildlebender Pflanzenarten in der Bundesrepublik Deutschland, Versuch einer okologischen Betrachtung. Forstarchiv 57, 127–133 (1983).
Google Scholar 
Deeb, M. et al. Using constructed soils for green infrastructure—challenges and limitations. Soil 6, 413–434 (2020).CAS 
Article 

Google Scholar 
BuGG-Marktreport Gebäudegrün 2020 Dach-, Fassaden-und Innenraumbegrünung Deutschland Neu begrünte Flächen Bestand und Potenziale Kommunale Förderung (BuGG, 2020).Reichard, S. H. & White, P. Horticulture as a pathway of invasive plant introductions in the United States: most invasive plants have been introduced for horticultural use by nurseries, botanical gardens, and individuals. Bioscience 51, 103–113 (2001).Article 

Google Scholar 
der Lippe, M. & Kowarik, I. Do cities export biodiversity? Traffic as dispersal vector across urban—rural gradients. Divers. Distrib. 14, 18–25 (2008).Article 

Google Scholar 
Razgour, O. et al. Considering adaptive genetic variation in climate change vulnerability assessment reduces species range loss projections. Proc. Natl Acad. Sci. USA 116, 10418–10423 (2019).CAS 
Article 

Google Scholar 
Goddard, M. A., Dougill, A. J. & Benton, T. G. Scaling up from gardens: biodiversity conservation in urban environments. Trends Ecol. Evol. 25, 90–98 (2010).Article 

Google Scholar 
Sharrock, S. Plant Conservation Report 2020: A Review of Progress Towards the Global Strategy for Plant Conservation 2011–2020 CBD Technical Series No. 95 (Convention on Biological Diversity, 2020).Ismail, S. A., Pouteau, R., van Kleunen, M., Maurel, N. & Kueffer, C. Horticultural plant use as a so-far neglected pillar of ex situ conservation. Conserv. Lett. 14, e12825 (2021).Article 

Google Scholar 
Wüstemann, H., Kalisch, D. & Kolbe, J. Access to urban green space and environmental inequalities in Germany. Landsc. Urban Plan. 164, 124–131 (2017).Article 

Google Scholar 
Kleingärten im Wandel—Innovationen für verdichtete Räume (BBSR, 2018).Rudd, H., Vala, J. & Schaefer, V. Importance of backyard habitat in a comprehensive biodiversity conservation strategy: a connectivity analysis of urban green spaces. Restor. Ecol. 10, 368–375 (2002).Article 

Google Scholar 
Kowarik, I. & von der Lippe, M. Plant population success across urban ecosystems: a framework to inform biodiversity conservation in cities. J. Appl. Ecol. 55, 2354–2361 (2018).Article 

Google Scholar 
du Toit, M. J., Shackleton, C. M., Cilliers, S. S. & Davoren, E. in Urban Ecology in the Global South (eds Shackleton, C. M. et al.) 433–461 (Springer, 2021).Sawyer, J. Saving threatened native plant species in cities—from traffic islands to real islands. In Greening the City: Bringing Biodiversity Back Into the Urban Environment: Proc. (ed Dawson, M. I.) 111–117 (Royal New Zealand Institute of Horticulture, 2005).Webb, E. L. A Guide to the Native Ornamental Trees of American Samoa (National Univ. Singapore, 2011).Pan, K. et al. Urban green spaces as potential habitats for introducing a native endangered plant, Calycanthus chinensis. Urban For. Urban Green. 46, 126444 (2019).Article 

Google Scholar 
Gardening sales value worldwide from 2015 to 2020, with a forecast up to 2024. Statista statista.com/statistics/1220222/global-gardening-sales-value/ (2021).Warenstromanalyse 2018: Blumen, Zierpflanzen & Gehölze (AMI, 2020).Nature Awareness Study (Bundesamt für Naturschutz (BfN), 2019).Abbandonato, H., Pedrini, S., Pritchard, H. W., De Vitis, M. & Bonomi, C. Native seed trade of herbaceous species for restoration: a European policy perspective with global implications. Restor. Ecol. 26, 820–826 (2018).Article 

Google Scholar 
Hancock, N., Gibson-Roy, P., Driver, M. & Broadhurst, L. The Australian Native Seed Survey Report (Australian Network for Plant Conservation, 2020).Wilkinson, D. M. Is local provenance important in habitat creation? J. Appl. Ecol. 38, 1371–1373 (2001).Article 

Google Scholar 
Pedrini, S. & Dixon, K. W. International principles and standards for native seeds in ecological restoration. Restor. Ecol. 28, S286–S303 (2020).
Google Scholar 
De Vitis, M. et al. The European native seed industry: characterization and perspectives in grassland restoration. Sustainability 9, 1682 (2017).Article 

Google Scholar 
Westwood, M., Cavender, N., Meyer, A. & Smith, P. Botanic garden solutions to the plant extinction crisis. Plants People Planet 3, 22–32 (2021).Article 

Google Scholar 
Mounce, R., Smith, P. & Brockington, S. Ex situ conservation of plant diversity in the world’s botanic gardens. Nat. Plants 3, 795–802 (2017).Article 

Google Scholar 
Pedrini, S. et al. Collection and production of native seeds for ecological restoration. Restor. Ecol. 28, S228–S238 (2020).
Google Scholar 
Groves, R. H. Can Australian native plants be weeds. Plant Prot. Q. 16, 114–117 (2001).
Google Scholar 
Brummitt, R. K., Pando, F., Hollis, S. & Brummitt, N. A. World Geographic Scheme for Recording Plant Distributions 2nd edn (Hunt Institute for Botanical Documentation, 2001).Davis, M. A. et al. Don’t judge species on their origins. Nature 474, 153–154 (2011).CAS 
Article 

Google Scholar 
Mumaw, L. & Bekessy, S. Wildlife gardening for collaborative public—private biodiversity conservation. Australas. J. Environ. Manag. 24, 242–260 (2017).Article 

Google Scholar 
Mumaw, L. Transforming urban gardeners into land stewards. J. Environ. Psychol. 52, 92–103 (2017).Article 

Google Scholar 
Abeli, T. et al. Ex situ collections and their potential for the restoration of extinct plants. Conserv. Biol. 34, 303–313 (2020).Article 

Google Scholar 
Ladouceur, E. et al. Native seed supply and the restoration species pool. Conserv. Lett. 11, e12381 (2018).Article 

Google Scholar 
Hyvärinen, M.-T. Rubus humulifolius rescued by narrowest possible margin, conserved ex situ, and reintroduced in the wild. J. Nat. Conserv. 55, 125819 (2020).Article 

Google Scholar 
Holz, H., Segar, J., Valdez, J. & Staude, I. R. Assessing extinction risk across the geographic ranges of plant species in Europe. Plants People Planet https://doi.org/10.1002/ppp3.10251 (2022).Brodie, J. F. et al. Global policy for assisted colonization of species. Science 372, 456–458 (2021).CAS 
Article 

Google Scholar 
Gann, G. D. et al. International principles and standards for the practice of ecological restoration. Restor. Ecol. 27, S1–S46 (2019).Article 

Google Scholar 
Bower, A. D., Clair, J. B. S. & Erickson, V. Generalized provisional seed zones for native plants. Ecol. Appl. 24, 913–919 (2014).Article 

Google Scholar 
Goddard, M. A., Dougill, A. J. & Benton, T. G. Why garden for wildlife? Social and ecological drivers, motivations and barriers for biodiversity management in residential landscapes. Ecol. Econ. 86, 258–273 (2013).Article 

Google Scholar 
Ignatieva, M. & Ahrné, K. Biodiverse green infrastructure for the 21st century: from “green desert” of lawns to biophilic cities. J. Archit. Urban. 37, 1–9 (2013).Article 

Google Scholar 
van Heezik, Y. M., Dickinson, K. J. M. & Freeman, C. Closing the gap: communicating to change gardening practices in support of native biodiversity in urban private gardens. Ecol. Soc. 17, 34 (2012).
Google Scholar 
Shaw, A. E. & Miller, K. K. Preaching to the converted? Designing wildlife gardening programs to engage the unengaged. Appl. Environ. Educ. Commun. 15, 214–224 (2016).Article 

Google Scholar 
Mumaw, L. M. & Raymond, C. M. A framework for catalysing the rapid scaling of urban biodiversity stewardship programs. J. Environ. Manag. 292, 112745 (2021).Article 

Google Scholar 
Niemiec, R., Jones, M. S., Lischka, S. & Champine, V. Efficacy-based and normative interventions for facilitating the diffusion of conservation behavior through social networks. Conserv. Biol. 35, 1073–1085 (2021).Article 

Google Scholar 
Haywood, B. K., Parrish, J. K. & Dolliver, J. Place-based and data-rich citizen science as a precursor for conservation action. Conserv. Biol. 30, 476–486 (2016).Article 

Google Scholar 
Lerman, S. B., Turner, V. K. & Bang, C. Homeowner associations as a vehicle for promoting native urban biodiversity. Ecol. Soc. 17, 45 (2012).Article 

Google Scholar 
Nassauer, J. I. Messy ecosystems, orderly frames. Landsc. J. 14, 161–170 (1995).Article 

Google Scholar 
Cavender, N., Smith, P. & Marfleet, K. BGCI Technical Review: The Role of Botanic Gardens in Urban Greening and Conserving Urban Biodiversity (BGCI, 2019). More

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

Study areaHarbin is located in the centre of Northeast Asia, between 44°04’–46° 40′ N and 125° 42′–130° 10′ E24,26. The site has a mid-temperate continental monsoon climate, with an average annual temperature of 3.6° C and an average annual precipitation is 569.1 mm. The main precipitation months being from June to September, accounting for about 60% of the annual precipitation, the main snow months are from November to January24,25. The overall topography is high in the east and low in the west, with mountains and hills predominating in the east and plains predominating in the west27. In this study, we identified the central district of Harbin, where urban construction activities are frequent and the population is dense, as the study area. According to the “Harbin City Urban Master Plan (2011–2020)” (revised draft in 2017), the specific scope includes Daoli District, Daowai District, Nangang District, Xiangfang District, Pingfang District, Songbei District’s administrative district, Hulan District, and Acheng District part of the area, with a total area of 4187 km2 (Fig. 2). The blue–green space in this study included woodland, grassland, cultivated land, wetland and water that permeate inside and outside the construction sites. They all have integrated functions such as ecology, supply, beautification, culture, and disaster prevention and avoidance, and have a decisive influence on the urban ecological environment.Figure 2Schematic of study area. The Figure is created using ArcGIS ver.10.2 (https://www.esri.com/).Full size imageData sourcesThe data used in this research included the following: land-cover date (30 m × 30 m) of two periods (2011, 2020) spported by the China Geographic National Conditions Data Cloud Platform (http://www.dsac.cn/), Meteorological datasets (1 km × 1 km) were obtained from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http:∥www.resdc.cn/), including air temperature, precipitation, and surface runoff. ASTER GDFM elevation data (30 m × 30 m) came from the Geospatial Data Cloud (http:∥www.gscloud.cn), from which the slope was extracted. Soil data (1 km × 1 km) were from the World Soil Database (HWSD) China Soil Data Set (v1.1). The normalized difference vegetation index (NDVI) and modified normalized difference water index (MNDWI) data (30 m × 30 m) came from the National Comprehensive Earth Observation Data Sharing Platform (http://www.chinageoss.org/), ET datasets (30 m × 30 m) were drawn from the NASA-USGS (https://lpdaac.usgs.gov/). Social and economic data were mainly obtained through the Harbin statistical yearbook and the Harbin social and economic bulletin.Framework of urban blue–green space LEH assessmentUrban blue–green space is a politically defined man-land coupling region composed of ecological, economic, and social systems, which is greatly disturbed by human activities11. The essence of urban blue–green space LEH is that the landscape ecological function sustainably meets human needs28,29. The landscape ecological function reflects the value orientation of human beings to blue–green space, and to a large extent affects the blue–green landscape ecological pattern and process. The interaction between the blue–green landscape ecological pattern and process drives the overall dynamics of blue–green space. Meanwhile, presenting certain landscape ecological function characteristics, which provide ecological support for various human activities30,31,32. While the pattern and process of blue–green space both profoundly influence and are influenced by human activities33,34. This influence is long-term, the standard of LEH should not be fixed in real-time health, but should fully consider the sustainability of the health state.In summary, the landscape ecological pattern, process, function, and sustainability are not separate, but a complex of mutual integration, and organic unity. In this study, we constructed an integrated assessment framework of blue–green space LEH that included four units: pattern, process, service, and sustainability (Fig. 3). In the assessment framework, the LEH of urban blue–green space involves two dimensions: the first is the health status of the urban blue–green space itself, emphasizing the maintenance of the ecological conditions, thereby potentially satisfying a series of diversity goals. The other is that urban blue–green space, as a part of social and economic development, could sustainably provide the ability to meet (subject) needs and goals.Figure 3Key units, interactions of urban blue–green space LEH.Full size imageLandscape ecological patternThe landscape ecological pattern of urban blue–green space is a spatial mosaic combination of landscape elements at different levels or the same level. Affected by human activities interference31, the landscape ecological pattern shows the changing trend of landscape structure complexity, landscape type diversification, and landscape fragmentation. The assessment of urban landscape ecological pattern should be a comprehensive reflection of this changing trend1. Landscape pattern indexes are the most frequently applied which could reflect the structural composition and spatial configuration characteristics of the landscape4,35. This study took landscape ecology as the entry point and selected the landscape pattern indexes that can quantitatively reflect the change characteristics of landscape structural composition and spatial configuration under the disturbance. In this way, the landscape disturbance index (U), landscape connectivity index (CON), and landscape adaptability index (LAI) were used as the indexes for the assessment of landscape ecological pattern health.

(1)

Landscape disturbance index (U)

There are two kinds of relationships between the landscape ecological pattern and the external disturbance: compatibility and conflict. As the landscape ecological pattern has accommodating characteristics, the disturbance beyond the accommodating capacity will degrade the landscape ecological pattern36,37. The landscape disturbance index (U) could characterize the degree of fragmentation, dispersion, and morphological changes in landscape pattern38. The index is a comprehensive index that can reflect the health of the landscape pattern by quantifying the ability of ecosystems to accommodate external disturbances. It consists of the landscape fragmentation index, the inverse of the fractional dimension, and the dominance index. They measure the response of the landscape pattern to external disturbance from the perspective of different landscape types, the same landscape type, and landscape diversity, respectively36,38, and their weights were determined by the entropy weight method. The formula is as follows:$$ U = alpha N_{{{Fi}}} + bD_{{{Fi}}} + cD_{{{Oi}}} $$
(1)
where NFi is the landscape fragmentation index, DFi is the inverse of the fractional dimension, DOi is the dominance index, and a, b, and c are the corresponding weights, which were 0.20, 0.5, and 0.3 in this study, respectively.

(2)

Landscape connectivity index (CON)

The most direct result of landscape ecological pattern degradation caused by external disturbance is that the flow of energy, material, and information among ecological patches is reduced or even blocked, ultimately the stability of the landscape pattern is decreased. The connectivity could characterize the ability of landscape ecological pattern to mitigate risk transmission, which is significant for the dynamic stability of landscape ecological pattern39,40. The landscape connectivity index (CON) could measure the connectivity between ecosystem components through the aggregation or dispersion trend of patches41. The better the connectivity, the stronger the stability of landscape ecological pattern. The formula is as follows:$$ CON = frac{{100sumlimits_{s = 1}^{q} {sumlimits_{h ne l}^{p} {C_{{{shl}}} } } }}{{sumlimits_{s = 1}^{s} {left[ {q_{{s}} (q_{{s}} – 1)/2} right]} }} $$
(2)
where qs is the number of plaques of patch type s, Cshl is the link between patch h and patch l in s within the delimited distance.

(3)

Landscape Restorability Index (LRI)

The ability to recover to its original structure when subjected to disturbances is an important criterion for the landscape ecological pattern42. Research confirmed that the restorability of the landscape ecological pattern is closely related to the structure, function, diversity, and uniformity of distribution. The landscape restorability index (LRI) combines the above landscape information and could indicate the restorability of the landscape ecological pattern in response to disturbance43. The index consists of the patch density, Shannon diversity index, and the landscape evenness, the patch density is the number of patches per square kilometer. The Shannon diversity index reflects the change in the proportion of landscape types. The landscape evenness index shows the distribution evenness of patches in terms of area. The larger the LRI index, the more complex and evenly distributed the structure is, and the more recovery ability of the landscape pattern against disturbance is. The formula is as follows:$$ LRI = PD times SHDI times SHEI $$
(3)
where PD is the patch density, SHDI is the Shannon diversity index, and SHEI is the landscape evenness index.Landscape ecological processThe landscape ecological process of urban blue–green space is extremely complex for it involves multiple factors such as natural ecology, economy, and culture. Landscape ecological process assessment is the measure of the self-organized capacity and the efficiency of ecological processes within and among patches44. A blue–green space with a healthy landscape ecological process should have the ability to adapt to conventional land use under human management and maintain physiological integrity while maintaining the balance of ecological components. Specifically, the landscape ecological process could quickly restore its balance after severe disturbances, with strong organization, suitability, recoverability, and low sensitivity45,46. A single model hardly to gets good research on landscape ecological process under the urban scale. The comprehensive application of multidisciplinary methods is effective means to solve the problem. Regarding this, we selected ecological indexes and models from four aspects: organization, suitability, restoration, and sensitivity to assess the landscape ecological process of urban blue–green space.

(1)

Organization index (O)

The organization of the landscape ecological process is the maintenance ability of stable and orderly material cycling and energy flow within and between landscapes47. The normalized vegetation index (NDVI) and the modified normalized difference water index (MNDWI) could reflect the efficiency and order of ecological processes. Such as accumulation of organic matter, fixation of solar energy, nutrient cycling, regeneration, and metabolism13. The indexes are the external performance of the internal dynamics and organizational capabilities of the ecological process. In recent years, it has been widely used in the assessment of related to landscape ecological process. The formulas are as follows:$$ NDVI = frac{NIR – R}{{NIR + R}} $$$$ MNDWI = frac{p(green) – p(MIR)}{{p(green) + p(MIR)}} $$
(4)
where (NDVI) is the normalized vegetation index, (MNDWI) is the modified water body index, (NIR) is the reflectance value in the near-infrared band, (R) is the reflectance value in the visible channel, (p(green)) and (p(MIR)) are the normalized values in the green and mid-infrared bands.

(2)

Suitability index (Q)

The suitability of the landscape ecological process is a measurement of the self-regulating ability of the landscape ecosystem. That is, to effectively maintain the ecological process in a state of being protected from disturbance during the occasional changes caused by the external environment2. The water conservation amount index (Q) can measure the operating capacity of ecosystems to maintain ecological balance, water conservation, climate regulation, and other ecological processes by integrating the water balance of rainfall, surface runoff, and evaporation41. It could reflect the suitability of landscape ecological process to regional environment and developmental conditions. The formula is as follows:$$ Q = R – J – ET $$
(5)
where Q is the water conservation amount, R is the annual rainfall, J is the surface runoff, ET is the evapotranspiration.

(3)

Recoverability index (ECO)

The recoverability of the landscape ecological process refers to the ability of an ecosystem to return to its original operating state after being subjected to external impacts. Land-use types play an essential role in landscape ecological recoverability48. The ecological recoverability index (ECO) uses the resilience coefficients of land-use types to reflect the level of ecosystem resilience38. Based on previous studies, the resilience coefficient of land-use types was assigned (Table 1).

(4)

Sensitivity index(A)

Table 1 Resilience coefficients of different land use types.Full size tableThe sensitivity index (A) could be used to indicate landscape ecological process formation, change, and vulnerability to disturbance31. We started from the physical effects of blue–green space on sand production, water confluence, and sediment transport, introduced the Soil Erosion Modulus to characterize the sensitivity of landscape ecological processes to disturbance. The index effectively combines landscape ecology, erosion mechanics, soil science, and sediment dynamics49. The formula is as follows:$$ begin{gathered} A = R_{{i}} cdot K cdot LS cdot C cdot P hfill \ L = (l/22.1)^{m} hfill \ S = left{ begin{gathered} 10.8sin theta + 0.03,theta < 5^{ circ } hfill \ 16.8sin theta - 0.50,5^{ circ } le theta < 10^{ circ } hfill \ 21.9sin theta - 0.96,theta ge 10^{ circ } hfill \ end{gathered} right. hfill \ C = left{ begin{gathered} 1,c = 0 hfill \ 0.6508 - 0.3436lg c,0 < c le 78.3% hfill \ 0,c > 78.3% hfill \ end{gathered} right. hfill \ end{gathered} $$
(6)
where A is the soil erosion modulus. Ri is the rainfall erosion factor, K is the soil erosion factor, L and S are slope the length factor and the slope factor respectively, C is the vegetation coverage and management factor, P is the soil and water conservation factor, l is the slope length value, m is the slope length index, and θ the is slope value.Landscape ecological functionThe landscape ecological function determines the ability of ecological service50,51,52, the ecological service of urban blue–green space depends on the human value orientation48. It includes four categories: supply, support, regulation, and culture. Based on Maslow’s Hierarchy of Needs and Alderfer’s ERG theory, scholars have summarized the three major needs of human beings for urban blue–green space. Namely, securing the living environment to meet the survival needs, improving social relationships to meet the interaction needs, and cultivating cultural cultivation to meet the development needs53. Specifically corresponding to the landscape ecological function of urban blue–green space, supply is not the main function, only plays a subsidiary role, support is the basic guarantee, regulation is the basic need for urban environmental construction, and culture is an important element of high-quality social life. Ecosystem service value (ESV) can realize the measurement of ecological service function by calculating the specific value of life support products and services produced by the ecosystem54,55,56. Considering the human value orientation of the urban blue–green space landscape ecological function, the weights were given by consulting 16 experts, with supply, regulation, support, and culture weights of 0.2, 0.3, 0.3, 0.2, respectively. The formula is as follows:$$ ESV = sumlimits_{k = 1}^{n} {S_{k} times V_{k}^{{}} } $$
(7)
where Sk is the area of landscape type k, Vk is the value coefficient of the ecosystem service function of landscape type k .Landscape ecological sustainabilityWu (2013) proposed a research framework for landscape sustainability based on a summary of related studies, stating that landscape ecological sustainability is the ability to provide ecosystem services in a long-term and stable manner34. The framework emphasized that landscape sustainability should focus on the analysis of ecosystem service trade-offs effect34,57. In the process of dynamic change of urban blue–green space ecosystem, there are complex trade-offs among various ecosystem services. This is important for promoting the optimal overall benefits of various ecosystem services and achieving sustainable development of urban ecology58. In addition, as a special type of human-centered ecosystem developed by humans based on nature, human well-being is also very important for the landscape ecological sustainability of urban blue–green space. For this reason, we introduced ecosystem service trade-offs (EST) and ecological construction input (IEC) as assessment indexes of landscape ecological sustainability.

(1)

Ecosystem service trade-offs (EST)

This study applied the root mean square deviation of ecological services to quantify ecosystem service trade-offs (EST). The index could effectively measure the average difference in standard deviation between individual ecosystem services and the average ecosystem services. It is a simple and effective way to evaluate the trade-offs among ecosystem services. The formula is as follows:$$ EST = sqrt {frac{1}{n – 1}sumnolimits_{i = 1}^{n} {(ES_{std} – overline{ES}_{std} } } )^{2} $$
(8)
where ESstd is the normalized ecosystem services, n is the number of ecosystem services , and (overline{ES}_{std}) is the mean value of normalized ecosystem services.

(2)

Ecological construction input (ECI)

Human well-being is a premise for the landscape ecological sustainability of urban blue–green spaces, it is closely related to government investment in ecological construction planning34. From the perspective of economics, this study assessed the human well-being obtained by urban blue–green space with the ratio of urban ecological construction investment to GDP, that is, the ecological construction input (ECI). The formula is as follows:$$ ECI = EI/G $$
(9)
where EI is the amount of ecological construction investment, and G is the gross regional product.Evaluation methodThe index weight determines its relative importance in the index system, and the selection of the weight calculation method in the decision-making of multi-attribute problems has an important impact on the assessment results21. Traditional weighting methods can be divided into two categories, subjective weighting method and objective weighting method21,38. The subjective weighting method is represented by the analytic hierarchy process (AHP), Delphi method, and so on. It has the advantage of simplicity, but the disadvantage is too subjective and randomness because it was completely dependent on the knowledge and experience of decision makers. The objective weighting method is represented by the entropy weighting method (EWM), principal component analysis, variation coefficient method, and so on. And it has been widely recognized for reflecting the variability of assessment results18, but the values of indexes have significant influence and the calculation results are not stable. Considering the limitations of the single weighting method, the weights of each assessment index in this study were determined by the combination of subjective weight and objective weight. Among them, the subjective weighting selected the AHP, and the objective weighting selected the EWM (Table 2). The formula is as follows:$$ w_{{j}} = alpha w_{{j}}^{{{AHP}}} + (1 – alpha )w_{{j}}^{{{EWM}}} $$
(10)
$$ w_{{j}}^{{{EWM}}} = d_{{j}} /sumlimits_{i = 1}^{m} {d_{{j}} } $$
(11)
$$ d_{{j}} = 1 – e_{{j}} $$
(12)
$$ e_{{j}} = – ksumlimits_{i = 1}^{n} {f_{{{ij}}} ln (f_{{{ij}}} )} ,;k = 1/ln (n) $$
(13)
$$ f_{{{ij}}} = X^{prime}_{{{ij}}} /sumlimits_{i = 1}^{n} {X^{prime}_{{{ij}}} } $$
(14)
where (W_{{j}}^{{}}) is the combined weight. (W_{{j}}^{{_{AHP} }}) is the weight of the j-th index of the AHP, (W_{{j}}^{{{EWM}}}) is the weight of the j-th index of the EWM, dj is the information entropy of the j-th index, ej is the entropy value of the j-th index, (f_{{{ij}}}) is the proportion of the index value of the j-th sample under the i-th indexm, (X^{prime}_{{{ij}}}) is the standardized value of the i-th sample of the j-th index, m is the number of index, n is the number of samples, and (alpha) was taken as 0.5.Table 2 Weight of assessment index.Full size tableSince the dimensions of indexes are different, it is necessary to unify the dimensions of the index to avoid the errors caused by direct calculation to make the evaluation results inaccurate. The range standardization was used to normalize the index data and bound its value in the interval [0, 1], the range standardization can be expressed as follows15,23:$$ {text{Positive indicator}}left( + right):A_{{{ij}}} = (X_{{{ij}}} – X_{{{jmin}}} )/(X_{{{jmax}}} – X_{{{jmin}}} ) $$
(15)
$$ {text{Negative indicator}}left( – right):A_{{{ij}}} = (X_{{{jmax}}} – X_{ij} )/(X_{{{jmax}}} – X_{{{jmin}}} ) $$
(16)
Additionally, we divided the LEH index into five levels from high to low using an equal-interval approach as follows40: [1–0.8) healthy, [0.8–0.6) sub-healthy, [0.6–0.4) moderately healthy, [0.4–0.2) unhealthy, [0.2–0] pathological, corresponding level I–V. And the level transfer of LEH in different periods was divided into three types: improvement type, degradation type, and stabilization type. For example, III-II means that the transfer from level III to level II is the improvement type.Spatial autocorrelation analysisSpatial autocorrelation analysis is one of the basic methods in theoretical geography. It could deeply investigate the spatial correlation characteristics of data, including global spatial autocorrelation and local spatial autocorrelation23. The global spatial autocorrelation uses global Moran’s I to evaluate the degree of their spatial agglomeration or differentiation of an attribute value in the study area. The local spatial autocorrelation is a decomposed form of the global spatial autocorrelation18,21, including four types: HH(High-High), LL(Low-Low), HL(High-Low), LH(Low–High). In this study, spatial autocorrelation analysis was applied to study the spatial correlation characteristics of blue–green space LEH. The calculation formulas are as follows:$$ I = frac{{Nsumlimits_{i} {sumlimits_{v} {W_{iv} (Y_{i} – overline{Y} )(Y_{v} – overline{Y} )} } }}{{(sumlimits_{i} {sumlimits_{v} {W_{iv} } } )sumlimits_{i} {(Y_{i} – overline{Y} )} }} $$
(17)
$$ I_{i} = frac{{Y_{i} – overline{Y} }}{{S_{x}^{2} }}sumlimits_{v} {left[ {W_{iv} (Y_{i} – overline{Y} )} right]} $$
(18)
where N is the number of space units, (W_{iv}) is the spatial weight, (Y_{i} ,Y_{v}) are the variable attribute values of the area (i,v), (overline{Y}) is the variable mean, (S_{x}^{2}) is the variance, (I) is the global Moran’s I index, and (I_{i}) is the local Moran’s I index. More

We selected eleven major wheat production provinces of China for the study area, which comprise the largest winter wheat-sowing fraction: Henan, Shandong, Anhui, Jiangsu, Hebei, Hubei, Shanxi, Shaanxi, Sichuan, Xinjiang, and Gansu (Fig. 1). The wheat planting area is about 22 million ha in these provinces, accounting for more than 93% of the total wheat planting area. The total wheat production in these regions contributes more than 96% of the total wheat production in China, with more than 128 million tons in 201933.We developed a methodological framework for high-resolution AGB mapping. It mainly includes three parts: (1) Data acquisition and processing. (2) The WOFOST model parameterization and calibration. (3) Data assimilation (Fig. 1). Each part is explained in more detail below.Data acquisition and processingMeteorological dataChina Meteorological Forcing Dataset34,35 is used as weather driving data for the WOFOST model. The dataset is based on the internationally existing Princeton reanalysis data, Global Land Data Assimilation System data, Global Energy and Water Cycle Experiment-Surface Radiation Budget radiation data, and Tropical Rainfall Measuring Mission precipitation data. It is made by fusing the conventional meteorological observation data of the China Meteorological Administration. It includes seven elements: near-surface air temperature, air pressure, near-surface total humidity, wind speed, ground downward shortwave radiation, ground downward longwave radiation, and ground precipitation rate. The meteorological drive elements required for WOFOST are daily radiation, minimum temperature, maximum temperature, water vapor pressure, average wind speed, and precipitation. Details of these variables that participated in the WOFOST model can be referred to in Table S1.Soil characteristics measurements and phenology observationsSoil and phenology data were collected at 177 agricultural meteorological stations (AMS) from 2007 to 2015 (Fig. 1). Soil characteristics include soil moisture content at wilting points, field capacity, and saturation. To be consistent with the corresponding units in the crop model, the original data in weight water content was converted into volume water content through the corresponding soil bulk density measurements. Winter wheat phenology observations include the date of emergence (more than 50% of the wheat seedlings in the field show the first green leaves and reached about 2 cm), anthesis (the inner and outer glumes of the middle and upper florets of more than 50% of the wheat ears in the whole field are open, and the anthers loose powder), and maturity (more than 80% of the wheat grains turn yellow, the glumes and stems turn yellow, and only the upper first and second nodes are still slightly green). In most cases, the phenological stage “anthesis” is missing. The anthesis date was calculated by adding seven days to the observed heading date (when more than 50% of the wheat in the whole field exposes the tip of the ear from the sheath of the flag leaf).County-level yield statistics dataThe county-level yield data was collected from city statistical yearbooks of the study area from 2007 to 2015. Since most statistical yearbooks do not directly record per-unit yield data, the county-level yield was obtained by dividing the total yield and planting area. It is worth noting that all yields were calculated in units of metric kilograms per cultivated hectares (kg·ha−1).The winter wheat land cover dataWe used a winter wheat land cover product from a 1 km resolution dataset named ChinaCropArea1km36. This data was derived from GLASS leaf area index products and crop phenology from 2000 to 2015. This dataset is the base map of our data production.The MODIS LAI dataWe used the improved 8-days MODIS LAI products (i.e., 1 km) generated by Yuan et al.32 to assimilate the WOFOST model. The products applied the modified temporal-spatial filter and Savitzky-Golay filter to overcome the spatial-temporal discontinuity and inconsistence of raw MODIS LAI products, which makes them more applicable for the realm of land surface and climate modeling. The products can be accessed via the Land-Atmosphere Interaction Research Group website at Sun Yat-sen University (http://globalchange.bnu.edu.cn/research/lai).The WOFOST model parameterization and calibrationThe WOFOST model introductionThe WOFOST model was initially developed as a crop growth simulation model to evaluate the yield potential of various crops in tropical countries37. In this study, we chose the WOFOST model because the model reaches a trade-off of the complexity of the crop model and is suitable for large-scale simulations3. The WOFOST model is a typical crop growth model that explains crop growth based on underlying processes such as photosynthesis and respiration and their response to changing environmental conditions38. The WOFOST model estimates phenology, LAI, aboveground biomass, and storage organ biomass (i.e., grain yield) at a daily time step39 (Fig. 2).Fig. 2Schematic overview of the major processes implemented in WOFOST. The Astronomical module calculates day length, some variables relating to solar elevation, and the fraction of diffuse radiation.Full size imageZonal parameterizationWe first divided the study area covered by AMS into seamless Thiessen polygon zones. Each Thiessen polygon contains only a single AMS. These zones represent the whole areas where any location is closer to its associated AMS point than any other AMS point. Then, we assigned parameters to the entire zone based on the AMS data, including crop calendar (date of emergence) and soil water retention parameters (soil moisture content at wilting point, field capacity, and saturation). Besides, we also optimized two main crop parameters for controlling phenological development stages, namely TSUM1 (accumulated temperature required from emergence to anthesis) and TSUM2 (accumulated temperature required from anthesis to maturity), by minimizing the cost function of the observational and simulated date corresponding to anthesis and maturity.Parameter calibration within a single zoneWe implemented the calibration of parameters within every single zone, as illustrated in Fig. 3. We calculated the average statistical yield of each county within every single zone from 2007 to 2015, then ranked the counties in descending order and divided them into three groups, namely high, medium, and low-level yield counties, by the 33% quantile and 67% quantile of the average statistical yield. The three counties corresponding to 17% quantile, 50% quantile, and 83% quantile would be used for subsequent calibration and represent the corresponding three yield level groups. We used the statistical yields (converted to dry matter mass based on the standard moisture content of 12.5%) of the three counties for multiple years and a harvest index for each province to convert the county-level yield to AGB for calibration. The harvest index of each province was mainly estimated from the AMS’s dynamic growth records on the biomass composition of the dominant winter wheat varieties of the province and a published literature40. Besides, we collected the maximum LAI observations on all agrometeorological stations in all years in the study area, according to its histogram. We found that the histogram follows a normal distribution with a mean of 6.5 and a standard deviation of 1.5. Finally, we calibrated three sets of parameters corresponding to three yield level groups in each single zone according to the three selected counties.Fig. 3Flow chart of parameter calibration within a single zone.Full size imageWe designed a three-step calibration strategy for a specific yield level group. Firstly, as winter wheat varieties did not change significantly according to information recorded by agrometeorological stations from 2007 to 2015, we assumed the crop parameters of winter wheat remain unchanged every three years to combine three years of observational data to calibrate the parameters of the WOFOST model better. We maximized a log-likelihood function based on the maximum LAI statistics and every three-year county-level yield and AGB data mentioned to optimize selected crop parameters (see Table S2 in the Supplement Materials).The log-likelihood function was constructed as follows:$$log;{{rm{L}}}_{{rm{LAI}}}=-frac{1}{2}left[dlogleft(2pi right)+logleft(left|{Sigma }_{{rm{LAI}}}right|right)+{rm{MD}}{left({{bf{x}}}_{{rm{LAI}}};{mu }_{{rm{LAI}}},{Sigma }_{{rm{LAI}}}right)}^{2}right]$$
(1)
$$log;{{rm{L}}}_{{rm{TWSO}}}=-frac{1}{2}left[dlog(2pi )+logleft(left|{{boldsymbol{Sigma }}}_{{rm{TWSO}}}right|right)+{rm{MD}}{left({{bf{x}}}_{{rm{TWSO}}};{{boldsymbol{mu }}}_{{rm{TWSO}}},{{boldsymbol{Sigma }}}_{{rm{TWSO}}}right)}^{2}right]$$
(2)
$$log;{{rm{L}}}_{{rm{AGB}}}=-frac{1}{2}left[dlog(2pi )+logleft(left|{{boldsymbol{Sigma }}}_{{rm{AGB}}}right|right)+{rm{MD}}{left({{bf{x}}}_{{rm{AGB}}};{{boldsymbol{mu }}}_{{rm{AGB}}},{{boldsymbol{Sigma }}}_{{rm{AGB}}}right)}^{2}right]$$
(3)
$$log;{rm{L}}=log;{L}_{{rm{LAI}}}+log;{L}_{{rm{TWSO}}}+log;{L}_{{rm{AGB}}}$$
(4)
Where log L is the natural logarithm of the likelihood function, d is the dimension, that is, the number of years of joint calibration, which is set to 3 in this study xLAI is the vector composed of the maximum value of the 3-year LAI simulated by WOFOST, μLAI and ΣLAI are the mean vector and error covariance matrix of maximum LAI based on observation statistics. The annual maximum LAI was assumed to be independent, and the mean and standard deviation for each year was set the same as the result of Fig. 3. Similarly, xTWSO and xAGB are the yield vector and AGB vector at maturity of 3 years simulated by WOFOST, and μTWSO, μAGB are their corresponding county-level statistic vector, ΣTWSO and ΣAGB are their corresponding error covariance matrix. In this study, we assumed that the annual yield or AGB was independent, and their corresponding standard deviation was 10% of their statistical value. |Σ| is the determinant of Σ. The expression ({rm{MD}}{({bf{x}};{boldsymbol{mu }},{boldsymbol{Sigma }})}^{2}={({bf{x}}-{boldsymbol{mu }})}^{{rm{T}}}{{boldsymbol{Sigma }}}^{-1}({bf{x}}-{boldsymbol{mu }})), where MD is the Mahalanobis distance between the point x and the mean vector μ.Secondly, we optimized the inter-annual irrigation. We optimized two parameters every year: the critical value of soil moisture (denoted as SMc) and the amount of irrigation (denoted as V). When the soil moisture simulated by WOFOST is lower than SMc, an irrigation event will be triggered, and the irrigation amount is V cm. In this study, we combined three-year data for calibration with six parameters for optimization. The optimization strategy is the same as the previous step by maximizing the log-likelihood function. Finally, we fixed the optimized irrigation parameters and repeated the first step to calibrate the selected crop parameters and obtain the final optimal parameters.Data assimilationConsidering that MODIS LAI is relatively low compared to the actual LAI of winter wheat41, we select a weak-constraint cost function based on the least square of normalized observational and simulated LAI as shown in Eq. (5), which is assimilating the trend information of MODIS LAI into the crop growth model.$$J={sum }_{{rm{t}}=1}^{{rm{n}}}{left(frac{{{rm{LAI}}}_{{rm{MODIS}}}^{{rm{t}}}-{{rm{LAI}}}_{{rm{MODIS}}}^{min}}{{{rm{LAI}}}_{{rm{MODIS}}}^{max}-{{rm{LAI}}}_{{rm{MODIS}}}^{min}}-frac{{{rm{LAI}}}_{{rm{WOFOS}}}^{{rm{t}}}-{{rm{LAI}}}_{{rm{WOFOS}}}^{min}}{{{rm{LAI}}}_{{rm{WOFOS}}}^{max}-{{rm{LAI}}}_{{rm{WOFOS}}}^{min}}right)}^{2}$$
(5)
Where ({{rm{LAI}}}_{{rm{MODIS}}}^{{rm{t}}}) and .. are MODIS LAI and WOFOST simulated LAI of time t. ({{rm{LAI}}}_{{rm{MODIS}}}^{max}) and ({{rm{LAI}}}_{{rm{WOFOS}}}^{max}) are maximum of MODIS LAI and WOFOST simulated LAI. ({{rm{LAI}}}_{{rm{MODIS}}}^{min}) and ({{rm{LAI}}}_{{rm{WOFOS}}}^{min}) are minimum of MODIS LAI and WOFOST simulated LAI. J is the value of the cost function.We reinitialize the day of emergence (IDEM), the life span of leaves growing at 35 °C (SPAN), and thermal time from emergence to anthesis (TSUM1) in the WOFOST model on each 1 km winter wheat pixel according to cost function between WOFOST LAI and MODIS LAI. Besides, we applied the Subplex algorithm from the NLOPT library (https://github.com/stevengj/nlopt) for parameter optimization. More