Philippa Kaur

These authors contributed equally: Yannick Garcin, Enno Schefuß, Greta C. Dargie, Simon L. LewisAix Marseille University, CNRS, IRD, INRAE, CEREGE, Aix-en-Provence, FranceYannick Garcin & Ghislain GassierInstitute of Geosciences, University of Potsdam, Potsdam, GermanyYannick GarcinMARUM—Center for Marine Environmental Sciences, University of Bremen, Bremen, GermanyEnno SchefußSchool of Geography, University of Leeds, Leeds, UKGreta C. Dargie, Bart Crezee, Dylan M. Young, Andy J. Baird, Paul J. Morris & Simon L. LewisSchool of Geography and Sustainable Development, University of St Andrews, St Andrews, UKDonna Hawthorne, Ian T. Lawson & George E. BiddulphIFP Energies Nouvelles, Earth Sciences and Environmental Technologies Division, Rueil-Malmaison, FranceDavid SebagInstitute of Earth Surface Dynamics, Geopolis, University of Lausanne, Lausanne, SwitzerlandDavid SebagFaculté des Sciences et Techniques, Université Marien Ngouabi, Brazzaville, Republic of the CongoYannick E. Bocko & Y. Emmanuel Mampouya WeninaÉcole Normale Supérieure, Université Marien Ngouabi, Brazzaville, Republic of the CongoSuspense A. IfoÉcole Normale Supérieure d’Agronomie et de Foresterie, Université Marien Ngouabi, Brazzaville, Republic of the CongoMackline MbembaFaculté de Gestion des Ressources Naturelles Renouvelables, Université de Kisangani, Kisangani, Democratic Republic of the CongoCorneille E. N. Ewango & Joseph Kanyama TabuFaculté des Sciences, Université de Kisangani, Kisangani, Democratic Republic of the CongoCorneille E. N. EwangoInstitut Supérieur Pédagogique de Mbandaka, Mbandaka, Democratic Republic of the CongoOvide Emba & Pierre BolaSchool of Geography, Geology and the Environment, University of Leicester, Leicester, UKGenevieve Tyrrell, Arnoud Boom & Susan E. PageSchool of Water, Energy and Environment, Cranfield University, Bedford, UKNicholas T. GirkinBritish Geological Survey, Centre for Environmental Geochemistry, Keyworth, UKChristopher H. VaneInstitute of Earth Sciences, University of Lausanne, Lausanne, SwitzerlandThierry AdatteNEIF Radiocarbon Laboratory, Scottish Universities Environmental Research Centre (SUERC), Glasgow, UKPauline GulliverSchool of Biosciences, University of Nottingham, Nottingham, UKSofie SjögerstenDepartment of Geography, University College London, London, UKSimon L. Lewis More

Sket, B. Can we agree on an ecological classification of subterranean animals?. J. Nat. Hist. 42, 1549–1563. https://doi.org/10.1080/00222930801995762 (2008).Article 

Google Scholar 
Mammola, S. et al. Scientists’ warning on the conservation of subterranean ecosystems. Bioscience 69, 641–650. https://doi.org/10.1093/biosci/biz064 (2019).Article 

Google Scholar 
Boulton, A. J., Fenwick, G. D., Hancock, P. J. & Harvey, M. S. Biodiversity, functional roles and ecosystem services of groundwater invertebrates. Invertebr. Syst. 22, 103–116. https://doi.org/10.1071/IS07024 (2008).Article 

Google Scholar 
Danielopol, D. L. & Griebler, C. Changing paradigms in groundwater ecology—From the ‘living fossils’ tradition to the ‘new groundwater ecology’. Int. Rev. Hydrobiol. 93, 565–577. https://doi.org/10.1002/iroh.200711045 (2008).Article 

Google Scholar 
Griebler, C., Malard, F. & Lefébure, T. Current developments in groundwater ecology—From biodiversity to ecosystem function and services. Curr. Opin. Biotechnol. 27, 159–167. https://doi.org/10.1016/j.copbio.2014.01.018 (2014).Article 
PubMed 

Google Scholar 
Fišer, C. Niphargus—A model system for evolution and ecology. In Encyclopedia of Caves (eds Culver, D. C. et al.) 746–755 (Academic Press, 2019).Chapter 

Google Scholar 
Riddle, M. R. et al. Insulin resistance in cavefish as an adaptation to a nutrient-limited environment. Nature 555, 647–651. https://doi.org/10.1038/nature26136 (2018).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Gibert, J. et al. Assessing and conserving groundwater biodiversity: Synthesis and perspectives. Freshw. Biol. 54, 930–941. https://doi.org/10.1111/j.1365-2427.2009.02201.x (2009).Article 

Google Scholar 
Trontelj, P. et al. A molecular test for cryptic diversity in ground water: How large are the ranges of macro stygobionts?. Freshw. Biol. 54, 727–744. https://doi.org/10.1111/j.1365-2427.2007.01877.x (2009).Article 

Google Scholar 
Cooper, J. E. Ecological and Behavioral Studies in Shelta Cave, Alabama, with Emphasis on Decapod Crustaceans (University of Kentucky, 1975).
Google Scholar 
Voituron, Y., de Fraipont, M., Issartel, J., Guillaume, O. & Clobert, J. Extreme lifespan of the human fish (Proteus anguinus): A challenge for ageing mechanisms. Biol. Lett. 7, 105–107. https://doi.org/10.1098/rsbl.2010.0539 (2011).Article 
PubMed 

Google Scholar 
Poulson, T. L. Cave adaptation in amblyopsid fishes. Am. Midl. Nat. 70, 257–290. https://doi.org/10.2307/2423056 (1963).Article 

Google Scholar 
Venarsky, M. P., Huryn, A. D. & Benstead, J. P. Re-examining extreme longevity of the cave crayfish Orconectes australis using new mark–recapture data: A lesson on the limitations of iterative size-at-age models. Freshw. Biol. 57, 1471–1481. https://doi.org/10.1111/j.1365-2427.2012.02812.x (2012).Article 

Google Scholar 
Culver, D. C., Kane, T. C. & Fong, D. W. Adaptation and Natural Selection in Caves: The Evolution of Gammarus minus (Harvard University Press, 1995).Book 

Google Scholar 
Niemiller, M. L. & Poulson, T. L. Subterranean fishes of North America: Amblyopsidae. In Biology of Subterranean Fishes (eds Trajano, E. et al.) 169–280 (CRC Press, 2010).Chapter 

Google Scholar 
Fišer, C., Zagmajster, M. & Zakšek, V. Coevolution of life history traits and morphology in female subterranean amphipods. Oikos 122, 770–778. https://doi.org/10.1111/j.1600-0706.2012.20644.x (2013).Article 

Google Scholar 
Purvis, A., Gittleman, J. L., Cowlishaw, G. & Mace, G. M. Predicting extinction risk in declining species. Proc. Biol. Sci. 267, 1947–1952. https://doi.org/10.1098/rspb.2000.1234 (2000).Article 
PubMed 
PubMed Central 

Google Scholar 
Pearson, R. G. et al. Life history and spatial traits predict extinction risk due to climate change. Nat. Clim. Change 4, 217–221. https://doi.org/10.1038/nclimate2113 (2014).Article 
ADS 

Google Scholar 
Niemiller, M. L., Bichuette, E. & Taylor, S. J. Conservation of cave fauna in Europe and the Americas. In Ecological Studies: Cave Ecology (eds Moldovan, O. T. et al.) 451–478 (Springer, 2018).Chapter 

Google Scholar 
Niemiller, M. L. & Taylor, S. J. Protecting cave life. In Encyclopedia of Caves (eds Culver, D. C. et al.) 822–829 (Academic Press, 2019).Chapter 

Google Scholar 
Niemiller, M. L., Taylor, S. J., Slay, M. E. & Hobbs, H. H. III. Biodiversity in the United States and Canada. In Encyclopedia of Caves (eds Culver, D. C. et al.) 163–176 (Academic Press, 2019).Chapter 

Google Scholar 
Hortal, J. et al. Seven shortfalls that beset large-scale knowledge of biodiversity. Annu. Rev. Ecol. Evol. Syst. 46, 523–529. https://doi.org/10.1146/annurev-ecolsys-112414-054400 (2015).Article 

Google Scholar 
MacKenzie, D. I. et al. Occupancy Estimation and Modeling: Inferring Patterns and Dynamics of Species Occurrence (Academic Press, 2018).MATH 

Google Scholar 
Roberto, P. & Pietro, B. Species rediscovery or lucky endemic? Looking for the supposed missing species Leistus punctatissimus through a biogeographer’s eye (Coleoptera, Carabidae). ZooKeys 740, 97–108. https://doi.org/10.3897/zookeys.740.23495 (2018).Article 

Google Scholar 
Chu, C., Mandrak, N. E. & Minns, C. K. Potential impacts of climate change on the distributions of several common and rare freshwater fishes in Canada. Divers. Distrib. 11, 299–310. https://doi.org/10.1111/j.1366-9516.2005.00153.x (2005).Article 

Google Scholar 
Larson, E. R. & Olden, J. D. Latent extinction and invasion risk of crayfishes in the southeastern United States. Conserv. Biol. 24, 1099–1110. https://doi.org/10.1111/j.1523-1739.2010.01462.x (2010).Article 
PubMed 

Google Scholar 
Filipe, A. F. et al. Selection of priority areas for fish conservation in Guadiana River Basin, Iberian Peninsula. Conserv. Biol. 18, 189–200. https://doi.org/10.1111/j.1523-1739.2004.00620.x (2004).Article 

Google Scholar 
Mammola, S. et al. Fundamental research questions in subterranean biology. Biol. Rev. 95, 1855–1872. https://doi.org/10.1111/brv.12642 (2020).Article 
PubMed 

Google Scholar 
Domínguez-Domínguez, O., Martínez-Meyer, E., Zombrano, L. & de León, G. P. Using ecological-niche modeling as a conservation tool for freshwater species: Live-bearing fishes in central Mexico. Conserv. Biol. 20, 1730–1739. https://doi.org/10.1111/j.1523-1739.2006.00588.x (2006).Article 
PubMed 

Google Scholar 
Mammola, S. & Leroy, B. Applying species distribution models to caves and other subterranean habitats. Ecography 41, 1194–1208. https://doi.org/10.1111/ecog.03464 (2018).Article 

Google Scholar 
Castellarini, F., Malard, F., Dole-Olivier, M. & Gibert, J. Modelling the distribution of stygobionts in the Jura Mountains (eastern France). Implications for the protection of ground waters. Divers. Distrib. 13, 213–224. https://doi.org/10.1111/j.1472-4642.2006.00317.x (2007).Article 

Google Scholar 
Foulquier, A., Malard, F., Lefébure, T., Douady, C. J. & Gibert, J. The imprint of Quaternary glaciers on the present-day distribution of the obligate groundwater amphipod Niphargus virei (Niphargidae). J. Biogeogr. 35, 552–564. https://doi.org/10.1111/j.1365-2699.2007.01795.x (2008).Article 

Google Scholar 
Johns, T. et al. Regional-scale drivers of groundwater faunal distributions. Freshw. Sci. 34, 316–328. https://doi.org/10.1086/678460 (2015).Article 

Google Scholar 
Camp, C. D., Wooten, J. A., Jensen, J. B. & Bartek, D. F. Role of temperature in determining relative abundance in cave twilight zones by two species of lungless salamander (family Plethodontidae). Can. J. Zool. 92, 119–127. https://doi.org/10.1139/cjz-2013-0178 (2014).Article 

Google Scholar 
Korbel, K. L., Hancock, P. J., Serov, P., Lim, R. P. & Hose, G. C. Groundwater ecosystems vary with land use across a mixed agricultural landscape. J. Environ. Qual. 42, 380–390. https://doi.org/10.2134/jeq2012.0018 (2013).Article 
PubMed 

Google Scholar 
Español, C. et al. Does land use impact on groundwater invertebrate diversity and functionality in floodplains?. Ecol. Eng. 103, 394–403. https://doi.org/10.1016/j.ecoleng.2016.11.061 (2017).Article 

Google Scholar 
Christman, M. C. et al. Predicting the occurrence of cave-inhabiting fauna based on features of the earth surface environment. PLoS One 11, e0160408. https://doi.org/10.1371/journal.pone.0160408 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Zagmajster, M. et al. Geographic variation in range size and beta diversity of groundwater crustaceans: Insights from habitats with low thermal seasonality. Glob. Ecol. Biogeogr. 23, 1135–1145. https://doi.org/10.1111/geb.12200 (2014).Article 

Google Scholar 
Poff, N. L. Landscape filters and species traits: Towards mechanistic understanding and prediction in stream ecology. J. North Am. Benthol. Soc. 16, 391–409. https://doi.org/10.2307/1468026 (1997).Article 

Google Scholar 
Stevenson, R. J. Scale-dependent determinants and consequences of benthic algal heterogeneity. J. North Am. Benthol. Soc. 16, 248–262. https://doi.org/10.2307/1468255 (1997).Article 
ADS 

