Philippa Kaur

1.Hoegh-Guldberg, O. et al. Coral reefs under rapid climate change and ocean acidification. Science 318, 1737–1742 (2007).ADS 
CAS 
Article 

Google Scholar 
2.Woodley, J. D. et al. Hurricane Allen’s impact on Jamaican coral reefs. Science 214, 749–755 (1981).ADS 
CAS 
Article 

Google Scholar 
3.Scoffin, T. P. The geological effects of hurricanes on coral reefs and the interpretation of storm deposits. Coral Reefs 12, 203–221 (1993).ADS 
Article 

Google Scholar 
4.Browning, T. N. et al. Widespread deposition in a Coastal Bay following three major 2017 Hurricanes (Irma, Jose, and Maria). Sci. Rep. 9, 7101. https://doi.org/10.1038/s41598-019-43062-4 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
5.Rogers, C. S. Immediate effects of hurricanes on a diverse coral/mangrove ecosystem in the US Virgin Islands and the potential for recovery. Diversity 11, 130 (2019).Article 

Google Scholar 
6.Graham, N. A. J. & Nash, K. L. The importance of structural complexity in coral reef ecosystems. Coral Reefs 32, 315–326 (2013).ADS 
Article 

Google Scholar 
7.Newman, S. P. et al. Reef flattening effects on total richness and species responses in the Caribbean. J. Anim Ecol. 84, 1678–1689 (2015).Article 

Google Scholar 
8.Burns, J. H. et al. Assessing the impact of acute disturbances on the structure and composition of a coral community using innovative 3D reconstruction techniques. Methods Oceanogr. 15, 49–59 (2016).Article 

Google Scholar 
9.Massel, S. R. & Done, T. J. Effects of cyclone waves on massive coral assemblages on the Great Barrier Reef: Meteorology, hydrodynamics and demography. Coral Reefs 12(3–4), 153–166 (1993).ADS 
Article 

Google Scholar 
10.Madin, J. S. & Connolly, S. R. Ecological consequences of major hydrodynamic disturbances on coral reefs. Nature 444(7118), 477–480 (2006).ADS 
CAS 
Article 

Google Scholar 
11.Madin, J. S., Baird, A. H., Dornelas, M. & Connolly, S. R. Mechanical vulnerability explains size-dependent mortality of reef corals. Ecol. Lett. 17(8), 1008–1015 (2014).Article 

Google Scholar 
12.Bellwood, D. R. & Hughes, T. P. Regional-scale assembly rules and biodiversity of coral reefs. Science 292, 1532–1535 (2001).ADS 
CAS 
Article 

Google Scholar 
13.Friedlander, A. et al. The state of coral reef ecosystems of the Northwestern Hawaiian Islands. The state of coral reef ecosystems of the United States and Pacific Freely Associated States. PLoS ONE 73, 263–306 (2005).
Google Scholar 
14.Houston, S. & Birchard, T. Central Pacific Hurricane Center Tropical Cyclone Report, Hurricane Walaka. https://www.nhc.noaa.gov/data/tcr/CP012018_Walaka.pdf (2020).15.Nugert, A. D. et al. Fire and rain: The legacy of Hurricane Lane in Hawai ‘i. Bull. Am. Meteorol. Soc. 101, 954–967 (2020).Article 

Google Scholar 
16.Fukunaga, A., Burns, J. H. R., Pascoe, K. H. & Kosaki, R. K. Associations between benthic cover and habitat complexity metrics obtained from 3D reconstruction of coral reefs at different resolutions. Remote Sens. 12, 1011 (2020).ADS 
Article 

Google Scholar 
17.Fukunaga, A. & Burns, J. H. R. Metrics of coral reef structural complexity extracted from 3D mesh models and digital elevation models. Remote Sens. 12, 2676 (2020).ADS 
Article 

Google Scholar 
18.Meyer, C. G., Holland, K. N. & Papastamatiou, Y. P. Seasonal and diel movements of giant trevally Caranx ignobilis at remote Hawaiian atolls: Implications for the design of marine protected areas. Mar. Ecol. Prog. Ser. 12(333), 13–25 (2007).ADS 
Article 

Google Scholar 
19.Fukunaga, A., Kosaki, R. K. & Hauk, B. B. Distribution and abundance of the introduced snapper Lutjanus kasmira (Forsskål, 1775) on shallow and mesophotic reefs of the Northwestern Hawaiian Islands. Bioinvasions Rec. 6, 259–268 (2017).Article 

Google Scholar 
20.Burns, J. H., Delparte, D., Gates, R. D. & Takabayashi, M. Integrating structure-from-motion photogrammetry with geospatial software as a novel technique for quantifying 3D ecological characteristics of coral reefs. PeerJ 3, e1077. https://doi.org/10.7717/peerj.1077 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
21.R Core Team. R: A Language and Environment for Statistical Computing, v. 3.5.3 (R Foundation for Statistical Computing, 2019). https://www.R-project.org (2019).22.Hijmans, R. J. raster: Geographic Data Analysis and Modeling, R package Version 2.9-5, https://CRAN.R-project.org/package=raster (2019)23.Bivand, R. & Rundel, C. rgeos: Interface to Geometry Engine—Open Source (GEOS), R package version 0.4-3. https://CRAN.R-project.org/package=rgeos (2019).24.Fukunaga, A., Burns, J. H. R., Craig, B. K. & Kosaki, R. K. Integrating three-dimensional benthic habitat characterization techniques into ecological monitoring of coral reefs. J. Mar. Sci. Eng. 7, 27 (2019).Article 

Google Scholar 
25.Risk, M. J. Fish diversity on a coral reef in the Virgin Islands. Atoll Res. Bull. 153, 1–4 (1972).Article 

Google Scholar 
26.Horn, B. K. P. Hill shading and the reflectance map. Proc. IEEE 69, 14–47 (1981).ADS 
Article 

Google Scholar 
27.Walbridge, S., Slocum, N., Pobuda, M. & Wright, D. J. Unified geomorphological analysis workflows with Benthic Terrain Modeler. Geosciences 8, 94 (2018).Article 

Google Scholar 
28.Sappington, J. M., Longshore, K. M. & Thompson, D. B. Quantifying landscape ruggedness for animal habitat analysis: A case study using bighorn sheep in the Mojave Desert. J. Wildl. Manage. 71, 1419–1426 (2007).Article 

Google Scholar 
29.Zevenbergen, L. W. & Thorne, C. R. Quantitative analysis of land surface topography. Earth Surf. Process Landf. 12, 47–56 (1987).ADS 
Article 

Google Scholar 
30.Couch, C. S. et al. Mass coral bleaching due to unprecedented marine heatwave in Papahānaumokuākea Marine National Monument (Northwestern Hawaiian Islands). PLoS ONE 12(9), e0185121. https://doi.org/10.1371/journal.pone.0185121 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
31.Magel, J. M. T., Burns, J. H., Gates, R. D. & Baum, J. K. Effects of bleaching-associated mass coral mortality on reef structural complexity across a gradient of local disturbance. Sci. Rep. 9, 2515. https://doi.org/10.1038/s41598-018-37713-1 (2019).ADS 
CAS 
Article 

Google Scholar 
32.Luckhurst, B. E. & Luckhurst, K. Analysis of the influence of substrate variables on coral reef fish communities. Mar. Biol. 49(4), 317–323 (1978).Article 

Google Scholar 
33.Graham, N. A. J. et al. Dynamic fragility of oceanic coral reef ecosystems. Proc. Natl. Acad. Sci. USA 103, 8425–8429 (2006).ADS 
CAS 
Article 

Google Scholar 
34.Mandelbrot, B. B. & Mandelbrot, B. B. The Fractal Geometry of Nature Vol. 1 (WH Freeman, 1982).MATH 

Google Scholar 
35.Young, G. C. et al. Cost and time-effective method for multi-scale measures of rugosity, fractal dimension, and vector dispersion from coral reef 3D models. PLoS ONE 12, e0175341. https://doi.org/10.1371/journal.pone.0175341 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar  More

