Philippa Kaur

Sunquist, M. E. & Sunquist, F. C. Family Felidae. In Handbook of the Mammals of the World Vol. 1 (eds Wilson, D. E. & Mittermeier, R. A.) 54–170 (Lynx Editions, 2009).
Google Scholar 
Breitenmoser, U. et al. Action plan for the conservation of the Eurasian Lynx (Lynx lynx) in Europe. Nat. Environ. 112, 1–70 (2000).
Google Scholar 
Linnell, J. D. C., Breitenmoser, U., Breitenmoser-Würsten, C., Odden, J. & von Arx, M. Recovery of Eurasian lynx in Europe: What part has reintroduction played? In Reintroduction of Top-Order Predators (eds Hayward, M. W. & Somers, M. J.) 72–91 (Blackwell Publishing, 2009).Chapter 

Google Scholar 
Schmidt, K., Ratkiewicz, M. & Konopiński, M. K. The importance of genetic variability and population differentiation in the Eurasian lynx Lynx lynx for conservation, in the context of habitat and climate change. Mammal Rev. 41, 112–124 (2011).Article 

Google Scholar 
von Arx, M. et al. Status and conservation of the Eurasian lynx (Lynx lynx) in Europe in 2001. KORA Bericht 19, 1–330 (2004).
Google Scholar 
Kaczensky, P. et al. Status, management and distribution of large carnivores—Bear, lynx, wolf and wolverine in Europe. Part 1 – Europe summaries. Report: 1–72. A Large Carnivore Initiative for Europe Report prepared for the European Commission (2013).Chapron, G. et al. Recovery of large carnivores in Europe’s modern human-dominated landscapes. Science 346, 1517–1519 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Franz, K. W. & Romanowski, J. Revisiting the reintroduced Eurasian lynx population in Kampinos National Park Poland. Eur. Zool. J. 88, 966–979. https://doi.org/10.1080/24750263.2021.1968046 (2021).Article 

Google Scholar 
Bieniek, M., Wolsan, M. & Okarma, H. Historical biogeography of the lynx in Poland. Acta Zool. Cracov. 41, 143–167 (1998).
Google Scholar 
Jędrzejewski, W., Nowak, S., Schmidt, K. & Jędrzejewska, B. Wilk i ryś w Polsce: Wyniki inwentaryzacji w 2001 roku. Kosmos 51, 491–499 (2002).
Google Scholar 
Mysłajek, R., Kwiatkowska, I., Diserens, T., Haidt, A. & Nowak, S. Occurrence of Eurasian lynx in western Poland after two decades of strict protection. CATnews 69, 12–13 (2019).
Google Scholar 
Schmidt, K. Program ochrony rysia Lynx lynx w Polsce – Project. Strategia ochrony Rysia Warunkująca Trwałość Populacji Gatunku w Polsce (Warsaw University of Life Sciences, 2011).
Google Scholar 
Kaczensky, P. et al. Status, management and distribution of large carnivores—Bear, lynx, wolf and wolverine in Europe. Part 2: Country Species Summaries. Report: 1–200. A Large Carnivore Initiative for Europe Report prepared for the European Commission (2013).Breitenmoser, U. et al. Lynx lynx (errata version published in 2017). The IUCN Red List of Threatened Species 2015: e.T12519A121707666. Accessed 30 Oct 2021 (2015).Vandel, J.-M., Stahl, P., Herrenschmidt, V. & Marboutin, E. Reintroduction of the lynx into the Vosges mountain massif: From animal survival and movements to population development. Biol. Conserv. 131, 370–385. https://doi.org/10.1016/j.biocon.2006.02.012 (2006).Article 

Google Scholar 
Zimmermann, F., Breitenmoser-Würsten, C. & Breitenmoser, U. Importance of dispersal for the expansion of a Eurasian lynx Lynx lynx population in a fragmented landscape. Oryx 41, 358–368. https://doi.org/10.1017/s0030605307000712 (2007).Article 

Google Scholar 
Schmidt, K., Kowalczyk, R., Ozolins, J., Mannil, P. & Fickel, J. Genetic structure of the Eurasian lynx population in north-eastern Poland and the Baltic states. Conserv. Genet. 10, 497–501. https://doi.org/10.1007/s10592-008-9795-7 (2009).Article 

Google Scholar 
Ratkiewicz, M. et al. Long-range gene flow and the effects of climatic and ecological factors on genetic structuring in a large, solitary carnivore: The Eurasian Lynx. PLoS ONE 9, e115160. https://doi.org/10.1371/journal.pone.0115160 (2014).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Port, M. et al. Rise and fall of a Eurasian lynx (Lynx lynx) stepping-stone population in central Germany. Mammal Res. https://doi.org/10.1007/s13364-020-00527-6 (2020).Article 

Google Scholar 
Drouilly, M. & O’Riain, J. M. Rewilding the world’s large carnivores without neglecting the human dimension: A response to reintroducing the Eurasian lynx to southern Scotland, England and Wales. Biodivers. Conserv. 30, 917–923. https://doi.org/10.1007/s10531-021-02112-y (2021).Article 

Google Scholar 
Böer, M., Smielowski, J. & Tyrala, P. Reintroduction of the European lynx (Lynx lynx) to the Kampinoski National Park/Poland field experiment with zooborn individuals. Part I: Selection, adaptation and training. Der Zool. Garten 70, 304–312 (1994).
Google Scholar 
Jakimiuk, S. (ed.). Aktywna ochrona populacji nizinnej rysia w Polsce. 1–144 (WWF, Poland, 2015).Huck, M. et al. Habitat suitability, corridors and dispersal barriers for large carnivores in Poland. Acta Theriol. 55, 177–192 (2010).Article 

Google Scholar 
Niedziałkowska, M. et al. Environmental correlates of Eurasian lynx occurrence in Poland: Large scale census and GIS mapping. Biol. Conserv. 133, 63–69. https://doi.org/10.1016/j.biocon.2006.05.022 (2006).Article 

Google Scholar 
Schmidt, K., Kowalczyk, R., Ozolins, J., Männil, P. & Fickel, J. Genetic structure of the Eurasian lynx population in north-eastern Poland and the Baltic states. Conserv. Genet. 10, 497–501. https://doi.org/10.1007/s10592-008-9795-7 (2009).Article 

Google Scholar 
Tracz, M. et al. The return of lynx to northwestern Poland. CATnews 14, 43–44 (2021).
Google Scholar 
The Return of Lynx to north-west Poland. http://www.rysie.org/en/rysie-strona-glowna. Accessed on 31 Oct 2021.IUCN/SSC. Guidelines for Reintroductions and Other Conservation Translocations. Version 1.0. 1–57 (IUCN Species Survival Commission, 2013).Rueda, C., Jiménez, J., Palacios, M. J. & Margalida, A. Exploratory and territorial behavior in a reintroduced population of Iberian lynx. Sci. Rep. 11, 14148. https://doi.org/10.1038/s41598-021-93673-z (2021).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Gusset, M. A framework for evaluating reintroduction success in carnivores: Lessons from African wild dogs. In Reintroduction of Top-Order Predators (eds Hayward, M. W. & Somers, M. J.) 307–320 (Blackwell Publishing, 2009).Chapter 

Google Scholar 
Breitenmoser, U. & Haller, H. Patterns of predation by reintroduced European Lynx in the Swiss Alps. J. Wildl. Manage. 57, 135–144 (1993).Article 

Google Scholar 
Drouilly, M. & O’Riain, M. J. Rewilding the world’s large carnivores without neglecting the human dimension. Biodivers. Conserv. 30, 917–923 (2021).Article 

Google Scholar 
Jędrzejewski, W. et al. Population dynamics (1869–1994), demography, and home ranges of the Lynx in Białowieza Primeval Forest (Poland and Belarus). Ecography 19, 122–138 (1996).Article 

Google Scholar 
Palmero, S. et al. Demography of a Eurasian lynx (Lynx lynx) population within a strictly protected area in Central Europe. Sci. Rep. 11, 19868. https://doi.org/10.1038/s41598-021-99337-2 (2021).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Maran, T., Põdra, M., Põlma, M. & Macdonald, D. The survival of captive-born animals in restoration programmes: Case study of the endangered European mink Mustela lutreola. Biol. Conserv. 142, 1685–1692 (2009).Article 

Google Scholar 
Moehrenschlager, A. & Macdonald, D. W. Movement and survival parameters of translocated and resident swift foxes Vulpes velox. Anim. Conserv. 6, 199–206 (2003).Article 

Google Scholar 
Böer, M., Reklewski, J., Śmiełowski, J. & Tyrała, P. Reintroduction of the European Lynx to the Kampinoski Nationalpark/Poland: A field experiment with zooborn individuals. Part III: Demographic development of the population from December 1993 until January 2000. Der Zool. Garten 70, 304–312 (2000).
Google Scholar 
Jule, K. R., Leaver, L. A. & Lea, E. G. L. The effects of captive experience on reintroduction survival in carnivores: A review and analysis. Biol. Conserv. 141, 355–363 (2008).Article 

Google Scholar 
Hellstedt, P. & Kallio, E. R. Survival and behaviour of captive-born weasels (Mustela nivalis nivalis) released in nature. J. Zool. 266, 37–44 (2005).Article 

Google Scholar 
Devineau, O. et al. Evaluating the Canada lynx reintroduction programme in Colorado: Patterns in mortality. J. Appl. Ecol. 47, 524–531 (2010).Article 

Google Scholar 
Lengger, J., Breitenmoser, U. & Sliwa, A. EAZA breeding programmes as sources for lynx reintroductions. CATnews 14, 76–77 (2021).
Google Scholar 
Reading, P. R. & Clark, T. W. Carnivore introductions: An interdisciplinary Examination. In Carnivore Behavior, Ecology and Evolution (ed. Gittleman, J. L.) 296–336 (Cornell University Press, 1996).
Google Scholar 
McCarthy, M. A., Armstrong, D. P. & Runge, M. C. Adaptive management of reintroduction. In Reintroduction Biology: Integrating Science and Management (eds Ewen, J. G. et al.) 256–289 (Wiley-Blackwell, 2012).Chapter 

Google Scholar 
Bremner-Harrison, S., Prodohl, P. A. & Elwood, R. W. Behavioural trait assessment as a release criterion: Boldness predicts early death in a reintroduction programme of captive-bred swift fox (Vulpes velox). Anim. Conserv. 7, 313–320 (2004).Article 

Google Scholar 
Harrington, L., Põdra, M., Macdonald, D. & Maran, T. Post-release movements of captive-born European mink Mustela lutreola. Endanger. Species Res. 24, 137–148 (2014).Article 