Google Scholar 
U.S. Geological Survey. NLCD 2011 land cover. Multi-Resolution Land Characteristics. https://www.mrlc.gov/data/nlcd-2011-land-cover-conus (2011).Adamski, J. C. Geochemistry of the Springfield Plateau Aquifer of the Ozark Plateaus Province in Arkansas, Kansas, Missouri and Oklahoma, USA. Hydrol. Process. 14, 849–866. https://doi.org/10.1002/(SICI)1099-1085(20000415)14:5%3c849::AID-HYP973%3e3.0.CO;2-7 (2000).Article 
ADS 

Google Scholar 
Woods, A. J. et al. Ecoregions of Oklahoma (Color Poster with Map, Descriptive Text, Summary Tables, and Photographs) (U.S. Geological Survey, 2005).
Google Scholar 
Unklesbay, A. G. & Vineyard, J. D. Missouri Geology: Three Billion Years of Volcanoes, Seas, Sediments, and Erosion (University of Missouri Press, 1992).
Google Scholar 
Eigenmann, C. H. A new blind fish. In Proceedings of the Indiana Academy of Science 1897 (ed Waldo, C. A.) 231 (1898).Graening, G. O., Fenolio, D. B., Niemiller, M. L., Brown, A. V. & Beard, J. B. The 30-year recovery effort for the Ozark cavefish (Amblyopsis rosae): Analysis of current distribution, population trends, and conservation status of this threatened species. Environ. Biol. Fish. 87, 55–88. https://doi.org/10.1007/s10641-009-9568-2 (2010).Article 

Google Scholar 
Niemiller, M. L., Near, T. J. & Fitzpatrick, B. M. Delimiting species using multilocus data: Diagnosing cryptic diversity in the southern cavefish, Typhlichthys subterraneus (Teleostei: Amblyopsidae). Evolution 66, 846–866. https://doi.org/10.1111/j.1558-5646.2011.01480.x (2012).Article 
PubMed 

Google Scholar 
Hobbs, H. H. Jr. & Brown, A. V. A new troglobitic crayfish from northwestern Arkansas (Decapoda: Cambaridae). Proc. Biol. Soc. Wash. 100, 1040–1048 (1987).
Google Scholar 
Graening, G. O., Slay, M. E., Brown, A. V. & Koppelman, J. B. Status and distribution of the endangered Benton cave crayfish, Cambarus aculabrum (Decapoda: Cambaridae). Southwest. Nat. 51, 376–381. https://doi.org/10.1894/0038-4909(2006)51[376:SADOTE]2.0.CO;2 (2006).Article 

Google Scholar 
Faxon, W. Cave animals from southwestern Missouri. Bull. Mus. Comp. Zool. 17, 225–240 (1889).
Google Scholar 
Graening, G. O., Hobbs, H. H. III., Slay, M. E., Elliott, W. R. & Brown, A. V. Status update for bristly cave crayfish, Cambarus setosus (Decapoda: Cambaridae), and range extension into Arkansas. Southwest. Nat. 51, 382–392. https://doi.org/10.1894/0038-4909(2006)51[382:SUFBCC]2.0.CO;2 (2006).Article 

Google Scholar 
Hobbs, H. H. III. Cambarus (Jugicambarus) subterraneus, a new cave crayfish (Decapoda: Cambaridae) from northeastern Oklahoma, with a key to the troglobitic members of the subgenus Jugicambarus. Proc. Biol. Soc. Wash. 106, 719–727 (1993).
Google Scholar 
Graening, G. O. & Fenolio, D. B. Status update of the Delaware County cave crayfish, Cambarus subterraneus (Decapoda: Cambaridae). Proc. Okla. Acad. Sci. 85, 85–89 (2005).
Google Scholar 
Hobbs, H. H. Jr. & Cooper, M. R. A new troglobitic crayfish from Oklahoma (Decapoda: Astacidae). Proc. Biol. Soc. Wash. 85, 49–56 (1972).
Google Scholar 
Graening, G. O. et al. Range extension and status update for the Oklahoma cave crayfish, Cambarus tartarus (Decapoda: Cambaridae). Southwest. Nat. 51, 94–99 (2006).Article 

Google Scholar 
Hobbs, H. H. III. A new cave crayfish of the genus Orconectes, subgenus Orconectes, from southcentral Missouri, USA, with a key to the stygobitic species of the genus (Decapoda, Cambaridae). Crustaceana 74, 635–646. https://doi.org/10.1163/156854001750377911 (2001).Article 

Google Scholar 
Miller, B. V. The Hydrology of the Carroll Cave-Toronto Springs System: Identifying and Examining Source Mixing Through Dye Tracing, Geochemical Monitoring, Seepage Runs, and Statistical Methods (Western Kentucky University, 2010).
Google Scholar 
Mouser, J. B., Brewer, S. K., Niemiller, M. L., Mollenhauer, R. & Van Den Bussche, R. Comparing visual and environmental DNA surveys for detection of stygobionts. Subterr. Biol. 39, 79–105. https://doi.org/10.3897/subtbiol.39.64279 (2021).Article 

Google Scholar 
Longmire, J. L., Maltbie, M. & Baker, R. J. Use of “Lysis Buffer” in DNA Isolation and Its Implication for Museum Collections (Museum of Texas Tech University, 1997).Book 

Google Scholar 
Mouser, J. B., Mollenhauer, R. & Brewer, S. K. Relationships between landscape constraints and a crayfish assemblage with consideration of competitor presence. Divers. Distrib. 25, 61–73. https://doi.org/10.1111/ddi.12840 (2019).Article 

Google Scholar 
U.S. Geological Survey. 1 Arc-second digital elevation models (DEMs)—USGS national map 3DEP downloadable data collection. https://data.usgs.gov/datacatalog/data/USGS:35f9c4d4-b113-4c8d-8691-47c428c29a5b (U.S. Geological Survey, 2017).Oak Ridge National Laboratory Distributed Active Archive Center. MODIS and VIIRS land products global subsetting and visualization tool. Oak Ridge National Laboratory Distributed Active Archive Center. https://doi.org/10.3334/ORNLDAAC/1379 (2018).Horton, J. D., San Juan, C. A. & Stoeser, D. B. The state geologic map compilation (SGMC) geodatabase of the conterminous United States. U.S. Geol. Surv. https://doi.org/10.3133/ds1052 (2017).Article 

Google Scholar 
MacKenzie, D. I. et al. Estimating site occupancy rates when detection probabilities are less than one. Ecology 83, 2248–2255. https://doi.org/10.1890/0012-9658(2002)083[2248:ESORWD]2.0.CO;2 (2002).Article 

Google Scholar 
Tyre, A. J. et al. Improving precision and reducing bias in biological surveys: Estimating false negative error rates. Ecol. Appl. 13, 1790–1801. https://doi.org/10.1890/02-5078 (2003).Article 

Google Scholar 
Gelman, A. & Hill, J. Data Analysis Using Regression and Multilevel Hierarchical Models (Cambridge University Press, 2007).
Google Scholar 
Kéry, M. & Royle, J. A. Applied Hierarchical Modeling in Ecology: Analysis of Distribution, Abundance and Species Richness in R and BUGS (Academic Press, 2016).MATH 

Google Scholar 
Plummer, M. JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. In Proceedings of the 3rd International Workshop on Distributed Statistical Computing (eds Hornik, K. et al.) 1–10 (Austrian Science Foundation, 2003).
Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2020).
Google Scholar 
Kellner, J. jagsUI: A wrapper around ‘rjags’ to streamline ‘JAGS’ analyses. https://CRAN.R-project.org/package=jagsUI (R-project, 2019).Brooks, S. P. & Gelman, A. General methods for monitoring convergence of iterative simulations. J. Comput. Graph. Stat. 7, 434–455. https://doi.org/10.1080/10618600.1998.10474787 (1998).Article 
MathSciNet 

Google Scholar 
Kruschke, J. K. Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan (Academic Press, 2015).MATH 

Google Scholar 
Hobbs, N. T. & Hooten, M. B. Bayesian Models (Princeton University Press, 2015). https://doi.org/10.1515/9781400866557.Book 

Google Scholar 
Conn, P. B., Johnson, D. S., Williams, P. J., Melin, S. R. & Hooten, M. B. A guide to Bayesian model checking for ecologists. Ecol. Monogr. 88, 526–542. https://doi.org/10.1002/ecm.1314 (2018).Article 

Google Scholar 
Allan, J. D. Landscapes and riverscapes: The influence of land use on stream ecosystems. Annu. Rev. Ecol. Evol. Syst. 35, 257–284. https://doi.org/10.1146/annurev.ecolsys.35.120202.110122 (2004).Article 

Google Scholar 
Paul, M. J. & Meyer, J. L. Streams in the urban landscape. Annu. Rev. Ecol. Evol. Syst. 32, 333–365. https://doi.org/10.1146/annurev.ecolsys.32.081501.114040 (2001).Article 

Google Scholar 
Wicks, C., Kelley, C. & Peterson, E. Estrogen in a karstic aquifer. Groundwater 42, 384–389. https://doi.org/10.1111/j.1745-6584.2004.tb02686.x (2004).Article 

Google Scholar 
Buřič, M., Kouba, A., Máchová, J., Mahovská, I. & Kozák, P. Toxicity of the organophosphate pesticide diazinon to crayfish of differing age. Int. J. Environ. Sci. Technol. 10, 607–610. https://doi.org/10.1007/s13762-013-0185-4 (2013).Article 

Google Scholar 
Sohn, L., Brodie, R. J., Couldwell, G., Demmons, E. & Sturve, J. Exposure to a nicotinoid pesticide reduces defensive behaviors in a non-target organism, the rusty crayfish Orconectes rusticus. Ecotoxicology 27, 900–907. https://doi.org/10.1007/s10646-018-1950-4 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Noltie, D. B. & Wicks, C. M. How hydrogeology has shaped the ecology of Missouri’s Ozark cavefish, Amblyopsis rosae, and southern cavefish, Typhlichthys subterraneus: Insights on the sightless from understanding the underground. Environ. Biol. Fish. 62, 171–194. https://doi.org/10.1023/A:1011815806589 (2001).Article 

Google Scholar 
Kuhajda, B. R. & Mayden, R. L. Status of the federally endangered Alabama cavefish, Speoplatyrhinus poulsoni (Amblyopsidae), in Key Cave and surrounding caves, Alabama. Environ. Biol. Fish. 62, 215–222. https://doi.org/10.1023/A:1011817023749 (2001).Article 

Google Scholar 
Hutchins, B. T. The conservation status of Texas groundwater invertebrates. Biodivers. Conserv. 27, 475–501. https://doi.org/10.1007/s10531-017-1447-0 (2018).Article 

Google Scholar 
Niemiller, M. L. et al. Discovery of a new population of the federally endangered Alabama cave shrimp, Palaemonias alabamae Smalley, 1961, in northern Alabama. Subterr. Biol. 32, 43–59. https://doi.org/10.3897/subtbiol.32.38280 (2019).Article 

Google Scholar 
Abell, R., Allan, J. D. & Lehner, B. Unlocking the potential of protected areas for freshwaters. Biol. Conserv. 134, 48–63. https://doi.org/10.1016/j.biocon.2006.08.017 (2007).Article 

Google Scholar 
Liu, Y. et al. A review on effectiveness of best management practices in improving hydrology and water quality: Needs and opportunities. Sci. Total Environ. 601–602, 580–593. https://doi.org/10.1016/j.scitotenv.2017.05.212 (2017).Article 
ADS 
PubMed 

Google Scholar  More