Departamento de Ciencias Forestales, Universidad Nacional de Colombia Sede Medellín, Medellín, ColombiaAlvaro Duque, Miguel A. Peña & Sebastián González-CaroGrupo de Investigación en Biodiversidad, Medio Ambiente y Salud -BIOMAS – Universidad de Las Américas (UDLA), Quito, EcuadorFrancisco Cuesta, Marco Calderón-Loor & Esteban PintoDepartment of Plant and Microbial Biology, University of Minnesota, Saint Paul, MN, USAPeter KennedySchool of Geography, University of Leeds, Leeds, UKOliver L. PhillipsCentre for Integrative Ecology, School of Life and Environmental Sciences, Deakin University, Melbourne, VIC, AustraliaMarco Calderón-LoorInstituto de Ecología Regional (IER), Universidad Nacional de Tucumán (UNT) – Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Tucumán, ArgentinaCecilia Blundo, Julieta Carilla, Ricardo Grau, Agustina Malizia & Oriana Osinaga-AcostaHerbario Nacional de Bolivia (LPB), La Paz, BoliviaLeslie Cayola, Alfredo Fuentes & María I. Loza-RiveraMissouri Botanical Garden, St. Louis, MO, USALeslie Cayola, Alfredo Fuentes & María I. Loza-RiveraCenter for Conservation and Sustainable Development, Missouri Botanical Garden, St. Louis, MO, USAWilliam Farfán-Ríos, María I. Loza-Rivera & J. Sebastián TelloLiving Earth Collaborative, Washington University in Saint Louis, St. Louis, MO, USAWilliam Farfán-RíosPlant Ecology and Ecosystems Research, University of Gottingen, Gottingen, GermanyJürgen HomeierCentre of Biodiversity and Sustainable Land Use (CBL), University of Gottingen, Gottingen, GermanyJürgen HomeierEnvironmental Change Institute, School of Geography and the Environment, University of Oxford, Oxford, UKYadvinder MalhiFacultad de Ciencias Agrarias, Universidad Nacional de Jujuy, Jujuy, ArgentinaLucio MaliziaUniversité du Quebec a Montreal, Montreal, QC, CanadaJohanna A. Martínez-VillaDepartment of Biology, Washington University in St. Louis, St. Louis, MO, USAJonathan A. MyersConsorcio para el Desarrollo Sostenible de la Ecorregión Andina (CONDESAN), Quito, EcuadorManuel PeralvoColumbus State University, University System of Georgia, Columbus, GA, USAEsteban PintoCarbon Cycle and Ecosystems, Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USASassan SaatchiCenter for Energy, Environment and Sustainability, Winston-Salem, NC, USAMiles SilmanCentro Jambatú de Investigación y Conservación de Anfibios, Quito, EcuadorAndrea Terán-ValdezBiology Department, University of Miami, Coral Gables, FL, USAKenneth J. Feeley More

1.Díaz et al. Summary for Policymakers of the Global Assessment.pdf.2.Kehoe, L. et al. Biodiversity at risk under future cropland expansion and intensification. Nat. Ecol. Evolut. 1, 1129–1135 (2017).Article 

Google Scholar 
3.Hendershot, J. N. et al. Intensive farming drives long-term shifts in avian community composition. Nature 579, 393–396 (2020).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
4.Bommarco, R., Kleijn, D. & Potts, S. G. Ecological intensification: Harnessing ecosystem services for food security. Trends Ecol. Evol. 28, 230–238 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
5.Michael, D. R., Wood, J. T., O’Loughlin, T. & Lindenmayer, D. B. Influence of land sharing and land sparing strategies on patterns of vegetation and terrestrial vertebrate richness and occurrence in Australian endangered eucalypt woodlands. Agr. Ecosyst. Environ. 227, 24–32 (2016).Article 

Google Scholar 
6.Tittonell, P. Ecological intensification of agriculture—Sustainable by nature. Curr. Opin. Environ. Sustain. 8, 53–61 (2014).Article 

Google Scholar 
7.Willer, E. H., Schlatter, B., Trávní, J., Kemper, L. & Lernoud, J. The World of Organic Agriculture Statistics and Emerging Trends 2020. 337.8.Reganold, J. P. & Wachter, J. M. Organic agriculture in the twenty-first century. Nat. Plants 2 (2016).9.Connor, D. J. Organic agriculture cannot feed the world. Field Crop Res. 106, 187–190 (2008).Article 

Google Scholar 
10.Seufert, V. & Ramankutty, N. Many shades of gray—The context-dependent performance of organic agriculture. Sci. Adv. 3, e1602638 (2017).11.Smith, O. M. et al. Landscape context affects the sustainability of organic farming systems. Proc. Natl. Acad. Sci. 117, 2870–2878 (2020).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
12.Bengtsson, J., Ahnström, J. & Weibull, A.-C. The effects of organic agriculture on biodiversity and abundance: A meta-analysis: Organic agriculture, biodiversity and abundance. J. Appl. Ecol. 42, 261–269 (2005).Article 

Google Scholar 
13.Tuck, S. L. et al. Land-use intensity and the effects of organic farming on biodiversity: A hierarchical meta-analysis. J. Appl. Ecol. 51, 746–755 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
14.Lichtenberg, E. M. et al. A global synthesis of the effects of diversified farming systems on arthropod diversity within fields and across agricultural landscapes. Glob. Change Biol. 23, 4946–4957 (2017).ADS 
Article 

Google Scholar 
15.Lori, M., Symnaczik, S., Mäder, P., De Deyn, G. & Gattinger, A. Organic farming enhances soil microbial abundance and activity—A meta-analysis and meta-regression. PLOS ONE 12, e0180442 (2017).16.Kleijn, D., Rundlöf, M., Scheper, J., Smith, H. G. & Tscharntke, T. Does conservation on farmland contribute to halting the biodiversity decline?. Trends Ecol. Evol. 26, 474–481 (2011).PubMed 
Article 
PubMed Central 

Google Scholar 
17.Birkhofer, K., Ekroos, J., Corlett, E. B. & Smith, H. G. Winners and losers of organic cereal farming in animal communities across Central and Northern Europe. Biol. Cons. 175, 25–33 (2014).Article 

Google Scholar 
18.Mackie, K. A., Müller, T., Zikeli, S. & Kandeler, E. Long-term copper application in an organic vineyard modifies spatial distribution of soil micro-organisms. Soil Biol. Biochem. 65, 245–253 (2013).CAS 
Article 

Google Scholar 
19.Buchholz, J. et al. Soil biota in vineyards are more influenced by plants and soil quality than by tillage intensity or the surrounding landscape. Sci. Rep. 7 (2017).20.Hole, D. G. et al. Does organic farming benefit biodiversity?. Biol. Cons. 122, 113–130 (2005).Article 

Google Scholar 
21.Power, A. G. Ecosystem services and agriculture: Tradeoffs and synergies. Philos. Trans. R. Soc. B Biol. Sci. 365, 2959–2971 (2010).Article 

Google Scholar 
22.Peigné, J. et al. Earthworm populations under different tillage systems in organic farming. Soil Tillage Res. 104, 207–214 (2009).Article 

Google Scholar 
23.Biondi, A., Desneux, N., Siscaro, G. & Zappalà, L. Using organic-certified rather than synthetic pesticides may not be safer for biological control agents: Selectivity and side effects of 14 pesticides on the predator Orius laevigatus. Chemosphere 87, 803–812 (2012).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
24.Mehrabi, Z., Seufert, V., Ramankutty, N. The conventional versus alternative agricultural divide: A response to Garibaldi et al. Trends Ecol. Evolut. 32, 720–721 (2017).25.Tscharntke, T., Klein, A. M., Kruess, A., Steffan-Dewenter, I. & Thies, C. Landscape perspectives on agricultural intensification and biodiversity – ecosystem service management. Ecol. Lett. 8, 857–874 (2005).Article 

Google Scholar 
26.Gámez-Virués, S. et al. Landscape simplification filters species traits and drives biotic homogenization. Nat. Commun. 6 (2015).27.Holzschuh, A., Steffan-Dewenter, I. & Tscharntke, T. Agricultural landscapes with organic crops support higher pollinator diversity. Oikos 117, 354–361 (2008).Article 

Google Scholar 
28.Muneret, L., Auriol, A., Thiéry, D. & Rusch, A. Organic farming at local and landscape scales fosters biological pest control in vineyards. Ecol. Appl. 29, e01818 (2019).29.Gabriel, D. et al. Scale matters: The impact of organic farming on biodiversity at different spatial scales: Scale matters in organic farming. Ecol. Lett. 13, 858–869 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
30.Agreste. Pratiques Phytosanitaires en Viticulture. Campagne 2016. (2020)31.Agreste. La Viticulture Bio en Nouvelle-Aquitaine: Un Dynamisme à Tous les Stades de la Filière. (2020).32.Gruber, S. & Claupein, W. Effect of tillage intensity on weed infestation in organic farming. Soil Tillage Res. 105, 104–111 (2009).Article 

Google Scholar 
33.Pfingstmann, A. et al. Contrasting effects of tillage and landscape structure on spiders and springtails in vineyards. Sustainability 11, 2095 (2019).Article 

Google Scholar 
34.Dittmer, S. & Schrader, S. Longterm effects of soil compaction and tillage on Collembola and straw decomposition in arable soil. Pedobiologia 44, 527–538 (2000).Article 

Google Scholar 
35.Kolb, S., Uzman, D., Leyer, I., Reineke, A. & Entling, M. H. Differential effects of semi-natural habitats and organic management on spiders in viticultural landscapes. Agric. Ecosyst. Environ. 287, 106695 (2020).36.Birkhofer, K. et al. Relationships between multiple biodiversity components and ecosystem services along a landscape complexity gradient. Biol. Cons. 218, 247–253 (2018).Article 