Google Scholar 
Andrén, H. et al. Survival rates and causes of mortality in Eurasian lynx (Lynx lynx) in multi-use landscapes. Biol. Conserv. 131, 23–32 (2006).Article 

Google Scholar 
Heurich, M. et al. Illegal hunting as a major driver of the source-sink dynamics of a reintroduced lynx population in Central Europe. Biol. Conserv. 224, 355–365 (2018).Article 

Google Scholar 
Schmidt-Posthaus, H., Breitenmoser, Ch., Posthaus, H., Bacciarini, L. & Breitenmoser, U. Causes of mortality in reintroduced Eurasian lynx in Switzerland. J. Wildl. Dis. 38, 84–92 (2002).PubMed 
Article 

Google Scholar 
Kołodziej-Sobocińska, M., Zalewski, A. & Kowalczyk, R. Sarcoptic mange vulnerability in carnivores of the Białowieża Primeval Forest, Poland: underlying determinant factors. Ecol. Res. 29, 237–244 (2014).Article 

Google Scholar 
Holt, G. & Berg, C. Sarcoptic mange in red fox and other wild carnivores in Norway. Nor Veterinaertidsskr 102, 427–432 (1990).
Google Scholar 
Mörner, T. Sarcoptic mange in Swedish wildlife. Rev. Sci. Tech. Off. Int. Epiz. 11, 1115–1121 (1992).Article 

Google Scholar 
Ryser-Degiorgis, M. P. et al. Notoedric and sarcoptic mange in free-ranging lynx from Switzerland. J. Wildl. Dis. 38, 228–232 (2002).PubMed 
Article 

Google Scholar 
Soulsbury, C. D. et al. The impact of sarcoptic mange Sarcoptes scabiei on the British fox Vulpes vulpes population. Mam. Rev. 37, 278–296 (2007).
Google Scholar 
Garrote, G., Fernández-López, J., López, G., Ruiz, G. & Simón, M. A. Prediction of Iberian lynx road–mortality in southern Spain: A new approach using the MaxEnt algorithm. Anim. Biodivers. Conserv. 41, 217–225 (2018).Article 

Google Scholar 
Bencin, H., Prange, S., Rose, Ch. & Popescu, V. Roadkill and space use data predict vehicle-strike hotspots and mortality rates in a recovering bobcat (Lynx rufus) population. Sci. Rep. 9, 15391 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Bouyer, Y. et al. Tolerance to anthropogenic disturbance by a large carnivore: The case of Eurasian lynx in south-eastern Norway. Anim. Conserv. https://doi.org/10.1111/acv.12168 (2014).Article 

Google Scholar 
López-Bao, J. V. et al. Eurasian lynx fitness shows little variation across Scandinavian human-dominated landscapes. Sci. Rep. 9, 8903 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Cimatti, M. et al. Large carnivore expansion in Europe is associated with human population density and land cover changes. Divers. Distrib. 27, 602–617 (2021).Article 

Google Scholar 
Wegner, M. (ed.). Statistical Yearbook of Zachodniopomorskie Voivodship. 1–213 (Statistical Office in Szczecin, 2020).Górny, M., Schmidt, K. & Kowalczyk, R. Analiza przydatności środowiska dla reintrodukcji rysia w północno-zachodniej Polsce oraz prognoza i perspektywy funkcjonowania populacji. Expert study under the project POIS.02.04.00–0143/16 “Return of the lynx to northwestern Poland”. 1–25.Woodford, M. H., Keet, D. F. & Bengis, R. G. Post-mortem Procedures for Wildlife Veterinarians and Field Biologists. 1–55 (IUCN Species Survival Commission (SSC) & Veterinary Specialist Group, Care for the Wild International, World Organisation for Animal Health (OIE), 2000).Fain, A. Ѐtude de la variabilitѐ de Sarcoptes scabiei avec une rѐvision des Sarcoptidae. Acta Zool. Pathol. Antverp 47, 1–196 (1968).
Google Scholar 
Kaplan, E. L. & Meier, P. Nonparametric estimation from incomplete observations. J. Am. Stat. Assoc. 53, 457–481 (2012).MathSciNet 
MATH 
Article 

Google Scholar 
Therneau, M., Lumley, T., Atkinson, E. & Crowson, C. Survival Analysis. R Package Version 3.2-13. http://CRAN.R-project.org/package=survival (2021).Kassambara, A., Kosinski, M., Biecek, P. & Scheipl, F. survminer. Drawing Survival Curves using ‘ggplot2’. R package version 0.4.9. http://CRAN.R-project.org/package=survminer (2021).Dardis, C. survMisc. Miscellaneous Functions for Survival Data. R package version 0.5.5. http://CRAN.R-project.org/package=survMisc (2018).R Core Team. R: A language and environment for statistical computing (R Foundation for Statistical Computing). https://www.R-project.org (2021).Snedecor, G. W. & Cochran, W. G. Statistical Methods 7th edn. (Iowa State University Press, 1980).MATH 

Google Scholar 
Cox, D. R. Regression models and life tables (with discussion). J. R. Stat. Soc. B. 34, 187–220 (1972).MATH 

Google Scholar 
Bradburn, M. J., Clark, T. G., Love, S. B. & Altman, D. G. Survival Analysis Part II: Multivariate data analysis: An introduction to concepts and methods. Br. J. Cancer. 89, 431–436 (2003).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Wald, A. Tests of statistical hypothesis concerning several parameters when the number of observations is large. Trans. Am. Math. Soc. 54, 426–482 (1943).MATH 
Article 

Google Scholar 
Aitchison, J. & Silvey, S. D. Maximum likelihood estimation of parameters subject to restraints. Ann. Math. Stat. 29, 813–828 (1958).MathSciNet 
MATH 
Article 

Google Scholar 
Mantel, N. Evaluation of survival data and two new rank order statistics arising in its consideration. Cancer Chemother. Rep. 50, 163–170 (1966).CAS 
PubMed 

Google Scholar  More

Díaz, S. et al. The global spectrum of plant form and function. Nature 529, 167–171 (2016).ADS 
Article 

Google Scholar 
Céréghino, R. et al. Constraints on the functional trait space of aquatic invertebrates in bromeliads. Funct. Ecol. 32, 2435–2447 (2018).Article 

Google Scholar 
Winemiller, K. O., Fitzgerald, D. B., Bowler, L. & Pianka, E. R. Functional traits, convergent evolution, and periodic tables of niches. Ecol. Lett. 18, 737–751 (2015).Article 

Google Scholar 
Díaz, S. et al. Functional traits, the phylogeny of function, and ecosystem service vulnerability. Ecol. Evol. 3, 2958–2975 (2013).Article 

Google Scholar 
Leimar, O. Evolutionary change and Darwinian demons. Selection 2, 65–72 (2001).Article 

Google Scholar 
Cummins, K. W. & Klug, M. J. Feeding ecology of stream invertebrates. Annu. Rev. Ecol. Syst. 10, 147–172 (1979).Article 

Google Scholar 
Pianka, E. R., Vitt, L. J., Pelegrin, N., Fitzgerald, D. B. & Winemiller, K. O. Toward a periodic table of niches, or exploring the lizard niche hypervolume. Am. Nat. 190, 601–616 (2017).Article 

Google Scholar 
Rosenberg, D. M. & Resh, V. H. Freshwater Biomonitoring and Benthic Macroinvertebrates (Chapman and Hall, 1993).
Google Scholar 
Allan, J. D. & Castillo, M. M. Stream Ecology. Structure and Function of Running Waters 2nd edn. (Springer, 2007).
Google Scholar 
Wallace, J. B. & Webster, J. R. The role of macroinvertebrates in stream ecosystem function. Annu. Rev. Entomol. 41, 115–139 (1996).CAS 
Article 

Google Scholar 
Southwood, T. R. E. Habitat, the templet for ecological strategies?. J. Anim. Ecol. 46, 336–365 (1977).Article 

Google Scholar 
Townsend, C. R. & Hildrew, A. G. Species traits in relation to a habitat templet for river systems. Freshw. Biol. 31, 265–275 (1994).Article 

Google Scholar 
Violle, C. et al. Let the concept of trait be functional!. Oikos 116, 882–892 (2007).Article 

Google Scholar 
Schmera, D., Heino, J., Podani, J., Erős, T. & Dolédec, S. Functional diversity: A review of methodology and current knowledge in freshwater macroinvertebrate research. Hydrobiologia 787, 27–44 (2017).Article 

Google Scholar 
Diehl, S. Fish predation and benthic community structure: The role of omnivory and habitat complexity. Ecology 73, 1646–1661 (1992).Article 

Google Scholar 
Brucet, S. et al. Effects of temperature, salinity and fish in structuring the macroinvertebrate community in shallow lakes: Implications for the effects of climate change. PLoS ONE 7, e30877 (2012).ADS 
CAS 
Article 

Google Scholar 
Usseglio-Polatera, P., Bournaud, M., Richoux, P. & Tachet, H. Biological and ecological traits of benthic freshwater macroinvertebrates: Relationships and definition of groups with similar traits. Freshw. Biol. 43, 175–205 (2000).Article 

Google Scholar 
Poff, N. L. et al. Functional trait nichees of North American lotic insects: Trait-based ecological applications in light of phylogenetic relationships. J. North Am. Soc. 25, 730–755 (2006).Article 

Google Scholar 
Bonada, N., Dolédec, S. & Statzner, B. Taxonomic and biological trait differences of stream macroinvertebrate communities between Mediterranean and temperate regions: Implications for future climatic scenarios. Glob. Change Biol. 13, 1658–1671 (2007).ADS 
Article 

Google Scholar 
Stazner, B., Bonada, N. & Dolédec, S. Biological attributes discriminating invasive from native European stream macroinvertebrates. Biol. Invasions 10, 517–530 (2008).Article 

Google Scholar 
Schmidt-Kloiber, A. & Hering, D. www.freshwaterecology.info—An online tool that unifies, standardises and codifies more than 20,000 European freshwater organisms and their ecological preferences. Ecol. Indic. 53, 271–282 (2015).Article 

Google Scholar 
Verberk, W. C. E. P., Siepel, H. & Esselink, H. Life-history strategies in freshwater macroinvertebrates. Freshw. Biol. 53, 1722–1738 (2008).Article 

Google Scholar 
Dolédec, S., Statzner, B. & Frainay, V. Accurate description of functional community structure: Identifying stream invertebrates to species-level?. Bull. North Am. Benthol. Soc. 15, 154–155 (1998).
Google Scholar 
Podani, J., Kalapos, T., Barta, B. & Schmera, D. Principal component analysis of incomplete data—A simple solution to an old problem. Ecol. Inform. 61, 101235 (2021).Article 