Automated tracking of four social insect speciesFifty queenright colonies were used in the tracking experiments (Table 1). Honeybee colonies (subspecies A. mellifera carnica) were housed in the campus apiary of the University of Lausanne. Colonies of L. niger were raised from single mated queens collected on campus. T. nylanderi colonies were collected from the University of Lausanne campus, and L. acervorum colonies collected from Anzeindaz, Switzerland. These four species were chosen because of their abundance and easy availability in Switzerland, and because they – or closely-related species – have previously been used as model systems for the study of spatial organisation in social insects20,21,23,27,44. The colony sizes used in our experiments (Table 1) fell within the natural range of sizes experienced by these species in nature, either as recently founded colonies (L. niger colonies are founded by a single queen and progressively grow from a few workers to mature sizes of up to 40,000 workers over the course of several years; new honeybee colonies are founded by swarms counting 2400–41,000 bees61) or as mature colonies (all colonies of T. nylanderi and L. acervorum used in our experiments were mature colonies collected whole from the field).In all species, a paper tag bearing a unique two-dimensional barcode was glued to the thorax of individuals to allow automated tracking of their movements (Fig. S1). In the ants, tagging of all individuals was performed in a single session two days before the beginning of the experiment, whilst in the bees, newly-emerged workers (one-day-old or less) were tagged every 3 days over the 21 days prior to the beginning of the experiment (Supplementary Note 1).Tagged colonies were kept in glass observation nests with a single entrance (internal nest dimensions, A. mellifera: 69 × 45 × 4 cm, L. niger: 70 × 40 × 8 mm; L. acervorum: 63 × 42 × 2 mm, T. nylanderi: 63 × 42 × 1.5 mm). The honeybee observation nests also included a 64 × 44 cm wooden frame enclosing a double-sided wax comb containing honey, pollen, and developing brood20. Bees were free to move between both sides of the comb. In all species, individuals were allowed to freely exit and enter the nest. Ants were provided with ad libitum food (Drosophila, sugar solution) and water in a foraging arena, while bees foraged on natural resources outside. Both the ant and honeybee observation nests were exposed to diurnal cycles of temperature and light (Supplementary Note 1).High resolution digital video cameras operating at two frames per second were used to identify the location and orientation of each tag across successive images22. All colonies were continuously tracked for three days, which corresponded to the inter-cohort time in the honeybee colonies. The trajectories of each worker, and the physical contacts between workers (Fig. S18 and Supplementary Note 14) were extracted using an existing software pipeline62.Building bipartite site-visit networksTo quantify the spatial preferences of individual ants and bees, the interior of the nest was discretised into a regular hexagonal lattice (Fig. 1a, b). Because the worker body lengths of our four study species span an order of magnitude (from ~ 1.5 mm for T. nylanderi to ~ 15 mm for A. mellifera), the width of the hexagonal bins were defined as 1/4 of the mean worker body-length.To characterise the spatial preferences of different individuals to different parts of the nest, we counted the number of times ({n}_{i}^{s}) that each individual i visited each hexagonal site s. A visit by individual i to site s began when i crossed the border into s, and was terminated when i crossed the border out of s, regardless of the amount of time spent inside. To prevent stationary individuals located on the border between two adjacent sites from rapidly accumulating many single-frame visits to the two sites, successive visits to a same site were only counted when at least 20s elapsed between the end of the previous visit and the start of the next.The site-visit data were used to construct a bipartite network, in which individuals (layer 1) were connected by undirected edges to the sites (layer 2) they visited (Fig. 1c, d). Because individuals typically made multiple visits to the same sites, each edge i–s was weighted according to the total number of times individual i visited site s, that is, ({n}_{i}^{s}).Partitioning site-visit networks into modulesThe extent to which the site-visit networks were partitioned into discrete ‘modules’ (i.e., set of workers with similar space-use patterns and the set of sites that they exhibit strong ties to) was assessed using the DIRTLPAwb+ algorithm for partitioning weighted bipartite networks39. This algorithm searches for the partition that maximises the number and strength of the links within modules, whilst minimizing connections between modules. The number of modules was not specified a priori by the user, but was identified by the algorithm. All site-visit networks had positive modularity (Fig. S3), indicating that they could be partitioned into a set of well-separated modules (Figs. 1e–h, S2, and S4–S5). The modules in each partition were then assigned functional labels according to the following rules. First, the module whose sites were on average closest to the nest entrance was labelled ‘forager’ module. Second, the module or modules with the greatest spatial overlap with the brood pile in the ant colonies or the broodnest(s) in the honeybee colonies were labelled ‘nurse’ module(s). After defining the forager and nurse modules, the remaining modules (if any) were labelled as follows. If there was only one module remaining after identifying the nurse and forager modules, as was typically the case in honeybee colonies, it was labelled ‘peripheral’. If there were two modules remaining, as was typically the case in ant colonies, then the module whose sites were on average closer to the nest borders (i.e., to the periphery of the nest) was labelled ‘peripheral’, and the remaining module labelled ‘intermediate’. In some cases, the DIRTLPAwb+ algorithm identified five or more modules (9.0% of all iterations across all species and colonies). In these cases, the supernumerary modules never contained more than 1 or 2 individuals, and as they could not be unequivocally assigned using our labelling scheme, they were left unclassified for these iterations.Validating network modulesAs a network constructed by a purely random process could exhibit apparent modular structure by chance, we tested whether the discovered modules represent statistically significant entities. To do so, we produced 1000 null model random networks for each observed network using an established permutation method for bipartite networks63 (Supplementary Note 2). Comparisons between the maximum modularity of the observed networks with that of the corresponding random networks showed that, in all four species, the observed modularity was significantly greater than expected by chance (Fig. S3).Constructing worker task profilesA unique labour profile for each ant and each honeybee was constructed by estimating the activity of each worker in the following four tasks:1. Entrance guarding: workers were classed as guarding when they were (i) within two body lengths of the entrance, (ii) roughly facing the entrance, i.e., with a body alignment diverging from the direct heading to the entrance by no more than π/2 radians, and (iii) ‘on station’ at the entrance, as defined by trajectory coordinates with an associated first passage time (ref. 64; time taken for the individual to pass beyond a circle centred on its current location with a radius of two body-lengths) in excess of 500s.2. Patrolling: workers were classed as patrolling65 when they were (i) active, and (ii) ‘roaming’, as defined by first passage times of More

Podos, J. Correlated evolution of morphology and vocal signal structure in Darwin’s finches. Nature 409, 185–188 (2001).ADS 
CAS 
PubMed 

Google Scholar 
Bradbury, J. W. & Vehrencamp, S. Principles of Animal Communication (Sinauer Associated, 1998).
Google Scholar 
Podos, J. & Warren, P. S. The evolution of geographic variation in birdsong. Adv. Study Behav. 37, 403–458 (2007).
Google Scholar 
Charlton, B. D., Owen, M. A. & Swaisgood, R. R. Coevolution of vocal signal characteristics and hearing sensitivity in forest mammals. Nat. Commun. 10, 1–7 (2019).CAS 

Google Scholar 
Soltis, J. Vocal communication in African elephants (Loxodonta africana). Zoo Biol. 29, 192–209 (2010).PubMed 

Google Scholar 
King, S. L. & Janik, V. M. Bottlenose dolphins can use learned vocal labels to address each other. Proc. Natl. Acad. Sci. 110, 13216–13221 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Ravignani, A. et al. Ontogeny of vocal rhythms in harbor seal pups: An exploratory study. Curr. Zool. 65, 107–120 (2019).PubMed 

Google Scholar 
Rauber, R. & Manser, M. B. Effect of group size and experience on the ontogeny of sentinel calling behaviour in meerkats. Anim. Behav. 171, 129–138 (2021).
Google Scholar 
Janik, V. M. & Slater, P. J. Vocal learning in mammals. Adv. Study Behav. 26, 59–100 (1997).
Google Scholar 
Fitch, W. T. Production of vocalizations in mammals. In Encyclopedia of Language and Linguistics 2nd edn, Vol. 1 115–121 (Elsevier, 2006).
Google Scholar 
Fletcher, N. H. A frequency scaling rule in mammalian vocalization. In Handbook of Behavioral Neuroscience Vol. 19 51–56 (Elsevier, 2010).
Google Scholar 
Fant, G. Acoustic Theory of Speech Production (Mouton, 1960).
Google Scholar 
Taylor, A. & Reby, D. The contribution of source–filter theory to mammal vocal communication research. J. Zool. 280, 221–236 (2010).
Google Scholar 
Fitch, W. T., Neubauer, J. & Herzel, H. Calls out of chaos: The adaptive significance of nonlinear phenomena in mammalian vocal production. Anim. Behav. 63, 407–418 (2002).
Google Scholar 
Domning, D. P. & Hayek, L. A. C. Interspecific and intraspecific morphological variation in manatees (Sirenia: Trichechus). Mar. Mamm. Sci. 2, 87–144 (1986).
Google Scholar 
Anderson, P. K. & Barclay, R. M. R. Acoustic signals of solitary Dugongs: Physical characteristics and behavioral correlates. J. Mamm. 76, 1226–1237 (1995).
Google Scholar 
Sousa-Lima, R., Paglia, A. P. & Da Fonseca, G. Signature information and individual recognition in the isolation calls of Amazonian manatees, Trichechus inunguis (Mammalia: Sirenia). Anim. Behav. 63, 301–310 (2002).
Google Scholar 
Sousa-Lima, R. S., Paglia, A. P. & Fonseca, G. A. B. Gender, age, and identity in the isolation calls of Antillean manatees (Trichechus manatus manatus). Aquat. Mamm. 34, 109–122 (2008).
Google Scholar 
O’Shea, T. J. & Poché, L. B. Aspects of underwater sound communication in Florida manatees (Trichechus manatus latirostris). J. Mamm. 87, 1061–1071 (2006).
Google Scholar 
Rosas, F. C. W. Biology, conservation and status of the Amazonian manatee Trichechus inunguis. Mamm. Rev. 24, 49–59 (1994).
Google Scholar 
Meirelles, A. C. O. & Carvalho, V. L. Peixe-boi marinho: biologia e conservação no Brasil. Aquasis, Bambu Editora e Artes Gráficas, São Paulo (2016).Alvarez-Alemán, A., Beck, C. A. & Powell, J. A. First report of a Florida manatee (Trichechus manatus latirostris) in Cuba. Aquat. Mamm. 36, 148 (2010).
Google Scholar 
Castelblanco-Martínez, D. N. et al. First documentation of long-distance travel by a Florida manatee to the Mexican Caribbean. Ethol. Ecol. Evol. 1–12 (2021).Packard, J. M. & Wetterqvist, O. F. Evaluation of manatee habitat systems on the northwestern Florida coast. Coast. Manag. 14, 279–310 (1986).
Google Scholar 
Luna, F. D. O. et al. Genetic connectivity of the West Indian manatee in the southern range and limited evidence of hybridization with Amazonian manatees. Front. Mar. Sci. 7, 1089 (2021).
Google Scholar 
Hartman, D. Ecology and behavior of the manatee (Trichechus manatus) in Florida. Spec. Publ. Am. Soc. Mammal. 5, 153 (1979).
Google Scholar 
D’AffonsecaNeto, J. A. & Vergara-Parente, J. E. Sirenia (peixe-boi-da-Amazônia, Peixe-boi-marinho). In Tratado de Animais Selvagens: medicina veterinária (eds Cubas, Z. S. et al.) 701–714 (Roca, 2006).
Google Scholar 
Deutsch, C. J., Reid, J. P., Bonde, R. K., Easton, D. E., Kochman, H. I. & O’Shea, T. J. Seasonal movements, migratory behavior, and site fidelity of West Indian manatees along the Atlantic coast of the United States. Wildl. Monogr. 1–77 (2003).Laist, D. W. & Reynolds, J. E. III. Influence of power plants and other warm-water refuges on Florida manatees. Mar. Mamm. Sci. 21, 739–764 (2005).
Google Scholar 
Reynolds, J. E. Aspects of the social behaviour and herd structure of a semi-isolated colony of West Indian manatees, Trichechus manatus. Mammalia 45, 431–452 (1981).
Google Scholar 
Dantas, G. A. Ontogenia do padrão vocal individual do peixe-boi da Amazônia Trichechus inunguis (Sirenia, trichechidae). Dissertação (Instituto Nacional de Pesquisas da Amazônia, 2009).
Google Scholar 
Brady, B., Moore, J. & Love, K. Behavior related vocalizations of the Florida manatee (Trichechus manatus latirostris). Mar. Mamm. Sci. 1–15 (2021).Nowacek, D. P., Casper, B. M., Wells, R. S., Nowacek, S. M. & Mann, D. A. Intraspecific and geographic variation of West Indian manatee (Trichechus manatus spp.) vocalizations. J. Acoust. Soc. Am. 114, 66–69 (2003).ADS 
PubMed 

Google Scholar 
Rycyk, A. M. et al. First characterization of vocalizations and passive acoustic monitoring of the vulnerable African manatee (Trichechus senegalensis). J. Acoust. Soc. Am. 150, 3028–3037 (2021).ADS 
PubMed 

Google Scholar 
Landrau-Giovannetti, N., Mignucci-Giannoni, A. A. & Reidenberg, J. S. Acoustical and anatomical determination of sound production and transmission in West Indian (Trichechus manatus) and Amazonian (T. inunguis) manatees. Anat. Rec. 297, 1896–1907 (2014).
Google Scholar 
Morton, E. On the occurrence and significance of motivation-structural rules in some bird and mammal sounds. Am. Nat. 111, 855–869 (1977).
Google Scholar 
Borges, J. C. et al. Growth pattern differences of captive born Antillean manatee (Trichechus manatus) calves and those rescued in the Brazilian northeastern coast. J. Zoo Wildl. Med. 43, 494–500 (2012).PubMed 

Google Scholar 
Lima, D. S., Vergara-Parente, J. E., Young, R. J. & Paszkiewicz, E. Training of Antillean manatee Trichechus manatus manatus Linnaeus, 1758 as a management technique for individual welfare. Lat. Am. J. Mar. Mamm. 4, 61–68 (2005).
Google Scholar 
K. Lisa Yang Center for Conservation Bioacoustics at the Cornell Lab of Ornithology. Raven Pro: Interactive Sound Analysis Software (Version 1.5) [Computer software]. https://ravensoundsoftware.com/ (The Cornell Lab of Ornithology, 2022).Zollinger, S. A., Podos, J., Nemeth, E., Goller, F. & Brumm, H. On the relationship between, and measurement of, amplitude and frequency in birdsong. Anim. Behav. 84, e1–e9 (2012).
Google Scholar 
Sokal, R. R. & Rohlf, F. J. Biometry (W. H. Freeman and Co., 1995).MATH 