Google Scholar 
37.Kratschmer, S. et al. Tillage intensity or landscape features: What matters most for wild bee diversity in vineyards?. Agric. Ecosyst. Environ. 266, 142–152 (2018).Article 

Google Scholar 
38.Ullmann, K. S., Meisner, M. H. & Williams, N. M. Impact of tillage on the crop pollinating, ground-nesting bee, Peponapis pruinosa in California. Agric. Ecosyst. Environ. 232, 240–246 (2016).Article 

Google Scholar 
39.Jiang, X., Wright, A. L., Wang, X. & Liang, F. Tillage-induced changes in fungal and bacterial biomass associated with soil aggregates: A long-term field study in a subtropical rice soil in China. Appl. Soil. Ecol. 48, 168–173 (2011).Article 

Google Scholar 
40.Zuber, S. M. & Villamil, M. B. Meta-analysis approach to assess effect of tillage on microbial biomass and enzyme activities. Soil Biol. Biochem. 97, 176–187 (2016).CAS 
Article 

Google Scholar 
41.Luff, M. L. The biology of the ground beetle Harpalus rufipes in a strawberry field in Northumberland. Ann. Appl. Biol. 94, 153–164 (1980).Article 

Google Scholar 
42.Shearin, A. F., Reberg-Horton, S. C. & Gallandt, E. R. Direct effects of tillage on the activity density of ground beetle (Coleoptera: Carabidae) weed seed predators. Environ. Entomol. 36, 1140–1146 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
43.Rundlöf, M. et al. Seed coating with a neonicotinoid insecticide negatively affects wild bees. Nature 521, 77–80 (2015).ADS 
PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
44.Martin, E. A. et al. The interplay of landscape composition and configuration: new pathways to manage functional biodiversity and agroecosystem services across Europe. Ecol. Lett. 22, 1083–1094 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
45.Goded, S., Ekroos, J., Azcárate, J. G., Guitián, J. A. & Smith, H. G. Effects of organic farming on plant and butterfly functional diversity in mosaic landscapes. Agric. Ecosyst. Environ. 284, 106600 (2019).46.Rusch, A., Valantin-Morison, M., Sarthou, J.-P. & Roger-Estrade, J. Multi-scale effects of landscape complexity and crop management on pollen beetle parasitism rate. Landsc. Ecol. 26, 473–486 (2011).Article 

Google Scholar 
47.Tamburini, G., De Simone, S., Sigura, M., Boscutti, F. & Marini, L. Conservation tillage mitigates the negative effect of landscape simplification on biological control. J. Appl. Ecol. 53, 233–241 (2016).Article 

Google Scholar 
48.Le Féon, V. et al. Intensification of agriculture, landscape composition and wild bee communities: A large scale study in four European countries. Agric. Ecosyst. Environ. 137, 143–150 (2010).Article 

Google Scholar 
49.Sousa, J. P. et al. Changes in Collembola richness and diversity along a gradient of land-use intensity: A pan European study. Pedobiologia 50, 147–156 (2006).Article 

Google Scholar 
50.Vanbergen, A. J. et al. Scale-specific correlations between habitat heterogeneity and soil fauna diversity along a landscape structure gradient. Oecologia 153, 713–725 (2007).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.Lehmitz, R., Russell, D., Hohberg, K., Christian, A. & Xylander, W. E. R. Active dispersal of oribatid mites into young soils. Appl. Soil. Ecol. 55, 10–19 (2012).Article 

Google Scholar 
52.Concepción, E. D., Díaz, M. & Baquero, R. A. Effects of landscape complexity on the ecological effectiveness of agri-environment schemes. Landsc. Ecol. 23, 135–148 (2008).Article 

Google Scholar 
53.Tscharntke, T. et al. Landscape moderation of biodiversity patterns and processes – Eight hypotheses. Biol. Rev. 87, 661–685 (2012).PubMed 
Article 
PubMed Central 

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

Google Scholar 
55.Naveed, M. et al. Simultaneous loss of soil biodiversity and functions along a copper contamination gradient: When soil goes to sleep. Soil Sci. Soc. Am. J. 78, 1239–1250 (2014).ADS 
Article 
CAS 

Google Scholar 
56.Eijsackers, H., Beneke, P., Maboeta, M., Louw, J. P. E. & Reinecke, A. J. The implications of copper fungicide usage in vineyards for earthworm activity and resulting sustainable soil quality. Ecotoxicol. Environ. Saf. 62, 99–111 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
57.Le Provost, G. et al. Land-use history impacts functional diversity across multiple trophic groups. Proc. Natl. Acad. Sci. 117, 1573–1579 (2020).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
58.Muneret, L. et al. Organic farming expansion drives natural enemy abundance but not diversity in vineyard-dominated landscapes. Ecol. Evol. https://doi.org/10.1002/ece3.5810 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
59.Lechenet, M., Dessaint, F., Py, G., Makowski, D. & Munier-Jolain, N. Reducing pesticide use while preserving crop productivity and profitability on arable farms. Nat. Plants 3 (2017).60.Le Féon, V. et al. Solitary bee abundance and species richness in dynamic agricultural landscapes. Agric. Ecosyst. Environ. 166, 94–101 (2013).Article 

Google Scholar 
61.McCravy, K. & Ruholl, J. Bee (Hymenoptera: Apoidea) diversity and sampling methodology in a midwestern USA deciduous forest. Insects 8, 81 (2017).PubMed Central 
Article 

Google Scholar 
62.Bano, R. & Roy, S. Extraction of Soil Microarthropods: A Low Cost Berlese-Tullgren Funnels Extractor. 4.63.Grueber, C. E., Nakagawa, S., Laws, R. J. & Jamieson, I. G. Multimodel inference in ecology and evolution: Challenges and solutions: Multimodel inference. J. Evol. Biol. 24, 699–711 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
64.Burnham, K. P. & Anderson, D. R. Multimodel inference: Understanding AIC and BIC in model selection. Sociol. Methods Res. 33, 261–304 (2004).MathSciNet 
Article 

Google Scholar 
65.Bartoń, K. MuMIn: Multi-Model Inference. R Package Version 1.43.17. https://CRAN.R-project.org/package=MuMIn (2020).66.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/ (2020).67.Hartig, F. DHARMa: Residual Diagnostics for Hierarchical (Multi-Level/Mixed) Regression Models (2020). R Package Version 0.3.3.0. https://CRAN.R-project.org/package=DHARMa68.Bates, D., Maechler, M., Bolker, B., Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01 (2015). More

1.Mutreja, A. et al. Evidence for several waves of global transmission in the seventh cholera pandemic. Nature 477, 462–465 (2011).CAS 
Article 

Google Scholar 
2.Ali, M., Nelson, A. R., Lopez, A. L. & Sack, D. A. Updated global burden of cholera in endemic countries. PLoS Negl. Trop. Dis. 9, e0003832 (2015).Article 

Google Scholar 
3.Bompangue, D. N. et al. Lakes as source of cholera outbreaks, Democratic Republic of Congo. Emerg. Infect. Dis. 14, 798–800 (2008).Article 

Google Scholar 
4.Weill, F. X. et al. Genomic history of the seventh pandemic of cholera in Africa. Science 358, 785–789 (2017).CAS 
Article 

Google Scholar 
5.Hounmanou, Y. M. G. et al. Genomic insights into Vibrio cholerae O1 responsible for cholera epidemics in Tanzania between 1993 and 2017. PLoS Neglect Trop. D. 13, e0007934 (2019).Article 

Google Scholar 
6.Colwell, R. R. Global climate and infectious disease: the cholera paradigm. Science 274, 2025–2031 (1996).CAS 
Article 

Google Scholar 
7.Singleton, F., Attwell, R., Jangi, M. & Colwell, R. Influence of salinity and organic nutrient concentration on survival and growth of Vibrio choleraein aquatic microcosms. Appl. Environ. Microbiol. 43, 1080–1085 (1982).CAS 
Article 

Google Scholar 
8.Kirschner, A. K. T. et al. Rapid growth of planktonic Vibrio cholerae non-O1/non-O139 strains in a large alkaline lake in Austria: Dependence on temperature and dissolved organic carbon quality. Appl. Environ. Microb. 74, 2004–2015 (2008).CAS 
Article 

Google Scholar 
9.Reid, P. C. et al. The continuous plankton recorder: concepts and history, from Plankton Indicator to undulating recorders. Prog. Oceanogr. 57, 117–173 (2003).Article 

Google Scholar 
10.Vezzulli, L. et al. Climate influence on Vibrio and associated human diseases during the past half-century in the coastal North Atlantic. Proc. Natl. Acad. Sci. USA. 113, E5062–E5071 (2016).CAS 
Article 