Google Scholar 
Podani, J. Introduction to the Exploration of Multivariate Biological Data (Backhuys Publishers, 2000).MATH 

Google Scholar 
Tachet, H., Richoux, P., Bournaud, M. & Usseglio-Polatera, P. Invertébrés d’eau douce—systématique, biologie, écologie 600 (CNRS Editions, 2010).
Google Scholar 
Chevenet, F., Dolédec, S. & Chessel, D. A fuzzy coding approach for the analysis of long-term ecological data. Freshw. Biol. 31, 295–309 (1994).Article 

Google Scholar 
Schmidt-Kloiber A. & Hering, D. www.freshwaterecology.info—The Taxa and Autecology Database for Freshwater Organisms, Version 7.0. (Accessed on 12.09.2019) (2019).Schmera, D., Podani, J., Heino, J., Erős, T. & Poff, N. L. A proposed unified terminology of species traits in stream ecology. Freshw. Sci. 34, 823–830 (2015).Article 

Google Scholar 
Schmera, D., Podani, J., Erős, T. & Heino, J. Combining taxon-by-trait and taxon-by-site matrices for analysing trait patterns of macroinvertebrate communities: A rejoinder to Monaghan & Soares (2014). Freshw. Biol. 59, 1551–1557 (2014).Article 

Google Scholar 
Bonada, N. et al. Do Mediterranean genera not included in Tachet et al. 2002 have Mediterranean characteristics?. Limnetica 30, 129–142 (2011).Article 

Google Scholar 
de Jong, Y. et al. Fauna Europaea—All European animal species on the web. Biodivers. Data J. 2, e4034 (2014).Article 

Google Scholar 
de Bello, F., Botta-Dukát, Z., Leps, J. & Fibich, P. Towards a more balanced combination of multiple traits when computing functional differences between species. Methods Ecol. Evol. 12, 443–448 (2021).Article 

Google Scholar 
Legendre, P. & Legendre, L. Numerical Ecology 2 English. (Elsevier, 1998).MATH 

Google Scholar 
Cornwell, W. K., Schwilk, D. W. & Ackerly, D. A trait-based test for habitat filtering: Convex hull volume. Ecology 87, 1465–1471 (2006).Article 

Google Scholar 
Anderson, M. J. A new method for non-parametric multivariate analysis of variance. Austral. Ecol. 26, 32–46 (2001).
Google Scholar 
Anderson, M. J. Distance-based tests for homogeneity of multivariate dispersions. Biometrics 62, 245–253 (2006).MathSciNet 
Article 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2020).
Google Scholar 
Kindt, R. & Coe, R. Tree Diversity Analysis. A Manual and Software for Common Statistical Methods for Ecological and Biodiversity Studies (World Agroforestry Centre (ICRAF), 2005).
Google Scholar 
Habel, K., Grasman, R., Gramacy, R. G., Mozharovskyi, P. & Sterratt, D. C. Geometry: Mesh Generation and Surface Tessellation. R package version 0.4.5. https://CRAN.R-project.org/package=geometry (2019).Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Minchin, P. R., O’Hara, R. G., Simpson, G. L., Solymos, P., Stevens, M. H. H., Szoecs, E. & Wagner, H. vegan: Community Ecology Package. R package version 2.5-6. https://CRAN.R-project.org/package=vegan (2019). More

After the successful protection of Himalayan areas on the border of China and Nepal, we propose that the two nations should create the world’s highest international peace park by combining the Qomolangma and Sagarmatha national parks. This would align with United Nations Sustainable Development Goal 17, to achieve sustainable development through international cooperation (see go.nature.com/3ixmini).
Competing Interests
The authors declare no competing interests. More

Fisher, M. C. et al. Emerging fungal threats to animal, plant and ecosystem health. Nature 484, 186–194 (2012).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Fisher, M. C., Gow, N. A. & Gurr, S. J. Tackling emerging fungal threats to animal health, food security and ecosystem resilience. Philos. Trans. R. Soc. B 371, 20160332. https://doi.org/10.1098/rstb.2016.0332 (2016).Article 

Google Scholar 
Lips, K. R. Overview of chytrid emergence and impacts on amphibians. Philos. Trans. R. Soc. B 371, 20150465. https://doi.org/10.1098/rstb.2015.0465 (2016).Article 

Google Scholar 
Lips, K. R., Diffendorfer, J., Mendelson, J. R. III. & Sears, M. W. Riding the wave: Reconciling the roles of disease and climate change in amphibian declines. PLoS Biol. 6, e72 (2008).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Caruso, N. M. & Lips, K. R. Truly enigmatic declines in terrestrial salamander populations in great smoky mountains national park. Divers. Distrib. 19, 38–48 (2013).Article 

Google Scholar 
Martel, A. et al. Recent introduction of a chytrid fungus endangers western palearctic salamanders. Science 346, 630–631 (2014).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
van der Spitzen Sluijs, A. et al. Rapid enigmatic decline drives the fire salamander (Salamandra salamandra) to the edge of extinction in the Netherlands. Amphib.-Reptil. 34, 233–239 (2013).Article 

Google Scholar 
Blehert, D. S. et al. Bat white-nose syndrome: An emerging fungal pathogen?. Science 323, 227–227 (2009).CAS 
PubMed 
Article 

Google Scholar 
Thogmartin, W. E., King, R. A., McKann, P. C., Szymanski, J. A. & Pruitt, L. Population-level impact of white-nose syndrome on the endangered Indiana bat. J. Mammal. 93, 1086–1098 (2012).Article 

Google Scholar 
Fisher, M. C., Garner, T. W. & Walker, S. F. Global emergence of Batrachochytrium dendrobatidis and amphibian chytridiomycosis in space, time, and host. Annu. Rev. Microbiol. 63, 291–310 (2009).CAS 
PubMed 
Article 

Google Scholar 
Martel, A. et al. Batrachochytrium salamandrivorans sp. Nov. causes lethal chytridiomycosis in amphibians. Proc. Natl. Acad. Sci. 110, 15325–15329 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Allender, M. C., Raudabaugh, D. B., Gleason, F. H. & Miller, A. N. The natural history, ecology, and epidemiology of Ophidiomyces ophiodiicola and its potential impact on free-ranging snake populations. Fungal Ecol. 17, 187–196. https://doi.org/10.1016/j.funeco.2015.05.003 (2015).Article 

Google Scholar 
Grioni, A. et al. Detection of Ophidiomyces ophidiicola in a wild Burmese python (Python bivittatus) in Hong Kong SAR, China. J. Herpetol. Med. Surg. 31, 283–291 (2021).Article 

Google Scholar 
Allender, M. C. et al. Chrysosporium sp. infection in eastern massasauga rattlesnakes. Emerg. Infect. Dis. 17, 2383–2384. https://doi.org/10.1136/vr.b4816 (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
Franklinos, L. H. V. et al. Emerging fungal pathogen Ophidiomyces ophiodiicola in wild European snakes. Sci. Rep. 7, 1–7. https://doi.org/10.1038/s41598-017-03352-1 (2017).CAS 
Article 

Google Scholar 
Lorch, J. M. et al. Experimental infection of snakes with Ophidiomyces ophiodiicola causes pathological changes that typify snake fungal disease. mBio 6, 1–9. https://doi.org/10.1128/mBio.01534-15 (2015).CAS 
Article 

Google Scholar 
Clark, R. W., Marchand, M. N., Clifford, B. J., Stechert, R. & Stephens, S. Decline of an isolated timber rattlesnake (Crotalus horridus) population: Interactions between climate change, disease, and loss of genetic diversity. Biol. Cons. 144, 886–891. https://doi.org/10.1016/j.biocon.2010.12.001 (2011).Article 

Google Scholar 
Chandler, H. C. et al. Ophidiomycosis prevalence in Georgia’s eastern indigo snake (Drymarchon couperi) populations. PLoS ONE 14, e0218351 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Guthrie, A. L., Knowles, S., Ballmann, A. E. & Lorch, J. M. Detection of snake fungal disease due to Ophidiomyces ophiodiicola in Virginia, USA. J. Wildl. Dis. 52, 143–149. https://doi.org/10.7589/2015-01-007 (2016).CAS 
Article 
PubMed 

Google Scholar 
Last, L. A., Fenton, H., Gonyor-McGuire, J., Moore, M. & Yabsley, M. J. Snake fungal disease caused by Ophidiomyces ophiodiicola in a free-ranging mud snake (Farancia abacura). J. Vet. Diagn. Invest. 28, 709–713. https://doi.org/10.1177/1040638716663250 (2016).CAS 
Article 
PubMed 

Google Scholar 
Lorch, J. M. et al. Snake fungal disease: An emerging threat to wild snakes. Philos. Trans. R. Soc. B 371, 20150457. https://doi.org/10.1098/rstb.2015.0457 (2016).Article 

Google Scholar 
Haynes, E. et al. First report of ophidiomycosis in a free-ranging California Kingsnake (Lampropeltis californiae) in California, USA. J. Wildl. Dis. 57, 246–249 (2021).CAS 
PubMed 
Article 

Google Scholar 
Burbrink, F. T., Lorch, J. M. & Lips, K. R. Host susceptibility to snake fungal disease is highly dispersed across phylogenetic and functional trait space. Sci. Adv. 3, 1–10. https://doi.org/10.1126/sciadv.1701387 (2017).Article 

Google Scholar 
Dixon, J. R. Amphibians and Reptiles of Texas: With Keys, Taxonomic synopses, Bibliography, and Distribution Maps 3rd edn. (Texas A&M University Press, 2000).
Google Scholar 
McKeown, S. A Field Guide to Reptiles and Amphibians in the Hawaiian Islands (Diamond Head Publishing, 1996).
Google Scholar 
Powell, R., Conant, R. & Collins, J. T. Peterson Field Guide to Reptiles and Amphibians of Eastern and Central NORTH AMERICA (Houghton Mifflin Harcourt, 2016).
Google Scholar 
Stebbins, R. C. & McGinnis, S. M. Peterson Field Guide to Western Reptiles and Amphibians (Houghton Mifflin Harcourt, 2018).
Google Scholar 
Texas Administrative Code. State‐listed threatened species in Texas. 31 TAC §65.175. (2020).Dixon, J. R., Werler, J. E. & Forstner, M. R. J. Texas Snakes: A Field Guide Revised. (University of Texas Press, 2020).Book 

Google Scholar 
Rodriguez, D., Forstner, M. R. J., McBride, D. L., Densmore, L. D. III. & Dixon, J. R. Low genetic diversity and evidence of population structure among subspecies of Nerodia harteri, a threatened water snake endemic to Texas. Conserv. Genet. 13, 977–986 (2012).Article 