Google Scholar 
Charrier, I. & Harcourt, R. G. Individual vocal identity in mother and pup Australian sea lions (Neophoca cinerea). J. Mamm. 87, 929–938 (2006).
Google Scholar 
Green, S. & Salkind, N. J. Using SPSS for Windows and Macintosh: Analyzing and Understanding Data 4th edn. (Prentice Hall, 2003).
Google Scholar 
Charlton, B. D. et al. Cues to body size in the formant spacing of male koala (Phascolarctos cinereus) bellows: Honesty in an exaggerated trait. J. Exp. Biol. 214(20), 3414–3422 (2011).PubMed 

Google Scholar 
IBM Corp. Released 2020. IBM SPSS Statistics for Windows, Version 27.0.Best, R. C. The aquatic mammals and reptiles of the Amazon. In The Amazon 371–412 (Springer, 1984).
Google Scholar 
Gerhardt, H. C. The evolution of vocalization in frogs and toads. Annu. Rev. Ecol. Syst. 25, 293–324 (1994).
Google Scholar 
Mendoza, P. et al. Growth curve of Amazonian manatee (Trichechus inunguis) in captivity. Aquat. Mamm. 45 (2019).Schwarz, L. K. Methods and models to determine perinatal status of Florida manatee carcasses. Mar. Mamm. Sci. 24, 881–898 (2008).
Google Scholar 
Hauser, M. D., Chomsky, N. & Fitch, W. T. The faculty of language: What is it, who has it, and how did it evolve?. Science 298, 1569–1579 (2002).ADS 
CAS 
PubMed 

Google Scholar 
Marler, P. & Peters, S. Developmental overproduction and selective attrition: New processes in the epigenesis of birdsong. Dev. Psychol. J. Int. Soc. Dev. Psychol. 15.4, 369–378 (1982).
Google Scholar 
Casey, C., Reichmuth, C., Costa, D. P. & Le Boeuf, B. The rise and fall of dialects in northern elephant seals. Proc. R. Soc. B 285, 2018–2176 (2018).
Google Scholar 
Hunter, M. E. et al. Puerto Rico and Florida manatees represent genetically distinct groups. Conserv. Genet. 13, 1623–1635 (2012).
Google Scholar 
Castelblanco-Martínez, D. N. et al. Analysis of body condition indices reveals different ecotypes of the Antillean manatee. Sci. Rep. 11, 1–14 (2021).
Google Scholar 
McCracken, K. G. & Sheldon, F. H. Avian vocalizations and phylogenetic signal. Proc. Natl. Acad. Sci. 94, 3833–3836 (1997).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Moron, J. R. et al. Whistle variability of Guiana dolphins in South America: Latitudinal variation or acoustic adaptation? Mar. Mamm. Sci. 1–32 (2018)Luís, A. R. et al. Vocal universals and geographic variations in the acoustic repertoire of the common bottlenose dolphin. Sci. Rep. 11, 1–9 (2021).
Google Scholar 
Ey, E. & Fisher, J. The ‘Acoustic adaptations hypothesis’ a review of the evidence from birds, anurans and mammals. Bioacoustics 19, 21–48 (2009).
Google Scholar 
Miksis-Olds, J. L. & Tyack, P. L. Manatee (Trichechus manatus) vocalization usage in relation to environmental noise levels. J. Acoust. Soc. Am. 125, 1806–1815 (2009).ADS 
PubMed 

Google Scholar 
Sun, W., Wang, Z., Jamalabdollahi, M. & Reza Zekavat, S. A. Experimental study on the difference between acoustic communication channels in freshwater rivers/lakes and in oceans. In 2014 48th Asilomar Conference on Signals, Systems and Computers, 333–337 (2004)Rivera Chavarría, M., Castro, J. & Camacho, A. The relationship between acoustic habitat, hearing and tonal vocalizations in the Antillean manatee (Trichechus manatus manatus, Linnaeus, 1758). Biol. Open 4, 1237–1242 (2015).PubMed 
PubMed Central 

Google Scholar 
Gaspard, J. C. III. et al. Audiogram and auditory critical ratios of two Florida manatees (Trichechus manatus latirostris). J. Exp. Biol. 215, 1442–1447 (2012).PubMed 

Google Scholar 
Gerstein, E. R., Gerstein, L., Forsythe, S. E. & Blue, J. E. The underwater audiogram of the West Indian manatee (Trichechus manatus). J. Acoust. Soc. Am. 105, 3575–3583 (1999).ADS 
CAS 
PubMed 

Google Scholar 
Klishin, V., Pezo, R., Popov, V., Ya, A. & Supin. Some characteristics of hearing of the Brazilian manatee, Trichechus inunguis. Aquat. Mamm. 16 (1990).Johnson, M., de Soto, N. A. & Madsen, P. T. Studying the behaviour and sensory ecology of marine mammals using acoustic recording tags: A review. Mar. Ecol. Prog. Ser. 395, 55–73 (2009).ADS 

Google Scholar  More

All methods were performed in accordance with the local relevant guidelines, regulations and legislation.InstrumentsLI-6400 photosynthesis system (LI-6400 Inc., Lincoln, NE, USA) and PAM-2500 portable chlorophyll fluorescence apparatus (PAM-2500, Walz, Germany) were used in the study.MaterialsAbout 90 living E. brevicornu plants were collected from Taihang Mountains in October 2018. The E. brevicornu was not in endangered or protected. The collection of these E. brevicornu plants was permitted by local government. These plants were averagely planted in nine plots of 2 m2. The roots of E. pubescens were planted 6–8 cm below ground. These plots were placed on farmland near Taihang Mountains and covered with sunshade net (about 70% light transmittance). These plants were timely irrigated after planting to ensure that they grew well but not fertilized.Determination of photosynthetic characteristicsThe photosynthetic characteristics of mature leaves on the E. brevicornu plants were determined between June 6–8, 2019 with the Li-6400 photosynthesis system. The diurnal variation of photosynthesis in three leaves of three plants was determined. When the light response curve was determined, the temperature of the leaf chamber was set at 28 °C, and the concentration of CO2 in the leaf chamber was set at 400 µmol mol−1. When determining the CO2 response curve, the light intensity in the leaf chamber was set at 1000 µmol m−2 s−1, and the temperature of the leaf chamber was set at 28 °C. The light response curve and CO2 response curve were determined three times in three leaves of three different plants.Determination of chlorophyll fluorescence characteristicsThe fluorescence characteristics of chlorophyll in E. brevicornu leaves were determined with PAM-2500 portable chlorophyll fluorescence apparatus between June 8–9, 2019. The leaves underwent dark adaptation for 30 min before determining slow kinetics of chlorophyll fluorescence. Then the light curves of chlorophyll fluorescence were determined. All of these determinations were repeated three times on three mature leaves of three plants.The data was analysed with SPSS (Statistical Product and Service Solutions, International Business Machines Corporation, USA). The light response curves were fitted with following modified rectangular hyperbola model11,12.$${text{Photo}}, = ,{text{E}}cdotleft( {{1} – {text{M}}cdot{text{PAR}}} right)cdotleft( {{text{PAR}} – {text{LCP}}} right)/({1}, + ,{text{N}}cdot{text{PAR}})$$PAR is the value of light intensity in leaf chamber of Li-6400 photosynthesis system. Photo is net photosynthetic rate. LCP is the light compensation point. E is the apparent quantum yield. M and N are parameters. The dark respiration rate under the LCP is calculated according to E·LCP. The light saturation point (LSP) is calculated according to (((M + N) ·(1 + N·LCP)/M)½)/−1)/N.The net photosynthetic rate under the light saturation point (LSP) can be calculated according to the above model.The CO2 response curves were fitted with below modified rectangular hyperbola model11,12.$${text{Photo}}, = ,{text{E}}cdotleft( {{1} – {text{M}}cdot{text{PAR}}} right)cdotleft( {{text{PAR}} – {text{CCP}}} right)/({1}, + ,{text{N}}cdot{text{PAR}})$$PAR is the value of light intensity in leaf chamber of Li-6400 photosynthesis system. Photo is net photosynthetic rate. CCP is CO2 compensation point. E is also the apparent quantum yield. M and N are parameters. The dark respiration rate under the CO2 calculated according to E·CCP. The CO2 saturation point (CSP) is calculated according to (((M + N) ·(1 + N·CCP)/M)½)/−1)/N.The net photosynthetic rate under the CO2 saturation point (CSP) can be alternatively calculated according to the above model.The light curves of chlorophyll fluorescence were fitted according to the below model of Eilers and Peeters12,13.$${text{ETR}}, = ,{text{PAR}}/({text{a}}cdot{text{PAR}}^{{2}} , + ,{text{b}}cdot{text{PAR}}, + ,{text{c}})$$ETR is the electron transport rate of photosynthetic system II. PAR is fluorescence intensity. The letters a, b and c are parameters. More

European Commission Directorate General for Research and Innovation. A sustainable Bioeconomy for Europe: Strengthening the Connection Between Economy, Society and the Environment: Updated Bioeconomy Strategy (Directorate General for Research and Innovation, 2018).
Google Scholar 
Teitelbaum, L., Boldt, C. & Patermann, C. Global Bioeconomy Policy Report (IV): A Decade of Bioeconomy policy (International Advisory Council on Global Bioeconomy, 2020).
Google Scholar 
European Parliament; European Council. Directive (EU) 2018/2001 of the European Parliament and of the Council of 11 December 2018 on the promotion of the use of energy from renewable sources (2018). (Online). http://data.europa.eu/eli/dir/2018/2001/oj.European Parliament; European Council. Directive 2009/28/EC on the Promotion of the Use of Energy from Renewable Sources (2009). (Online). http://data.europa.eu/eli/dir/2009/28/oj.Glasenapp, S., & McCusker, A. Wood energy data: the joint wood, in Wood Energy in the ECE Region: Data, Trends and Outlook in Europe, the Commonwealth of Independent States and North America, Geneva, United Nations’ Economic Commission for Europe: ECE/TIM/SP/42, 12–29 (2018).Eurostat. Wood Products—Production and Trade (2021). (Online). https://ec.europa.eu/eurostat/statistics-explained/index.php?title=Wood_products_-_production_and_trade#Wood-based_industries. Accessed 10 9 2021.Food and Agriculture Organization of the United Nations. FAOSTAT: Forestry Production and Trade (2021). (Online). http://www.fao.org/faostat/en/#data. Accessed 13 September 2021.The Intergovernmental Panel on Climate Change. Refinement to the 2006 IPCC Guidelines for National Greenhouse Gas Inventories (PCC Task Force on National Greenhouse Gas Inventories, 2019).
Google Scholar 
European Parliament; European Council. Commission Delegated Regulation (EU) 2019/807 of 13 March 2019 Supplementing Directive (EU) 2018/2001 of the European Parliament and of the Council as Regards the Determination of High Indirect Land-Use Change-Risk (2018) (Online). fttps://eur-lex.europa.eu/eli/reg_del/2019/807/oj.de Oliveira Garcia, W., Amann, T. & Hartmann, J. Increasing biomass demand enlarges negative forest nutrient budget areas in wood export regions. Sci. Rep. 8, 5280 (2018).ADS 
PubMed 
PubMed Central 

Google Scholar 
Searchinger, T. et al. Europe’s renewable energy directive poised to harm global forests. Nat. Commun. 9, 3741 (2018).ADS 
PubMed 
PubMed Central 

Google Scholar 
Galik, C. S. & Abt, R. C. Sustainability guidelines and forest market response: An assessment of European Union pellet demand in the southeastern United States. GCB Bioenergy 8, 658–669 (2016).
Google Scholar 
Favero, A. D. & Sohngen, B. Forests: Carbon sequestration, biomass energy, or both?. Sci. Adv. 6(13), eaay6792 (2020).ADS 
PubMed 
PubMed Central 

Google Scholar 
Cowie, A. et al. Applying a science-based systems perspective to dispel misconceptions about climate effects of forest bioenergy. GCB-Bioenergy 13, 1210–1231 (2021).
Google Scholar 
Camia, A, Jonsson, G. J. R., Robert, N., Cazzaniga, N., Jasinevičius, G., Avitabile, V., Grassi, G., Barredo, J., & Mubareka, S. The Use of Woody Biomass for Energy Production in the EU (European Commission, Joint Research Center, 2021).Aguilar, F. X., Mirzaee, A., McGarvey, R., Shifley, S. & Burtraw, D. Expansion of US wood pellet industry points to positive trends but the need for continued monitoring. Sci. Rep. 10, 18607 (2020).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Dale, V., Parish, E., Kline, K. & Tobin, E. How is wood-based pellet production affecting forest conditions in the southeastern United States?. For Ecol Manag 396, 143–14 (2017).
Google Scholar 
Ceccherini, G. et al. Abrupt increase in harvested forest area over Europe after 2015. Nature 583, 72–77 (2020).ADS 
CAS 
PubMed 

Google Scholar 
FORISK Consulting. U.S. Wood Bioenergy Database (2020). (Online). https://forisk.com/. Accessed 2020.Domke, G. et al. Toward inventory-based estimates of soil organic carbon in forests of the United States. Ecol. Appl. 27(4), 1223–1235 (2017).CAS 
PubMed 