Google Scholar 
11.Huq, A. et al. Detection, isolation, and identification of Vibrio cholerae from the environment. Curr. Protoc. Microbiol. https://doi.org/10.1002/9780471729259.mc06a05s26 (2012).12.Thompson, F. L. et al. Phylogeny and molecular identification of vibrios on the basis of multilocus sequence analysis. Appl. Environ. Microb. 71, 5107–5115 (2005).CAS 
Article 

Google Scholar 
13.Vezzulli, L. et al. GbpA as a novel qPCR target for the species-specific detection of Vibrio cholerae O1, O139, Non-O1/Non-O139 in environmental, stool, and historical continuous plankton recorder samples. s. PLoS ONE 10, e0123983 (2015).Article 

Google Scholar 
14.Alam, M. et al. Toxigenic Vibrio cholerae in the aquatic environment of Mathbaria, Bangladesh. Appl. Environ. Microb. 72, 2849–2855 (2006).CAS 
Article 

Google Scholar 
15.Senoh, M. et al. Isolation of viable but nonculturable Vibrio cholerae O1 from environmental water samples in Kolkata, India, in a culturable state. Microbiol. Open 3, 239–246 (2014).CAS 
Article 

Google Scholar 
16.Vezzulli, L. et al. Whole-genome enrichment provides deep insights into vibrio cholerae metagenome from an African river. Microb. Ecol. 73, 734–738 (2017).CAS 
Article 

Google Scholar 
17.Kaboré, S. et al. Occurrence of Vibrio cholerae in water reservoirs of Burkina Faso. Res Microbiol 169, 1–10 (2018).Article 

Google Scholar 
18.Bwire, G. et al. Environmental surveillance of Vibrio cholerae O1/O139 in the five African Great Lakes and other major surface water sources in Uganda. Front. Microbiol. 9, 1560 (2018).Article 

Google Scholar 
19.Vezzulli, L., Baker-Austin, C., Kirschner, A., Pruzzo, C. & Martinez-Urtaza, J. Global emergence of environmental non-O1/O139 Vibrio cholerae infections linked with climate change: a neglected research field? Environ. Microbiol. 22, 4342–4355 (2020).CAS 
Article 

Google Scholar  More

The data needed for estimating the SR relationship consist of spawning biomass (S) and recruitment (R) observed over time. A lognormal distribution is frequently used as the distribution of errors for SR relationships13. We therefore assume that the residuals from a regression model having (r=log (R)) as a response variable and the logarithm of the latent SR relationship as the mean will have a normal distribution. In addition, we assume that the latent SR relationship is likely to be contaminated by some outliers given that fish populations often suffer from nonnegligible contamination, such as sporadic strong cohorts5.Figure 3Parameter estimates of the density-independent parameter (a), density-dependent parameter (b), and autocorrelation ((rho)) for the simulation using the HS SR function with autocorrelation (true (rho = 0.8)) in the residuals.Full size imageFigure 4Application of the robust SR model to fish population data from Japan. (Top) Estimates of ((b-min (S))/(max (S)-min (S))) using the LS and RSR methods. (Bottom) Examples of fitted SR curves using the LS (black line) and RSR (red line) methods (left, walleye pollock in the Sea of Japan; right, round herring in the Tsushima warm current).Full size imageA robust regression approachSuppose that the logarithm of recruitment ((r_t = log (R_t), (t = 1, ldots , T))) has the following autocorrelated normal distribution,$$begin{aligned} r_t = f(S_t|{varvec{theta }})+varepsilon _t, end{aligned}$$
(1)
where (varepsilon _t) is a scaled autoregressive error of order one, that is, (sqrt{lambda _t}(varepsilon _t-rho sqrt{lambda _{t-1}} varepsilon _{t-1})= e_t) with a gaussian noise (e_t) of mean zero and variance (sigma ^2), (S_t) is the spawning biomass, (f(S_t|{varvec{theta }})) is the logarithm of a density-dependent population growth model (spawner-recruitment (SR) curve), ({varvec{theta }}) is the parameter (vector) of the SR curve, (rho) is the autocorrelation, and (sigma ^2) is the base variance of the normal distribution. (lambda _t , (in (0,1])) is the weight for a datum in year t. Rearranging the equation for (varepsilon _t), we have (varepsilon _t sim N(rho sqrt{lambda _{t-1}} varepsilon _{t-1}, sigma ^2/lambda _t)) (Appendix A). We define (lambda _t) to be related to the magnitude of the residual (varepsilon _t),$$begin{aligned} lambda _t = exp left( – phi varepsilon _t^2 right) , end{aligned}$$where (phi , ( >0)) is the parameter that adjusts the influence of outliers. Given that the base variance (sigma ^2) is divided by (lambda _t), the variance is inflated when the difference between the datum and the SR curve is large. The model is equivalent to the AR(1) model when (lambda _t equiv 1) (i.e., (phi =0)) for any t. (sqrt{lambda _t}) is interpreted as the probability of the datum being generated from an uncontaminated normal distribution. When changing the (phi) parameter with (rho =0), the shapes of the probability density function and its derivative are similar to the Tukey’s biweight (also called bisquare) function14, which is close to the gaussian function near zero but decays swiftly as the datum becomes farther from zero (Fig. 1).By solving the equation at equilibrium, the mean deviance residual at (t=1) is zero and the variance at (t=1) is given by ({text{var}}({varepsilon_{1}} ){ = }{sigma ^{2}} {{/}}left[ {lambda _{1}} left( {1} – {rho ^{2}} {tilde{lambda }} right) right]), where ({tilde{lambda }}) is calculated by substituting the sample mean of (lambda _t), (tilde{lambda } = (1/T) sum _{t=1}^T lambda _t) (Appendix B). Incorporating the initial status, the log-likelihood function to be maximized is given by$$begin{aligned} log (L) = sum _{t=1}^T log left( N(r_t|f(S_t|{varvec{theta }})+delta _t, nu _t sigma ^2 lambda _t^{-1}) right) , end{aligned}$$
(2)
where (delta _{t} = 0) and (nu _{t} = (1-rho ^2 tilde{lambda })^{-1}) if (t = 1), and (delta _{t} = rho sqrt{lambda _{t-1}} varepsilon _{t-1}) and (nu _{t} = 1) if (t > 1). Because (varepsilon _{t-1}) increases and (lambda _{t-1}) decreases when there is an outlier at (t-1), the multiplication of (rho) and (sqrt{lambda _{t-1}}) mitigates the influence of an extreme outlier on autocorrelation and contributes to the restoration of the original autocorrelation.We need to estimate the parameters (sigma), (rho), and (phi) in addition to the SR relationship parameters ({varvec{theta }}). The parameter (phi) determines the mixing proportion of contamination and governs the predictive ability of the model. We use time series cross-validation15, which is also called retrospective forecasting16 (RF), to stably determine the value of (phi). First we delete the last datum. Then we use the SR relationship estimated from the data excluding the last datum to forecast recruitment and calculate its error assuming that the deleted recruitment for the last year is true. Next, we delete the two last data, forecast the second-to-last recruitment, and calculate the error assuming that the deleted second-to-last year’s recruitment is true. After the procedure is repeated on a rolling basis, the (phi) parameter having the smallest average error is finally selected. The optimum (phi) is determined by minimizing the following RF error:$$begin{aligned} RF_R = exp left( frac{1}{P} sum _{t=1}^P log left[ left( r_{T-(t-1)} -hat{r}_{T-(t-1)}^{1:(T-t)} right) ^2 right] right) . end{aligned}$$
(3)
This is the geometric mean of predicted errors, which stabilizes the performance of retrospective forecasting. (r_{T-(t-1)}) is the logarithm of observed recruitment in year (T-(t-1)) and (hat{r}_{T-(t-1)}^{1:(T-t)}) is the predicted value estimated using the data from years 1 to (T-t), which is given by$$begin{aligned} hat{r}_{T-(t-1)}^{1:(T-t)} = f(S_{T-(t-1)}|hat{varvec{theta }})+hat{rho } sqrt{hat{lambda }_{T-t}} hat{varepsilon }_{T-t}, end{aligned}$$where (t = 1, ldots , P). We adopt (P=10) for stable estimation in this paper, though we commonly take 5 as the minimum P17.All subsequent analyses are performed using R18 and its package TMB19 (Template Model Builder).SimulationWe generate the simulated data ((left{ (R_t, S_t) ; t = 1, ldots , T right})) with some outliers and autocorrelated errors and test the performance of our robust SR (RSR) method in comparison with the LS and LAD methods. LAD was chosen because it is a typical robust method and is generally superior to the least median squares method used in Chen & Paloheimo (1995)11. The average recruitment data are generated from the Hockey–Stick (HS) SR function12, (f(S_t|{varvec{theta }}) = log left( a min (S_t, b) right)), where ({varvec{theta }} = (a, b) = (1.2, 500)). Stochastic normal errors are added to the log recruitment data with or without autocorrelation. When there is an autocorrelation in the residuals of log recruitment, the autocorrelation is set to (rho = 0.8). To examine the effect of outliers, we add the outliers that occur at the expected frequency of twice per 10 years ((p=0.2)) to the residuals of log recruitment. The patterns of outlier occurrence are threefold: evenly occurring positive and negative outliers ((q=0.5)), all positive outliers ((q=1.0)), and all negative outliers ((q=0.0)) (see Appendix C for the definition of q). We then have eight types of simulated data (no outliers, positive and negative outliers, all positive outliers, and all negative outliers for autocorrelation in the normal residual (rho = 0) and (rho =0.8), respectively). The simulations are replicated 1,000 times for each of the eight types. The length of each SR data time series (T) is set to 30 years which is typical for SR time series data9,12. The performance of the methods is evaluated by two indicators that represent long-term and short-term predictive abilities ((hat{R}_0 – R_0)/R_0) and ((hat{R}_{T+1} – R_{T+1})/R_{T+1}), respectively, where the former is the asymptotic maximum recruitment ((R_0 = ab) for the HS SR function) and the latter is recruitment in the ensuing year (T+1), which is given by (R_{T+1} = exp (f(S_{T+1}|{varvec{theta }}) + rho omega _{T} + eta _{T+1})), where (omega _T) and (eta _{T+1}) are independent gaussian noises (Appendix C). Note that the true recruitment at (T+1) does not include any outliers. The mathematical details of the simulation are given in Appendix C. Autocorrelation is always estimated such that (rho) is set to zero when an estimate of (rho) is equal to or less than zero because a negative autocorrelation is usually impractical20. The parameter (log (phi )) in RSR is chosen from the grid values from (-3.0) to 3.0 in increments of 0.5. The best (phi) is a minimizer of the RF error (RF_R) (Eq. 3).For sensitivity tests, we conduct the following additional simulations: (S1) same as the above base case scenario (S0) except that (a = 1.8); (S2) same as S0 except that (p = 0.1) (the expected frequency of outliers is once every 10 years) in place of (p=0.2); (S3) same as S0 except that (p = 0.3) (the expected frequency of outliers is three times every 10 years) in place of (p=0.2); (S4) same as S0 except that (f(S_t|{varvec{theta }})) is the logarithm of the Beverton–Holt function; (S5) same as S0 except that (f(S_t|{varvec{theta }})) is the logarithm of the Ricker function; S6) same as S0 except for the spawner-abundance dependent p, in which the expected frequency of outliers is higher for lower spawner abundances than for higher spawner abundances.Finally, we calculate biological reference points related to maximum sustainable yield (MSY), i.e., fishing rate at MSY ((F_{rm {msy}})) and spawning biomass at MSY ((S_{rm {msy}})), for each scenario and evaluate their relative biases. To calculate (F_{rm {msy}}) and (S_{rm {msy}}), we require additional information on survival and growth as well as an assumption about population dynamics. For simplicity, we use the delay-difference model as the population dynamics model5. The mathematical details are given in Appendix D.Real data analysisIchinokawa, Okamura & Kurota (2017) fitted the SR curves to fish population data from Japan which comprise 26 SR datasets (Appendix E), demonstrating that some populations showed strong density dependence but others had weak or low density dependence. We fit the HS SR curves to the same 26 SR datasets used in Ichinokawa, Okamura & Kurota (2017). Because Ichinokawa, Okamura & Kurota (2017) used LS as the fitting method, we use LS and RSR to compare the density-independent parameter (log (hat{a})), standardized density-dependent parameter (( hat{b}-min (S) )/( max (S) – min (S) )), autocorrelation in the residuals (hat{rho }), and predictability (hat{RF}_R) in the HS SR curves. More