Google Scholar 
Scott, N. J., Maxwell, T. C., Thornton, O. W., Fitzgerald, L. A. & Flury, J. W. Distribution, habitat, and future of Harter’s water snake, Nerodia harteri Texas. J. Herpetol. 23, 373–389 (1989).Article 

Google Scholar 
Whiting, M. J., Dixon, J. R. & Greene, B. D. Spatial ecology of the Concho water snake (Nerodia harteri paucimaculata) in a large lake system. J. Herpetol. 31, 327–335 (1997).Article 

Google Scholar 
McBride, D. L. Distribution and status of the Brazos water snake (Nerodia harteri harteri) Master of Science thesis, Tarleton State University (2009).United States Office of the Federal Register. Endangered and threatened wildlife and plants; determination of Nerodia harteri paucimaculata (Concho water snake) to be a threatened species Final rule. Fed. Regist. 51, 31412–31422 (1986).
Google Scholar 
United States Office of the Federal Register. Endangered and threatened wildlife and plants; removal of the Concho water snake from the federallist of endangered and threatened wildlife and removal of designated critical habitat. Fed. Reg. 76, 66779–66804 (2011).
Google Scholar 
United States Office of the Federal Register. Endangered and threatened wildlife and plants; findings on petitions and initiation of status review. Fed. Reg. 50, 29238–29239 (1985).
Google Scholar 
United States Office of the Federal Register. Endangered and threatened wildlife and plants; animal candidate review for listing as endangered or threatened species. Fed. Reg. 59, 58982–59028 (1994).
Google Scholar 
Gibbons, J. W. & Dorcas, M. E. North American Watersnakes: A Natural History (University of Oklahoma Press, 2004).
Google Scholar 
Werler, J. E. & Dixon, J. R. Texas Snakes: Identification, Distribution, and Natural History (University of Texas Press, 2000).
Google Scholar 
Lind, C. M., McCoy, C. M. & Farrell, T. M. Tracking outcomes of snake fungal disease in free-ranging pigmy rattlesnakes (Sistrurus miliarius). J. Wildl. Dis. 54, 352–356. https://doi.org/10.7589/2017-05-109 (2018).Article 
PubMed 

Google Scholar 
McBride, M. P. et al. Ophidiomyces ophiodiicola dermatitis in eight free-ranging timber rattlesnakes (Crotalus horridus) from Massachusetts. J. Zoo Wildl. Med. 46, 86–94. https://doi.org/10.1638/2012-0248R2.1 (2015).Article 
PubMed 

Google Scholar 
McCoy, C. M., Lind, C. M. & Farrell, T. M. Environmental and physiological correlates of the severity of clinical signs of snake fungal disease in a population of pigmy rattlesnakes Sistrurus miliarius. Conserv. Physiol. 5, cow077. https://doi.org/10.1093/conphys/cow077 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Haynes, E. et al. Ophidiomycosis surveillance of snakes in Georgia, USA reveals new host species and taxonomic associations with disease. Sci. Rep. 10, 1–15 (2020).Article 
CAS 

Google Scholar 
Stengle, A. G. et al. Evidence of vertical transmission of the snake fungal pathogen Ophidiomyces ophiodiicola. J. Wildl. Dis. 55, 961–964 (2019).PubMed 
Article 

Google Scholar 
Britton, M., Allender, M. C., Hsiao, S.-H. & Baker, S. J. Postnatal mortality in neonate rattlesnakes associated with Ophidiomyces ophiodiicola. J. Zoo Wildl. Med. 50, 672–677 (2019).PubMed 
Article 

Google Scholar 
Allender, M. C., Hileman, E., Moore, J. & Tetzlaff, S. Detection of Ophidiomyces, the caustive agent of snake fungal disease, in the eastern massasauga (Sistrurus catenatus) in Michigan, USA, 2014. J. Wildl. Dis. 52, 694–698. https://doi.org/10.7589/2015-12-333 (2016).CAS 
Article 
PubMed 

Google Scholar 
Hileman, E. T. et al. Estimation of Ophidiomyces prevalence to evaluate snake fungal disease risk. J. Wildl. Manag. 82, 173–181. https://doi.org/10.1002/jwmg.21345 (2018).Article 

Google Scholar 
McKenzie, J. M. et al. Field diagnostics and seasonality of Ophidiomyces ophiodiicola in wild snake populations. EcoHealth 16, 141–150 (2019).PubMed 
Article 

Google Scholar 
Snyder, S. D., Sutton, W. B. & Walker, D. M. Prevalence of Ophidiomyces ophiodiicola, the causative agent of Snake Fungal Disease, in the Interior Plateau Ecoregion of Tennessee, USA. J. Wildl. Dis. 56, 907–911 (2020).PubMed 
Article 

Google Scholar 
Tetzlaff, S. J. et al. Snake fungal disease affects behavior of free-ranging massasauga rattlesnakes (Sistrurus catenatus). Herpetol. Conserv. Biol. 12, 624–634 (2017).
Google Scholar 
Aldridge, R. D., Flanagan, W. P. & Swarthout, J. T. Reproductive biology of the water snake Nerodia rhombifer from Veracruz, Mexico, with comparisons of tropical and temperate snakes. Herpetologica 51, 182–192 (1995).
Google Scholar 
Greene, B. D., Dixon, J. R., Whiting, M. J. & Mueller, J. M. Reproductive ecology of the Concho water snake Nerodia harteri paucimaculata. Copeia 1999, 701–709 (1999).Article 

Google Scholar 
Kofron, C. P. Reproduction of aquatic snakes in south-central Louisiana. Herpetologica 35, 44–50 (1979).
Google Scholar 
Green, B. D. Life History and Ecology of the Concho Water Snake, Nerodia harteri paucimaculata. Dissertation (Texas A&M University, 1993).
Google Scholar 
McKenzie, C. M. et al. Ophidiomycosis in red cornsnakes (Pantherophis guttatus): potential roles of brumation and temperature on pathogenesis and transmission. Vet. Pathol. 57, 825–837 (2020).CAS 
PubMed 
Article 

Google Scholar 
Gregoire, D. R. Nerodia rhombifer (Hallowell, 1852): U.S. geological survey, nonindigenous aquatic species database, Gainesville, FL, Retrieved from 27 Oct 2009 https://nas.er.usgs.gov/queries/FactSheet.aspx?SpeciesID=2577.Janecka, M. J., Janecka, J. E., Haines, A. M., Michaels, A. & Criscione, C. D. Post-delisting genetic monitoring reveals population subdivision along river and reservoir localities of the endemic Concho water snake (Nerodia harteri paucimaculata). Conserv. Genet. 22, 1005–1021 (2021).CAS 
Article 

Google Scholar 
Madsen, T., Stille, B. & Shine, R. Inbreeding depression in an isolated population of adders Vipera berus. Biol. Cons. 75, 113–118 (1996).Article 

Google Scholar 
Carter, J. et al. Variation in pathogenicity associated with the genetic diversity of Fusarium graminearum. Eur. J. Plant Pathol. 18, 573–583 (2002).Article 

Google Scholar 
Charlesworth, D. & Willis, J. H. The genetics of inbreeding depression. Nat. Rev. Genet. 10, 783 (2009).CAS 
PubMed 
Article 

Google Scholar 
Keller, L. F. & Waller, D. M. Inbreeding effects in wild populations. Trends Ecol. Evol. 17, 230–241 (2002).Article 

Google Scholar 
Nieminen, M., Singer, M. C., Fortelius, W., Schöps, K. & Hanski, I. Experimental confirmation that inbreeding depression increases extinction risk in butterfly populations. Am. Nat. 157, 237–244 (2001).CAS 
PubMed 
Article 

Google Scholar 
Roelke, M. E., Martenson, J. S. & O’Brien, S. J. The consequences of demographic reduction and genetic depletion in the endangered Florida panther. Curr. Biol. 3, 340–350. https://doi.org/10.1016/0960-9822(93)90197-v (1993).CAS 
Article 
PubMed 

Google Scholar 
United States Office of the Federal Register. Endangered and threatened wildlife and plants: findings on petitions involving the Yacare Caiman and Harter’s water snake. Fed. Reg. 49, 21089–21090 (1984).
Google Scholar 
NatureServe. NatureServe Explorer: An Online Encyclopedia of Life [web application]. Version 7.0. NatureServe, Arlington, Virginia., http://www.natureserve.org/explorer (2020).Hammerson, G. A. Nerodia harteri (Trapido, 1941). The IUCN red list of threatened species 2007. https://doi.org/10.2305/IUCN.UK.2007.RLTS.T62238A12583490.en (2007).Allender, M. C. et al. Hematology in an eastern massasauga (Sistrurus catenatus) population and the emergence of Ophidiomyces in Illinois, USA. J. Wildl. Dis. 52, 258–269 (2016).CAS 
PubMed 
Article 

Google Scholar 
Becker, C. G., Rodriguez, D., Lambertini, C., Toledo, L. F. & Haddad, C. F. Historical dynamics of Batrachochytrium dendrobatidis in Amazonia. Ecography 39, 954–960 (2016).Article 

Google Scholar 
Rodriguez, D., Becker, C. G., Pupin, N. C., Haddad, C. F. B. & Zamudio, K. R. Long-term endemism of two highly divergent lineages of the amphibian-killing fungus in the Atlantic forest of Brazil. Mol. Ecol. 23, 774–787. https://doi.org/10.1111/mec.12615 (2014).CAS 
Article 
PubMed 

Google Scholar 
Fitch, H. S. Collecting and Life-History Techniques. In Snakes: Ecology and Evolutionary Biology (eds Seigel, Richard A. et al.) 143–164 (Macmillan, 1987).
Google Scholar 
Winne, C. T., Willson, J. D., Andrews, K. M. & Reed, R. N. Efficacy of marking snakes with disposable medical cautery units. Herpetol. Rev. 37, 52–54 (2006).
Google Scholar 
Greene, B. D., Dixon, J. R., Mueller, J. M., Whiting, M. J. & Thornton, O. W. Jr. Feeding ecology of the Concho water snake, Nerodia harteri paucimaculata. J. Herpetol. 28, 165–172 (1994).Article 

Google Scholar 
Lacki, M. J., Hummer, J. W. & Fitzgerald, J. L. Population patterns of copperbelly water snakes (Nerodia erythrogaster neglecta) in a riparian corridor impacted by mining and reclamation. Am. Midl. Nat. 153, 357–369 (2005).Article 

Google Scholar 
Hyatt, A. D. et al. Diagnostic assays and sampling protocols for the detection of Batrachochytrium dendrobatidis. Dis. Aquat. Org. 73, 175–192 (2007).CAS 
Article 