Google Scholar 
Python Org. Python Programming Language (2022) (Online). https://www.python.org/. Accessed 1 January 2018.STATA. Stata: statistical software for data science (2022) (Online). https://www.stata.com/. Accessed 1 January 2018.QGIS. Free and Open Source Geographic Information System (2021). (Online). https://qgis.org/en/site/.US Department of Agriculture, Forest Service. Forest Inventory and Analysis National Program (2020). (Online). https://www.fia.fs.fed.us/.Burrill, E. A., Wilson, A. M., Turner, J. A., Pugh, S. A., Menlove, J., Christiansen, G., Conkling, B., & David, W. The Forest Inventory and Analysis Database: Database Description and User Guide Version 8.0 for Phase 2 (US Department of Agriculture, US Forest Service, 2018).Ahmed, M. et al. Spatially-explicit modeling of multi-scale drivers of aboveground forest biomass and water yield in watersheds of the Southeastern United States. J. Environ. Manag. 199, 158–171 (2017).
Google Scholar 
Timilsina, N. et al. A framework for identifying carbon hotspots and forest management drivers. J. Environ. Manag. 114, 293–302 (2012).
Google Scholar 
Coulston, J., Ritters, K., McRoberts, R., Reams, G. & Smith, W. True versus perturbed forest inventory plot locations for modeling: A simulation study. Can. J. For. Res. 36, 801–807 (2006).
Google Scholar 
Anselin, L. Spatial effects in econometric practice in environmental and resource economics. Am. J. Agric. Econ. 83(3), 705–710 (2001).MathSciNet 

Google Scholar 
Strange-Olesen, A., Bager, S., Kittler, B., Price, W., & Aguilar, F. Environmental Implications of Increased Reliance of the EU on Biomass from the South East US (European Commission Report ENV.B.1/ETU/2014/0043, 2015).Spelter, H., & Toth, D. North America’s Wood Pellet Sector (U.S. Department of Agriculture, Forest Service, Forest Products Laboratory, 2009).Goerndt, M., Aguilar, F. & Skog, K. Drivers of biomass co-firing in US coal-fired power plants. Biomass Bioenerg. 58, 158–167 (2013).
Google Scholar 
US Department of Agriculture, Forest Service. Forest Inventory and Analysis National Program: Timber Products Output Studies (2022). (Online). https://www.fia.fs.fed.us/program-features/tpo/. Accessed 2022.Sonter, L. et al. Mining drives extensive deforestation in the Brazilian Amazon. Nat. Commun. 8(1013), 66. https://doi.org/10.1038/s41467-017-00557-w (2017).CAS 

Google Scholar 
Mirzaee, A., McGarvey, R., Aguilar, F. & Schliep, E. Impact of biopower generation on eastern US forests. Environ. Dev. Sustain. https://doi.org/10.1007/s10668-022-02235-4 (2022).
Google Scholar 
Brandeis, C., Taylor, M., Abt, K., & Alderman, D. Status and Trends for the U.S. Forest Products Sector: A Technical Document Supporting the Forest Service 2020 RPA Assessment (US Department of Agriculture, Forest Service Southern Research Station, Forest Inventory and Analysis, 2021).US Environmental Protection Agency. Emissions & Generation Resource Integrated Database (eGRID) (2021) (Online). https://www.epa.gov/egrid.US Department of Transportation. Ports: ArcGIS Online (2021) (Online). https://data-usdot.opendata.arcgis.com/datasets/usdot::ports/about.US Census Bureau. TIGER/Line Shapefiles (2021) (Online). https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html.US Census Bureau. Population and Housing Units Estimates Datasets (2021) (Online). https://www.census.gov/programs-surveys/popest/data/data-sets.html.McCann, P. The Economics of Industrial Location: A Logistics-Costs Approach (Springer, 1998).Singh, D., Cubbage, F., Gonzalez, R. & Abt, R. Locational determinants for wood pellet plants: A review and case study of North and South America. BioResources 11(3), 7928–7952 (2016).
Google Scholar 
Boukherroub, T., LeBel, L. & Lemieux, S. An integrated wood pellet supply chain development: Selecting among feedstock sources and a range of operating scales. Appl. Energy 198, 385–400 (2017).
Google Scholar 
Heckman, J., Ichimura, H. & Todd, P. Matching as an econometric evaluation estimator: Evidence from evaluating a JobTraining Programme. Rev. Econ. Stud. 64(4), 605–654 (1997).MATH 

Google Scholar 
Caliendo, M. & Kopeinig, S. Some practical guidance for the implementation of propensity score matching. J. Econ. Surv. 22(1), 31–72 (2008).
Google Scholar 
Woo, H., Eskelson, B. & Monleon, V. Matching methods to quantify wildfire effects on forest carbon mass in the U.S. Pacific Northwest. Ecol. Appl. 31(3), e02283 (2021).PubMed 

Google Scholar 
Morreale, L., Thompson, J., Tang, X., Reinmann, A. & Hutyra, L. Elevated growth and biomass along temperate forest edges. Nat. Commun. 12(7181), 66 (2021).
Google Scholar 
Isard, W. The general theory of location and space-economy. Q. J. Econ. 63(4), 476–506 (1949).
Google Scholar 
Aguilar, F. X. Spatial econometric analysis of location drivers in a renewable resource-based industry: The U.S. South Lumber Industry. For. Policy Econ. 11(3), 184–193 (2009).
Google Scholar 
Aguilar, F. X. Conjoint analysis of industry location preferences: evidence from the softwood lumber industry in the US. Appl. Econ. 66, 3265–3274 (2010).
Google Scholar 
Aguilar, F. X., Goerndt, M., Song, N. & Shifley, S. Internal, external and location factors influencing cofiring of biomass with coal in the US northern region. Energy Econ. 34, 1790–1798 (2012).
Google Scholar 
Ferraro, P. J. et al. Estimating the impacts of conservation on ecosystem services and poverty by integrating modeling and evaluation. Proc. Natl. Acad. Sci. 112(24), 7420–7425 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Zhang, D. & Pearse, P. Forest Economics 412 (UBC Press, 2011).
Google Scholar 
Villalobos, L., Coria, J. & Nordén, L. Has forest certification reduced forest degradation in Sweden?. Land Econ. 94, 220–238 (2018).
Google Scholar 
Wooldridge, J. Econometric Analysis of Cross Section and Panel Data (MIT Press, 2010).Blackman, A., Corral, L., Lima, E. & Asner, G. Titling indigenous communities protects forests in the Peruvian Amazon. PNAS 114(16), 4123–4128 (2016).ADS 

Google Scholar 
Abt, K. L., Abt, R. C., Galik, C. S., & Skog, K. E. Effect of Policies on Pellet Production and Forests in the U.S. South: A Technical Document Supporting the Forest Service Update of the 2010 RPA Assessment USDA (Forest Service GTR Srs-202, 2014).Hardie, P. Parks, P. Gottleib and D. Wear, “Responsiveness of rural and urban land uses to land rent determinants in the U.S. South,” Land Economics, vol. 76, no. 4, pp. 659–673, 2000.Parish, E., Herzberger, A., Phifer, C. & Dale, V. Transatlantic wood pellet trade demonstrates telecoupled benefits. Ecol. Soc. 23(1), 28 (2018).
Google Scholar 
Titus, B. et al. Sustainable forest biomass: A review of current residue harvesting guidelines. Energy Sustain. Soc. 11, 66. https://doi.org/10.1186/s13705-021-00281-w (2021).
Google Scholar 
Jandl, R. et al. How strongly can forest management influence soil carbon sequestration?. Geoderma 137(3), 253–268 (2007).ADS 
CAS 

Google Scholar 
Nave, L., Vance, E., Swanston, C. & Cepas, P. S. Harvest impacts on soil carbon storage in temperate forests. For. Ecol. Manag. 259, 857–866 (2010).
Google Scholar 
Mayer, M. et al. Tamm review: Influence of forest management activities on soil organic carbon stocks: A knowledge synthesis. For. Ecol. Manag. 466, 118127 (2020).
Google Scholar 
Berryman, E., Hatten, J., Page-Dumroese, D. S., Heckman, K. A., D’Amore, D. V., Puttere, J., & Domke, G. M. Soil carbon in Forest and Rangeland Soils of the United States Under Changing Conditions 9–31 (Springer, 2020).Nave, L. E. et al. Land use and management effects on soil carbon in US Lake States, with emphasis on forestry, fire, and reforestation. Ecol. Appl. 66, 2356 (2021).
Google Scholar 
Cao, B., Domke, G. M., Russell, M. B. & Walters, B. Spatial modeling of litter and soil carbon stocks on forest land in the conterminous United States. Sci. Total Environ. 654, 94–106 (2019).ADS 
CAS 
PubMed 

Google Scholar 
Coulston, J. & Wear, D. From sink to source: Regional variation in U.S. forest carbon futures. Sci. Rep. 5, 66. https://doi.org/10.1038/srep16518 (2015).
Google Scholar 
Röder, M., Whittaker, C. & Thornley, P. How certain are greenhouse gas reductions from bioenergy? Life cycle assessment and uncertainty analysis of wood pellet-to-electricity supply chains from forest residues. Biomass Bioenerg. 79, 50–63 (2015).
Google Scholar 
Hanssen, S., Duden, A., Junginger, M., Dale, D. & D. vander Hilst,. Wood pellets, what else? Greenhouse gas parity times of European electricity from wood pellets produced in the south-eastern United States using different softwood feedstocks. GC-Bioenergy 9(9), 1406–1422 (2017).CAS 

Google Scholar 
Picciano, P., Aguilar, F., Burtraw, D. & Mirzaee, A. Environmental and socio-economic implications of woody biomass co-firing at coal-fired power plants. Resour. Energy Econ. 6, 66 (2022).
Google Scholar 
Hetchner, S., Schelhas, J., & Brosius, J. Forests as Fuel: Energy, Landscape, Climate, and Race in the U.S. South (Lexington Books, 2022).Coulston, J., Wear, D. & Vose, J. Complex forest dynamics indicate potential for slowing carbon accumulation in the southeastern United States. Sci. Rep. 5, 8002 (2015).ADS 
PubMed 
PubMed Central 

Google Scholar 
Palahí, M. et al. Concerns about reported harvests in European forests. Nature 592, E15–E17 (2021).PubMed 

Google Scholar  More

Hug, L. A. et al. A new view of the tree of life. Nat. Microbiol. 1, 16048 (2016).Article 
PubMed 
CAS 

Google Scholar 
Spang, A. et al. Complex archaea that bridge the gap between prokaryotes and eukaryotes. Nature 521, 173–179 (2015).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Tyson, G. W. et al. Community structure and metabolism through reconstruction of microbial genomes from the environment. Nature 428, 37–43 (2004).Article 
PubMed 
CAS 

Google Scholar 
Anantharaman, K. et al. Thousands of microbial genomes shed light on interconnected biogeochemical processes in an aquifer system. Nat. Commun. 7, 13219 (2016).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Parks, D. H. et al. Recovery of nearly 8,000 metagenome-assembled genomes substantially expands the tree of life. Nat. Microbiol. 2, 1533–1542 (2017).Article 
PubMed 
CAS 

Google Scholar 
Tully, B. J. & Graham, E. D. & Heidelberg, J. F. The reconstruction of 2,631 draft metagenome-assembled genomes from the global oceans. Sci. Data 5, 170203 (2018).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Stewart, R. D. et al. Assembly of 913 microbial genomes from metagenomic sequencing of the cow rumen. Nat. Commun. 9, 870 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Pasolli, E. et al. Extensive unexplored human microbiome diversity revealed by over 150,000 genomes from metagenomes spanning age, geography and lifestyle. Cell 176, 649–662 (2019).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Nayfach, S. et al. A genomic catalog of Earth’s microbiomes. Nat. Biotechnol. 39, 499–509, https://doi.org/10.1038/s41587-020-0718-6 (2021).Article 
PubMed 
CAS 

Google Scholar 
Gilbert, J. A., Jansson, J. K. & Knight, R. The Earth Microbiome project: successes and aspirations. BMC Biol 12, 69 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
Saheb Kashaf, S., Almeida, A., Segre, J. A. & Finn, R. D. Recovering prokaryotic genomes from host-associated, short-read shotgun metagenomic sequencing data. Nat. Protoc. 16, 2520–2541 (2021).Article 
PubMed 
CAS 

Google Scholar 
Chong, J., Liu, P., Zhou, G. & Xia, J. Using MicrobiomeAnalyst for comprehensive statistical, functional, and meta-analysis of microbiome data. Nat. Protoc. 15, 799–821 (2020).Article 
PubMed 
CAS 

Google Scholar 
Arkin, A. P. et al. KBase: The United States Department of Energy Systems Biology Knowledgebase. Nat. Biotechnol. 36, 566–569 (2018).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Sayers, E. W. et al. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res 49, D10–D17 (2021).Article 
PubMed 
CAS 