1.Climate Change: Vital Signs of the Planet (NASA, 2021); https://climate.nasa.gov/vital-signs/carbon-dioxide2.Kimball, B. A. Crop responses to elevated CO2 and interactions with H2O, N, and temperature. Curr. Opin. Plant Biol. 31, 36–43 (2016).CAS 
Article 

Google Scholar 
3.Leakey, A. D. B. et al. Elevated CO2 effects on plant carbon, nitrogen, and water relations: six important lessons from FACE. J. Exp. Bot. 60, 2859–2876 (2009).CAS 
Article 

Google Scholar 
4.Lee, T. D., Barrott, S. H. & Reich, P. B. Photosynthetic responses of 13 grassland species across 11 years of free-air CO2 enrichment is modest, consistent and independent of N supply. Glob. Change Biol. 17, 2893–2904 (2011).Article 

Google Scholar 
5.Ainsworth, E. A. & Long, S. P. 30 years of free-air carbon dioxide enrichment (FACE): what have we learned about future crop productivity and its potential for adaptation? Glob. Change Biol. 27, 27–49 (2021).Article 

Google Scholar 
6.Dusenge, M. E., Duarte, A. G. & Way, D. A. Plant carbon metabolism and climate change: elevated CO2 and temperature impacts on photosynthesis, photorespiration and respiration. New Phytol. 221, 32–49 (2019).CAS 
Article 

Google Scholar 
7.Sardans, J. et al. Ecometabolomics for a better understanding of plant responses and acclimation to abiotic factors linked to global change. Metabolites 10, 239 (2020).CAS 
Article 

Google Scholar 
8.Poorter, H. & Navas, M.-L. Plant growth and competition at elevated CO2: on winners, losers and functional groups. New Phytol. 157, 175–198 (2003).Article 

Google Scholar 
9.Parmesan, C. & Hanley, M. E. Plants and climate change: complexities and surprises. Ann. Bot. 116, 849–864 (2015).Article 

Google Scholar 
10.Ode, P. J., Johnson, S. N. & Moore, B. D. Atmospheric change and induced plant secondary metabolites—are we reshaping the building blocks of multi-trophic interactions? Curr. Opin. Insect Sci. 5, 57–65 (2014).Article 

Google Scholar 
11.Robinson, E. A., Ryan, G. D. & Newman, J. A. A meta-analytical review of the effects of elevated CO2 on plant–arthropod interactions highlights the importance of interacting environmental and biological variables. New Phytol. 194, 321–336 (2012).CAS 
Article 

Google Scholar 
12.Willeit, M., Ganopolski, A., Calov, R. & Brovkin, V. Mid-Pleistocene transition in glacial cycles explained by declining CO2 and regolith removal. Sci. Adv. 5, eaav7337 (2019).CAS 
Article 

Google Scholar 
13.Busch, F. A. & Sage, R. F. The sensitivity of photosynthesis to O2 and CO2 concentration identifies strong Rubisco control above the thermal optimum. New Phytol. 213, 1036–1051 (2017).CAS 
Article 

Google Scholar 
14.Drake, B. G., Gonzàlez-Meler, M. A. & Long, S. P. More efficient plants: a consequence of rising atmospheric CO2? Annu. Rev. Plant Physiol. Plant Mol. Biol. 48, 609–639 (1997).CAS 
Article 

Google Scholar 
15.Ziska, L. H., Sicher, R. C., George, K. & Mohan, J. E. Rising atmospheric carbon dioxide and potential impacts on the growth and toxicity of poison ivy (Toxicodendron radicans). Weed Sci. 55, 288–292 (2007).CAS 
Article 

Google Scholar 
16.Ziska, L. H. & Caulfield, F. A. Rising CO2 and pollen production of common ragweed (Ambrosia artemisiifolia L.), a known allergy-inducing species: implications for public health. Funct. Plant Biol. 27, 893 (2000).Article 

Google Scholar 
17.Ziska, L. H., Panicker, S. & Wojno, H. L. Recent and projected increases in atmospheric carbon dioxide and the potential impacts on growth and alkaloid production in wild poppy (Papaver setigerum DC.). Clim. Change 91, 395 (2008).CAS 
Article 

Google Scholar 
18.Del Fabbro, C. & Prati, D. The relative importance of immediate allelopathy and allelopathic legacy in invasive plant species. Basic Appl. Ecol. 16, 28–35 (2015).Article 

Google Scholar 
19.Ni, G. et al. Exploring the novel weapons hypothesis with invasive plant species in China. Allelopath. J. 29, 199–214 (2012).
Google Scholar 
20.Peñuelas, J. et al. Higher allocation to low cost chemical defenses in Iinvasive species of Hawaii. J. Chem. Ecol. 36, 1255–1270 (2010).Article 

Google Scholar 
21.Bajwa, A. A., McClay, A. & Adkins, S. W. in Parthenium Weed: Biology, Ecology and Management (eds Adkins, S., Shabbir, A. et al.) 7–39 (CABI, 2019).22.Adkins, S. & Shabbir, A. Biology, ecology and management of the invasive parthenium weed (Parthenium hysterophorus L.): management of parthenium weed. Pest Manag. Sci. 70, 1023–1029 (2014).CAS 
Article 

Google Scholar 
23.Niranjan, A. et al. Identification and quantification of heterologous compounds parthenin and organic acids in Parthenium hysterophorus L. using HPLC-PDA-MS-MS. Anal. Lett. 46, 48–59 (2013).Article 