Google Scholar 
Allender, M. C., Bunick, D., Dzhaman, E., Burrus, L. & Maddox, C. Development and use of a real-time polymerase chain reaction assay for the detection of Ophidiomyces ophiodiicola in snakes. J. Vet. Diagn. Invest. 27, 217–220. https://doi.org/10.1177/1040638715573983 (2015).CAS 
Article 
PubMed 

Google Scholar 
Bohuski, E., Lorch, J. M., Griffin, K. M. & Blehert, D. S. TaqMan real-time polymerase chain reaction for detection of Ophidiomyces ophiodiicola, the fungus associated with snake fungal disease. BMC Vet. Res. 11, 1–10. https://doi.org/10.1186/s12917-015-0407-8 (2015).CAS 
Article 

Google Scholar 
Longo, A. V. et al. ITS1 copy number varies among Batrachochytrium dendrobatidis strains: Implications for qPCR estimates of infection intensity from field-collected amphibian skin swabs. PLoS ONE 8, e59499 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ohkura, M. et al. Genome sequence of Ophidiomyces ophiodiicola, an emerging fungal pathogen of snakes. Genome Announc. 5, 1–2 (2017).Article 

Google Scholar 
Falk, B. G., Snow, R. W. & Reed, R. N. A validation of 11 body-condition indices in a giant snake species that exhibits positive allometry. PLoS ONE 12, e0180791 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Garrow, J. S. & Webster, J. Quetelet’s index (W/H2) as a measure of fatness. Int. J. Obes. 9, 147–153 (1985).CAS 
PubMed 

Google Scholar 
Dorai-Raj, S. binom: Binomial confidence intervals for several parameterizations. https://CRAN.R-project.org/package=binom (2014).Thiele, C. & Hirschfeld, G. cutpointr: Improved estimation and validation of optimal cutpoints in R. J. Stat. Softw. 98, 1–27 (2021).Article 

Google Scholar 
Diggle, P. J. Estimating prevalence using an imperfect test. Epidemiol. Res. Int. 2011, 1–5 (2011).Article 

Google Scholar 
Bender, R. & Lange, S. Adjusting for multiple testing: When and how?. J. Clin. Epidemiol. 54, 343–349 (2001).CAS 
PubMed 
Article 

Google Scholar 
Brooks, M. et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 9, 378–400 (2017).Article 

Google Scholar 
Barton, K. MuMIn: Multi-model inference. R package version 1.43.6 (2019).Lenth, R., Singmann, H., Love, J., Buerkner, P. & Herve, M. emmeans: Estimated marginal means, aka least-squares means, R package version 1.4.8. https://CRAN.R-project.org/package=emmeans (2020). More

LED trapsWe used a commercially available portable light trap (Eco-chu trap, Konan Shisetsu Kanri, Okinawa, Japan) to modify the light source. A prototype trap equipped with 12 UV-LED bulbs was developed to catch the green chafer Anomala albopilosa (Hope)35, but it was not sufficiently attractive to stink bugs. Light sources with different numbers of LEDs, from 12 to 84, were used. Either or both bullet-type UV-LED bulbs (NS395L-ERLO; 395 nm, 20 mA, Nitride Semiconductors, Tokushima, Japan) and green LED bulbs (NEPG510S; 525 nm, 20 mA, Nichia, Tokushima, Japan) were used. LEDs were arranged vertically on a stainless-steel cylinder (4.8 cm in diameter, 20 cm in height). Light sources with 12 LEDs were arranged in six rows around the circumference. Each row was arranged as two LEDs at 7.8 cm intervals. Adjacent LEDs were arranged in a left-handed spiral (depression angle of 53°, at approximately 2.5 cm intervals). Light sources with 21 LEDs were arranged in eight rows around the circumference. Each row was arranged as two or three LEDs at 7.2 cm intervals. Adjacent LEDs were arranged in a left-handed spiral (depression angle of 63°, at approximately 2.0 cm intervals). Light sources with 42 LEDs were arranged in eight rows around the circumference. Each row was arranged as five or six LEDs at 3.6 cm intervals. Adjacent LEDs were arranged in a left-handed spiral (elevation angle 63°, at approximately 2.0 cm intervals). Light sources with 84 LEDs were arranged in eight rows around the circumference. Each row was arranged as 10 or 11 LEDs at 1.8 cm intervals. Adjacent LEDs were arranged in a left-handed spiral (elevation angle 63°, at approximately 2.0 cm interval). When both UV and green LEDs were used, both LEDs were arranged alternately in a row (Fig. 4). The cylinder with the LEDs was covered with a transparent acrylic cylinder (9.8 cm in diameter, 20 cm in height).Figure 4Photograph of combined UV and green LED trap used in the experiments.Full size imageThe light source was mounted on a funnel (31 cm in diameter, 24 cm in height), and the lower part of the light source was approximately 100 cm above the ground. A cylindrical chamber (23 cm in diameter, 20 cm in height) was placed under the funnel so that insects that were attracted to the light fell into the funnel and were trapped. The legs of the trap were anchored to the ground using steel stakes. A dimethyl-dichloro-vinyl-phosphate (DDVP) plate containing 10.7 g dichlorvos (Bapona, Earth Chemical, Tokyo, Japan) was placed inside the chamber to kill the insects. The lights were turned on at 18:00 and turned off at 6:00 the next day. The power for the lights was supplied by rechargeable car batteries (N-40B19R/SB; DC 12 V, 28 Ah, Panasonic, Osaka, Japan) or domestic electricity power supplies (AC100V).Emission spectra of combined UV and green lightThe spectral intensity of combined UV and green light was measured using a high-speed spectrometer (HSU-100S, Asahi Spectra, Tokyo, Japan) in a dark room. An attached sensor fiber was placed 50 cm in front of the light source. The measurement was performed five times, the light source was rotated for each measurement to minimize the angle effect, and the average was used as a representative value. The UV- and green-LED emission spectra showed single peaks at wavelengths of 400 and 526 nm, respectively (Fig. 5). Calculated light intensities of UV (350–450 nm) and green (451–600 nm) regions were 2.12 × 1017 and 2.03 × 1017 photons m−2 s−1, respectively; that is, the light intensities of UV- and green-LEDs were almost equal.Figure 5Emission spectra of light source with UV- and green-LEDs. The light source was composed of alternating 42 UV-LEDs and 42 green-LEDs. The intensity of light was measured using a high-speed spectrometer (HSU-100S). An attached sensor fiber was placed 50 cm in front of the light source.Full size imageField evaluation of attractiveness to light sourcesField experiments were conducted at three locations in Japan: Central Region Agricultural Research Center (CARC), Hokuriku Research Station (37° 07′ 00″ N, 138° 16′ 23″ E) in Niigata; Yamaguchi Prefectural Agriculture & Forestry General Technology Center (YPATC) (34° 09′ 37″ N, 131° 29′ 47″ E) in Yamaguchi; and Okinawa Prefectural Agricultural Research Center (OPARC) (26° 06′ 18″ N, 127° 40′ 53″ E) in Okinawa. The distribution of Nezara spp. varies among the regions in Japan. Only N. antennata is distributed in Niigata, and only N. viridula is distributed in Okinawa. Both N. antennata and N. viridula were found in Yamaguchi.Experiment 1: Attractiveness of UV light at different intensitiesField experiments to evaluate the attractiveness of UV light at different intensities were conducted from August 2 to 29, 2017, around a soybean field at the CARC in Niigata and from July 12 to September 9, 2019, in grassland at the OPARC in Okinawa. Light traps with different numbers of UV-LEDs (12, 21, 42, and 84) were used as light sources. Each of the four LED traps was spaced more than 30 m apart and placed randomly around the soybean field or grassland. Due to time constraints, the numbers of N. viridula and N. antennata captured in traps were counted every 3–4 days at Niigata (total eight replicates) and every 7 days in Okinawa (total eight replicates). The traps were randomly repositioned every week to minimize the effect of trap location. The raw capture data for each trap are listed in Supplementary Table S1.Experiment 2: Attractiveness of green light at different intensitiesField experiment to evaluate the attractiveness of green light at different intensities was conducted from July 5 to August 5, 2019, around a soybean field at the YPATC in Yamaguchi. Light traps with different numbers of green LEDs (12, 21, 42, and 84) were used as light sources. Light trap with 84 UV-LEDs was used as the positive control. Each of the five LED traps was spaced more than 30 m apart and placed randomly around the soybean field. The numbers of Nezara bugs captured in traps were counted every 3–4 days (total nine replicates). The traps were randomly repositioned every week. The raw capture data for each trap are listed in Supplementary Table S2.Experiment 3: Attractiveness of combined-UV and green lightField experiments to evaluate the attractiveness of combinations of UV- and green-LEDs were conducted from June 13 to September 4, 2017, in the grassland at the OPARC in Okinawa, and from July 15 to September 1, 2017, around a soybean field at the YPATC in Yamaguchi. Light traps with 84 UV-LEDs, 84 green-LEDs, and a combination of 42 UV-LEDs and 42 green-LEDs were used as light sources. Each of the three LED traps was spaced more than 30 m apart and placed randomly around the soybean field or grassland. Although insects other than Nezara bugs (mainly coleopteran species) were captured in the light traps, for soybean pests, the funnel-type light traps are intended for monitoring large coleopteran and heteropteran insects ( > 1 cm). Therefore, we targeted and counted insects that meet these conditions. Statistical analysis was performed on species with a total capture number of more than 20 individuals in the three traps. The species were as follows: in addition to Nezara bugs, heteropteran bugs, Piezodorus hybneri (Gmelin), Glaucias subpunctatus (Walker), Halyomorpha halys (Stål), and Plautia stali Scott, as well as coleopteran beetles, Anomala albopilosa (Hope), A. cuprea Hope, A. rufocuprea Motschulsky, and Holotrichia parallela Motschulsky. The insects captured in traps were counted for each species every 7 days at Okinawa (total 12 replicates) and every 3–4 days at Yamaguchi (total 14 replicates). The traps were randomly repositioned every week. The raw capture data for each trap are listed in Supplementary Table S3.Data analysisIn Experiment 1, the effect of UV light intensities for trap catches were analyzed using a nonparametric one-tailed Shirley–Williams test under an assumption that higher light intensity attracts larger amounts of insects. In Experiment 2, the effect of green light intensities for trap catches were analyzed using the Shirley–Williams test. Subsequently, the attractiveness of each green light was compared to that of UV light using Wilcoxon matched pairs signed-rank test. In Experiment 3, the effect of light sources for trap catches was analyzed using the Friedman test, followed by the Wilcoxon signed-rank test, with Bonferroni correction for multiple comparisons. Statistical analyses were performed using R version 4.2.0 (R Core Team, 2022). More