Google Scholar 
Kluyver, T., et al. Jupyter Notebooks – a publishing format for reproducible computational workflows. In: Loizides F, Schmidt B, editors. Positioning and Power in Academic Publishing: Players, Agents and Agendas. p. 87–90 (2016).Banfield, J. Development of a Knowledgebase to Integrate, Analyze, Distribute, and Visualize Microbial Community Systems Biology Data. (2015). Report number: DOE-UCB-4918, OSTI ID: 1167269.Chen, I.-M. A. et al. IMG/M v.5.0: an integrated data management and comparative analysis system for microbial genomes and microbiomes. Nucleic Acids Res 47, D666–D677 (2019).Article 
PubMed 
CAS 

Google Scholar 
Afgan, E. et al. The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2016 update. Nucleic Acids Res 44, W3–W10 (2016).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Devisetty, U. K., Kennedy, K., Sarando, P., Merchant, N. & Lyons, E. Bringing your tools to CyVerse discovery environment using Docker. F1000Res. 5, 1442 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Wang, L., Lu, Z., Van Buren, P. & Ware, D. SciApps: a bioinformatics workflow platform powered by XSEDE and CyVerse. in Proceedings of the Practice and Experience on Advanced Research Computing 1–5 (Association for Computing Machinery, 2018).Eren, A. M. et al. Community-led, integrated, reproducible multi-omics with anvi’o. Nat. Microbiol. 6, 3–6 (2021).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Wattam, A. R. et al. Improvements to PATRIC, the all-bacterial bioinformatics database and analysis resource center. Nucleic Acids Res 45, D535–D542 (2017).Article 
PubMed 
CAS 

Google Scholar 
Mitchell, A. L. et al. MGnify: the microbiome analysis resource in 2020. Nucleic Acids Res. 48, D570–D578 (2020).PubMed 
CAS 

Google Scholar 
Wu, Y.-W. et al. Ionic liquids impact the bioenergy feedstock-degrading microbiome and transcription of enzymes relevant to polysaccharide hydrolysis. mSystems 1, e00120–16 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Rajeev, L. et al. Dynamic cyanobacterial response to hydration and dehydration in a desert biological soil crust. ISME J 7, 2178–2191 (2013).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Foster, I. Globus Online: accelerating and democratizing science through cloud-based services. IEEE Internet Comput 15, 70–73 (2011).Article 

Google Scholar 
Nurk, S., Meleshko, D., Korobeynikov, A. & Pevzner, P. A. metaSPAdes: a new versatile metagenomic assembler. Genome Res 27, 824–834 (2017).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Zhang, H. et al. dbCAN2: a meta server for automated carbohydrate-active enzyme annotation. Nucleic Acids Res 46, W95–W101 (2018).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Chaumeil, P.-A., Mussig, A. J., Hugenholtz, P. & Parks, D. H. GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database. Bioinformatics 36, 1925–1927 (2019).PubMed Central 

Google Scholar 
Camacho, C. et al. BLAST+: architecture and applications. BMC Bioinforma 10, 421 (2009).Article 

Google Scholar 
Nordberg, H. et al. The genome portal of the Department of Energy Joint Genome Institute: 2014 updates. Nucleic Acids Res 42, D26–D31 (2014).Article 
PubMed 
CAS 

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

Google Scholar 
Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. J. 17, 10–12 (2011).Article 

Google Scholar 
Menzel, P., Ng, K. L. & Krogh, A. Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nat. Commun. 7, 11257 (2016).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Freitas, T. A. K., Li, P.-E., Scholz, M. B. & Chain, P. S. G. Accurate read-based metagenome characterization using a hierarchical suite of unique signatures. Nucleic Acids Res 43, e69 (2015).Article 
PubMed 
PubMed Central 

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

Google Scholar 
Truong, D. T. et al. MetaPhlAn2 for enhanced metagenomic taxonomic profiling. Nat. Methods 12, 902–903 (2015).Article 
PubMed 
CAS 

Google Scholar 
Milanese, A. et al. Microbial abundance, activity and population genomic profiling with mOTUs2. Nat. Commun. 10, 2014 (2019).Article 

Google Scholar 
Youngblut, N. D. & Ley, R. E. Struo2: efficient metagenome profiling database construction for ever-expanding microbial genome datasets. Peer J 9, e12198 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Ondov, B. D., Bergman, N. H. & Phillippy, A. M. Interactive metagenomic visualization in a Web browser. BMC Bioinform 12, 385 (2011).Article 

Google Scholar 
Li, D., Liu, C.-M., Luo, R., Sadakane, K. & Lam, T.-W. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics 31, 1674–1676 (2015).Article 
PubMed 
CAS 

Google Scholar 
Peng, Y., Leung, H. C. M., Yiu, S. M. & Chin, F. Y. L. IDBA-UD: a de novo assembler for single-cell and metagenomic sequencing data with highly uneven depth. Bioinformatics 28, 1420–1428 (2012).Article 
PubMed 
CAS 

Google Scholar 
Orakov, A. et al. GUNC: detection of chimerism and contamination in prokaryotic genomes. Genome Biol 22, 178 (2021).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Gurevich, A., Saveliev, V., Vyahhi, N. & Tesler, G. QUAST: quality assessment tool for genome assemblies. Bioinformatics 29, 1072–1075 (2013).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Wu, Y.-W., Simmons, B. A. & Singer, S. W. MaxBin 2.0: an automated binning algorithm to recover genomes from multiple metagenomic datasets. Bioinformatics 32, 605–607 (2016).Article 
PubMed 
CAS 

Google Scholar 
Kang, D. D. et al. MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ 7, e7359 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Alneberg, J. et al. Binning metagenomic contigs by coverage and composition. Nat. Methods 11, 1144–1146 (2014).Article 
PubMed 
CAS 

Google Scholar 
Sieber, C. M. K. et al. Recovery of genomes from metagenomes via a dereplication, aggregation and scoring strategy. Nat. Microbiol. 3, 836–843 (2018).Article 
PubMed 
PubMed Central 
CAS 

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

Google Scholar 
Delcher, A. L., Salzberg, S. L. & Phillippy, A. M. Using MUMmer to identify similar regions in large sequence sets. Curr. Protoc. Bioinform. Chapter 10, Unit 10.3 (2003).
Google Scholar 
Darling, A. C. E., Mau, B., Blattner, F. R. & Perna, N. T. Mauve: multiple alignment of conserved genomic sequence with rearrangements. Genome Res 14, 1394–1403 (2004).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Parks, D. H. et al. GTDB: an ongoing census of bacterial and archaeal diversity through a phylogenetically consistent, rank normalized and complete genome-based taxonomy. Nucleic Acids Res 50, D785–D794 (2022).Article 
PubMed 
CAS 

Google Scholar 
Bowers, R. M. et al. Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea. Nat. Biotechnol. 35, 725–731 (2017).Article 
PubMed 
PubMed Central 
CAS 

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

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

Google Scholar 
Seemann, T. Prokka: rapid prokaryotic genome annotation. Bioinformatics 30, 2068–2069 (2014).Article 
PubMed 
CAS 

Google Scholar 
Hyatt, D. et al. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinform 11, 119 (2010).Article 

Google Scholar 
Parks, D. H. et al. A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life. Nat. Biotechnol. 36, 996–1004 (2018).Article 
PubMed 
CAS 

Google Scholar 
Rinke, C. et al. A standardized archaeal taxonomy for the Genome Taxonomy Database. Nat. Microbiol. 6, 946–959 (2021).Article 
PubMed 
CAS 

Google Scholar 
Haft, D. H. et al. RefSeq: an update on prokaryotic genome annotation and curation. Nucleic Acids Res 46, D851–D860 (2018).Article 
PubMed 
CAS 

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

Google Scholar 
Shaffer, M. et al. DRAM for distilling microbial metabolism to automate the curation of microbiome function. Nucleic Acids Res 48, 8883–8900 (2020).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Galperin, M. Y., Makarova, K. S., Wolf, Y. I. & Koonin, E. V. Expanded microbial genome coverage and improved protein family annotation in the COG database. Nucleic Acids Res 43, D261–D269 (2015). (Database Issue).Article 
PubMed 
CAS 

Google Scholar 
El-Gebali, S. et al. The Pfam protein families database in 2019. Nucleic Acids Res 47, D427–D432 (2019).Article 
PubMed 
CAS 

Google Scholar 
Haft, D. H. et al. TIGRFAMs and Genome Properties in 2013. Nucleic Acids Res 41, D387–D395 (2013). (Database issue).Article 
PubMed 
CAS 

Google Scholar 
Eddy, S. R. Accelerated Profile HMM Searches. PLoS Comput. Biol. 7, e1002195 (2011).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Lombard, V., Golaconda Ramulu, H., Drula, E., Coutinho, P. M. & Henrissat, B. The carbohydrate-active enzymes database (CAZy) in 2013. Nucleic Acids Res 42, D490–D495 (2014).Article 
PubMed 
CAS 

Google Scholar 
Chivian, D., Dehal, P. S., Keller, K. & Arkin, A. P. MetaMicrobesOnline: phylogenomic analysis of microbial communities. Nucleic Acids Res 41, D648–D654 (2013).Article 
PubMed 
CAS 

Google Scholar 
Karaoz, U. & Brodie, E. L. microTrait: a toolset for a trait-based representation of microbial genomes. Front. Bioinform. https://doi.org/10.3389/fbinf.2022.918853 (2022).Article 
PubMed 
PubMed Central 

Google Scholar 
Wood-Charlson, E. M. et al. The National Microbiome Data Collaborative: enabling microbiome science. Nat. Rev. Microbiol. 18, 313–314 (2020).Article 
PubMed 
CAS 

Google Scholar 
Hofmeyr, S. et al. Terabase-scale metagenome coassembly with MetaHipMer. Sci. Rep. 10, 10689 (2020).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Kolmogorov, M. et al. metaFlye: scalable long-read metagenome assembly using repeat graphs. Nat. Methods 17, 1103–1110 (2020).Article 
PubMed 
CAS 

Google Scholar 
Koren, S. et al. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res 27, 722–736 (2017).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Bertrand, D. et al. Hybrid metagenomic assembly enables high-resolution analysis of resistance determinants and mobile elements in human microbiomes. Nat. Biotechnol. 37, 937–944 (2019).Article 
PubMed 
CAS 

Google Scholar 
Chen, L.-X. et al. Accurate and complete genomes from metagenomes. Genome Res 30, 315–333 (2020).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Lui, L. M., Nielsen, T. N. & Arkin, A. P. A method for achieving complete microbial genomes and improving bins from metagenomics data. PLoS Comput Biol 17, e1008972 (2021).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Miller, C. S., Baker, B. J., Thomas, B. C., Singer, S. W. & Banfield, J. F. EMIRGE: reconstruction of full-length ribosomal genes from microbial community short read sequencing data. Genome Biol 12, R44 (2011).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Chivian, D. et al. Genome extraction from shotgun metagenome sequence data. KBase n/33233/628 https://doi.org/10.25982/33233.606/1831502 (2022).Article 

Google Scholar 
Chivian, D., et al. Moab desert crust – sample 4E. KBase n/62384/334 (2022). https://doi.org/10.25982/62384.253/1831503Jain, C., Rodriguez-R, L. M., Phillippy, A. M., Konstantinidis, K. T. & Aluru, S. High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat. Commun. 9, 5114 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Matsen, F. A., Kodner, R. B. & Armbrust, E. V. pplacer: linear time maximum-likelihood and Bayesian phylogenetic placement of sequences onto a fixed reference tree. BMC Bioinform 11, 538 (2010).Article 

Google Scholar 
Benson, D. A. et al. GenBank. Nucleic Acids Res 46, D41–D47 (2018).Article 
PubMed 
CAS 

Google Scholar 
Ewing, B. & Green, P. Base-calling of automated sequencer traces using phred. II. Error probabilities. Genome Res. 8, 186–194 (1998).Article 
PubMed 
CAS 

Google Scholar 
Teiling, C. BaseSpace: Simplifying metagenomic analysis. 26th European Congress of Clinical Microbiology and Infectious Diseases (2016) 10.26226/morressier.56d5ba2ed462b80296c9509dReich, M. et al. The GenePattern notebook environment. Cell Syst 5, 149–151.e1 (2017).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Uritskiy, G. V., DiRuggiero, J. & Taylor, J. MetaWRAP-a flexible pipeline for genome-resolved metagenomic data analysis. Microbiome 6, 158 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Karp, P. D. et al. A comparison of microbial genome web portals. Front. Microbiol. 10, 208 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Yue, Y. et al. Evaluating metagenomics tools for genome binning with real metagenomic datasets and CAMI datasets. BMC Bioinform 21, 334 (2020).Article 
CAS 

Google Scholar 
Nelson, W. C., Tully, B. J. & Mobberley, J. M. Biases in genome reconstruction from metagenomic data. PeerJ 8, e10119 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Olm, M. R., Brown, C. T., Brooks, B. & Banfield, J. F. dRep: a tool for fast and accurate genomic comparisons that enables improved genome recovery from metagenomes through de-replication. ISME J 11, 2864–2868 (2017).Article 
PubMed 
PubMed Central 
CAS 

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

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