Google Scholar 
24.Belz, R. G., van der Laan, M., Reinhardt, C. F. & Hurle, K. Soil degradation of parthenin—does it contradict the role of allelopathy in the invasive weed Parthenium hysterophorus L.? J. Chem. Ecol. 35, 1137–1150 (2009).CAS 
Article 

Google Scholar 
25.Hanif, Z., Adkins, S. W., Prentis, P. J., Navie, S. C. & O’Donnell, C. J. Characterization of the reproductive behaviour and invasive potential of parthenium weed in Australia. Pak. J. Weed Sci. Res. 18, 767–774 (2012).
Google Scholar 
26.Bajwa, A. A., Chauhan, B. S. & Adkins, S. Morphological, physiological and biochemical responses of two Australian biotypes of Parthenium hysterophorus to different soil moisture regimes. Environ. Sci. Pollut. Res. 24, 16186–16194 (2017).CAS 
Article 

Google Scholar 
27.Nguyen, T., Bajwa, A. A., Navie, S., O’Donnell, C. & Adkins, S. Parthenium weed (Parthenium hysterophorus L.) and climate change: the effect of CO2 concentration, temperature, and water deficit on growth and reproduction of two biotypes. Environ. Sci. Pollut. Res. 24, 10727–10739 (2017).Article 

Google Scholar 
28.Chadwick, M., Trewin, H., Gawthrop, F. & Wagstaff, C. Sesquiterpenoids lactones: benefits to plants and people. Int. J. Mol. Sci. 14, 12780–12805 (2013).Article 

Google Scholar 
29.Ojija, F., Arnold, S. E. J. & Treydte, A. C. Impacts of alien invasive Parthenium hysterophorus on flower visitation by insects to co-flowering plants. Arthropod Plant Interact. 13, 719–734 (2019).Article 

Google Scholar 
30.Bajwa, A. A., Chauhan, B. S. & Adkins, S. W. Germination ecology of two Australian biotypes of ragweed parthenium (Parthenium hysterophorus) relates to their invasiveness. Weed Sci. 66, 62–70 (2018).Article 

Google Scholar 
31.Bajwa, A. A. et al. Toxic potential and metabolic profiling of two Australian biotypes of the invasive plant parthenium weed (Parthenium hysterophorus L.). Toxins 12, 447 (2020).CAS 
Article 

Google Scholar 
32.Grime, J. P. Plant Strategies, Vegetation Processes, and Ecosystem Properties (Wiley, 2001).33.Grime, J. P. in Plant Evolutionary Biology (eds Gottlieb, L. D. & Jain, S. K.) 371–393 (Springer, 1988).34.Craine, J. M. Reconciling plant strategy theories of Grime and Tilman. J. Ecol. 93, 1041–1052 (2005).Article 

Google Scholar 
35.Bae, J. et al. Effect of elevated atmospheric carbon dioxide on the allelopathic potential of common ragweed. J. Ecol. Environ. 43, 21 (2019).Article 

Google Scholar 
36.Wang, R.-L. et al. Responses of Mikania micrantha, an invasive weed to elevated CO2: induction of β-caryophyllene synthase, changes in emission capability and allelopathic potential of β-caryophyllene. J. Chem. Ecol. 36, 1076–1082 (2010).CAS 
Article 

Google Scholar 
37.Robinson, J. M. Photosynthetic carbon metabolism in leaves and isolated chloroplasts from spinach plants grown under short and intermediate photosynthetic periods. Plant Physiol. 75, 397–409 (1984).CAS 
Article 

Google Scholar 
38.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2019).39.Filion, M., Dutilleul, P. & Potvin, C. Optimum experimental design for Free-Air Carbon dioxide Enrichment (FACE) studies. Glob. Change Biol. 6, 843–854 (2000).Article 

Google Scholar 
40.Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. https://www.jstatsoft.org/article/view/v067i01 (2015).41.Searle, S. R., Speed, F. M. & Milliken, G. A. Population marginal means in the linear model: an alternative to least squares means. Am. Stat. 34, 216–221 (1980).
Google Scholar  More

Our findings illustrate both the range and extent of wildlife exploitation in Wuhan markets, prior to new trading bans linked to the COVID-19 outbreak, along with the poor conditions under which these animals were kept prior to sale. Circumstantially, the absence of pangolins (and bats, not typically eaten in Central China; media footage generally depicts Indonesia) from our comprehensive survey data corroborates that pangolins are unlikely implicated as spill-over hosts in the COVID-19 outbreak. This is unsurprising because live pangolin trading has largely ceased in China13.We should therefore not be complacent, because the original source of COVID-19 does not seem to have been established. This is doubly important because false attribution can lead to extreme and irresponsible animal persecution. For instance, civets were killed en masse following the SARS-CoV outbreak5, and any unwarranted vilification or persecution of pangolins and bats in relation to COVID-19 would risk undermining otherwise very successful efforts to better protect and conserve wildlife in China.Regarding our insights into broader IWT issues in Wuhan, the animals sold were relatively expensive, representing luxury food items, not cheap bushmeat (Table 1). We thus make an ethical distinction here between the subsistence consumption of bush meat in poorer nations, versus the sort of cachet attached to wild animal consumption in parts of the developed world, notably China14, but also Japan15. While c. 30% of mammals were clearly wild-caught, indicated by trapping and shooting wounds, the captive breeding of other species is commonplace in China. Raccoon dog fur farming is legal in China; however, due to a drop in fur prices, raccoon dogs are now frequently sold off in live animal markets, augmented by wild-caught individuals. Similarly, all American mink (Neovison vison) originated from fur farms—noting that SARS-CoV-2 has been reported in mink farms in Europe and North America16, 17. In contrast, the captive breeding and sale of Siberian weasels (Mustela sibirica), is totally illegal in China, yet they are easy to breed, and sold openly, without attracting law enforcement. Indeed, prior to COVID-19 reforms, although enforcement officers from the Wuhan Forestry Bureau issued permits to market vendors, they were broadly disinterested in what species were sold. Furthermore, although animals were required to have an origin certificate and be quarantined to ensure they did not exhibit overt disease symptoms, no clear policy was enforced on these conditions. This is important because the species that were traded are capable of hosting a wide range of infectious zoonotic diseases or disease-baring parasites (Supplementary Table S1), aside from COVID-19. These range from potentially lethal viruses, for example, rabies, SFTS, H5N1, to common bacterial infections that, nevertheless, represent a risk to human health (e.g., Streptococcus). Indeed, globally, wildlife is thought to be the source of at least 70% of all emerging diseases18.Legislative reform is also vital to clarify unequivocally which species are considered ‘wild’ and cannot be traded legally and safely. Another problem, as encountered by the WHO report is that, retrospectively, it proved difficult to ascertain which species were on sale, even to the genus level, relying solely on the responsible market authority’s official sales records and disclosures1. As we19, 20, and others21, have proposed previously, China’s LFSSP and LESS must be updated to apply proper binomials, and to align with recent taxonomic revisions; for instance, cobra snakes (Nada atra) can be farmed legally for food with permits, but wild caught species, such as water snakes and wolf snakes were also sold in Wuhan, labelled simply as ‘snakes’. Such an application of clear species names would allow for more effective prosecutions19. Furthermore, the WHO reports that market authorities claimed all live and frozen animals sold in the Huanan market were acquired from farms officially licensed for breeding and quarantine, and as such no illegal wildlife trade was identified1. In reality, however, because China has no regulatory authority regulating animal trading conducted by small-scale vendors or individuals it is impossible to make this determination1, 21. Similar discrepancies concerning species identification and origins afflict investigations around the world22.Another important animal trade that requires attention, outside of exploitation as food, is the supply of pets, like the squirrels and crested myna birds sold in Wuhan’s market. Our previous research found annual trade volumes equivalent to c. 17,000 parrots and c. 160,000 turtles (many turtles being invasive if escaping to the wild) sold online as pets via Taobao.com between 2016–2017, in contravention of China’s WACL and/or the Animal Epidemic Prevention Law23,24,25. While not currently the vector of any major viral epidemics, it would be naive to imagine that unconventional pets do not still also pose a serious concern for public health26. This potential for disease is likely exacerbated by poor sanitary and welfare conditions (Fig. 2). More