All figures were produced using the R package ggplot2 v3.3.534.eBird and covariate dataeBird data are structured as follows. Birders submit observations as species checklists with counts of each species they identify. They report associated metadata, such as location, date and time, duration of the observation period, number of observers, and sampling protocol25,26,31. The birder indicates whether their checklist is “complete”; complete checklists yield inferred zeroes for all species not reported on a checklist.We retrieved the eBird Basic Dataset containing all eBird observations and sampling metadata. We extracted all complete checklists that occurred within the U.S. state of California between April 1 and June 30, 2019. Four survey-level covariates were retrieved from eBird checklist metadata as detection covariates: number of observers, checklist duration, date of year, and time of day; any checklist that failed to report one or more of these variables was dropped. Corresponding to best practices for use of eBird data, we filtered the data for quality according to the following criteria: we discarded checklists other than those following the “Stationary” survey protocol (observations made at a single spatial location) with duration shorter than 4 hours and at most 10 observers in the group31,35.We selected twenty circular regions of high sampling intensity with 10 km radii across California (Supplemental Fig. 1). These spanned the state’s many habitats including coastal, agricultural, wetland, and mountain areas, and contained active birding areas such as parks and human population centers. In each subregion, we selected 10 species with the highest reporting rate (proportion of checklists including that species) and 10 representing an intermediate reporting rate. An additional 10 species were selected that were detected in many regions to enable cross-region comparisons, yielding 407 species-subregion (SSR) datasets (with overlaps between the two species selection protocols; see Supplemental Section 2 for the full algorithm). Across 20 subregions, we accepted 6094 eBird checklists for analysis, each with an associated count (potentially zero) for each species. Observations were aggregated to sampling sites defined by a 50 m spatial grid. The 50 m grid was chosen to conservatively identify related surveys and was not motivated by biological processes, nor does it represent the sampling area of each survey. In this context, the concept of “closure” in the latent state is already suspect due to the fact that eBird checklist sampling areas are inconsistent. Data were processed in R using the ‘auk’ package36,37.An elevation surface for the state of California was retrieved from WorldClim at (8.3 times 10^{-3}) decimal degrees resolution using the R package raster38,39. This commonly used covariate was included as a baseline spatial covariate to enable comparison of estimation properties across sites, but its biological relevance to abundance is not crucial to our analysis31. Land cover data were retrieved from the LandFire GIS database’s Existing Vegetation Type layer40. For each unique survey location, a 500 m buffer was calculated around the reported location, and the percent of the buffer which was water, tree cover, agriculture or other vegetation (shrub or grassland) was calculated. We used the following five site-level covariates: elevation, and percent of the landscape within a 500 m buffer of the site that was water, trees, agricultural land, or other vegetation. We included six checklist-level covariates: duration, number of observers, time of day, time of day squared, Julian date, and Julian date squared. Covariates were dropped in datasets where only a single unique value was observed for that covariate.Model implementation and selectionWe considered four variants of the N-mixture model and two variants of the GLMM comprising a total of 6 distinct models, defined by the distributions used in the model or sub-model.The GLMM for count data that we considered is defined as$$begin{array}{*{20}l} {y_{{ij}} sim D(mu _{{ij}} ,[theta ])} hfill \ {log (mu _{{ij}} ) = beta _{0} + {mathbf{x}}_{{ij}}^{T} user2{beta } + alpha _{i} } hfill \ {alpha _{i} sim {mathcal{N}}(0,sigma _{alpha } )} hfill \ end{array}$$where (y_{ij}) is th jth observation at site i, D is a probability distribution (which may contain an extra parameter (theta) to account for overdispersion), (mu _{ij}) represents the mean expected count and is a logit-linear combination of observed site- and observation-level covariates (x_{ij}), (beta) are coefficients representing the effect of those covariates, (beta _0) is a log-scale intercept corresponding to the expected log count at the mean site (i.e. with all centered covariates set to 0), and (alpha _i) is the random effect of site i following a normal distribution. Due to the right skew of (exp (y_{ij})), by log-normal distribution theory the log of the expected count at the mean site is (beta _0 + 0.5 sigma _{alpha }^2). We considered two forms of this model, where D was either a Poisson or a negative binomial distribution, in the latter case with the extra parameter (theta).The N-mixture model is defined as$$begin{array}{*{20}l} {y_{{ij}} sim D_{w} (N_{i} ,p_{{ij}} ,[theta _{w} ])} hfill \ {N_{i} sim D_{b} (lambda _{i} ,[theta _{b} ])} hfill \ {{text{logit}}(p_{{ij}} ) = {text{}}{text{logit}}(p_{0} ) + {mathbf{x}}_{{ij(w)}} {mathbf{beta }}_{w} } hfill \ {log (lambda _{i} ) = log (lambda _{0} ) + {mathbf{x}}_{{i(b)}} {mathbf{beta }}_{b} } hfill \ {p_{0} = e^{{frac{{phi _{1} + phi _{2} }}{2}}} } hfill \ {lambda _{0} = e^{{frac{{phi _{1} – phi _{2} }}{2}}} } hfill \ end{array}$$where (D_b) and (D_w) are probability distributions representing between- and within-site variation, respectively; (N_i) is a site-level latent variable normally representing the “true” abundance at site i; (p_{ij}) is the detection probability of each individual on the jth observation event at site i; (lambda _i) is the mean abundance at site i; and (x_{(w)}) and (x_{(b)}) are covariate vectors for detection and abundance, respectively, with corresponding coefficients (beta _w) and (beta _b). For reasons described below, we reparameterize the intercept parameters of the N-mixture submodels, (log (lambda _0)) and (text{ logit }(p_0)), in terms of two orthogonal parameters (phi _1 = log (lambda _0 p_0)) and (phi _2=log (p_0 / lambda _0)). Now (phi _1) and (phi _2) represent the expected log count and the contrast between detection and abundance, respectively, at the mean site. This parameterization allows us to investigate stability of parameter estimation. The log-scale expected count of the N-mixture model is (phi _1 = log (lambda _0 p_0)), analogous to (beta _0 + 0.5 sigma _{alpha }^2) in the GLMM (see Supplemental Section 6). Each submodel distribution D could include or not include an overdispersion parameter ((theta _w) and (theta _b)), yielding four possible N-mixture model variants: binomial-Poisson (B-P), binomial-negative binomial (B-NB), beta-binomial-Poisson (BB-P), and beta-binomial-negative binomial (BB-NB)8,11.We chose to fit models with maximum likelihood estimation (MLE) for computational feasibility and because key diagnostic tools, such as AIC and methods for checking goodness-of-fit and autocorrelation, were best suited to MLE estimation15. We fit N-mixture models with the nimble and nimbleEcology R packages starting with a conservatively large choice of K, the truncation value of the infinite sum in the N-mixture likelihood calculation33,41 (see Supplemental Section 4 for a discussion of maximum likelihood estimation with NIMBLE). We fit GLMMs with the R package glmmTMB42. We applied forward AIC selection to choose the best covariates for each model with each dataset (illustrated in Fig. S1). One spatial covariate (elevation) and two checklist metadata covariates (duration and number of observers) were treated as a priori important and were included in all models. In the N-mixture model, checklist-specific sampling metadata were only allowed in the detection submodel, while land cover covariates and the interactions between them were allowed in both the detection and abundance submodels. Interactions were dropped in datasets when interaction values showed a correlation of > 0.8 with one of their first-order terms. In N-mixture models, additions to both submodels were considered simultaneously during forward AIC selection.For comparisons between models, we selected a heuristic threshold of (Delta text {AIC} > 2) to say that one model is supported over another30.Fit, estimation, and computationGoodness-of-fitWe used the Kolmogorov-Smirnov (KS) test, a p-value based metric, to evaluate goodness-of-fit on each selected model. For GLMMs, residuals were obtained using the DHARMa R package’s ‘simulateResiduals’ and the KS test was applied using the ‘testUniformity’ function43. For N-mixture models, we considered the site-sum randomized quantile (SSRQ) residuals described by Knape et al.15, computing these for each N-mixture model and running a KS test against the normal CDF. We assumed that covariate effects did not vary by space within subregions and chose not to use spatially explicit models31,44. To test this assumption, we applied Moran’s I test to the SSRQ or DHARMa-generated residuals for each site or observation.Parameter estimationWe compared two abundance parameters of interest across models: coefficients for elevation and log expected count at a standard site (in the GLMM, (beta _0 + 0.5 sigma _alpha ^2); in the N-mixture model, (log (lambda _0 p_0))). We examined absolute differences in point estimates and the log-scale ratios between their standard errors.Stability of estimated parametersAttempting to decompose the expected value of observed data into within- and between-site components can lead to ridged likelihood surfaces with difficult-to-estimate optima. Kéry found that instability of model estimates with increasing K occurred when there was a likelihood tradeoff between detection and abundance, resulting in a tendency in abundance toward positive infinity restrained only by K10. Dennis et al. showed that N-mixture models could in fact yield estimates of absolute abundance at infinity18. We interpreted this as a case of a boundary parameter estimate rather than non-identifiability and explored it by reparametrizing as follows. We estimated the intercepts for detection and abundance with two orthogonal parameters (rotated in log space) (phi _1 = log (lambda _0 p_0)) and (phi _2 = log (p_0 / lambda _0)), where (lambda _0) and (p_0) are real-scale abundance and detection probability at the mean site. We hypothesized that in unstable cases, (phi _1), log expected count, is well-informed by the data, but (phi _2), the contrast between abundance and detection, is not well-informed, corresponding to a likelihood ridge as (phi _2 rightarrow -infty) due to detection probability approaching 0 and abundance approaching infinity. This reparameterization isolates the likelihood ridge to one parameter direction, similar to a boundary estimate as (exp (phi _2) rightarrow 0). Boundary estimates occur in many models and are distinct from non-identifiability in that they result from particular datasets. Confidence regions extending from a boundary estimate may include reasonable parameters, reflecting that there is information in the data. We defined a practical lower bound for (phi _2). When (phi _2) was estimated very near that bound, we conditioned on that boundary for (phi _2) when estimating confidence regions for other parameters.In the N-mixture case, diagnosing a boundary estimate for (phi _2) is made more difficult by the need to increase K for large negative (phi _2) to calculate the likelihood accurately. We used an approach like that of Dennis et al.18 to numerically diagnose unstable cases. For each N-mixture variant in each SSR, the final model was refitted twice, using values of K 2000 and 4000 greater than the initial choice. Estimates were considered unstable if the absolute value of the difference in AIC between these two large-K refits was above a tolerance of 0.1. We monitored whether MLE estimates of (phi _1) and (phi _2) also varied with increasing K.Evaluating the fast N-mixture calculationWe extended previous work by Meehan et al. to drastically improve the efficiency of N-mixture models using negative binomial or beta-binomial distributions in submodels45 (see Supplemental Section 3).We ran benchmarks of this likelihood calculation for a single site against the traditional algorithm, which involves iterating over values of N to compute a truncated infinite sum. We calculated the N-mixture likelihood at 5,000 sites and compared the computation time between the two methods for all four N-mixture model variations. We ran benchmarks along gradients of (text {length}(y_i)) (number of replicate observations at the simulated site) and K (the upper bound of the truncated infinite sum) for each variant. More