Google Scholar 
Kim, D., Paggi, J. M., Park, C., Bennett, C. & Salzberg, S. L. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37, 907–915 (2019).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Pertea, M. et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295 (2015).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
Kumari, S. et al. A KBase case study on genome-wide transcriptomics and plant primary metabolism in response to drought stress in sorghum. Curr. Plant Biol. 28, 100229 (2021).Article 
CAS 

Google Scholar 
Seaver, S. M. D. et al. The ModelSEED biochemistry database for the integration of metabolic annotations and the reconstruction, comparison and analysis of metabolic models for plants, fungi and microbes. Nucleic Acids Res 49, D575–D588 (2021).Article 
PubMed 
CAS 

Google Scholar 
Schloss, P. D. et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75, 7537–7541 (2009).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Caporaso, J. G. et al. QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336 (2010).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar  More

Sample collection and isolation of V. cholerae O1 possessing the CT geneTwenty-three patients (patient numbers 9 to 31) who were diagnosed with cholera were examined in this study. The diagnosis was confirmed by the isolation of V. cholerae O1 from the stool of each patient. The age of patients, date of hospital admission, stool sampling date, pathogen isolated and medicines administered to the patients as treatments are described in Supplementary Table S1. Twenty-one of the stool samples were collected on the first day of hospitalization, while the remaining two stool samples were collected on the second day (patient number 29) and fourth day (patient number 10) of hospitalization. All patients had not been given any antibiotics and the samples of diarrheal stool were taken during severe diarrhea.To confirm the presence of the CT gene (ctx) in these 23 isolates, we examined the presence of ctxA in these isolates by PCR. The PCR to detect ctxA was performed as reported by Keasler and Hall6. In this PCR, amplification was performed in 30 cycles. The size of the amplified ctxA fragment was 302 bp. The target fragment was amplified from each of the V. cholerae O1 isolates. This indicated that all of the V. cholerae O1 isolates from the 23 cholera patients possessed ctxA.CT production from the isolatesThe production of CT from these 23 isolates was examined by detecting secreted CT in the medium. The 23 isolates were cultured statically in AKI medium7, and the secreted CT in the culture supernatants was measured using the GM1-ganglioside enzyme-linked immunosorbent assay (ELISA) method8. The detection limit of CT by the ELISA method used is 1.0 ng ml−1. All the samples examined were found to have CT above this concentration (Fig. 1). This shows that all isolates examined are toxigenic V. cholerae O1.Figure 1Amount of cholera toxin produced by V. cholerae O1 isolated from patients with diarrhea. Twenty-three strains of V. cholerae O1 were isolated from 23 patients with diarrhea. These isolates were cultured statically in AKI medium7 at 37 °C for 24 h. After removing the cells by centrifugation, the CT in the culture supernatants was measured using a GM1-ganglioside ELISA method8. The samples indicated by blue circle are isolates obtained by bacterial culture from two patients (patient 12 and patient 18), who are focused on in this study.Full size imageAnalysis of the stool samples of patients diagnosed with cholera diseaseMetagenomic sequencing analysisThe primary objective of this metagenomic analysis is to show the proportion of V. cholerae living in the diarrhea stool. Subsequently, if the number of V. cholerae infected in the intestinal tract is small, it is required to clarify the etiological microorganisms that cause diarrhea in that patient. For this analysis, it is necessary to investigate the presence of pathogenic microorganisms other than V. cholerae in the stool. To do this, we need to analyze the gene reads obtained by metagenomic analysis with a comprehensive manner. Therefore, we planned to obtain reads with the Burrows-Wheeler Alignment tool (BWA) with default parameters, a matching software with the ability to fulfill these objectives9 (http://bio-bwa.sourceforge.net).However, we were concerned that the genes derived from organisms other than V. cholerae in the stool were counted mistakenly as genes derived from V. cholerae in the analysis using BWA. We therefore first examined genes in stool from people unrelated to cholera disease to ensure that the analysis method we planned to use in this study would correctly detect genes from V. cholerae in stool. For this analysis we used DNA sequences reported by the NIH Human Microbiome Project (https://www.hmpdacc.org/hmp/hmp/hmasm2/). The genes we have analyzed are DNA derived from feces of 20 healthy individuals (10 males and 10 females). The results are shown in Supplementary Table S2.The number of reads analyzed in this analysis varied from sample to sample. The largest number obtained after quality filtering was 60,975,797. The lowest number was 10,301,809. However, the number of reads detected as originating from V. cholerae was very small (12 reads or less) in all samples, and none of them were detected in 7 samples. This very small number shows that the analytical method used is suitable for detecting the genes from V. cholerae in these samples.Therefore, we analyzed DNA and RNA samples from prepared diarrheal stool by the method using BWA. All raw sequencing data obtained were deposited into the DDBJ Sequence Read Archive under the accession code PRJDB10675. This number can be searched not only from DDBJ but also from EMBL and GenBank.Diarrheal stools are mostly composed of liquid, and their properties are very different from those of normal stools. The origin of the nucleic acids in diarrheal stools varies from patient to patient and is not constant. One sample may contain many genes derived from human cells, while another sample may contain many genes derived from microorganisms. To clarify the nature of the reads we obtained, we determined the proportion of reads of bacterial origin to the total number of reads in the samples analyzed, and presented this proportion in order of patient age (Fig. 2a). The ratios were not consistent, indicating that the cells of eukaryotic origin and microorganisms existing in the stool of patients with diarrhea varied from person to person.Figure 2Age of patient and the ratio of the number of read detected by metagenomic sequencing analysis of their stools. The DNA in the stool samples from 23 patients who were diagnosed with cholera disease were extracted using a commercially available kit. Patient ages are listed in Supplementary Table S1. The extracted DNA were investigated by a metagenomic sequencing analysis to clarify the origin of individual DNA. The origin of the DNA sequences was assigned by mapping to a database that included human and microorganism sequences. The obtained numbers of total reads, total bacterial reads, reads originating from V. cholerae, reads originating from ctxA in each sample are shown in Supplementary Table S3 (the data from DNA sample). The age of each patient and the ratio of the number of reads from all bacteria to the total number of reads after filtering (a) and the ratio of the number of reads from V. cholerae to the number of reads from all bacteria (b) were calculated. The horizontal axis of these figures shows the age of each patient and is the same arrangement in both (a) and (b). The numbers in parentheses indicate the sample numbers. This sample number is also the patient’s number.Full size imageThis result implied that it was difficult to detect V. cholerae in a sample with a small number of read derived from bacteria. Therefore, it was unclear whether the data obtained by the analysis was suitable for the detection of V. cholerae. In order to examine whether the data shown in Fig. 2a can be used to clarify the infection status of V. cholerae, the ratio of the reads from V. cholerae to the reads of all bacteria in the sample was calculated (Fig. 2b). As a result, the reads from V. cholerae were detected even in samples with a low ratio of bacterial genes, as seen in patients 13, 25, and 29. Conversely, some patients, such as patients 10, 18, and 17, had a high proportion of bacterial genes but a low detection rate of the read from V. cholerae (Fig. 2a,b). From these results, we thought that the data obtained are useful for analyzing the infection status of V. cholerae in the intestinal tract of the examined patients. The data also showed that patient age did not affect the intestinal retention of V. cholerae.In order to more clearly illustrate the presence of V. cholerae in the diarrheal stools of the patients examined, the ratio of reads from V. cholerae to total reads for each sample which was determined in Fig. 2a was sorted in descending order. The results are shown in Fig. 3a. The ratio (percentage) in each patient is indicated by the blue bar in the figure. The numbers in parentheses after the sample number, with D as the first letter, indicate the order from lowest to highest percentage obtained. As shown in Fig. 3a, the percentage of V. cholerae that the patients carried in their stools varied from 0.003% (sample 12(D1)) to 38.337% (sample 28(D23)).Figure 3The ratios of DNA and RNA derived from V. cholerae in stool samples. The DNA and RNA in the stool samples from 23 patients who were diagnosed with cholera disease were extracted using a commercially available kit. Subsequently, the RNA samples were treated with DNase I to remove DNA from the samples. Reverse-transcribed DNA was prepared from these RNA samples using random primers and reverse transcriptase. The extracted DNA and reverse-transcribed DNA were investigated by a metagenomic sequencing analysis to clarify the origin of individual DNA and RNA. The origin of the DNA sequences was assigned by mapping to a database that included human and microorganism sequences. The obtained numbers of total reads, total bacterial reads, reads originating from V. cholerae, reads originating from ctxA in each sample are shown in Supplementary Tables S3 (the data from DNA) and S4 (the data from RNA). The percentages of reads of DNA from V. cholerae and from ctxA relative to the total reads are presented by blue bar and red bar in panel a, respectively. The percentages of reads of DNA from V. cholerae relative to the total bacterial reads are presented in panels b. Samples are arranged in ascending order of the ratio of reads from V. cholerae to the total reads in the DNA analysis. The ranking of each sample is presented by the numbers in parentheses starting with the letter D. The samples in these panels are arranged in the order of the D number. Similarly, the ratio of reads from V. cholerae to the total RNA reads and the total bacterial RNA are presented in panels c and d, respectively. The samples indicated by red circle are samples from the diarrheal stools of a patients who are focused on in this analysis.Full size imageHowever, what we want to reveal in this study is the presence of toxigenic V. cholerae producing CT. The genes presented by blue bar in Fig. 3a appear to contain the genes derived from toxigenic V. cholerae, but it cannot be concluded that they are. It is highly possible that other genes derived from such as V. cholerae not possessing ctx or bacteria having the same gene sequence as V. cholerae, are included. So, in order to examine the existence of V. cholerae possessing ctx, we examined the number of reads derived from ctxA (Supplementary Table S3). In the samples with D number 8 or more, the gene derived from ctxA was detected in all the samples except one sample (D9). The ratio of read from ctxA to the number of reads from total DNA is shown by the red bar in Fig. 3a. The ratio of the number of reads derived from ctxA to the total DNA was correlated with the ratio of the number of reads derived from the V. cholerae gene to the total DNA (Fig. 3a). From these results, it seems that most of the genes of V. cholerae detected in Fig. 3a are derived from V. cholerae possessing ctx.Furthermore, the ratio of the number of reads of V. cholerae to the number of reads derived from total bacteria which was obtained in Fig. 2a, was arranged in the order used for the array in Fig. 3a (the order of the ratio of the number of reads from V. cholerae to the number of reads from total DNA) (Fig. 3b). From this arrangement of Fig. 3b, it can be seen that the sample with a large D head number has a large proportion of V. cholerae in the bacteria. The highest value was obtained from sample 24 (D22). The sample showed that 95.917% of the bacteria was V. cholerae.On the other hand, in many samples with small D numbers, this ratio is small, but there are exceptions. For example, in samples 25 (D3), 29 (D8) and 13 (D11), the presence of V. cholerae is clear. Although not as clear as these three samples, the presence of V. cholerae in other samples such as 22 (D4), 21 (D5), 9 (D10) and 11(D12) is evident, although in small quantities (Fig. 3b). Therefore, it was considered that these patients were infected with V. cholerae. These results seem to accurately reflect the actual state of V. cholerae in the stool. Therefore, it was considered that the infection status of V. cholerae in the patient could be inferred from the obtained data.As shown in Fig. 3b, in the samples of 18 (D2), 12 (D1), 17 (D7), 10 (D9)) and 23 (D6), the ratio of the read from V. cholerae to the read from total bacteria is very low at 0.032%, 0.118%, 0.225%, 0.244% and 0.285%, respectively. It was unknown whether these patients were infected with V. cholerae and developed diarrhea due to the infection with V. cholerae. Therefore, further examination was needed to determine if these patients were infected with V. cholerae. These five samples are marked by red circles in Fig. 3a,b.Subsequently, we examined the ratio of the reads of RNA of V. cholerae to clarify the expression of the genes of V. cholerae in the intestinal lumen of these patients. RNA samples were prepared by different methods from the patient’s stool and the RNA in these samples was analyzed by metagenomic sequencing analysis. The ratio of the number of reads derived from the RNA of V. cholerae to the number of reads derived from total RNA and to the number of reads derived from total bacterial RNA in the sample was determined. The results are shown in Fig. 3c,d, respectively. Samples that had fewer reads for genes derived from V. cholerae in the previous analysis of DNA reads (Fig. 3a,b)were also indicated with a red circle in Figs. 3c,d. These samples also had low amounts of RNA read from V. cholerae. In particular, the ratio of RNA read from V. cholerae to total bacterial RNA in samples 12 (D1) and 18 (D2) was low, 0.038% and 0.236%, respectively (Supplementary Table S4, Fig. 3d). Judging from these low values, it is doubtful that these two patients, patients 12 and 18, had diarrhea due to infection with V. cholerae.Detection of ctxA by PCRSubsequently, we amplified ctxA in the DNA samples extracted from the stool samples by PCR, in order to reconfirm the presence of ctx in stool samples. The PCR was performed using the same conditions used for the detection of ctxA in the isolates as described above in the “Sample collection and isolation of V. cholerae O1 possessing the CT gene” section of the “Results”. Amplification in this PCR was also done for 30 cycles.From the results of metagenomic sequencing shown in Fig. 3, we found that the samples from patient 12 (D1) and patient 18 (D2) contained few genes derived from V. cholerae O1. The results obtained by PCR are shown in Fig. 4. The samples from the two patients, 12 (D1) and 18 (D2), are indicated by blue circle. No distinct bands corresponding to ctxA were detected in the lanes analyzed sample 12(D1). Meanwhile, a very faint band was visible in the lane where the sample from 18(D2) was analyzed. However, it often happens that small amounts of sample are mixed into adjacent lanes when adding the sample to be analyzed in agar electrophoresis. Hence, we concluded that the amount of ctxA in these two samples amplified by PCR was very low. This supports our inference that the diarrhea in these two patients was not caused by the infection with V. cholerae O1.Figure 4PCR to detect ctxA in the stool samples of diarrhea patients. DNA was extracted from the stool samples of 23 patients who were diagnosed with cholera disease. PCR to amplify ctxA in these DNA samples was performed using the specific primers ctcagacgggatttgttaggcacg and tctatctctgtagcccctattacg6, and the products were analyzed by agarose gel electrophoresis. The sample numbers are the same as the numbers shown in the footnotes of Fig. 3. Numbers beginning with D in parentheses show the order of the content of DNA from V. cholerae among these samples. The samples indicated by blue circle are samples from the diarrheal stools of patients (patients 12 and 18), who are focused on in this study. S: the size marker for gel electrophoresis; N: the negative control in which DNA was not added to the reaction mixture; P: the positive control in which DNA prepared from V. cholerae O1 N1696128 was added to the reaction mixture.Full size imageSimilarly, clear bands were not detected in samples 9(D10), 10(D9), 13(D11), 22(D4), and 25(D3). The results of metagenomic analysis of these samples showed that the number of read from V. cholerae was low and ctxA was either not detected (samples 10(D9), 22(D4) and 25(D3)) or was detected but in small amounts (samples 9(D10) and 13(D11) (Supplementary Table S3, Fig. 3a).The amount of sample added to the reaction solution in the PCR reaction was as small as 5 µl, and it is not clear whether this small volume of solution contained the necessary amount of ctxA for the amplification in PCR. It is also possible that the sample contained substances that would inhibit amplification by PCR. For these reasons, we believe that no clear band corresponding to ctxA appeared in this PCR. However, it is clear from the results of Fig. 3b,d that these samples, (9(D10), 10(D9), 13(D11), 22(D4), and 25(D3)) contain the gene derived from V. cholerae (ctx). Therefore, we considered these four patients to be patients infected with V. cholerae.The levels of CT and proteolytic activity in the stool samplesFrom the genetic studies in Figs. 3 and 4, it was inferred that V. cholerae O1 was not involved in the onset of the diarrhea in two patients (12(D1) and 18(D2)). However, this inference was based on amplification and analysis of genetic sample prepared from diarrhea stool of patients. There is no proof that the sample procurement and the analysis of sample was done reliably with high probability. Hence, we thought that it was necessary to analyze samples adjusted from different perspectives by different means.Then, we challenged to measure the amount of CT. CT is the toxin responsible for the diarrhea caused by V. cholerae O1. CT is released into the intestinal lumen, where it acts on the intestinal cells of patients to induce diarrhea. Thus, we measured the CT content in the stool samples. In addition, we also measured the proteolytic activity in the stool samples, because CT is sensitive to proteolytic activity, and we were concerned that the CT would be degraded by proteases during storage outside of the body.The CT content and the proteolytic activity in the stool samples of the 23 cholera patients were measured by the GM1-ganglioside ELISA method and the lysis of casein, respectively8,10, and the results are presented in Fig. 5a,b, respectively.Figure 5The levels of CT and proteolytic activity in the stool samples. Twenty-three stool samples of patients who were diagnosed with cholera disease were centrifuged at 10,000×g for 10 min. The CT content of the supernatants was determined using a GM1-ganglioside ELISA method8 (a). The proteolytic activity of the supernatants was determined by the lysis of casein10 (b). The sample numbers are the same as the numbers shown in the footnotes of Fig. 3. Numbers beginning with D in parentheses show the order of the content of DNA from V. cholerae among these samples. The samples in this figure are arranged in the order of the D numbers. A bar indicating the amount of CT is not drawn in the figure for the sample whose CT amount was below the detection limit. From the tests shown in Figs. 2, 3 and 4, samples of two patients who are unlikely to have diarrhea caused by the infection with V. cholerae are marked with a blue circle. O.D.: optical density.Full size imageProteolytic activity was detected in all samples, although there were differences in the strengths of the activity. It was also found that high protease activity was not associated with decreased levels of CT in the samples, e.g., sample 11(D12) showed the highest protease activity among the samples examined, and the amount of CT in that sample was also high. Therefore, we considered that the proteolytic activity had almost no influence on the amount of CT in this study. Furthermore, the fact that protease activity was found in all samples indicated that these samples were collected and stored without any significant denaturation.The ELISA method used in this assay can accurately detect CT at concentrations above 1.0 ng ml−1, but it is impossible to accurately determine the concentration of CT at concentrations below 1.0 ng ml−1. Therefore, we treated samples containing less than 1.0 ng ml−1 of CT as containing no CT.As described above, we considered that the diarrhea in the two patients (12(D1) and 18(D2)) was not due to the infection with V. cholerae O1 from the genetic analysis. The analysis of CT in stool samples showed that the CT concentrations of these two samples were below the detection limit (Fig. 5a). This indicates that the number of V. cholerae O1 in the intestinal lumen of these patients, (12(D1) and 18(D2)), was extremely low at the time of sampling.Investigation of diarrheagenic microorganisms in diarrheal stoolIt was shown that diarrhea in patients 12 (D1) and 18 (D2) may have been caused by infection with microorganisms other than V. cholerae. Then we examined the data of metagenomic sequencing of these two patients to reveal the infected diarrhea-causing microorganisms (DDBJ Sequence Read Archive under the accession code PRJDB10675). As a result, we found that that DNA from the two bacteria, Streptococcus pneumoniae and Salmonella enterica was abundant in the stools of patients 12(D1) and 18(D2), respectively.The ratios of DNA read of St. pneumoniae in DNA samples of patient 12(D1) to the total DNA and to the total bacterial DNA are 0.095% and 3.988%, respectively. These ratios of V. cholerae in this patient, 12 (D1), are 0.003% and 0.118%, respectively. And those of S. enterica in the stools of patients 18(D2) are 0.536% and 1.118%, respectively. And these ratios of V. cholerae in this patient, 18 (D2), are 0.015% and 0.032%, respectively (Supplementary Table 2).These two bacteria, St. pneumoniae and S. enterica, are bacteria that are not detected as normal intestinal bacteria. As shown, these ratios of DNA of each bacteria in diarrheal stool are much higher than these of V. cholerae. Therefore, these two bacteria are considered to be related to these patients’ symptom, respectively.Nonetheless, toxigenic V. cholerae O1 was also isolated from these two patients in laboratory bacteriology tests. It is likely that some of the very few V. cholerae O1 in the intestinal tract were extruded with the diarrhea and were subsequently detected by the enrichment culture for V. cholerae. This indicated that V. cholerae O1 may cause subclinical infections in residents of the Kolkata region of India. With this subclinical infection, the number of V. cholerae O1 inhabiting the intestinal tract might be small.Surveillance of patient samples where no diarrhea-causing microorganisms were detectedTo detect people with a subclinical infection of V. cholerae O1, we further analyzed the specific-pathogen-free stool samples of diarrhea patients. “Specific-pathogen-free stool sample” refers to the stool samples in which no etiological agent of diarrhea, including V. cholerae, was detected by our bacterial examination in the laboratory.The number of samples examined in this analysis was 22 (samples number 1001 to 1022). All 22 diarrhea patients examined were inpatients at ID hospital, Kolkata. From the 22 patients, 20 patient stool samples were collected on the 1st day of hospitalization, and the stools of the remaining two patients (patients 1004 and 2022) were collected on the 2nd day of hospitalization. Antibiotics were used in a limited manner in these patients. Ofloxacin was the only antibiotic administered, and only four patients (patients 1001, 1011, 1012, and 1021) were administered with it (Supplementary Table S1).DNA and RNA were extracted from the stool samples, and the DNA and RNA were analyzed by a metagenomic sequencing analysis using the same method used in the analysis of diarrheal stools from cholera patients.Reads of the genes from V. cholerae were detected in every sample, although the value varied from sample to sample (Supplementary Tables S5 and S6). Although reads of the genes from V. cholerae were detected in every sample, we do not believe that every stool sample examined contained V. cholerae. In the metagenomic analysis, if the base sequence of a read was common to multiple bacteria, the read was recognized as being derived from those multiple bacteria. Therefore, even if a bacterium is not present in the sample, the reads in common with other bacteria are counted as the reads of those bacteria, i.e., if a read from bacteria other than V. cholerae is homologous to a corresponding gene of V. cholerae, its detection indicates that one gene derived from V. cholerae was found in the sample. The total number of such reads is finally counted as the number of reads of V. cholerae. Therefore, it is unclear whether bacteria presenting a low read count are present in the sample. In order to solve these problems, not only the reads derived from V. cholerae but also the reads derived from ctxA were searched for in the sample.In addition, as described above, other DNA present in diarrheal stool, such as food-derived DNA, might hinder the analysis of the bacteria in the stool. As such, we determined four relative values of the number of reads from the genes of V. cholerae: the ratio of DNA reads of V. cholerae to the total DNA; the ratio of the DNA reads of V. cholerae to the total bacterial DNA; the ratio of the RNA reads of V. cholerae to the total RNA; and the ratio of the RNA reads of V. cholerae to the total bacterial RNA. Furthermore, we determined the relative value of the number of reads from ctxA to the total DNA (Supplementary Tables S5 and S6). These ratios are also shown in Fig. 6a–d.Figure 6The ratio of DNA and RNA derived from V. cholerae in stool samples of the specific-pathogen-free patients. The stool samples from 22 diarrheal patients in which no etiological agent of diarrhea, including V. cholerae, was detected by our bacterial examination in the laboratory were analyzed in this examination. The extraction of DNA and RNA, and the preparation of reverse-transcribed DNA samples from the RNA samples were performed in the same manner as in Fig. 2. The origin of the reads obtained in this analysis was assigned by mapping to a database that included human and microorganism sequences. The obtained numbers of total reads, total bacterial reads, reads originating from V. cholerae, reads from ctxA in each sample are shown in Supplementary Tables S5 (the data from DNA) and S6 (the data from RNA). The percentages of reads of DNA from V. cholerae (blue bar in a) and of reads of DNA from ctxA (red bar in a) relative to the total DNA reads, and the percentages of reads of DNA from V. cholerae relative to the total bacterial DNA reads (b) are presented. Similarly, the results obtained from the RNA samples are presented in (c) and (d). The (c) and (d) show the percentages of reads of RNA from V. cholerae relative to the total RNA reads and to the reads of total bacterial RNA, respectively. The samples indicated by green circles are the samples of interest in this manuscript, as described in the text.Full size imageThe ratios of the number of reads derived from DNA of V. cholerae and the number of reads derived from ctxA to the number of reads of total DNA genes in these samples are shown by the blue and red bars in Fig. 6a, respectively. Reads from ctxA were detected in samples 1004, 1006, 1010, 1017 and 1018. This indicates that V. cholerae possessing ctx were alive in these samples; 1004, 1006, 1010, 1017 and 1018.The ratio of V. cholerae to total bacterial DNA in these samples was examined. The results are shown in Fig. 6b. The proportion of DNA of V. cholerae to total bacteria DNA in the stool of patients 1004, 1006, 1010, 1017, and 1018 is 28.633%, 0.234%, 73.068%, 2.282%, and 2.774%, respectively (Fig. 6b).In addition, the read of RNA from V. cholerae was examined. The ratio of the RNA to total RNA and to total bacterial RNA was calculated. RNA derived from V. cholerae was reliably detected in 4 of the 5 samples (1004, 1010, 1017, 1018). The ratio of the remaining one sample (1006) were low (Fig. 6c,d). However, it has been shown that the sample (1006) contains the read from DNA of ctxA (Supplementary Table S5). Therefore, we considered these five samples to be those containing toxigenic V. cholerae.As antibiotics were not administered to these five patients, the effects of antibacterial agents could be disregarded in our examination of the bacterial species in the stools. Among these 5 samples, the ratio of samples 1004 and 1010 examined in this examination was high and comparable to those of the samples of the cholera patients (Figs. 3 and 6). We considered that the diarrhea of the patients 1004 and 1010, might have been caused by the infection with V. cholerae O1.On the other hand, the samples of patients 1006, 1017 and 1018 did not show high values that could indicate that the diarrhea was caused by the infection with V. cholerae. It is probable that the diarrhea of these three patients (1006, 1017 and 1018) was caused by the actions of factors other than V. cholerae O1, and that a small number of V. cholerae inhabits the intestinal tract as a form of subclinical infection; this would explain why a gene derived from V. cholerae was detected by the metagenomic sequencing analysis. These results support the hypothesis that subclinical infections of V. cholerae occur in Kolkata. More