Sample collection and incubationThree replicates of soil samples were collected from the top 10 cm in of plant-free patches in four ecosystems along the C. Hart Merriam elevation gradient in Northern Arizona25 beginning at high desert grassland (1760 m), and followed at higher elevations by piñon-pine juniper woodland (2020 m), ponderosa pine forest (2344 m), and mixed conifer forest (2620 m). Soils were air-dried for 24 h at room temperature, homogenized, and passed through a 2 mm sieve before being stored at 4 °C for another 24 h. Soil incubations were performed on soils with mass of 20 g of dry soil for measurements of CO2 and microbial biomass carbon (MBC), while 2 g of dry soil aliquots were incubated separately (but under equivalent conditions) for quantitative stable isotope probing (qSIP). We applied three treatments to these soils through the addition of water (up to 70% water-holding capacity): water alone (control), with glucose (C treatment; 1000 µg C g−1 dry soil), or with glucose and nitrogen (C + N treatment; [NH4]2SO4 at 100 µg N g−1 dry soil). All samples for qSIP were incubated with 18O-enriched water (97 atom%) and matching controls necessary to calculate the change in 18O enrichment across the microbial community. We applied water at natural abundance (i.e., no 18O-enriched water) to the larger soil samples prepared for measurement of carbon flux. All soils were incubated in the dark for one week. Following incubation, soils were frozen at −80 °C for 1 week prior to DNA extraction.Soil, CO2, and microbial biomass measurementsWe analyzed headspace gas of soils for CO2 concentration and δ13CO2 three times during the week-long incubation using a LI-Cor 6262 (LI-Cor Biosciences Inc. Lincoln, NE, USA) and a Picarro G2201 (Picarro Inc., Sunnyvale, CA, USA), respectively. Prior to incubation we analyzed soil MBC using the chloroform-fumigation extraction method on 10 g of soil. One sub-sample was immediately extracted with 25 ml of a 0.05 M K2SO4 solution, while a second sub-sample was first fumigated with chloroform (for 5 days), after which it was similarly extracted. Following K2SO4 addition, we agitated soils for 1 h, filtered the extract through a Whatman #3 filter paper, and dried the filtered solution (60 °C, 4 days). Salts with extracted C were ground and analyzed for total C using an elemental analyzer coupled to a mass spectrometer. MBC was calculated as the difference between the fumigated and immediately extracted samples’ soil C using an extraction efficiency of 0.45 (as per Liu et al.26).Quantitative stable isotope probingWe performed DNA extraction and 16S amplicon sequencing on 18O-incubated qSIP soils11,12,13. The procedures targeted the V4 region of the 16S gene as specified by the Earth Microbiome Project (EMP, http://www.earthmicrobiome.org) standard protocols27,28. We used PowerSoil DNA extraction kits following manufacture instructions to isolate DNA from soil (MoBio laboratories, Carlsbad, CA, USA). We quantified extracted DNA using the Qubit dsDNA High-Sensitivity assay kit and a Qubit 2.0 Fluorometer (Invitrogen, Eugene, OR, USA). To quantify the degree of 18O isotope incorporation into bacterial DNA, we performed density fractionation and sequenced 15–18 fractions separately following methods modified from the canonical publication7. We added 1 µg of DNA to 2.6 mL of saturated CsCl solution in combination with a gradient buffer (200 mM Tris, 200 mM KCL, 2 mM EDTA) in a 3.3 mL OptiSeal ultracentrifuge tube (Beckman Coulter, Fullerton, CA, USA). The solution was centrifuged to produce a gradient of increasingly labeled (heavier) DNA in an Optima Max bench top ultracentrifuge (Beckman Coulter, Brea, CA, USA) with a Beckman TLN-100 rotor (127,000 × g for 72 h) at 18 °C. We separated each sample from the continuous gradient into approximately 20 fractions (150 µL) using a modified fraction recovery system (Beckman Coulter). We then measured the density of each separate fraction with a Reichart AR200 digital refractometer (Reichert Analytical Instruments, Depew, NY, USA) and retained fractions with densities between 1.640 and 1.735 g cm−3. We cleaned and purified DNA in these fractions using isopropanol precipitation, quantified DNA using the Quant-IT PicoGreen dsDNA assay (Invitrogen) and a BioTek Synergy HT plate reader (BioTek Instruments Inc., Winooski, VT, USA), and quantified bacterial 16S gene copies using qPCR (primers: Supplementary Table 1) in triplicate. We used 8 µL reactions consisting of 0.2 mM of each primer, 0.01 U µL−1 Phusion HotStart II Polymerase (Thermo Fisher Scientific, Waltham, MA), 1× Phusion HF buffer (Thermo Fisher Scientific), 3.0 mM MgCl2, 6% glycerol, and 200 µL of dNTPs. We amplified DNA using a Bio-Rad CFX384 Touch real-time PCR detection system (Bio-Rad, Hercules, CA, USA) with the following cycling conditions: 95 °C at 1 min and 44 cycles of 95 °C (30 s), 64.5 °C (30 s), and 72 °C (1 min).We sequenced the 16S V4 region (primers: EMP standard 515F—806R; Supplementary Table 1) on an Illumina MiSeq (Illumina, Inc., San Diego, CA, USA). Sequences were amplified using the same reaction mix as qPCR amplification but cycling at 95 °C for 2 min followed by 15 cycles of 95 °C (30 s), 55 °C (30 s), and 60 °C (4 min). In addition to post-incubation soils, we extracted, amplified, and sequenced DNA of the bacterial community at the start of the incubation.Sequence processing and qSIP analysisThe raw sequence data of forward and reverse reads (FASTQ) were processed within the QIIME 2 environment (release 2018.6)29,30, denoising sequences with the available DADA2 pipeline31. We clustered the remaining sequences into amplicon sequence variants or ASVs (at 100% sequence identity) against the SILVA 132 database32 using an open-reference Naïve Bayes feature classifier33. We removed global singletons and doubleton ASVs, non-bacterial lineages, and samples with less than 4000 sequence reads. Removal of global singletons and doubletons resulted in the removal of 2241 unique ASVs from the feature table yielding 115,647 out of 117,888 (a retention of 98% of all ASVs) as well as the loss of 4018 sequences leaving 37,765,678 (a retention >99% of all sequences). We combined taxonomic information and ASV sequence counts with per-fraction qPCR and density measurements using the phyloseq package (version 1.24.2), in R (version 3.5.1)34. Because high-throughput sequencing produces relativized measures of abundance, we converted ASV sequencing abundances in each fraction to the number of 16S rRNA gene copies per g dry soil based on the known amount of dry soil added and the amount of DNA in each soil sample. All data and analytical code have been made publicly accessible35.To perform qSIP analysis and calculate per-capita growth rates of each ASV, we used our in-house qsip package (https://github.com/bramstone/qsip) based on previously published research7,10. Because rare and infrequent taxa are more likely to be lost in samples with poor sequencing depth with their absences affecting DNA density changes, we invoked a presence or absence-based filtering criteria on ASVs prior to calculation of per-capita growth rates. Within each ecosystem, we kept only ASVs that appeared in two of the three replicates of a treatment (18O, C, and C + N) and at that appeared in at least five of the fractions within each of those two replicates. ASVs filtered out of one treatment were allowed to appear in another if they met the frequency threshold.For all remaining ASVs (1081 representing less than 1% of all ASVs but 58% of all sequence reads), we calculated per-capita gross growth (i.e., cell division) rates observed in each replicate using an exponential growth model10. We applied these per-capita rates to the number of 16S rRNA gene copies to estimate the production of new 16S rRNA gene copies of each ASV per g dry soil per week using the following equation:$$frac{{rm{d}}{N}_{{rm{i}}}}{{{rm{d}}t}}={N}_{{rm{i,t}}}-{N}_{{rm{i,t}}}{e}^{-{g}_{{rm{i}}}t},$$
(1)
Where Ni,t is the number of 16S rRNA gene copies of taxon i at time t (here after 7 days) and gi represents the per-capita growth rate (calculated as a daily rate). See Supplementary Fig. 3 for results on the production of 16S gene copies.Calculation of 16S rRNA gene copy numbers and cell massIn parallel to taxonomic assignment, we compared quality-filtered 16S sequences against a database of 12,415 complete prokaryote genomes obtained from GenBank. From these genomes, we extracted data on 16S rRNA gene copy number, total genome size, and 16S gene sequence. We used BLAST to find matches against this database to the ASVs generated from QIIME 2 to make per-taxon assignments of 16S rRNA gene copy number and total genome size13. For ASVs that did not find an exact match, we assigned 16S rRNA gene copy number values and genome sizes based on the median values observed in the most specific possible taxonomic rank. We estimated the mass of individual cells for each population using published allometric scaling relationships between genome length and cellular mass from West and Brown:36$${{{log }}}_{10}({M}_{{rm{i}}})=frac{{{{log }}}_{10}left({G}_{{rm{i}}}right)-9.4}{0.24},$$
(2)
where Mi indicates cellular mass (g) and Gi indicates genome length (bp) for taxon i. We obtained this relationship by digitizing Fig. 436 using DataThief III and re-fitting the trend line in log–log space. We estimated that 20% of the cellular mass was carbon37. To validate this approach, cellular mass estimates and initial 16S copy number measurements were used to estimate population-level biomass C values which were summed and compared to initial community-level MBC. We found that these values overestimated initial MBC by an order of magnitude. As such, cellular carbon mass was divided by 10 in our final calculations. We applied cellular mass and 16S copy number estimates to the production of 16S copies to estimate the production of biomass carbon for each taxon during the incubation period (t):$${P}_{{rm{i}}}=frac{{rm{d}}{N}_{{rm{i}}}/{{rm{d}}t}}{C_{{rm{i}}}}cdot {M}_{{rm{i}}}cdot 0.2,$$
(3)
where Pi indicates production of biomass carbon (µg C g dry soil−1 week−1) and Ci indicates 16S copy number per cell for taxon i. The 0.2 coefficient represents an estimate that 20% of cellular mass is composed of carbon.Efficiency and respiration modelingWe estimated rates of respiration using qSIP-informed growth rates and community-level carbon use efficiency (CUE). CUE estimates were based on the incorporation of 18O-water into DNA as a measure of gross biomass production38,39 and measured CO2 in headspace gas from soil incubations. We calculated the production of 18O-labeled biomass carbon (18P) at the community-level for each sample by summing the products of per-taxon 18O enrichment (excess atom fraction, EAF) and relative abundance:$${, }^{18}{P}=mathop{sum }limits_{i=1}^{n}({,}^{18}{{{rm{EAF}}}}_{{rm{i}}}cdot {y}_{{rm{i}}})cdot {rm{DN}}{rm{A}}_{0}cdot fleft({{rm{MB}}}{rm{C}}_{0} sim {rm{DN}}{rm{A}}_{0}right),$$
(4)
where 18P indicates the gross production of 18O-labeled microbial biomass carbon per gram of dry soil per week, 18EAFi indicates the enrichment of DNA of taxon i and yi indicates its relative abundance, DNA0 indicates the concentration of DNA per gram of dry soil prior to incubation, and MBC0 indicates the microbial biomass carbon per gram of dry soil prior to incubation. Here, the MBC0 ~ DNA0 function indicates the linear relationship between MBC and DNA concentration. We used the output from Eq. 4 to calculate community CUE for each sample:$${{rm{CUE}}}=frac{{,}^{18}{{P}}}{(!{,}^{18}P+R)},$$
(5)
where R indicates the total CO2 respired per gram dry soil per week.We used the community CUE values from each sample (Eq. 5) to constrain/as upper and lower limits our estimates of per-taxon CUE. For a group of three replicates from a given ecosystem and treatment, we used the minimum and maximum observed community-level CUE values as the acceptable range of per-taxon CUE values. These constraints were used to control the shape of the function of per-taxon CUE and growth rate, though functions were modeled both with and without constraints (i.e., per-taxon CUE values were bounded only by 0 and 0.7). The range of community-level CUE values for each treatment were 0.18–0.53 for control soils, 0.04–0.13 for carbon amended soils and 0.03–0.08 for carbon and nitrogen amended soils and did not vary much between ecosystems. As a result of uncertainty in the literature about the relationship between growth rate and CUE14, several different relationships were postulated to model per-taxon CUE as a function of per-taxon growth rate: linear increase, linear decrease, exponential decrease, unimodal with peak CUE at growth rate of 0.5, and unimodal with peak CUE at a growth rate of 0.05 (the median of all per-taxon growth rates in the data). Comparisons between functions were made by calculating AIC values from per-taxon respiration, summed, and regressing against measured respiration values. Likewise, for each function, we tested how well per-taxon CUE estimates reconstructed community-level CUE by weighting the CUE value of each taxon by its relative abundance, summing, and regressing against community-level CUE. To select the best per-taxon CUE function, AIC values from both scaling efforts were combined. To make AIC values comparable, all respiration and CUE terms were z-transformed prior to regression scaling. To reflect our priority of estimating per-taxon respiration, AIC values from the respiration scaling regression models were multiplied by two and summed with AIC values from CUE scaling such that AICTotal = 2(AICResp) + AICCUE. Across these comparisons, the best estimate of per-taxon CUE was the unimodal function of growth rate, constrained by community-level CUE and peaking at growth rates of 0.5 (Table 1), such that:$${{rm{CUE}}}_{{rm{i}}}=-4({{rm{CUE}}}_{{rm{E}}{rm{:}}{rm{T}}{rm{:}}{{rm{range}}}})cdot {left({g}_{{rm{i}}}-0.5right)}^{2}+({{rm{CUE}}}_{{rm{E}}{rm{:}}{rm{T}}{rm{:}}{max }}),$$
(6)
where CUEi indicates per-taxon CUE, CUEE:T:max indicates the maximum CUE values observed for a group of replicates within a given ecosystem and treatment (E:T). With this function, higher per-capita growth rate values were parameterized to produce higher CUE values initially and then decrease reflecting a growth-CUE tradeoff14, here bound by the difference in maximum and minimum CUE values. We applied per-taxon CUE estimates from Eq. 6 to per-taxon growth rates to yield estimates of per-taxon respiration:$${r}_{{rm{i}}}={r}_{{rm{g,i}}}+{r}_{{rm{m,i}}}=left(frac{{g}_{{rm{i}}}}{{{rm{CUE}}}_{{rm{i}}}}-{g}_{{rm{i}}}right)+left(frac{{g}_{{rm{i}}}}{{{rm{CUE}}}_{{rm{i}}}}-{g}_{{rm{i}}}right)cdot beta,$$
(7)
where ri indicates per-capita respiration for taxon i, rg,i indicates growth-related respiration, rm,i indicates maintenance-related respiration, and β is a constant of 0.01 that represents the maintenance requirements as a proportion of total energy use40. We used these values of per-taxon, per-capita respiration rates to estimate per-taxon respiration per gram of dry soil per week:$${R}_{{rm{i}}}={P}_{{rm{i}}}cdot {r}_{{{rm{g,i}}}}+{P}_{{rm{i}}}cdot {r}_{{{rm{m,i}}}},$$
(8)
where Ri indicates respiration of CO2–C (µg C g dry soil−1 week−1) for taxon i.In addition to per-taxon respiration estimates based on 18O enrichment, we used another model for comparison. Here, respiration was calculated based on 16S abundance alone:$${R}_{{rm{i}}}={N}_{{rm{i}}}cdot f(R sim N+0),$$
(9)
where Ni indicates final 16S abundance for taxon i, R indicates microbial respiration of CO2-C (µg C g dry soil−1 week−1) and N indicates total 16S abundance at the end of the incubation. Here, the R ~ N function indicates the linear relationship, with an intercept of 0, between CO2 respiration and 16S gene concentration across all samples.Diversity, compositional, and statistical analysisFor patterns of evenness in bacterial carbon use and relative abundance, we used Pielou’s evenness which is the quotient of Shannon’s diversity and the observed richness. For each sample, we applied Pielou’s evenness to bacterial abundances as well as bacterial carbon use (relativized to sum to one, in both cases).We created a linear mixed model to test the relationship between the carbon use (the sum of biomass production and respiration) and relative abundance of bacterial genera from the dominant phyla, which accounted for >90% of all C flux. Here, we averaged carbon use and relative abundance for all replicates in a given ecosystem and treatment. We used the lme4 R package (version 1.1-20)41 and obtained p-values using the Satterthwaite method in the lmerTest R package (version 3.1-0)42. To limit pseudo-replication, we accounted for differences in carbon use across ecosystems and due to bacterial Genus by implementing random intercepts. We selected for the optimal random and fixed components by dropping individual terms and comparing models with likelihood ratio tests, disregarding models that failed to converge. Our final model fit was:$${{{log }}}_{10}({C}_{{rm{i}}}) sim {{{log }}}_{10}left({y}_{{rm{i}}}right)ast T+left(1|Eright)+(1|{{rm{Genus}}}),$$
(10)
where Ci indicates the relativized carbon use for taxon i (averaged across all three replicates in a given ecosystem and treatment), yi indicates the relative abundance of taxon i (averaged across all three replicates), T indicates soil treatment, and E indicates ecosystem.For differences in composition, we created species abundance tables using the traditional abundances, as well as measures of carbon use (growth and maintenance respiration) of each ASV in each sample. To account for differences in absolute abundances and flux rates between sites, we relativized all abundance tables. We summarized compositional differences using Bray–Curtis dissimilarities then identified multivariate centroids for all replicates in a site by treatment group. We tested the effect of site and nutrient amendment on the resulting group centroids using PERMANOVA tests implemented with the adonis function in the vegan package (version 2.5-3)43. We related compositional shifts in relative abundance to those in relativized growth and maintenance using Mantel tests with the mantel function in vegan.To test for changes in the type of soil C preferred by microbial genera (either 13C-labeled glucose or 12C soil carbon) in response to nitrogen addition, we used Levene’s test with the car package (version 3.0-10)44. Specifically, we analyzed the relationship between 13C use and 12C use (both relativized) on bacterial genera across all replicates and in C and C + N treatments using a linear model. We then extracted model residuals and tested whether variance was significantly different across treatments by focusing on the interaction between individual replicates and treatment. This produced a significance test describing treatment-level differences in 13C–12C use.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More