Massive drift of pumice along the northeastern coast of Okinawa IslandA large amount of pumice stones reached and was deposited along the northeastern coast of Okinawa Island, that were brought by strong seasonal northeasterly winds (Supplementary Video 1). The pumice was thought to be brought by the Kuroshio countercurrent from sites near the Ogasawara Archipelago 1300 km away. Because the Kuroshio countercurrent is composed of various medium-sized eddies in the ocean, the current does not always flow in one direction and as a continuous flow27,28. The pumice drift was more strongly controlled by the seasonal northwesterly winds to be transported to Okinawa across the Philippine Sea (Fig. 1a). The pumice raft reached the northern part of Okinawa approximately 2 months after the eruption (Figs. 2, 3 and 4). According to a very recent report, the pumice clasts were drifting ashore in Thailand (traveling 4000 km-long distance) across the South China Sea within half a year of this eruption29. Most pumice stones were gray, but some pumice was banded, and others were black reflecting some compositional variation25,29 (Fig. 2d,e). The Kuroshio Current is faster than the Kuroshio countercurrent27, so some pumice clasts have already reached the main island of Japan25. Tracking the dispersal of the pumice will allow a better forecasting model based on observed raft trajectories by considering exact wind effects in the Philippine Sea30.Figure 2An example of a natural beach on Okinawa Island where pumice has washed ashore. (a) Appearance of natural sandy beaches on the northern part of Okinawa Island (Ibu beach, Kunigami Village, 26°75′57.88″ N, 128°32′23.32″ E). Photo was taken on 24 October 2021. Pumice drifted onto the sandy beach and formed a striped pattern. The white-capped waves indicate on the place where the reef edge exist. The white arrow points to the mangrove river estuary corresponding to Fig. 9. (b) Estimation of the pumice sedimentation depth on the original sand beach surface. (c) The high tide zone of the natural sandy beach is covered with pumice pebbles and stones. Yellow arrows indicate black pumice stones. Scale bar: 10 cm. (d, e) Front and back of examples of relatively large pumice stones from the same beach. The left image is mostly light brown, whereas the right image is almost black. Scale bars: 5 cm.Full size imageFigure 3Short-term migration of pumice from beaches as revealed by stationary observations. These four photos were taken at two sites on northern Okinawa Island on two consecutive days, 23 and 24 October 2021. (a, b) A sandy beach along the Sate Coast (26°78′84.56″ N, 128°22′30.57″ E). It was windy on the first day, and pumice stones were washed up with the waves. Almost all the pumice stones were removed from the beach and transported offshore on the following day. The black arrow in photo (a) indicates Cape Hedo, the northernmost tip of Okinawa Island. (c, d) At this gravelly beach (26°80′83.25″ N, 128°23′38.56″ E), pumice fully covers the seawall on the first day, but all of the pumice stones washed away, leaving the original gravels, on the following day. The white arrow in each photo indicates an identical marker stone placed on the beach. Weather data of northern Okinawa (https://www.data.jma.go.jp/obd/stats/etrn/view/daily_a1.php?prec_no=91&block_no=0901&year=2021&month=10&day=23&view=g_wsp) and tidal data (Naha: 26°13′ N, 127°40′ E) (https://www.data.jma.go.jp/gmd/kaiyou/db/tide/genbo/genbo.php) are provided by Japan Meteorological Agency.Full size imageFigure 4Pumice stones settled by marine organisms. (a) Pumice collected from Ibu beach on 31 October 2021. Two marine benthos coexist close together on a pumice stone. Scale bar: 1 cm. (b) Enlarged image of the Lepas barnacle. Scale bar: 3 mm. (c) Enlarged image of the bryozoan. Scale bar: 3 mm. (d) Stereo microscopic image of pumice pebbles of a few millimeters in diameter collected from Ibu beach on 15 January 2022. The light brown coloration indicates some algal/cyanobacterial growth on the pumice. Scale bar 1 mm. (e) Red autofluorescence was detected from pumice pebbles. Image corresponds to (d). Autofluorescence from microalgae was confirmed by Supplementary Fig. 2. Scale bar 1 mm. (f) Enlarged image of the center of the figure of (e) shows red autofluorescent signals with a diameter of 10–30 µm. Scale bar: 200 µm.Full size imageChanges in the coastal landscape: natural beaches and estuariesMarine calcifiers, including corals, calcareous algae, and foraminifers, produce white sandy beaches on Okinawa Island. However, the gray pumice drifting ashore changed the white sand beach, especially along the northeastern coastline. We observed several lines of pumice aggregations, suggesting that pumice was brought ashore by wavefronts several times produced by a strong north wind at the tide lines (Supplementary Video 1; Fig. 2a). At the same sampling site, the thickest depth of beached pumice was more than 30 cm (Fig. 2b; Supplementary Video 2). Most of the pumice stones were from 0.5 cm to 3 cm in diameter, with a few black pumice stones included (Fig. 2c: yellow arrow). Pumice stones arrived at the estuaries of some brackish rivers (Fig. 8, Supplementary Fig. 1a) and mangrove forests in northwest Okinawa (Fig. 9).Pumice stones and pumice rafts show dynamic behavior in a short period. We captured photographs 24 h apart at two positions on the shore of Okinawa, which allowed us to compare the pumice dynamics during this period (Fig. 3). Within that time frame, there were two high tides, and the tide level changed by up to 170 cm. As seen in Fig. 3a, on the first day, the coast was covered with pumice, and floating pumice could be seen on the seafront. The north wind was strong that day, as shown by the relatively high waves near the shore as well as white‐crested waves near the reef edge. By the following day, most of the pumice had been moved offshore by tides and winds (Fig. 3b), indicating that newly beached pumice raft deposits were removed quickly from open beach areas. At another site on a gravelly beach, pumice fully covered the seawall on the first day, but almost all of the pumice stones were washed away, leaving the original gravels, on the following day (Fig. 3c,d). Japan Meteorological Agency (Oku station: 232 m above sea level, latitude 26°50.1, longitude 128°16.3′) reported that northerly winds were blowing (mean wind speed: 3.4 m/s) on 23rd October in northern Okinawa. The following day, the wind direction changed to the east-southeast; blowing offshore (mean wind speed: 2.9 m/s), resulting in the dramatic removal of pumice form the coast (Fig. 3). These observations indicate that surface winds rather than ocean currents had a strong influence on the raft trajectory and residence time on beaches, and are consistent with past research5. These observations lead us to expect that the pumice rafts will disappear from the coast of Okinawa fairly quickly, but in fact, there have been many cases where they have come back again in a few days. Although the overall amount of pumice drifting has been decreasing, a small amount of pumice has been drifting in coastal area of Okinawa in May, 202231. It is unlikely that large amounts of pumice will drift repeatedly throughout Okinawa Prefecture as reported in this report, but it should be noted that detached pumice material remains in beach and river runoff.Biofouling of sessile organisms on pumice arriving to OkinawaIt is noteworthy that the pumice rafts traveled over the deep Philippine Sea for over 2 months, and on arrival in Okinawa there was little to no biofouling of the pumice (Fig. 2). Some stranded pumices observed on Okinawa beaches had become habitats for sessile organisms (Fig. 4), as reported in previous studies1,2,3,4,5,6,29. Goose barnacles (Lepas sp.) without external damage to the shell were the most abundant species observed on the pumice (Fig. 4b). Lepas is a common biofouling taxon distributed globally and plays a role in biofouling as a foundation organism. The shell growth rate is more than 1 mm/day in some Lepas species32 suggesting that the Lepas had been growing on the pumice for about two weeks. Measurements of the shell size of Lepas attached to the pumice collections conducted in the same area (Supplementary Video 2) showed a bias toward larger sizes in the second collection (5.92 ± 3.86 mm (average ± S.D.), n = 75, 13 November 2021) than in the first one (3.43 ± 1.08 mm, n = 21, 31 October 2021), and significant differences were detected between the measurement periods (Mann–Whitney U test, p  More

Temporal differentiation of wild populations of T. hispidus and climatic parametersComparison of morphometric features of T. hispidus shells collected in different years in two geographic regions, i.e., Wrocław and Lubawka, showed significant differences depending on the year of collection. The largest number of differences was revealed in shells from Wrocław (Figs. 1 and 2A; Additional file 2: Table S1). Out of 210 comparisons (15 pairs of collection years × 14 features), 84 were statistically significant (Additional file 2: Table S2), e.g., shell diameter (D) was significantly different in 11 cases, shell height (H) and shell width (W) in 10 cases, body whorl height (bwH), the number of whorls (whl), umbilicus major (U) and minor (u) diameters in 9 cases and aperture height/width ratio (h/w) in 7 cases. Nine features obtained more than 10% difference between shells in at least one comparison of mean values, e.g., U 24%, u 19%, H 16% and D 15% (Additional file 2: Table S2). Umbilicus major (U) and minor (u) diameters showed the largest average percentage difference, i.e., 12% and 10%, respectively, in comparisons of all years.Figure 1Shells of Trochulus hispidus collected in different years in Wrocław.Full size imageFigure 2Changes in: mean values of selected morphometric features of shells collected in various years in Wrocław (A) as well as the mean temperature (B) and the relative humidity (C) recorded in four seasons in Wrocław in eight-year period. Abbreviations: D—shell diameter (in mm), H—shell height (in mm), h/w—aperture height/width ratio, whl—number of whorls. The summary statistics for A is included in Table S1 and original data in Table S10 in Additional file 2.Full size imageFor snails from Lubawka, out of 84 comparisons (6 pairs of collection years × 14 features) only 8 were statistically significant (Additional file 2: Table S3). The shells differed significantly in their aperture height (h) and width (w) in 3 comparisons. The h feature showed the percentage difference up to 9% (Additional file 2: Table S3) and the largest average difference was 4.5%.Besides the phenotypic variation, climatic parameters also showed high fluctuations in the studied period (Fig. 2B,C, Additional file 2: Table S4). The maximum difference reported between temperature parameters in some years prior to sample collection in Wrocław was up to 3.7 °C for the maximum winter temperature, while the maximum difference in the relative humidity was up to 11% for autumn. The maximum temperature difference in Jelenia Góra close to Lubawka was up to 3.5 °C for the minimum winter temperature, while the relative humidity differed at most by up to 8% in summer.Differences in shell morphometry under various climatic conditionsThe distinction between shells collected in individual years and changes in climatic parameters along the same period suggest that these differences can be associated with the climate. Therefore, we calculated the average value of a given climatic parameter for each season and studied region and next divided the collected shell data into two groups according to this value. The first group included the shells that developed in conditions above this average and the second below this average (Additional file 2: Table S5). The differences between these groups were statistically significant for 15 out of 16 considered climatic parameters for at least two shell features (Fig. 3). Similarly, each of 14 features significantly separated the groups based on at least two climatic conditions. The results demonstrated that the mean winter temperature substantially influenced nine morphometric shell features, whereas eight characters were changed due to the maximum winter temperature as well as the mean and minimum temperatures in spring, summer and autumn. Umbilicus major (U) and minor (u) diameters as well as umbilicus relative diameter (U/D) were significantly different in 14 pairs of groups characterized by various climatic parameters. In 11 pairs, the height/width ratio (H/W) was significantly different and shell height (H) in 10 pairs.Figure 3Mean percentage differences in morphometric features between shells that were grown in different conditions. The shells were divided into two groups according to the average value of a given climatic parameter for each season and studied region. The first group included the shells that developed in conditions above this average and the second below this average. Positive values indicate that the given feature was greater in the first group, whereas negative values indicate that this feature was greater in the second group. Dendrograms cluster the features and the parameters according to their similarity in the percentage differences. Values marked in bold indicate statistically significant differences between the compared groups of shells. Values at the dendrogram nodes indicate significance assessed according to approximately unbiased test (au) and bootstrap resampling (bp).Full size imageThe umbilicus diameters (u and U) as well as umbilicus relative diameter (U/D) clustered together in the dendrogram based on the mean percentage difference, which indicates that they similarly responded to climatic conditions (Fig. 3). The features u and U revealed the strongest average increase of all features, from 4.1 to 10.5% in shells developed in higher temperatures in all seasons. The largest percentage difference exceeding 10% was recorded for groups separated according to the mean summer and autumn temperatures as well as the maximum summer and minimum autumn temperatures. The U/D ratio was also significantly greater with the mean percentage difference of 2.8–7.6% in shells grown under high temperatures for all seasons and almost all temperature types. On the other hand, the u and U diameters as well as the U/D ratio were on average by 3.7–6.0% significantly smaller in shells developed under higher humidity in summer and winter.The height/width shell ratio (H/W) was grouped with H and bwH features in the dendrogram and was on average by up to 3.6% significantly smaller in shells grown under higher temperatures in all seasons for almost all types of parameters. The maximum winter temperature caused a significant increase, on average by ca. 3%, in shell height (H) and body whorl height (bwH), whereas higher temperatures in other seasons led to their decrease by up to 3.4% (Fig. 3).The shells that were grown in autumn with a relatively high maximum temperature were characterized by ca. 3% significantly smaller aperture height (h) and aperture height/width ratio (h/w), which were clustered together in the dendrogram (Fig. 3).Other four features, shell diameter (D), number of whorls (whl) as well as shell (W) and aperture width (w), formed an additional cluster in the dendrogram (Fig. 3). All of them were on average significantly greater in shells collected one year after winter that was characterized by relatively higher mean and maximum temperatures. The percentage difference was greater, with 3.6–3.9% for W and D.In the dendrogram, the climatic parameters were clustered in several groups indicating their similar influence on the morphometric features of shells (Fig. 3). There are separate clusters for temperature and humidity parameters with the exception of the autumn maximum temperature and autumn humidity, which are grouped together. The other temperature parameters for warmer seasons are separated from those for winter, which indicates that they differently influenced the shell morphometry.Correlations between morphometric shell features and climatic parametersThe influence of climatic conditions on the shells collected in individual years was also assessed using Spearman’s correlation coefficient between the morphometric features and climatic parameters (Fig. 4). Of 224 potential relationships 113 were statistically significant. The spring mean temperature was significantly correlated with 10 morphometric features. Summer humidity and six temperature parameters, i.e., the minimum temperatures as well as the spring and winter maximum temperatures, significantly correlated with eight shell features. Minor umbilicus diameter (u) and umbilicus relative diameter (U/D) were significantly correlated with almost all climatic parameters, i.e., 15, umbilicus major diameter (U) and height/width ratio (H/W) with 13 and the ratio of umbilicus minor to its major diameter (u/U) with 11.Figure 4Spearman’s correlation coefficients between morphometric features of shells with climatic parameters under which the snails were grown. Dendrograms cluster the features and the parameters according to their similarity in the coefficients. Values marked in bold are statistically significant. Values at the dendrogram nodes indicate significance assessed according to approximately unbiased test (au) and bootstrap resampling (bp).Full size imageAs in the case of percentage difference, we can also recognize groups of morphometric features that were similarly correlated with climatic parameters (Fig. 4). Features U/D, U and u were significantly positively correlated with all or almost all temperature parameters for four seasons with the coefficients up to 0.34, 0.30 and 0.36, respectively. On the other hand, the significant correlation coefficients between these features and the humidity in spring, summer and winter were negative and reached − 0.34.Another group of features included shell height/width ratio (H/W), shell height (H) and body whorl height (bwH) (Fig. 4). All of them showed significant negative correlations with all temperature parameters for spring and summer as well as the minimum autumn temperature, and H/W also with the mean and maximum autumn temperatures as well as the mean and minimum winter temperatures. The correlation coefficients reached − 0.28, − 0.27 and − 0.28, respectively. These three features significantly correlated with summer and spring humidity, at up to 0.23.The number of whorls (whl), shell width (W), shell diameter (D), demonstrated a similar correlation with climatic parameters (Fig. 4). They showed the largest and significant correlation coefficients with winter temperatures: up to 0.24, 0.22 and 0.22, respectively. The ratio of umbilicus minor to its major diameter (u/U) showed significant positive correlation up to 0.22 with temperature of warmer seasons.The climatic parameters were grouped into several clusters indicating their similar relationships with morphometric features (Fig. 4). Humidity parameters of warmer seasons formed a separate cluster and temperature parameters were grouped according to seasons. The winter parameters were connected with autumn humidity and separated from temperatures for warmer seasons.Modelling relationships between morphometric shell features and climatic parametersThe joint influence of many climatic parameters on morphometry of shells collected in individual years was studied using a linear mixed-effects (LME) model after exclusion of correlated parameters and a linear ridge regression (LRR) model including all climatic parameters. The latter allows for the inclusion of correlated variables. We separately investigated the seasonal maximum, mean and minimum temperature parameters in combination with seasonal humidity parameters (Additional file 2: Table S6) because they are obviously correlated.Umbilicus minor (u) and major (U) diameters as well as umbilicus relative diameter (U/D) turned out best explained by the climatic parameters (with R2  > 0.15) in two models (Additional file 2: Table S6). Moreover, u, U and U/D were described in LME models by the largest number of significant climatic parameters, i.e., 15. The features u and U had also the largest number of significant parameters in LRR models, i.e., 18 out of 24 possibilities. The largest average values of temperature coefficients for the LRR models were 0.66 for D, 0.58 for W, 0.32 for H, 0.26 for U and 0.22 for u. Thus, all the above-mentioned features were under the strongest influence of the climatic conditions.In the case of LRR models, the coefficients at the winter mean temperature were most often selected as significant, in 12 out of 14 possibilities (Additional file 2: Table S6). The humidity coefficients for autumn were significant in 30 cases of 42 possibilities. The highest average absolute values of coefficients in climatic variables were those for the summer (0.63), spring (0.31) and autumn (0.24) minimum temperatures as well as the summer mean temperature (0.31). Thus, the temperatures of warmer seasons were more important for developing shell morphology. Seasonal humidity coefficients showed similar values compared to each other.Comparison of shell morphometry of T. hispidus and T. sericeus kept under various conditionsIn order to verify the influence of different climatic parameters on Trochulus shell morphometry in selected conditions, we compared shells from three groups of T. hispidus, which represented several subsequent generations: (1) parental snails collected in the wild in Wrocław-Jarnołtów, (2) their offspring bred in the laboratory for two generations and (3) offspring of the second laboratory-bred generation transplanted again into a garden in Wrocław (Fig. 5A–C). The comparison of the group 2 and 1 was to verify if laboratory conditions with controlled temperature and humidity can influence the shell morphometry within only one generation, whereas including the group 3 in the comparison, we wanted to check if snails raised in wild garden conditions can recover the original phenotype. Furthermore, we transplanted into the same garden conditions T. sericeus, which was collected in the wild in Muszkowice (Fig. 5D,E). In this case, we verified if two originally different ecophenotypes T. hispidus and T. sericeus, develop the same shell morphometry under the same conditions.Figure 5Shells of two Trochulus ecophenotypes: parental T. hispidus from wild habitat in Wrocław (A), the first generation of T. hispidus raised in laboratory (B); T. hispidus reared in garden in Wrocław (C); T. sericeus from wild habitat in Muszkowice (D); T. sericeus reared in garden in Wrocław (E).Full size imageConditions in which these snails developed were different. According to WorldClim, the wild environment of T. hispidus in Wrocław was generally warmer than that of T. sericeus in Muszkowice (Additional file 2: Table S7). The largest difference was 1.4 °C for the maximum summer temperature. Relative humidity was lower in Wrocław by up to 2% for warmer seasons but was higher in winter by 1.6%. The difference between the wild and garden localities in Wrocław was much smaller and did not exceed 0.41 °C. The garden conditions were less humid, by up to 2%. However, data from WorldClim are generalized over a longer period and wider regions, so may not well reflect local conditions in the studied places. Actually, the Wrocław site was an open habitat covered with a nettle community like a garden patch, while the Muszkowice site was overgrown by a beech forest, which most likely maintained a higher humidity and a more stable temperature.Laboratory temperatures were substantially different from those in the field, especially for winter (by 18–19.7 °C) as well as for spring and autumn (by 8.2–12 °C). Laboratory humidity was by up to 4.5% lower compared to winter and 5.9–9.9% higher than in spring and summer.A discriminant function analysis (DFA) for the defined groups of snails provided their interesting grouping and separation (Fig. 6). The analysis identified three significant discriminant functions (p  More