Ecology

Bay, R. A. et al. Genetic coupling of female mate choice with polygenic ecological divergence facilitates stickleback speciation. Curr. Biol. 27, 3344–3349 (2017).CAS 

Google Scholar 
Johannesson, K., Butlin, R. K., Panova, M. & Westram, A. M. Population Genomics: Marine Organisms (eds. Oleksiak, M. F. & Rajora, O. P.) 277–301 (Springer, 2017).Riesch, R. et al. Transitions between phases of genomic differentiation during stick-insect speciation. Nat. Ecol. Evol. 1, 1–13 (2017).
Google Scholar 
Feder, J. L. & Nosil, P. The efficacy of divergence hitchhiking in generating genomic islands during ecological speciation. Evolution 64, 1729–1747 (2010).
Google Scholar 
Rosenblum, E. B., Hickerson, M. J. & Moritz, C. A multilocus perspective on colonization accompanied by selection and gene flow. Evolution 61, 2971–2985 (2007).CAS 

Google Scholar 
Nosil, P., Egan, S. P. & Funk, D. J. Heterogeneous genomic differentiation between walking‐stick ecotypes: “isolation by adaptation” and multiple roles for divergent selection. Evolution 62, 316–336 (2008).
Google Scholar 
Orsini, L., Vanoverbeke, J., Swillen, I., Mergeay, J. & Meester, L. Drivers of population genetic differentiation in the wild: isolation by dispersal limitation, isolation by adaptation and isolation by colonization. Mol. Ecol. 22, 5983–5999 (2013).
Google Scholar 
Sexton, J. P., Hangartner, S. B. & Hoffmann, A. A. Genetic isolation by environment or distance: which pattern of gene flow is most common? Evolution 68, 1–15 (2014).CAS 

Google Scholar 
Roderick, G. K. & Gillespie, R. G. Speciation and phylogeography of Hawaiian terrestrial arthropods. Mol. Ecol. 7, 519–531 (1998).CAS 

Google Scholar 
Juan, C., Emerson, B. C., Oromı́, P. & Hewitt, G. M. Colonization and diversification: towards a phylogeographic synthesis for the Canary Islands. Trends Ecol. Evol. 15, 104–109 (2000).CAS 

Google Scholar 
Brown, R. P., Hoskisson, P. A., Welton, J. H. & Báez, M. Geological history and within‐island diversity: a debris avalanche and the Tenerife lizard Gallotia galloti. Mol. Ecol. 15, 3631–3640 (2006).CAS 

Google Scholar 
O’Connell, K. A., Prates, I., Scheinberg, L. A., Mulder, K. P. & Bell, R. C. Speciation and secondary contact in a fossorial island endemic, the São Tomé caecilian. Mol. Ecol. 30, 2859–2871 (2021).
Google Scholar 
Malhotra, A. & Thorpe, R. S. The dynamics of natural selection and vicariance in the Dominican anole: patterns of within‐island molecular and morphological divergence. Evolution 54, 245–258 (2000).CAS 

Google Scholar 
Brown, R. P., Woods, M. & Thorpe, R. S. Historical volcanism and within-island genetic divergence in the Tenerife skink (Chalcides viridanus). Biol. J. Linnean Soc. 122, 166–175 (2017).
Google Scholar 
Losos, J. Lizards in an Evolutionary Tree: Ecology and Adaptive Radiation of Anoles (University of California Press, 2009).Mahler, D. L., Revell, L. J., Glor, R. E. & Losos, J. B. Ecological opportunity and the rate of morphological evolution in the diversification of Greater Antillean anoles. Evolution 64, 2731–2745 (2010).
Google Scholar 
Wang, I. J., Glor, R. E. & Losos, J. B. Quantifying the roles of ecology and geography in spatial genetic divergence. Ecol. Lett. 16, 175–182 (2013).
Google Scholar 
Beerli, P. & Felsenstein, J. Maximum-likelihood estimation of migration rates and effective population numbers in two populations using a coalescent approach. Genetics 152, 763–773 (1999).CAS 

Google Scholar 
Hey, J. & Nielsen, R. Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis. Genetics 167, 747–760 (2004).CAS 

Google Scholar 
Hey, J. Recent advances in assessing gene flow between diverging populations and species. Curr. Opin. Genet. Dev. 16, 592–596 (2006).CAS 

Google Scholar 
Excoffier, L., Dupanloup, I., Huerta-Sánchez, E., Sousa, V. C. & Foll, M. Robust demographic inference from genomic and SNP data. PLoS Genet. 9, 1003905 (2013).
Google Scholar 
Butlin, R. K. et al. Parallel evolution of local adaptation and reproductive isolation in the face of gene flow. Evolution 68, 935–949 (2014).
Google Scholar 
Rosenblum, E. B., Hoekstra, H. E. & Nachman, M. W. Adaptive reptile color variation and the evolution of the MCIR gene. Evolution 58, 1794–1808 (2004).CAS 

Google Scholar 
Rosenblum, E. B. Convergent evolution and divergent selection: lizards at the White Sands ecotone. Am. Nat. 167, 1–15 (2006).
Google Scholar 
Sumner, F. B. An analysis of geographic variation in mice of the Peromyscus polionotus group from Florida and Alabama. J. Mammal. 7, 149–184 (1926).
Google Scholar 
Davenport, J., & Dellinger, T. Melanism and foraging behaviour in an intertidal population of the Madeiran lizard Podarcis (= Lacerta) dugesii (Milne-Edwards, 1829). Herpetol. J. 5, 200–203 (1995).
Google Scholar 
Galán, P. Demography and population dynamics of the lacertid lizard Podarcis bocagei in north-west Spain. J. Zool. 249, 203–218 (1999).
Google Scholar 
Censky, E. J., Hodge, K. & Dudley, J. Over-water dispersal of lizards due to hurricanes. Nature 395, 556 (1998).CAS 

Google Scholar 
Rolán‐Alvarez, E., Erlandsson, J., Johannesson, K. & Cruz, R. Mechanisms of incomplete prezygotic reproductive isolation in an intertidal snail: testing behavioural models in wild populations. J. Evol. Biol. 12, 879–890 (1999).
Google Scholar 
Ludt, W. B. & Rocha, L. A. Shifting seas: the impacts of Pleistocene sea‐level fluctuations on the evolution of tropical marine taxa. J. Biogeogr. 42, 25–38 (2015).
Google Scholar 
Lambeck, K. Late Pleistocene, Holocene and present sea-levels: constraints on future change. Glob. Planet Change 3, 205–217 (1990). & J.
Google Scholar 
Rosenblum, E. B. The role of phenotypic plasticity in color variation of Tularosa Basin lizards. Copeia 2005, 586–596 (2005).
Google Scholar 
Jin, Y. et al. Dorsal pigmentation and its association with functional variation in MC1R in a lizard from different elevations on the Qinghai–Tibetan plateau. Genome Biol. Evol. 12, 2303–2313 (2020).CAS 

Google Scholar 
Corl, A. et al. The genetic basis of adaptation following plastic changes in coloration in a novel environment. Curr. Biol. 28, 2970–2977 (2018).CAS 

Google Scholar 
Sacchi, R. et al. Genetic and phenotypic component in head shape of common wall lizard Podarcis muralis. Amphib.-Reptilia 37, 301–310 (2016).
Google Scholar 
Dice, L. R. Variation of the deer-mouse (Peromyscus maniculatus) on the Sand Hills of Nebraska and adjacent areas. Contrib. Lab Vertebrate Biol. Univ. Mich. 15, 1–19 (1941).
Google Scholar 
Vitt, L. J., Caldwell, J. P., Zani, P. A. & Titus, T. A. The role of habitat shift in the evolution of lizard morphology: evidence from tropical Tropidurus. Proc. Natl Acad. Sci. USA 94, 3828–3832 (1997).CAS 

Google Scholar 
Pfeifer, S. P. et al. The evolutionary history of Nebraska deer mice: local adaptation in the face of strong gene flow. Mol. Biol. Evol. 35, 792–806 (2018).CAS 

Google Scholar 
Scherrer, R., Donihue, C. M., Reynolds, R. G., Losos, J. B. & Geneva, A. J. Dewlap colour variation in Anolis sagrei is maintained among habitats within islands of the West Indies. J. Evol. Biol. 35, 680–692 (2022).
Google Scholar 
Janson, K. Selection and migration in two distinct phenotypes of Littorina saxatilis in Sweden. Oecologia 59, 58–61 (1983).CAS 

Google Scholar 
Richardson, J. L., Urban, M. C., Bolnick, D. I. & Skelly, D. K. Microgeographic adaptation and the spatial scale of evolution. Trends Ecol. Evol. 29, 165–176 (2014).
Google Scholar 
Engelstoft, C., Robinson, J., Fraser, D. & Hanke, G. Recent rapid expansion of common wall lizards (Podarcis muralis) in British Columbia, Canada. Northwest. Naturalist 101, 50–55 (2020).
Google Scholar 
Cascio, P. L. & Pasta, S. Preliminary data on the biometry and the diet of a microinsular population of Podarcis wagleriana (Reptilia: Lacertidae). Acta Herpetol. 1, 147–152 (2006).
Google Scholar 
Janssen, J., Towns, D. R., Duxbury, M. & Heitkönig, I. M. Surviving in a semi-marine habitat: dietary salt exposure and salt secretion of a New Zealand intertidal skink. Comp. Biochem Physiol. A Mol. Integr. Physiol. 189, 21–29 (2015).CAS 

Google Scholar 
Grismer, L. L. Three new species of intertidal side-blotched lizards (genus Uta) from the Gulf of California, Mexico. Herpetologica 50, 451–474 (1994).
Google Scholar 
Sepúlveda, M., Sabat, P., Porter, W. P. & Fariña, J. M. One solution for two challenges: the lizard Microlophus atacamensis avoids overheating by foraging in intertidal shores. PLoS One 9, 97735 (2014).
Google Scholar 
Hobson, E. S. Observations on diving in the Galapagos marine iguana, Amblyrhynchus cristatus (Bell). Copeia 1965, 249–250 (1965).Balakrishna, S., Amdekar, M. S. & Thaker, M. Morphological divergence, tail loss, and predation risk in urban lizards. Urban Ecosyst. 24, 1391–1398 (2021).
Google Scholar 
Falvey, C. H., Aviles-Rodriguez, K. J., Hagey, T. J. & Winchell, K. M. The finer points of urban adaptation: intraspecific variation in lizard claw morphology. Biol. J. Linn. Soc. 131, 304–318 (2020).
Google Scholar 
Marnocha, E., Pollinger, J. & Smith, T. B. Human‐induced morphological shifts in an island lizard. Evol. Appl 4, 388–396 (2011).
Google Scholar 
Rocha, R., Paixão, M. & Gouveia, R. Predation note: Anthus berthelotii madeirensis (Passeriformes: Motacillidae) catches Teira dugesii mauli (Squamata: Lacertidae) in Deserta Grande, Madeira Archipel. Herpetol. Notes 3, 77–78 (2010).
Google Scholar 
Völkl, W. & Brandl, R. Tail break rate in the Madeiran lizard (Podarcis dugesii). Amphibia-Reptilia 9, 213–218 (1988).Malhotra, A. & Thorpe, R. S. Microgeographic variation in Anolis oculatus, on the island of Dominica, West Indies. J. Evol. Biol. 4, 321–335 (1991).
Google Scholar 
Thorpe, R. S. & Brown, R. P. Microgeographic variation in the colour pattern of the lizard Gallotia galloti within the island of Tenerife: distribution, pattern and hypothesis testing. Biol. J. Linn. Soc. 38, 303–322 (1989).
Google Scholar 
Brown, R. P., Thorpe, R. S. & Báez, M. Parallel within-island microevolution of lizards on neighbouring islands. Nature 352, 60–62 (1991).
Google Scholar 
Báez, M. & Brown, R. P. Testing multivariate patterns of within‐island differentiation in Podarcis dugesii from Madeira. J. Evol. Biol. 10, 575–587 (1997).
Google Scholar 
Cook, L. M. Density of lizards in Madeira. Bocagiana (Funchal) 66, 1–3 (1983).
Google Scholar 
Sadek, R. A. The diet of the Madeiran lizard Lacerta dugesii. Zool. J. Linn. Soc. 73, 313–341 (1981).
Google Scholar 
Brehm, A. et al. Phylogeography of the Madeiran endemic lizard Lacerta dugesii inferred from mtDNA sequences. Mol. Phylogenetics Evol. 26, 222–230 (2003).CAS 

Google Scholar 
Suárez, N. M., Pestano, J. & Brown, R. P. Ecological divergence combined with ancient allopatry in lizard populations from a small volcanic island. Mol. Ecol. 23, 4799–4812 (2014).
Google Scholar 
Towns, D. R. Ecology of the black shore skink, Leiolopisma suteri (Lacertilia: Scincidae), in boulder beach habitats. N. Z. J. Zool. 2, 389–407 (1975).
Google Scholar 
Cook, L. M. Variation in the Madeiran lizard Lacerta dugesii. J. Zool. 187, 327–340 (1979).
Google Scholar 
Troscianko, J. & Stevens, M. Image calibration and analysis toolbox–a free software suite for objectively measuring reflectance, colour, and pattern. Methods Ecol. Evol. 6, 1320–1331 (2015).
Google Scholar 
Schneider, C. A., Rasband, W. S. & Eliceiri, K. W. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671–675 (2012).CAS 

Google Scholar 
Rohlf, F. J. The tps series of software. Hystrix, Ital. J. Mammal. 26, 9–12 (2015).
Google Scholar 
Bookstein, F. L. Morphometric Tools for Landmark Data: Geometry and Biology (Cambridge University Press, 1991).Klingenberg, C. P. MorphoJ: an integrated software package for geometric morphometrics. Mol. Ecol. Resour. 11, 353–357 (2011).
Google Scholar 
Rohlf, F. J. & Slice, D. Extensions of the Procrustes method for the optimal superimposition of landmarks. Syst. Biol. 39, 40–59 (1990).
Google Scholar 
Klingenberg, C. P., Barluenga, M. & Meyer, A. Shape analysis of symmetric structures: quantifying variation among individuals and asymmetry. Evolution 56, 1909–1920 (2002).
Google Scholar 
Andrews, S. FastQC: a Quality Control Tool for High Throughput Sequence Data. Babraham Bioinformatics version 0.11.7. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (2010).Melo, A. T., Bartaula, R. & Hale, I. GBS-SNP-CROP: a reference-optional pipeline for SNP discovery and plant germplasm characterization using variable length, paired-end genotyping-by-sequencing data. BMC Bioinform. 17, 1–15 (2016).
Google Scholar 
Sabadin, F., Carvalho, H. F., Galli, G. & Fritsche-Neto, R. Population-tailored mock genome enables genomic studies in species without a reference genome. Mol. Genet. Genom. 297, 33–46 (2022).CAS 

Google Scholar 
Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, 2156–2158 (2011).CAS 

Google Scholar 
Pfeifer, B., Wittelsbürger, U., Ramos-Onsins, S. E. & Lercher, M. J. PopGenome: an efficient swiss army knife for population genomic analyses in R. Mol. Biol. Evol. 31, 1929–1936 (2014).CAS 

Google Scholar 
Team, R. C. R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/ (2022).Jombart, T. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24, 1403–1405 (2008).CAS 

Google Scholar 
Luu, K., Bazin, E. & Blum, M. G. pcadapt: an R package to perform genome scans for selection based on principal component analysis. Mol. Ecol. Resour. 17, 67–77 (2017).CAS 

Google Scholar 
Günther, T. & Coop, G. Robust identification of local adaptation from allele frequencies. Genetics 195, 205–220 (2013).
Google Scholar 
Dray, S. et al. Package ‘adespatial.’ Available from: https://cran.r-project.org/package=adespatial (2018).Montano, V. & Jombart, T. An eigenvalue test for spatial principal component analysis. BMC Bioinform. 18, 1–7 (2017).
Google Scholar 
Pickrell, J. K. & Pritchard, J. K. Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet. 8, e1002967 (2012).CAS 

Google Scholar  More

Muller, Z. et al. Giraffa camelopardalis. The IUCN red list of threatened species 2016:e.T9194A109326950 (2018).Oconnor, D. et al. Updated geographic range maps for giraffe, Giraffa spp., throughout sub-Saharan Africa, and implications of changing distributions for conservation. Mamm. Rev. 49, 285–299. https://doi.org/10.1111/mam.12165 (2019).Article 

Google Scholar 
Brown, M. B. et al. Conservation status of giraffe: Evaluating contemporary distribution and abundance with evolving taxonomic perspectives. Ref. Module Earth Syst. Environ. Sci. https://doi.org/10.1016/B978-0-12-821139-7.00139-2 (2021).Article 

Google Scholar 
Dunn, M. E. et al. Investigating the international and pan-African trade in giraffe parts and derivatives. Conserv. Sci. Pract. 3, e390. https://doi.org/10.1111/csp2.390 (2021).Article 

Google Scholar 
Hassanin, A. et al. Mitochondrial DNA variability in Giraffa camelopardalis: Consequences for taxonomy, phylogeography and conservation of giraffes in West and Central Africa. C.R. Biol. 330, 265–274. https://doi.org/10.1016/j.crvi.2007.02.008 (2007).Article 
CAS 

Google Scholar 
Groves, C. & Grubb, P. Ungulate Taxonomy (Johns Hopkins University Press, 2011).Book 

Google Scholar 
Fennessy, J. et al. Multi-locus analyses reveal four giraffe species instead of one. Curr. Biol. 26, 1–7. https://doi.org/10.1016/j.cub.2016.07.036 (2016).Article 
CAS 

Google Scholar 
Winter, S., Fennessy, J. & Janke, A. Limited introgression supports division of giraffe into four species. Ecol. Evol. 8, 10156–10165. https://doi.org/10.1002/ece3.4490 (2018).Article 

Google Scholar 
Bercovitch, F. B. Giraffe taxonomy, geographic distribution, and conservation. Afr. J. Ecol. 58, 150–158. https://doi.org/10.1111/aje.12741 (2020).Article 

Google Scholar 
Petzold, A. & Hassanin, A. A comparative approach for species delimitation based on multiple methods of multi-locus DNA sequence analysis: A case study of the genus Giraffa (Mammalia, Cetartiodactyla). PLoS ONE 15, e0217956. https://doi.org/10.1371/journal.pone.0217956 (2020).Article 
CAS 

Google Scholar 
Petzold, A. et al. First insights into past biodiversity of giraffes based on mitochondrial sequences from museum specimens. Eur. J. Taxon. 703, L57-63. https://doi.org/10.1371/journal.pone.0217956 (2020).Article 
CAS 

Google Scholar 
Coimbra, R. T. F. et al. Whole-genome analysis of giraffe supports four distinct species. Curr. Biol. 31, 2929-2938.e5. https://doi.org/10.1016/j.cub.2021.04.033 (2021).Article 
CAS 

Google Scholar 
Muneza, A. B. et al. Giraffa camelopardalis ssp. reticulata. The IUCN Red List of Threatened Species 2018:e.T88420717A88420720 (2018).Miller, M. F. Dispersal of Acacia seeds by ungulates and ostriches in an African Savanna. J. Trop. Ecol. 12, 345–356. https://doi.org/10.1017/S0266467400009548 (1996).Article 

Google Scholar 
Palmer, T. M. et al. Breakdown of an ant-plant mutualism follows the loss of large herbivores from an African savanna. Science 319, 192–195. https://doi.org/10.1126/science.1151579 (2008).Article 
ADS 
CAS 

Google Scholar 
Kalema, G. Investigation of a skin disease in giraffe in Murchison Falls National Park. Report Submitted to Uganda National Park. Uganda National Parks. Kampala, Uganda (1996).Muneza, A. B. et al. Regional variation of the manifestation, prevalence, and severity of giraffe skin disease: A review of an emerging disease in wild and captive giraffe populations. Biol. Conserv. 198, 145–156. https://doi.org/10.1016/j.biocon.2016.04.014 (2016).Article 

Google Scholar 
Epaphras, A. M., Karimuribo, E. D., Mpanduji, D. G. & Meing’ataki, G. E. Prevalence, disease description and epidemiological factors of a novel skin disease in giraffes (Giraffa camelopardalis) in Ruaha National Park, Tanzania. Res. Opin. Anim. Vet. Sci. 2, 60–65 (2012).
Google Scholar 
Lee, D. E. & Bond, M. L. The occurrence and prevalence of giraffe skin disease in protected areas of northern Tanzania. J. Wildl. Dis. 52, 753–755. https://doi.org/10.7589/2015-09-24 (2016).Article 

Google Scholar 
Muneza, A. B. et al. Examining disease prevalence for species of conservation concern using non-invasive spatial capture–recapture techniques. J. Appl. Ecol. 54, 709–717. https://doi.org/10.1111/1365-2664.12796 (2017).Article 

Google Scholar 
Brown, M. Murchison falls giraffe project: Field report. Giraffid 9, 5–10 (2015).
Google Scholar 
Muneza, A. B. et al. Quantifying the severity of an emerging skin disease affecting giraffe populations using photogrammetry analysis of camera trap data. J. Wildl. Dis. 55, 770–781. https://doi.org/10.7589/2018-06-149 (2019).Article 

Google Scholar 
Han, S. et al. Giraffe skin disease: Clinicopathologic characterization of cutaneous filariasis in the critically endangered Nubian giraffe (Giraffa camelopardalis camelopardalis). Vet. Pathol. https://doi.org/10.1177/03009858221082606 (2022).Article 

Google Scholar 
Whittier, C. A. et al. Cutaneous filariasis in free-ranging Rothschild’s giraffes (Giraffa Camelopardalis rothschildi) in Uganda. J. Wildl. Dis. 56, 1–5. https://doi.org/10.7589/2018-09-212 (2020).Article 

Google Scholar 
Pellew, R. Food consumption and energy budgets of the giraffe. J. Appl. Ecol. 21, 141–159. https://doi.org/10.2307/2403043 (1984).Article 

Google Scholar 
Strauss, M. K. L. & Packer, C. Using claw marks to study lion predation on giraffes of the Serengeti. J. Zool. 289, 134–142. https://doi.org/10.1111/j.1469-7998.2012.00972.x (2013).Article 

Google Scholar 
Muneza, A. B. et al. Exploring the connections between giraffe skin disease and lion predation. J. Zool. https://doi.org/10.1111/jzo.12930 (2021).Article 

Google Scholar 
Lindsey, P. A. et al. The bushmeat trade in African savannas: Impacts, drivers, and possible solutions. Biol. Conserv. 160, 80–96. https://doi.org/10.1016/j.biocon.2012.12.020 (2013).Article 

Google Scholar 
Becker, M. et al. Evaluating wire-snare poaching trends and the impacts of by-catch on elephants and large carnivores. Biol. Conserv. 158, 26–36. https://doi.org/10.1016/j.biocon.2012.08.017 (2013).Article 

Google Scholar 
Mudumba, T., Jingo, S., Heit, D. & Montgomery, R. A. The landscape configuration and lethality of snare poaching of sympatric guilds of large carnivores and ungulates. Afr. J. Ecol. 59, 51–62. https://doi.org/10.1111/aje.12781 (2020).Article 

Google Scholar 
Strauss, M. K. L., Kilewo, M., Rentsch, D. & Packer, C. Food supply and poaching limit giraffe abundance in the Serengeti. Popul. Ecol. 57, 505–516. https://doi.org/10.1007/s10144-015-0499-9 (2015).Article 

Google Scholar 
Munn, J. Effects of injury on the locomotion of free-ranging chimpanzees in the Budongo Forest Reserve, Uganda. In Primates of Western Uganda: Developments in Primatology: Progress and Prospects (eds. Newton-Fisher, N. E., Notman, H., Paterson, J. D., & Reynolds, V.) 259–280 (Springer, 2006).Yersin, H., Asiimwe, C., Voordouw, M. J. & Zuberbühler, K. Impact of snare injuries on parasite prevalence in wild chimpanzees (Pan troglodytes). Int. J. Primatol. 38, 21–30. https://doi.org/10.1007/s10764-016-9941-x (2017).Article 

Google Scholar 
Dagg, A. I. Gaits of the giraffe and okapi. J. Mammal. 41, 282–282. https://doi.org/10.2307/1376381 (1960).Article 

Google Scholar 
Dagg, A. I. The role of the neck in the movements of the giraffe. J. Mammal. 43, 88–97. https://doi.org/10.2307/1376883 (1962).Article 

Google Scholar 
Dagg, A. I. & Vos, A. D. The walking gaits of some species of Pecora. J. Zool. 155, 103–110. https://doi.org/10.1111/j.1469-7998.1968.tb03031.x (1968).Article 

Google Scholar 
Alexander, R. M. N., Langman, V. A. & Jayes, A. S. Fast locomotion of some African ungulates. J. Zool. 183, 291–300. https://doi.org/10.1111/j.1469-7998.1977.tb04188.x (1977).Article 

Google Scholar 
Basu, C., Deacon, F., Hutchinson, J. R. & Wilson, A. M. The running kinematics of free-roaming giraffes, measured using a low cost unmanned aerial vehicle (UAV). PeerJ 7, e6312. https://doi.org/10.7717/peerj.6312 (2019).Article 

Google Scholar 
Basu, C., Wilson, A. M. & Hutchinson, J. R. The locomotor kinematics and ground reaction forces of walking giraffes. J. Exp. Biol. 222, jeb159277. https://doi.org/10.1242/jeb.159277 (2019).Article 

Google Scholar 
Hildebrand, M. The adaptive significance of tetrapod gait selection. Am. Zool. 20, 255–267. https://doi.org/10.1093/icb/20.1.255 (1980).Article 

Google Scholar 
Flower, F. C., Sanderson, D. J. & Weary, D. M. Hoof pathologies influence kinematic measures of dairy cow gait. J. Dairy Sci. 88, 3166–3173. https://doi.org/10.3168/jds.s0022-0302(05)73000-9 (2005).Article 
CAS 

Google Scholar 
Brown, M. B., Bolger, D. T. & Fennessy, J. All the eggs in one basket: A countrywide assessment of current and historical giraffe population distribution in Uganda. Glob. Ecol. Conserv. 19, e00612. https://doi.org/10.1016/j.gecco.2019.e00612 (2019).Article 

Google Scholar 
Foster, J. B. The giraffe of Nairobi National Park: Home range, sex ratios, the herd, and food. Afr. J. Ecol. 4, 139–148. https://doi.org/10.1111/j.1365-2028.1966.tb00889.x (1966).Article 

Google Scholar 
Bond, M. L., Strauss, M. K. L. & Lee, D. E. Soil correlates and mortality from giraffe skin disease in Tanzania. J. Wildl. Dis. 52, 953–958. https://doi.org/10.7589/2016-02-047 (2016).Article 

Google Scholar 
Dunham, N. T., McNamara, A., Shapiro, L., Hieronymus, T. & Young, J. W. A user’s guide for the quantitative analysis of substrate characteristics and locomotor kinematics in free-ranging primates. Am. J. Phys. Anthropol. 167, 569–584. https://doi.org/10.1002/ajpa.23686 (2018).Article 

Google Scholar 
Rueden, C. T. et al. Imagej 2: Imagej for the next generation of scientific image data. BMC Bioinform. 18, 529. https://doi.org/10.1186/s12859-017-1934-z (2017).Article 

Google Scholar 
Cartmill, M., Lemelin, P. & Schmitt, D. Support polygons and symmetrical gaits in mammals. Zool. J. Linn. Soc. 136, 401–420. https://doi.org/10.1046/j.1096-3642.2002.00038.x (2002).Article 

Google Scholar 
Hildebrand, M. Analysis of the symmetrical gaits of tetrapods. Folia Biotheoretica 6, 1–22. https://doi.org/10.2307/1379571 (1966).Article 

Google Scholar 
Shapiro, L. J. & Young, J. W. Kinematics of quadrupedal locomotion in sugar gliders (Petaurus breviceps): Effects of age and substrate size. J. Exp. Biol. 215, 480–496. https://doi.org/10.1242/jeb.062588 (2012).Article 

Google Scholar 
Shapiro, L. J., Young, J. W. & VandeBerg, J. L. Body size and the small branch niche: Using marsupial ontogeny to model primate locomotor evolution. J. Hum. Evol. 68, 14–31. https://doi.org/10.1016/j.jhevol.2013.12.006 (2014).Article 

Google Scholar 
Dunham, N. T., McNamara, A., Shapiro, L., Phelps, T. & Young, J. W. Asymmetrical gait kinematics of free-ranging callitrichines in response to changes in substrate diameter, orientation, and displacement. J. Exp. Biol. 223, jeb217562. https://doi.org/10.1242/jeb.217562 (2020).Article 

Google Scholar 
Robinson, R., Herzog, W. & Nigg, B. Use of force platform variables to quantify the effects of chiropractic manipulation on gait symmetry. J. Manipulative Physiol. Ther. 10, 172–176 (1987).CAS 

Google Scholar 
Vanden Hole, C. et al. How innate is locomotion in precocial animals? A study on the early development of spatiotemporal gait variables and gait symmetry in piglets. J. Exp. Biol. 220, 2706–2716. https://doi.org/10.1242/jeb.157693 (2017).Article 

Google Scholar 
Jacobs, B. Y., Kloefkorn, H. E. & Allen, K. D. Gait analysis methods for rodent models of osteoarthritis. Curr. Pain Headache Rep. 18, 456–475. https://doi.org/10.1007/s11916-014-0456-x (2014).Article 

Google Scholar 
Pfau, T., Spence, A., Starke, S., Ferrari, M. & Wilson, A. Modern riding style improves horse racing times. Science 325, 289–289. https://doi.org/10.1126/science.1174605 (2009).Article 
ADS 
CAS 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2019). http://www.R-project.org/.Kuznetsova, A., Brockhoff, P. B. & Christensen, R. H. B. LmerTest package: Tests in linear mixed effects models. J. Stat. Softw. https://doi.org/10.18637/jss.v082.i13 (2017).Article 

Google Scholar 
Length, R. emmeans: Estimated marginal means, aka least‐squares means. R package version 0.9. https://CRAN.R-project.org/package=emmeans (2017).Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B (Methodol.) 57, 289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x (1995).Article 
MathSciNet 
MATH 

Google Scholar 
Merkens, H. W. & Schamhardt, H. C. Evaluation of equine locomotion during different degrees of experimentally induced lameness I: Lameness model and quantification of ground reaction force patterns of the limbs. Equine Vet. J. 20, 99–106. https://doi.org/10.1111/j.2042-3306.1988.tb04655.x (1988).Article 

Google Scholar 
Fanchon, L. & Grandjean, D. Accuracy of asymmetry indices of ground reaction forces for diagnosis of hind limb lameness in dogs. Am. J. Vet. Res. 68, 1089–1094. https://doi.org/10.2460/ajvr.68.10.1089 (2007).Article 

Google Scholar 
Bragança, F. M. S., Rhodin, M. & van Weeren, P. R. On the brink of daily clinical application of objective gait analysis: What evidence do we have so far from studies using an induced lameness model?. Vet. J. 234, 11–23. https://doi.org/10.1016/j.tvjl.2018.01.006 (2018).Article 

Google Scholar 
Brown, M. B. & Bolger, D. T. Male-biased partial migration in a giraffe population. Front. Ecol. Evol. 7, 524. https://doi.org/10.3389/fevo.2019.00524 (2020).Article 

Google Scholar 
Dagg, A. I. Giraffe: Biology, Behaviour and Conservation (Cambridge University Press, 2014).Book 

Google Scholar 
Castles, M. P. et al. Relationships between male giraffes’ colour, age and sociability. Anim. Behav. 157, 13–25. https://doi.org/10.1016/j.anbehav.2019.08.003 (2019).Article 

Google Scholar  More

Anae, J. et al. Recent advances in biochar engineering for soil contaminated with complex chemical mixtures: Remediation strategies and future perspectives. Sci. Total Environ. 767, 144351. https://doi.org/10.1016/j.scitotenv.2020.144351 (2021).Article 
ADS 
CAS 

Google Scholar 
Kiran, B. R. & Prasad, M. N. V. Biochar and rice husk ash assisted phytoremediation potentials of Ricinus communis L. for lead-spiked soils. Ecotoxicol Environ Saf 183, 109574. https://doi.org/10.1016/j.ecoenv.2019.109574 (2019).Article 
CAS 

Google Scholar 
Bolan, N. et al. Remediation of heavy metal(loid)s contaminated soils – To mobilize or to immobilize?. J. Hazard. Mater. 266, 141–166. https://doi.org/10.1016/j.jhazmat.2013.12.018 (2014).Article 
CAS 

Google Scholar 
Burachevskaya, M. et al. The effect of granular activated carbon and biochar on the availability of Cu and Zn to Hordeum sativum distichum in contaminated soil. Plants https://doi.org/10.3390/plants10050841 (2021).Article 

Google Scholar 
Cao, P. et al. Mercapto propyltrimethoxysilane- and ferrous sulfate-modified nano-silica for immobilization of lead and cadmium as well as arsenic in heavy metal-contaminated soil. Environ. Pollut. 266, 115152. https://doi.org/10.1016/j.envpol.2020.115152 (2020).Article 
CAS 

Google Scholar 
Ok, Y. S. et al. Ameliorants to immobilize Cd in rice paddy soils contaminated by abandoned metal mines in Korea. Environ. Geochem. Health 33(Suppl 1), 23–30. https://doi.org/10.1007/s10653-010-9364-0 (2011).Article 
CAS 

Google Scholar 
Qin, Y. et al. Dual-wastes derived biochar with tailored surface features for highly efficient p-nitrophenol adsorption. J. Clean. Prod. 353, 131571. https://doi.org/10.1016/j.jclepro.2022.131571 (2022).Article 
CAS 

Google Scholar 
Rajput, V. D. et al. Nano-biochar: A novel solution for sustainable agriculture and environmental remediation. Environ. Res. 210, 112891. https://doi.org/10.1016/j.envres.2022.112891 (2022).Article 
CAS 

Google Scholar 
Ding, Y. et al. Biochar to improve soil fertility. A review. Agron. Sustain. Dev. 36, 36. https://doi.org/10.1007/s13593-016-0372-z (2016).Article 
CAS 

Google Scholar 
Oni, B. A., Oziegbe, O. & Olawole, O. O. Significance of biochar application to the environment and economy. Ann. Agric. Sci. 64, 222–236. https://doi.org/10.1016/j.aoas.2019.12.006 (2019).Article 

Google Scholar 
He, E. et al. Two years of aging influences the distribution and lability of metal(loid)s in a contaminated soil amended with different biochars. Sci. Total Environ. 673, 245–253. https://doi.org/10.1016/j.scitotenv.2019.04.037 (2019).Article 
ADS 
CAS 

Google Scholar 
Netherway, P. et al. Phosphorus-rich biochars can transform lead in an urban contaminated soil. J. Environ. Qual. 48, 1091–1099. https://doi.org/10.2134/jeq2018.09.0324 (2019).Article 
CAS 

Google Scholar 
O’Connor, D. et al. Biochar application for the remediation of heavy metal polluted land: A review of in situ field trials. Sci. Total Environ. 619–620, 815–826. https://doi.org/10.1016/j.scitotenv.2017.11.132 (2018).Article 
ADS 
CAS 

Google Scholar 
Xu, X. et al. Effect of physicochemical properties of biochar from different feedstock on remediation of heavy metal contaminated soil in mining area. Surf. Interfaces 32, 102058. https://doi.org/10.1016/j.surfin.2022.102058 (2022).Article 
CAS 

Google Scholar 
Melo, L. C. A. et al. Sorption and desorption of cadmium and zinc in two tropical soils amended with sugarcane-straw-derived biochar. J. Soils Sediments 16, 226–234. https://doi.org/10.1007/s11368-015-1199-y (2016).Article 

Google Scholar 
Uchimiya, M., Chang, S. & Klasson, K. T. Screening biochars for heavy metal retention in soil: Role of oxygen functional groups. J. Hazard. Mater. 190, 432–441. https://doi.org/10.1016/j.jhazmat.2011.03.063 (2011).Article 
CAS 

Google Scholar 
Jatav, H. S. et al. Sustainable approach and safe use of biochar and its possible consequences. Sustainability https://doi.org/10.3390/su131810362 (2021).Article 

Google Scholar 
Varalta, F. & Sorvari, J. In Organic Waste Composting through Nexus Thinking: Practices, Policies, and Trends (eds Hettiarachchi, H. et al.) 213–232 (Springer International Publishing, 2020).Chapter 

Google Scholar 
Pinotti, L. et al. Recycling food leftovers in feed as opportunity to increase the sustainability of livestock production. J. Clean. Prod. 294, 126290. https://doi.org/10.1016/j.jclepro.2021.126290 (2021).Article 

Google Scholar 
Jafri, N., Wong, W. Y., Doshi, V., Yoon, L. W. & Cheah, K. H. A review on production and characterization of biochars for application in direct carbon fuel cells. Process Saf. Environ. Prot. 118, 152–166. https://doi.org/10.1016/j.psep.2018.06.036 (2018).Article 
CAS 

Google Scholar 
Jin, Y. et al. Characterization of biochars derived from various spent mushroom substrates and evaluation of their adsorption performance of Cu(II) ions from aqueous solution. Environ. Res. 196, 110323. https://doi.org/10.1016/j.envres.2020.110323 (2021).Article 
CAS 

Google Scholar 
Tomczyk, A., Sokołowska, Z. & Boguta, P. Biomass type effect on biochar surface characteristic and adsorption capacity relative to silver and copper. Fuel 278, 118168. https://doi.org/10.1016/j.fuel.2020.118168 (2020).Article 
CAS 

Google Scholar 
FAO. Food Outlook – Biannual Report on Global Food Markets: November 2020. Rome. Phytoremediation of copper-contaminated soil by Artemisia absinthium: comparative effect of chelating agents. Environmental Geochemistry and Health. (2020). https://doi.org/10.4060/cb1993enRussian-Statistical-Year-Book. Statistical handbook. P76 M., 2020 – 700 p. ISBN 978-5-89476-497-9 (2020).Cheng, C.-H., Lehmann, J., Thies, J. E. & Burton, S. D. Stability of black carbon in soils across a climatic gradient. J. Geophys. Res. Biogeosci. 113, 55. https://doi.org/10.1029/2007JG000642 (2008).Article 
CAS 

Google Scholar 
Singh, B. P., Cowie, A. L. & Smernik, R. J. Biochar carbon stability in a clayey soil as a function of feedstock and pyrolysis temperature. Environ. Sci. Technol. 46, 11770–11778. https://doi.org/10.1021/es302545b (2012).Article 
ADS 
CAS 

Google Scholar 
He, Y. et al. Effects of biochar application on soil greenhouse gas fluxes: A meta-analysis. GCB Bioenergy 9, 743–755. https://doi.org/10.1111/gcbb.12376 (2017).Article 
CAS 

Google Scholar 
Janu, R. et al. Biochar surface functional groups as affected by biomass feedstock, biochar composition and pyrolysis temperature. Carbon Resour. Convers. 4, 36–46. https://doi.org/10.1016/j.crcon.2021.01.003 (2021).Article 
CAS 

Google Scholar 
Tan, X. et al. Application of biochar for the removal of pollutants from aqueous solutions. Chemosphere 125, 70–85. https://doi.org/10.1016/j.chemosphere.2014.12.058 (2015).Article 
ADS 
CAS 

Google Scholar 
Ni, B.-J. et al. Competitive adsorption of heavy metals in aqueous solution onto biochar derived from anaerobically digested sludge. Chemosphere 219, 351–357. https://doi.org/10.1016/j.chemosphere.2018.12.053 (2019).Article 
ADS 
CAS 

Google Scholar 
Park, J.-H. et al. Competitive adsorption of heavy metals onto sesame straw biochar in aqueous solutions. Chemosphere 142, 77–83. https://doi.org/10.1016/j.chemosphere.2015.05.093 (2016).Article 
ADS 
CAS 

Google Scholar 
Methodological-Guidelines. Methodological guidelines for the determination of heavy metals in the soils of agricultural land and crop production – M., TSINAO, 61 (1992)Zhang, A., Li, X., Xing, J. & Xu, G. Adsorption of potentially toxic elements in water by modified biochar: A review. J. Environ. Chem. Eng. 8, 104196. https://doi.org/10.1016/j.jece.2020.104196 (2020).Article 
CAS 

Google Scholar 
Avramiotis, E., Frontistis, Z., Manariotis, I. D., Vakros, J. & Mantzavinos, D. On the performance of a sustainable rice husk biochar for the activation of persulfate and the degradation of antibiotics. Catalysts 11, 1303 (2021).Article 
CAS 

Google Scholar 
Maiti, S., Dey, S., Purakayastha, S. & Ghosh, B. Physical and thermochemical characterization of rice husk char as a potential biomass energy source. Biores. Technol. 97, 2065–2070. https://doi.org/10.1016/j.biortech.2005.10.005 (2006).Article 
CAS 

Google Scholar 
Herrera, K., Morales, L. F., Tarazona, N. A., Aguado, R. & Saldarriaga, J. F. Use of biochar from rice husk pyrolysis: Part A: Recovery as an adsorbent in the removal of emerging compounds. ACS Omega 7, 7625–7637. https://doi.org/10.1021/acsomega.1c06147 (2022).Article 
CAS 

Google Scholar 
Szewczuk-Karpisz, K., Tomczyk, A., Grygorczuk-Płaneta, K. & Naveed, S. Rhizobium leguminosarum bv. trifolii exopolysaccharide and sunflower husk biochar as factors affecting immobilization of both tetracycline and Cd2+ ions on soil solid phase. J. Soils Sediments 22, 2620–2639. https://doi.org/10.1007/s11368-022-03255-3 (2022).Article 
CAS 

Google Scholar 
Hubetska, T. S., Kobylinska, N. G. & García, J. R. Sunflower biomass power plant by-products: Properties and its potential for water purification of organic pollutants. J. Anal. Appl. Pyrolysis 157, 105237. https://doi.org/10.1016/j.jaap.2021.105237 (2021).Article 
CAS 

Google Scholar 
Braghiroli, F. L. et al. The influence of pilot-scale pyro-gasification and activation conditions on porosity development in activated biochars. Biomass Bioenerg. 118, 105–114. https://doi.org/10.1016/j.biombioe.2018.08.016 (2018).Article 
CAS 

Google Scholar 
Braghiroli, F. L. et al. The conversion of wood residues, using pilot-scale technologies, into porous activated biochars for supercapacitors. J. Porous Mater. 27, 537–548. https://doi.org/10.1007/s10934-019-00823-w (2020).Article 
CAS 

Google Scholar 
Boraah, N., Chakma, S. & Kaushal, P. Attributes of wood biochar as an efficient adsorbent for remediating heavy metals and emerging contaminants from water: A critical review and bibliometric analysis. J. Environ. Chem. Eng. 10, 107825. https://doi.org/10.1016/j.jece.2022.107825 (2022).Article 
CAS 

Google Scholar 
Phillips, C. L. et al. Towards predicting biochar impacts on plant-available soil nitrogen content. Biochar 4, 9. https://doi.org/10.1007/s42773-022-00137-2 (2022).Article 
CAS 

Google Scholar 
Sun, L. & Gong, K. Silicon-based materials from rice husks and their applications. Ind. Eng. Chem. Res. 40, 5861–5877. https://doi.org/10.1021/ie010284b (2001).Article 
CAS 

Google Scholar 
Islam, T. et al. Synthesis of rice husk-derived magnetic biochar through liquefaction to adsorb anionic and cationic dyes from aqueous solutions. Arab. J. Sci. Eng. 46, 233–246. https://doi.org/10.1007/s13369-020-04537-z (2021).Article 
CAS 

Google Scholar 
Mohan, D. et al. Biochar production and applications in soil fertility and carbon sequestration – a sustainable solution to crop-residue burning in India. RSC Adv. 8, 508–520. https://doi.org/10.1039/C7RA10353K (2018).Article 
ADS 
CAS 

Google Scholar 
Li, F. et al. Preparation and characterization of biochars from Eichornia crassipes for cadmium removal in aqueous solutions. PLoS ONE 11, e0148132. https://doi.org/10.1371/journal.pone.0148132 (2016).Article 
CAS 

Google Scholar 
Song, H. et al. Potential of novel biochars produced from invasive aquatic species outside food chain in removing ammonium nitrogen: Comparison with conventional biochars and clinoptilolite. Sustainability https://doi.org/10.3390/su11247136 (2019).Article 

Google Scholar 
Yang, G. et al. Effects of pyrolysis temperature on the physicochemical properties of biochar derived from vermicompost and its potential use as an environmental amendment. RSC Adv. 5, 40117–40125. https://doi.org/10.1039/C5RA02836A (2015).Article 
ADS 
CAS 

Google Scholar 
Enders, A., Hanley, K., Whitman, T., Joseph, S. & Lehmann, J. Characterization of biochars to evaluate recalcitrance and agronomic performance. Bioresour. Technol. 114, 644–653. https://doi.org/10.1016/j.biortech.2012.03.022 (2012).Article 
CAS 

Google Scholar 
Zhang, Y., Wang, J. & Feng, Y. The effects of biochar addition on soil physicochemical properties: A review. CATENA 202, 105284. https://doi.org/10.1016/j.catena.2021.105284 (2021).Article 
CAS 

Google Scholar 
Özçimen, D. & Ersoy-Meriçboyu, A. Characterization of biochar and bio-oil samples obtained from carbonization of various biomass materials. Renew. Energy 35, 1319–1324. https://doi.org/10.1016/j.renene.2009.11.042 (2010).Article 
CAS 

Google Scholar 
Lin, Q. et al. Effects of biochar-based materials on the bioavailability of soil organic pollutants and their biological impacts. Sci. Total Environ. 826, 153956. https://doi.org/10.1016/j.scitotenv.2022.153956 (2022).Article 
ADS 
CAS 

Google Scholar 
Yang, H. et al. Thermogravimetric analysis−fourier transform infrared analysis of palm oil waste pyrolysis. Energy Fuels 18, 1814–1821. https://doi.org/10.1021/ef030193m (2004).Article 
CAS 

Google Scholar 
Pasangulapati, V. et al. Effects of cellulose, hemicellulose and lignin on thermochemical conversion characteristics of the selected biomass. Biores. Technol. 114, 663–669. https://doi.org/10.1016/j.biortech.2012.03.036 (2012).Article 
CAS 

Google Scholar 
Kim, P. et al. Surface functionality and carbon structures in lignocellulosic-derived biochars produced by fast pyrolysis. Energy Fuels 25, 4693–4703. https://doi.org/10.1021/ef200915s (2011).Article 
CAS 

Google Scholar 
Keiluweit, M., Nico, P. S., Johnson, M. G. & Kleber, M. dynamic molecular structure of plant biomass-derived black carbon (biochar). Environ. Sci. Technol. 44, 1247–1253. https://doi.org/10.1021/es9031419 (2010).Article 
ADS 
CAS 

Google Scholar 
Wijeyawardana, P. et al. Removal of Cu, Pb and Zn from stormwater using an industrially manufactured sawdust and paddy husk derived biochar. Environ. Technol. Innov. 28, 102640. https://doi.org/10.1016/j.eti.2022.102640 (2022).Article 
CAS 

Google Scholar 
Kołodyńska, D., Krukowska, J. & Thomas, P. Comparison of sorption and desorption studies of heavy metal ions from biochar and commercial active carbon. Chem. Eng. J. 307, 353–363. https://doi.org/10.1016/j.cej.2016.08.088 (2017).Article 
CAS 

Google Scholar 
Uchimiya, M. et al. Immobilization of heavy metal ions (CuII, CdII, NiII, and PbII) by broiler litter-derived biochars in water and soil. J. Agric. Food Chem. 58, 5538–5544. https://doi.org/10.1021/jf9044217 (2010).Article 
CAS 

Google Scholar 
Misono, M., Ochiai, E. I., Saito, Y. & Yoneda, Y. A new dual parameter scale for the strength of lewis acids and bases with the evaluation of their softness. J. Inorg. Nucl. Chem. 29, 2685–2691. https://doi.org/10.1016/0022-1902(67)80006-X (1967).Article 
CAS 

Google Scholar 
McBride, M. B. Environmental Chemistry of Soils (Oxford University Press, 1994).
Google Scholar 
Basta, N. T. & Tabatabai, M. A. Effect of cropping systems on adsorption of metals by soils: III. Competitive adsorption1. Soil Sci. 153, 331–337 (1992).Article 
ADS 
CAS 

Google Scholar 
Sposito, G. The Chemistry of Soils (Oxford University Press, 2016).Bauer, T. V. et al. Application of XAFS and XRD methods for describing the copper and zinc adsorption characteristics in hydromorphic soils. Environ. Geochem. Health 44, 335–347. https://doi.org/10.1007/s10653-020-00773-2 (2022).Article 
CAS 

Google Scholar 
Abd-Elfattah, A. L. Y. & Wada, K. Adsorption of lead, copper, zinc, cobalt, and cadmium by soils that differ in cation-exchange materials. J. Soil Sci. 32, 271–283. https://doi.org/10.1111/j.1365-2389.1981.tb01706.x (1981).Article 
CAS 

Google Scholar 
Etesami, H., Fatemi, H. & Rizwan, M. Interactions of nanoparticles and salinity stress at physiological, biochemical and molecular levels in plants: A review. Ecotoxicol. Environ. Saf. 225, 112769. https://doi.org/10.1016/j.ecoenv.2021.112769 (2021).Article 
CAS 

Google Scholar 
Soria, R. I., Rolfe, S. A., Betancourth, M. P. & Thornton, S. F. The relationship between properties of plant-based biochars and sorption of Cd(II), Pb(II) and Zn(II) in soil model systems. Heliyon 6, e05388. https://doi.org/10.1016/j.heliyon.2020.e05388 (2020).Article 

Google Scholar 
Alfarra, A., Frackowiak, E. & Béguin, F. The HSAB concept as a means to interpret the adsorption of metal ions onto activated carbons. Appl. Surf. Sci. 228, 84–92. https://doi.org/10.1016/j.apsusc.2003.12.033 (2004).Article 
ADS 
CAS 

Google Scholar 
Hu, J., Zhou, X., Shi, Y., Wang, X. & Li, H. Enhancing biochar sorption properties through self-templating strategy and ultrasonic fore-modified pre-treatment: Characteristic, kinetic and mechanism studies. Sci. Total Environ. 769, 144574. https://doi.org/10.1016/j.scitotenv.2020.144574 (2021).Article 
ADS 
CAS 

Google Scholar 
Ward, J., Rasul, M. G. & Bhuiya, M. M. K. Energy recovery from biomass by fast pyrolysis. Proced. Eng. 90, 669–674. https://doi.org/10.1016/j.proeng.2014.11.791 (2014).Article 
CAS 

Google Scholar 
Al-Wabel, M. I., Al-Omran, A., El-Naggar, A. H., Nadeem, M. & Usman, A. R. A. Pyrolysis temperature induced changes in characteristics and chemical composition of biochar produced from conocarpus wastes. Biores. Technol. 131, 374–379. https://doi.org/10.1016/j.biortech.2012.12.165 (2013).Article 
CAS 

Google Scholar 
Calvelo Pereira, R. et al. Contribution to characterisation of biochar to estimate the labile fraction of carbon. Org. Geochem. 42, 1331–1342. https://doi.org/10.1016/j.orggeochem.2011.09.002 (2011).Article 
CAS 

Google Scholar 
Vorob’eva, L. A. Theory and Practice Chemical Analysis of Soils (GEOS Press, Moscow, 2006).
Google Scholar 
Pinskii, D. L. et al. Copper adsorption by chernozem soils and parent rocks in Southern Russia. Geochem. Int. 56, 266–275. https://doi.org/10.1134/S0016702918030072 (2018).Article 
CAS 

Google Scholar 
Wang, Q., Wang, B., Lee, X., Lehmann, J. & Gao, B. Sorption and desorption of Pb(II) to biochar as affected by oxidation and pH. Sci. Total Environ. 634, 188–194. https://doi.org/10.1016/j.scitotenv.2018.03.189 (2018).Article 
ADS 
CAS 

Google Scholar 
Pourret, O. & Houben, D. Characterization of metal binding sites onto biochar using rare earth elements as a fingerprint. Heliyon 4, e00543. https://doi.org/10.1016/j.heliyon.2018.e00543 (2018).Article 

Google Scholar 
Huang, L. et al. High-resolution insight into the competitive adsorption of heavy metals on natural sediment by site energy distribution. Chemosphere 197, 411–419. https://doi.org/10.1016/j.chemosphere.2018.01.056 (2018).Article 
ADS 
MathSciNet 
CAS 

Google Scholar 
Ming, H. et al. Competitive sorption of cadmium and zinc in contrasting soils. Geoderma 268, 60–68. https://doi.org/10.1016/j.geoderma.2016.01.021 (2016).Article 
ADS 
CAS 

Google Scholar 
Musso, T. B., Parolo, M. E., Pettinari, G. & Francisca, F. M. Cu(II) and Zn(II) adsorption capacity of three different clay liner materials. J. Environ. Manag. 146, 50–58. https://doi.org/10.1016/j.jenvman.2014.07.026 (2014).Article 
CAS 

Google Scholar 
Cui, H. et al. Immobilization of Cu and Cd in a contaminated soil: One- and four-year field effects. J. Soils Sediments 14, 1397–1406. https://doi.org/10.1007/s11368-014-0882-8 (2014).Article 
CAS 

Google Scholar 
Elbana, T. A. et al. Freundlich sorption parameters for cadmium, copper, nickel, lead, and zinc for different soils: Influence of kinetics. Geoderma 324, 80–88. https://doi.org/10.1016/j.geoderma.2018.03.019 (2018).Article 
ADS 
CAS 

Google Scholar  More

Gillooly, J. F., Brown, J. H., West, G. B., Savage, V. M. & Charnov, E. L. Effects of size and temperature on metabolic rate. Science 293, 2248–2251 (2001).Article 
CAS 

Google Scholar 
Davidson, E. A. & Janssens, I. A. Temperature sensitivity of soil carbon decomposition and feedbacks to climate change. Nature 440, 165–173 (2006).Article 
CAS 

Google Scholar 
Lopez-Urrutia, A., San Martin, E., Harris, R. P. & Irigoien, X. Scaling the metabolic balance of the oceans. Proc. Natl Acad. Sci. USA 103, 8739–8744 (2006).Article 
CAS 

Google Scholar 
Yvon-Durocher, G. et al. Reconciling the temperature dependence of respiration across timescales and ecosystem types. Nature 487, 472–476 (2012).Article 
CAS 

Google Scholar 
Crowther, T. W. et al. Quantifying global soil carbon losses in response to warming. Nature 540, 104–108 (2016).Article 
CAS 

Google Scholar 
Bar-On, Y. M., Phillips, R. & Milo, R. The biomass distribution on Earth. Proc. Natl Acad. Sci. USA 115, 6506–6511 (2018).Article 
CAS 

Google Scholar 
Rivkin, R. B. & Legendre, L. Biogenic carbon cycling in the upper ocean: effects of microbial respiration. Science 291, 2398–2400 (2001).Article 
CAS 

Google Scholar 
Friedlingstein, P. et al. Uncertainties in CMIP5 climate projections due to carbon cycle feedbacks. J. Clim. 27, 511–526 (2014).Article 

Google Scholar 
Smith, T. P. et al. Community-level respiration of prokaryotic microbes may rise with global warming. Nat. Commun. 10, 5124 (2019).Article 

Google Scholar 
Antwis, R. E. et al. Fifty important research questions in microbial ecology. FEMS Microbiol. Ecol. 93, fix044 (2017).Bardgett, R. D., Freeman, C. & Ostle, N. J. Microbial contributions to climate change through carbon cycle feedbacks. ISME J. 2, 805–814 (2008).Article 
CAS 

Google Scholar 
Enquist, B. J. et al. Scaling from traits to ecosystems: developing a general trait driver theory via integrating trait-based and metabolic scaling theories. Adv. Ecol. Res. 52, 249–318 (2015).Article 

Google Scholar 
Allen, A. P., Gillooly, J. F. & Brown, J. H. Linking the global carbon cycle to individual metabolism. Funct. Ecol. 19, 202–213 (2005).Article 

Google Scholar 
Schramski, J. R., Dell, A. I., Grady, J. M., Sibly, R. M. & Brown, J. H. Metabolic theory predicts whole-ecosystem properties. Proc. Natl Acad. Sci. USA 112, 2617–2622 (2015).Article 
CAS 

Google Scholar 
Alster, C. J., Koyama, A., Johnson, N. G., Wallenstein, M. D. & von Fischer, J. C. Temperature sensitivity of soil microbial communities: an application of macromolecular rate theory to microbial respiration. J. Geophys. Res. Biogeosci. 121, 1420–1433 (2016).Article 

Google Scholar 
Yvon-Durocher, G. et al. Five years of experimental warming increases the biodiversity and productivity of phytoplankton. PLoS Biol. 13, e1002324 (2015).Article 

Google Scholar 
Garzke, J., Connor, S. J., Sommer, U. & O’Connor, M. I. Trophic interactions modify the temperature dependence of community biomass and ecosystem function. PLoS Biol. 17, e2006806 (2019).Foster, K. R. & Bell, T. Competition, not cooperation, dominates interactions among culturable microbial species. Curr. Biol. 22, 1845–1850 (2012).Article 
CAS 

Google Scholar 
Coyte, K. Z., Schluter, J. & Foster, K. R. The ecology of the microbiome: networks, competition, and stability. Science 350, 663–666 (2015).Article 
CAS 

Google Scholar 
Machado, D. et al. Polarization of microbial communities between competitive and cooperative metabolism. Nat. Ecol. Evol. 5, 195–203 (2021).Article 

Google Scholar 
Bradford, M. A. et al. Cross-biome patterns in soil microbial respiration predictable from evolutionary theory on thermal adaptation. Nat. Ecol. Evol. 3, 223–231 (2019).Article 

Google Scholar 
Garcia-Martin, E. E., McNeill, S., Serret, P. & Leakey, R. J. G. Plankton metabolism and bacterial growth efficiency in offshore waters along a latitudinal transect between the UK and Svalbard. Deep Sea Res. I 92, 141–151 (2014).Article 
CAS 

Google Scholar 
Davidson, E. A., Richardson, A. D., Savage, K. E. & Hollinger, D. Y. A distinct seasonal pattern of the ratio of soil respiration to total ecosystem respiration in a spruce-dominated forest. Glob. Change Biol. 12, 230–239 (2006).Article 

Google Scholar 
Dutkiewicz, S., Follows, M. J. & Bragg, J. G. Modeling the coupling of ocean ecology and biogeochemistry. Glob. Biogeochem. Cycles 23, GB4017 (2009).Article 

Google Scholar 
Follows, M. J., Dutkiewicz, S., Ward, B. & Follett, C. in Microbial Ecology of the Oceans 3rd edn (eds Gasol, J. & Kirchman, D.) Ch. 12 (John Wiley, 2018).Letten, A. D. & Stouffer, D. B. The mechanistic basis for higher-order interactions and non-additivity in competitive communities. Ecol. Lett. 22, 423–436 (2019).Article 

Google Scholar 
Grilli, J., Barabás, G., Michalska-Smith, M. J. & Allesina, S. Higher-order interactions stabilize dynamics in competitive network models. Nature 548, 210–213 (2017).Article 
CAS 

Google Scholar 
Maynard, D. S., Crowther, T. W. & Bradford, M. A. Competitive network determines the direction of the diversity–function relationship. Proc. Natl Acad. Sci. USA 114, 11464–11469 (2017).Article 
CAS 

Google Scholar 
Fiegna, F., Moreno-Letelier, A., Bell, T. & Barraclough, T. G. Evolution of species interactions determines microbial community productivity in new environments. ISME J. 9, 1235–1245 (2015).Article 

Google Scholar 
Lawrence, D. et al. Species interactions alter evolutionary responses to a novel environment. PLoS Biol. 10, e1001330 (2012).Article 
CAS 

Google Scholar 
Harcombe, W. R., Chacón, J. M., Adamowicz, E. M., Chubiz, L. M. & Marx, C. J. Evolution of bidirectional costly mutualism from byproduct consumption. Proc. Natl Acad. Sci. USA 115, 12000–12004 (2018).Article 
CAS 

Google Scholar 
Goldford, J. E. et al. Emergent simplicity in microbial community assembly. Science 361, 469–474 (2018).Article 
CAS 

Google Scholar 
Yvon-Durocher, G. et al. Methane fluxes show consistent temperature dependence across microbial to ecosystem scales. Nature 507, 488–491 (2014).Article 
CAS 

Google Scholar 
Fox, J. W. & Harpole, W. S. Revealing how species loss affects ecosystem function: the trait-based price equation partition. Ecology 89, 269–279 (2008).Article 

Google Scholar 
Kontopoulos, D., Smith, T. P., Barraclough, T. G. & Pawar, S. Adaptive evolution shapes the present-day distribution of the thermal sensitivity of population growth rate. PLoS Biol. 18, e3000894 (2020).Article 
CAS 

Google Scholar 
Wilson, W. G. & Lundberg, P. Biodiversity and the Lotka–Volterra theory of species interactions: open systems and the distribution of logarithmic densities. Proc. R. Soc. Lond. B 271, 1977–1984 (2004).Article 

Google Scholar 
Rossberg, A. G. in Food Webs and Biodiversity 181–191 (John Wiley & Sons, 2013).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 
CAS 

Google Scholar 
Garcia, F. C., Bestion, E., Warfield, R. & Yvon-Durocher, G. Changes in temperature alter the relationship between biodiversity and ecosystem functioning. Proc. Natl Acad. Sci. USA 115, 10989–10999 (2018).Article 
CAS 

Google Scholar 
Padfield, D., O’Sullivan, H. & Pawar, S. rTPC and nls.multstart: a new pipeline to fit thermal performance curves in R. Methods Ecol. Evol. 12, 1138–1143 (2021).Article 

Google Scholar  More

Study areaThe Qinghai–Tibet Plateau (26°00′-39°47′N, 73°19′-104°47′E), one of the most important pastoral areas in the world, straddles the southwest regions of China, and it includes 244 counties, which belong to six provinces: Tibet, Qinghai, Xinjiang, Gansu, Sichuan, and Yunnan. It is characterized by rich natural grassland resources, including desert steppes, alpine steppes, and alpine meadows (Fig. 1a). The grassland areas account for over 56% of this region34. The grassland plays a vital role in providing regional and national animal husbandry products and fodder35, which enables the local herders to obtain almost all of the resources required for survival36. The grazing density distribution is extremely unbalanced (Fig. 1a) owing to the high spatial heterogeneity of economic development (Fig. 1b-1) and grassland production (Fig. 1b-2), resulting from the differences in resources and environmental factors37. Over the past few decades, there has been a significant change in the number of livestock animals, and the number of sheep exceeded 160 million by 2020. Therefore, it is urgent to obtain a high-resolution gridded grazing dataset for its evaluating spatiotemporal changes and coordinating the relationship between human beings and the grassland ecosystem.Fig. 1Location of the Qinghai–Tibet Plateau: (a) grassland type and distribution, and grazing density (GD) in 244 counties; (b) spatial heterogeneity of economic development (ED) and grassland production (GP) in 244 counties. GD, ED, and GP are represented by sheep unit per grassland area per county (SU/hm2), human footprint index per pixel (HF/pixel) per county, and net primary production per grassland area per county (gC/m2), respectively.Full size imageFig. 2Methodological framework for grazing spatialization.Full size imageMethodological frameworkWe developed a methodological framework for high-resolution gridded grazing dataset mapping. The framework mainly includes four parts: (i) identifying features affecting grazing, (ii) extracting theoretical suitable grazing areas, (iii) building grazing spatialization model, and (iv) correcting the grazing spatialization dataset. Each step is explained in more detail below (Fig. 2).Step 1: Identifying features affecting grazingGrazing activities are affected by the spatial heterogeneity of resources and environmental factors, regulated by the grazing behavior of herders and the foraging behavior of herds, and restricted by ecological protection policies. Therefore, the specific implications of the 14 influencing factors from the above four aspects are presented in Table 1. These factors are necessary for spatializing the county-level grazing data.Table 1 The identified features affecting grazing.Full size tableStep 2: Extracting theoretical suitable grazing areasThe decision tree approach38 was adopted to extract the theoretical suitable grazing areas for further grazing spatialization (step 2 in Fig. 2). First, the potential grazing area was identified according to the boundary of the grassland ecosystem, because grazing behavior only occurs in the grassland. Then, the unsuitable areas for grazing, i.e., extremely-high-altitude areas and areas adjacent to towns, were removed from the potential grazing area stepwise. The areas strictly prohibited for grazing, i.e., the core areas of national nature reserves39 within grassland areas, were also deemed unsuitable for grazing. Finally, the extracted areas were the theoretically suitable grazing areas.Step 3: Building grazing spatialization model(i) Extracting cross-scale feature (CSFs)In the traditional method, the spatial resolution of the training data (i.e., the average value at the administrative level) differs from that of the predicting data (i.e., the value at the pixel level), and the trained model can only capture the characteristics within the training data. However, the extreme value of the predicting data inevitably exceeds the range of the training data, which can result in underestimation in these parts40. To reduce these mismatches, we built an improved method for CSFs extraction (Fig. 2, first part of step 3).First, the census grazing data are simply distributed from county level to pixel level using the weight of the absolute disturbance (AD) index as Eq. (1). The AD index is measured by Mahalanobis distance using Eq. (2), which is calculated according to the deviation between the potential and observed normalized difference vegetation index (NDVI) values22. Second, the distributed grazing data are graded via the hierarchical clustering method, and the optimal number of the group can be determined using the Davies–Bouldin index (DBI)41 as Eq. (3), an index for evaluating the quality of clustering algorithm. The smaller the DBI, the smaller the distance within each group. Therefore, the DBI can be used to select the best similar values to minimize the deviation within each group. Finally, we can obtain the scope of the groups within each county using the above two steps and obtain the average value of all independent variables and the dependent variable accordingly. As expected, we can decompose the average value at the county level (traditional features in Fig. 2) into the average value at the group level (improved features in Fig. 2).$$S{U}_{i}=S{U}_{j}^{C}frac{{w}_{A{D}_{i}}}{{w}_{A{D}_{j}}}$$
(1)
where SUi and (S{U}_{j}^{C}) are the grazing value for pixel i and the census grazing value for county j; ({w}_{A{D}_{i}}) is the weight of the AD index for pixel i and ({w}_{A{D}_{j}}) represents the summed weight of the AD index values for all pixels in county j.$$begin{array}{cll}A{D}_{i} & = & sqrt{{({D}_{i}-u)}^{T}co{v}^{-1}({D}_{i}-u)}\ {D}_{i} & = & NDV{I}_{i}^{T}-NDV{I}_{i}^{P}end{array}$$
(2)
where ADi is the AD index for pixel i; the vector composed of its observed NDVI (left(NDV{I}_{i}^{T}right)) and potential NDVI (left(NDV{I}_{i}^{P}right)) time-series data could be considered as two points in the feature space for pixel i, and Di and u are the difference and the mean value of the vector, respectively; cov is the covariance matrix.$$DB{I}_{k}=frac{1}{k}{sum }_{x=1}^{k}ma{x}_{yne x}left(frac{overline{{a}_{x}}+overline{{a}_{y}}}{left|{delta }_{x}-{delta }_{y}right|}right)$$
(3)
where DBIk is the DBI coefficient when the cluster number is k; (overline{{a}_{x}}) and (overline{{a}_{y}}) are the average distances of the group xth and the group yth, respectively; δx and δy are the center distance of the group xth and the group yth, respectively.Different from the traditional method, our method can decompose features into multiple features using the grading AD index. The differences among counties will not be easily averaged out. Moreover, our method is less affected by scale mismatch and can be transferred to cross-scale modeling26.(ii) Building RF model with partitioningA single model cannot accurately obtain the variation information of the Qinghai–Tibet Plateau with high spatial heterogeneity. The partition model, a widely used method for estimating population distribution and others42,43, can be incorporated into the proposed model to improve its performance. The thresholds (0.43, 0.35 and 0.21 SU/hm2), determined according to the theoretical livestock carrying capacity (equation S1), were calculated and used to separate independent variables and dependent variable for each grassland types: alpine meadow, alpine steppe and alpine desert steppe (see Section 6.1 for details). Then, the RF models were established, and the training and testing samples were randomly divided in the proportion of 3:1. It is notable that transforming the response variable using natural log prior to RF model fitting is necessary to achieve higher prediction accuracies44. Finally, the independent variables at the pixel level were inputted into the two trained RF models, and the corresponding grid grazing dataset was output by combining the two results (Fig. 2, second part of step 3).(iii) Validating the accuracy of the methodsThe performance of the grazing spatialization model was evaluated through a comparison of the predicted value with census value26. Accuracy validation indexes, including coefficients of determination (R2), root mean square error (RMSE), and mean absolute error (MAE), were used to evaluate the performances of the proposed RF-based models (Table 2), as presented in Eq. (4).$$begin{array}{ccc}{R}^{2} & = & 1-frac{{sum }_{j=1}^{N}{left(S{U}_{j}^{C}-S{U}_{j}^{P}right)}^{2}}{{sum }_{j=1}^{N}{left(S{U}_{j}^{C}-overline{S{U}^{C}}right)}^{2}}\ RMSE & = & sqrt{frac{{sum }_{j=1}^{N}{left(S{U}_{j}^{C}-S{U}_{j}^{P}right)}^{2}}{N}}\ MAE & = & frac{{sum }_{j=1}^{N}| S{U}_{j}^{C}-S{U}_{j}^{P}| }{N}end{array}$$
(4)
where (S{U}_{j}^{C}) and (S{U}_{j}^{P}) are the census grazing value and the predicted grazing value for county j, respectively; (overline{S{U}^{C}}) is the average census data for all counties; and N is the number of all counties.Table 2 The proposed methods and their descriptions.Full size tableStep 4: Correcting grazing spatialization dataset(i) Correcting residuals of datasetCorrecting residuals is necessary to obtain datasets with higher accuracy45,46, because propagating the cross-scale relationship in the RF models will inevitably generate errors47. The residuals, calculated by the difference between the average census grazing and predicted grazing values at the administrative level, were used to calibrate the errors related to all pixels within this county. The revised dataset after residual correction is the final product provided in this study. The residual correction method is expressed by Eq. (5), and the process is shown in the fourth step in Fig. 2.$$S{U}_{i}^{RP}=S{U}_{i}^{P}+{R}_{j}$$
(5)
where (S{U}_{i}^{RP}) denotes the predicted grazing value revised by the residuals for pixel i, (S{U}_{i}^{P}) denotes the predicted grazing for pixel i, and Rj denotes the residuals calculated from the difference between census grazing and predicted grazing data for county j.(ii) Validating the accuracy of datasetTwo goodness-of-fit indexes were used to validate the consistency of spatial distribution and the temporal trend between predicted grazing data and census grazing data. Generally, the coefficient of determination (R2), defined in Eq. (4), is used to verify the consistency of spatial distribution, and the Nash–Sutcliffe efficiency (NSE, Eq. (6)) is used to verify the consistency of temporal trend. An index value closer to 1 corresponds to a more accurate dataset. Meanwhile, we also collected field grazing data from 56 sites to further validate the spatial accuracy of the dataset, and it measured using the R2 in Eq. (4).$$NSE=1-frac{{sum }_{t=1}^{T}{left(S{U}_{t}^{RP}-S{U}_{t}^{C}right)}^{2}}{{sum }_{t=1}^{T}{left(S{U}_{t}^{C}-overline{S{U}^{{C}^{{prime} }}}right)}^{2}}$$
(6)
where (S{U}_{t}^{RP}) and (S{U}_{t}^{C}) are the predicted grazing value revised by residuals and the census grazing value of all counties in year t, respectively; (overline{S{U}^{{C}^{{prime} }}}) is the average census grazing value of all years; and T is the number of time steps.(iii) Identifying uncertainties associated with datasetThe uncertainties associated with the dataset originate from the following two aspects: First, the unreasonableness of our method, owing to the errors related to cross-scale modeling or the inappropriate selection of influencing factors, is an important source of uncertainties. Second, the incompleteness of auxiliary variables also introduces uncertainties. In this instance, grassland-free areas are not accurately identified in some counties, but livestock animals are raised in these counties. These counties have no effective value for livestock density prediction. Overall, the uncertainties can be identified in terms of the mean relative error (MRE) in Eq. (7).$$MRE=frac{{sum }_{j=1}^{N}left|frac{S{U}_{j}^{C}-S{U}_{j}^{RP}}{S{U}_{j}^{C}}right|}{N}ast 100 % $$
(7)
where (S{U}_{j}^{C}) is the census grazing value for county j, (S{U}_{j}^{RP}) is the predicted grazing value revised by residuals for county j, and N is the number of counties.Data sourceCensus grazing data at county levelEight types of livestock, namely cattle, yaks, horses, donkeys, mules, camels, goats, and sheep, were considered according to the regional characteristics, and livestock stocking quantity at the end of year for each county can be determined from statistical yearbooks. However, the numbers of livestock at the county level for some years between 1982 and 2015 were not recorded. The missing data were indirectly approximated from city- or provincial-level data (e.g., interpolation using their temporal trends). Each type of livestock stocking quantity was converted into standard sheep unit (SU) according to the national standards using Eq. (8)48, namely the calculation of rangeland carrying capacity (NY/T 635-2015). Of the 244 counties of the Qinghai–Tibet Plateau, only 242 counties were considered, as the census grazing data for the other 2 counties were unavailable. The unit of grazing statistics data at the county level is defined as SU per county per year (SU·county−1·year−1).$$begin{array}{l}SU={N}_{sheep}+0.8times {N}_{goats}+5times {N}_{cattle}+5times {N}_{yaks+}+\ 6times {N}_{horses}+3times {N}_{donkeys}+6times {N}_{mules}+7times {N}_{camels}end{array}$$
(8)
where Nsheep, Ngoats, Ncattle, Nyaks, Nhorses, Ndonkeys, Nmules, Ncamels are the number of sheep, goats, cattle, yaks, horses, donkeys, mules, and camels at the year-end, respectively. SU denotes the standard sheep unit (SU·county−1·year−1).Data of grazing influencing factors at pixel levelThe types of features affecting grazing were obtained from the first step described in Methods, and the detailed information, such as original spatiotemporal resolution, format, and source, is shown in Table 3. The format (i.e., GeoTIFF), spatial resolution (i.e., 0.083°), and the number of rows and columns of the gridded features were leveraged to further produce a high-resolution grazing dataset.Table 3 Data source of grazing influence factors.Full size table More

Now consider a generalized model in which the species interactions are heterogeneous. A natural way of introducing heterogeneity in the system is by having a species diversify into subpopulations with different interaction strengths12,13,14,15. This way of modeling heterogeneity is useful as it can describe different kinds of heterogeneity. For example, the subpopulations could represent polymorphic traits that are genetically determined or result from plastic response to heterogeneous environments. A population could also be divided into local subpopulations in different spatial patches, which can migrate between patches and may face different local predators. We can also model different behavioral modes as subpopulations that, for instance, spend more time foraging for food or hiding from predators. We study several kinds of heterogeneity after we introduce a common mathematical framework. By studying these different scenarios using variants of the model, we show that our main results are not sensitive to the details of the model.We focus on the simple case where only the prey species splits into two types, (C_1) and (C_2), as illustrated in Fig. 1b. The situation is interesting when predator A consumes (C_1) more readily than predator B and B consumes (C_2) more readily than A (i.e., (a_1 / a_0 > b_1 / b_0) and (b_2 / b_0 > a_2 / a_0), which is equivalent to the condition that the nullclines of A and B cross, see section “Resources competition and nullcline analysis”). The arrows between (C_1) and (C_2) in Fig. 1b represent the exchange of individuals between the two subpopulations, which can happen by various mechanisms considered below. Such exchange as well as intraspecific competition between (C_1) and (C_2) result from the fact that the two prey types remain a single species.The system is now described by an enlarged Lotka-Volterra system with four variables, A, B, (C_1), and (C_2): $$begin{aligned} dot{A}&= varepsilon _A ,alpha _{A1} , A , C_1 + alpha _{A2} , A , C_2 – beta _A , A end{aligned}$$
(3a)
$$begin{aligned} dot{B}&= varepsilon _B , alpha _{B1} , B , C_1 + alpha _{B2} , B , C_2 – beta _B , B end{aligned}$$
(3b)
$$begin{aligned} dot{C_1}&= C_1 , (beta _C – alpha _{CC} , C)-alpha _{A1} , C_1 A-alpha _{B1} , C_1 B – sigma _1 , C_1 + sigma _2 , C_2 end{aligned}$$
(3c)
$$begin{aligned} dot{C_2}&= C_2 , (beta _C – alpha _{CC} , C) -alpha _{A2} , C_2 A -alpha _{B2} , C_2 B + sigma _1 , C_1 – sigma _2 , C_2 end{aligned}$$
(3d)
The parameters in these equations and their meanings are listed in Table 1. Here we assume that the prey types (C_1) and (C_2) have the same birth rate and intraspecific competition strength, but different interaction strengths with A and B. Note that (C_1) and (C_2) are connected by the (sigma _i) terms, which represent the exchange of individuals between these subpopulations through mechanisms studied below; these terms indicate a major difference between our model of a prey with intraspecific heterogeneity and other models of two prey species. For the convenience of analysis, we transform the variables (C_1) and (C_2) to another pair of variables C and (lambda), where (C equiv C_1 + C_2) is the total population of C as before, and (lambda equiv C_2 / (C_1 + C_2)) represents the composition of the population (Fig. 1c). After this transformation and rescaling of variables (described in “Methods”), the new dynamical system can be written as: $$begin{aligned} dot{A}&= A , big ( C , (a_1 (1-lambda ) + a_2 lambda ) – a_0 big ) end{aligned}$$
(4a)
$$begin{aligned} dot{B}&= B , big ( C , (b_1 (1-lambda ) + b_2 lambda ) – b_0 big ) end{aligned}$$
(4b)
$$begin{aligned} dot{C}&= C , big ( 1 – C – A (a_1 (1-lambda ) + a_2 lambda ) – B (b_1 (1-lambda ) + b_2 lambda ) big ) end{aligned}$$
(4c)
$$begin{aligned} dot{lambda }&= lambda (1-lambda ) , big ( A (a_1 – a_2) + B (b_1 – b_2) big ) + eta _1 (1-lambda ) – eta _2 lambda end{aligned}$$
(4d)
Here, (a_i) and (b_i) are the (rescaled) feeding rates of the predators on the prey type (C_i); (a_0) and (b_0) are the death rates of the predators as before; (eta _1) and (eta _2) are the exchange rates of the prey types (Table 1). The latter can be functions of other variables, representing different kinds of heterogeneous interactions that we study below. Notice that Eqs. (4a–4c) are equivalent to the homogeneous Eqs. (2a–2c) but with effective interaction strengths (a_text {eff} = (1-lambda ) , a_1 + lambda , a_2) and (b_text {eff} = (1-lambda ) , b_1 + lambda , b_2) that both depend on the prey composition (lambda) (Fig. 1c).Table 1 Model parameters (before/after rescaling) and their meanings.Full size tableThe variable (lambda) can be considered an internal degree of freedom within the C population. In all of the models we study below, (lambda) dynamically stabilizes to a special value (lambda ^*) (a bifurcation point), as shown in Fig. 3a. Accordingly, a new equilibrium point (P_N) appears (on the line (mathscr {L}) in Fig. 2), at which all three species coexist. For comparison, Fig. 3b shows the equilibrium points if (lambda) is held fixed at any other values, which all result in the exclusion of one of the predators. Thus, heterogeneous interactions give rise to a new coexistence phase (see Fig. 4 below) by bringing the prey composition (lambda) to the value (lambda ^*), instead of having to fine-tune the interaction strengths. The exact conditions for this new equilibrium to be stable are detailed in “Methods”.Figure 3(a) Time series of (lambda) for systems with each kind of heterogeneity. All three systems stabilize at the same (lambda ^*) value, which is the bifurcation point in panel (b). (b) Equilibrium population of each species (X = A), B, or C, with (lambda) held fixed at different values. Solid curves represent stable equilibria and dashed curves represent unstable equilibria (see Eq. (9) in “Methods”). The vertical dashed line is where (lambda = lambda ^*), which is also the bifurcation point. Notice that the equilibrium population of C is maximized at this point (for (a_1 > a_2) and (b_2 > b_1)). Parameters used here are ((a_0, a_1, a_2, b_0, b_1, b_2, rho , eta _1, eta _2, kappa ) = (0.25, 0.5, 0.2, 0.4, 0.2, 0.6, 0.5, 0.05, 0.05, 50)).Full size imageInherent heterogeneityWe first consider a scenario where individuals of the prey species are born as one of two types with a fixed ratio, such that a fraction (rho) of the newborns are (C_2) and ((1-rho )) are (C_1). This could describe dimorphic traits, such as the winged and wingless morphs in aphids12 or the horned and hornless morphs in beetles13. We call this “inherent” heterogeneity, because individuals are born with a certain type and cannot change in later stages of life. The prey type given at birth determines the individual’s interaction strength with the predators. This kind of heterogeneity can be described by Eq. (4d) with (eta _1 = rho (1-C)) and (eta _2 = (1-rho ) (1-C)) (see “Methods”).Figure 4Phase diagrams showing regions of parameter space identified by the stable equilibrium points. Yellow region represents (P_C) (predators A, B both extinct), red represents (P_A) (A excludes B), blue represents (P_B) (B excludes A), and green represents (P_N) (A, B coexist). The middle point (black dot) is where the preferences of the two predators are identical, (a_2/a_0=b_2/b_0) and (b_1/b_0=a_1/a_0). The coexistence phase appears in all three kinds of heterogeneity modeled here. (a–d) Inherent heterogeneity: Individuals of the prey population are born in two types with a fixed composition (rho). In the extreme cases of (rho = 0) and 1, the prey is homogeneous and there is no coexistence of the predators. (e–h) Reversible heterogeneity: Individual prey can switch types with fixed switching rates (eta _1) and (eta _2). As the switching rates increase, the coexistence region shrinks because the prey population becomes effectively homogeneous (the occasional green spots are numerical artifacts because the time to reach the equilibrium becomes long in this limit). (i–l) Adaptive heterogeneity: The switching rates (eta _i) dynamically adapt to the predator densities, so as to maximize the growth rate of the prey. As the sharpness (kappa) of the sigmoidal decision function is increased, the prey adapts more optimally and the region of coexistence expands. Parameters used here are ((a_0, a_1, b_0, b_2) = (0.3, 0.5, 0.4, 0.6)).Full size imageThe stable equilibrium of the system can be represented by phase diagrams that show the identities of the species at equilibrium. We plot these phase diagrams by varying the parameters (a_2) and (b_1) while keeping (a_1) and (b_2) constant. As shown in Fig. 4a–d, the equilibrium state depends on the parameter (rho). In the limit (rho = 0) or 1, we recover the homogeneous case because only one type of C is produced. The corresponding phase diagrams (Fig. 4a, d) contain only two phases where either of the predators is excluded, illustrating the competitive exclusion principle. For intermediate values of (rho), however, there is a new phase of coexistence that separates the two exclusion phases (Fig. 4b, c). There are two such regions of coexistence, which touch at a middle point and open toward the bottom left and upper right, respectively. The middle point is at ((a_2/a_0 = b_2/b_0, b_1/b_0 = a_1/a_0)), where the feeding preferences of the two predators are identical (hence their niches fully overlap). Towards the origin and the far upper right, the predators consume one type of C each (hence their niches separate). The coexistence region in the bottom left is where the feeding rates of the predators are the lowest overall. There can be a region (yellow) where both predators go extinct, if their primary prey type alone is not enough to sustain each predator. Increasing the productivity of the system by increasing the birth rate ((beta _C)) of the prey eliminates this extinction region, whereas lowering productivity causes the extinction region to take over the lower coexistence region. Because the existence and identity of the phases is determined by the configuration of the equilibrium points (Fig. 2, see also section “Mathematical methods”), the qualitative shape of the phase diagram is not sensitive to changes of parameter values.The new equilibrium is not only where the predators A and B can coexist, but also where the prey species C grows to a larger density than what is possible for a homogeneous population. This is illustrated in Fig. 3b, which shows the equilibrium population of C if we hold (lambda) fixed at different values. The point (lambda = lambda ^*) is where the system with a dynamic (lambda) is stable, and also where the population of C is maximized (when A and B prefer different prey types). That means the population automatically stabilizes at the optimal composition of prey types. Moreover, the value of (C^*) at this coexistence point can even be larger than the equilibrium population of C when there is only one predator A or B. This is discussed further in section “Multiple-predator effects and emergent promotion of prey”. These results suggest that heterogeneity in interaction strengths can potentially be a strategy for the prey population to leverage the effects of multiple predators against each other to improve survival.Reversible heterogeneityWe next consider a scenario where individual prey can switch their types. This kind of heterogeneity can model reversible changes of phenotypes, i.e., trait changes that affect the prey’s interaction with predators but are not permanent. For example, changes in coat color or camouflage14,16,17, physiological changes such as defense15, and biomass allocation among tissues18,19. One could also think of the prey types as subpopulations within different spatial patches, if each predator hunts at a preferred patch and the prey migrate between the patches20,21. With some generalization, one could even consider heterogeneity in resources, such as nutrients located in different places, that can be reached by primary consumers, such as swimming phytoplankton22. We can model this “reversible” kind of heterogeneity by introducing switching rates from one prey type to the other. In Eq. (4d), (eta _1) and (eta _2) now represent the switching rates per capita from (C_1) to (C_2) and from (C_2) to (C_1), respectively. Here we study the simplest case where both rates are fixed.In the absence of the predators, the natural composition of the prey species given by the switching rates would be (rho equiv eta _1 / (eta _1 + eta _2)), and the rate at which (lambda) relaxes to this natural composition is (gamma equiv eta _1 + eta _2). Compared to the previous scenario where we had only one parameter (rho), here we have an additional parameter (gamma) that modifies the behavior of the system. Fig. 4e–h shows phase diagrams for the system as (rho) is fixed and (gamma) varies. We again find the new equilibrium (P_N) where all three species coexist. When (gamma) is small, the system has a large region of coexistence. As (gamma) is increased, this region is squeezed into a border between the two regions of exclusion, where the slope of the border is given by (eta _1/eta _2) as determined by the parameter (rho). However, this is different from the exclusion we see in the case of inherent heterogeneity, which happens only for (rho rightarrow 0) or 1, where the borders are horizontal or vertical (Fig. 4a,d). Here the predators exclude each other despite having a mixture of prey types in the population.This special limit can be understood as follows. For a large (gamma), (lambda) is effectively set to a constant value equal to (rho), because it has a very fast relaxation rate. In other words, the prey types exchange so often that the population always maintains a constant composition. In this limit, the system behaves as if it were a homogeneous system with effective interaction strengths (a_text {eff} = (1-rho ) , a_1 + rho , a_2) and (b_text {eff} = (1-rho ) , b_1 + rho , b_2). As in a homogeneous system, there is competitive exclusion between the predators instead of coexistence. This demonstrates that having a constant level of heterogeneity is not sufficient to cause coexistence. The overall composition of the population must be able to change dynamically as a result of population growth and consumption by predators.An interesting behavior is seen when we examine a point inside the shrinking coexistence region as (gamma) is increased. Typical trajectories of the system for such parameter values are shown in Fig. 5. As (gamma) increases, the system relaxes to the line (mathscr {L}) quickly, then slowly crawls along it towards the final equilibrium point (P_N). This is because increasing (gamma) increases the speed that (lambda) relaxes to (lambda ^*), and when (lambda rightarrow lambda ^*), (mathscr {L}) becomes marginally stable. Therefore, the attraction to (mathscr {L}) in the perpendicular direction is strong, but the attraction towards the equilibrium point along (mathscr {L}) is weak. This leads to a long transient behavior that makes the system appear to reach no equilibrium in a limited time23,24. It is especially true when there is noise in the dynamics, which causes the system to diffuse along (mathscr {L}) with only a weak drift towards the final equilibrium (Fig. 5). Thus, the introduction of a fast timescale (quick relaxation of (lambda) due to a large (gamma)) actually results in a long transient.Figure 5Trajectories of the system projected in the A-B plane, with parameters inside the coexistence region (by holding the position of (P_N) fixed). As (gamma) increases, the system tends to approach the line (mathscr {L}) quickly and then crawl along it. The grey trajectory is with independent Gaussian white noise ((sim mathscr {N}(0,0.5))) added to each variable’s dynamics. Noise causes the system to diffuse along (mathscr {L}) for a long transient period before coming to the equilibrium point (P_N). Parameters used here are ((a_0, a_1, a_2, b_0, b_1, b_2) = (0.2, 0.8, 0.5, 0.2, 0.6, 0.9)), chosen to place (P_N) away from the middle of (mathscr {L}) to show the trajectory drifting toward the equilibrium.Full size imageAdaptive heterogeneityA third kind of heterogeneity we consider is the change of interactions in time. By this we mean an individual can actively change its interaction strength with others in response to certain conditions. This kind of response is often invoked in models of adaptive foraging behavior, where individuals choose appropriate actions to maximize some form of fitness25,26. For example, we may consider two behaviors, resting and foraging, as our prey types. Different predators may prefer to strike when the prey is doing different things. In response, the prey may choose to do one thing or the other depending on the current abundances of different predators. Such behavioral modulation is seen, for example, in systems of predatory spiders and grasshoppers27. Phenotypic plasticity is also seen in plant tissues in response to consumers28,29,30.This kind of “adaptive” heterogeneity can be modeled by having switching rates (eta _1) and (eta _2) that are time-dependent. Let us assume that the prey species tries to maximize its population growth rate by switching to the more favorable type. From Eq. (4c), we see that the growth rate of C depends linearly on the composition (lambda) with a coefficient (u(A,B) equiv (a_1 – a_2) A + (b_1 – b_2) B). Therefore, when this coefficient is positive, it is favorable for C to increase (lambda) by switching to type (C_2). This can be achieved by having a positive switching rate (eta _2) whenever (u(A,B) > 0). Similarly, whenever (u(A,B) < 0), it is favorable for C to switch to type (C_1) by having a positive (eta _1). In this way, the heterogeneity of the prey population constantly adapts to the predator densities. We model such adaptive switching by making (eta _1) and (eta _2) functions of the coefficient u(A, B), e.g., (eta _1(u) = 1/(1+mathrm {e}^{kappa u})) and (eta _2(u) = 1/(1+mathrm {e}^{-kappa u})). The sigmoidal form of the functions means that the switching rate in the favorable direction for C is turned on quickly, while the other direction is turned off. The parameter (kappa) controls the sharpness of this transition.Phase diagrams for the system with different values of (kappa) are shown in Fig. 4i–l. A larger (kappa) means the prey adapts its composition faster and more optimally, which causes the coexistence region to expand. In the extreme limit, the system changes its dynamics instantaneously whenever it crosses the boundary where (u(A,B) = 0), like in a hybrid system31. Such a system can still reach a stable equilibrium that lies on the boundary, if the flow on each side of the boundary points towards the other side32. This is what happens in our system and, interestingly, the equilibrium is the same three-species coexistence point (P_N) as in the previous scenarios. The region of coexistence turns out to be largest in this limit (Fig. 4l).Our results suggest that the coexistence of the predators can be viewed as a by-product of the prey’s strategy to maximize its own benefit. The time-dependent case studied here represents a strategy that involves the prey evaluating the risk posed by different predators. This is in contrast to the scenarios studied above, where the prey population passively creates phenotypic heterogeneity regardless of the presence of the predators. These two types of behavior are analogous to the two strategies studied for adaptation in varying environments, i.e., sensing and bet-hedging33,34. The former requires accessing information about the current environment to make optimal decisions, whereas the latter relies on maintaining a diverse population to reduce detrimental effects caused by environmental changes. Here the varying abundances of the predators play a similar role as the varying environment. From this point of view, the heterogeneous interactions studied here can be a strategy of the prey species that is evolutionarily favorable. More

Ren, W. et al. Changes of periphyton abundance and biomass driven by factors specific to flooding inflow in a river inlet area in Erhai Lake, China. Front. Environ. Sci. 9, 680718. https://doi.org/10.3389/fenvs.2021.680718 (2021).Article 

Google Scholar 
Woodruff, S. L. et al. The effects of a developing biofilm on chemical changes across the sediment-water interface in a freshwater environment. Int. Rev. Hydrobiol. 84(5), 509–532 (1999).CAS 

Google Scholar 
Muñoz, I., Real, M., Guasch, H., Navarro, E. & Sabater, S. Effects of atrazine on periphyton under grazing pressure. Aquat. Toxicol. 55(3–4), 239–249 (2001).
Google Scholar 
Hoagland, K. D., Roemer, S. C. & Rosowski, J. R. Colonization and community structure of two periphyton assemblages, with emphasis on the diatoms (Bacillariophyceae). Am. J. Bot. 69, 188–213. https://doi.org/10.2307/2443006 (1982).Article 

Google Scholar 
Steinman, A. D. & McIntire, C. D. Effects of current velocity and light energy on the structure of periphyton assemblages in laboratory streams. J. Phycol. 22, 352–361. https://doi.org/10.1111/J.1529-8817.1986.TB00035.X (1986).Article 

Google Scholar 
Tonkin, J. D., Death, R. G. & Barquín, J. Periphyton control on stream invertebrate diversity: Is periphyton architecture more important than biomass?. Mar. Freshw. Res. 65(9), 818–829 (2014).
Google Scholar 
Beck, W. S., Markman, D. W., Oleksy, I. A., Lafferty, M. H. & Poff, N. L. Seasonal shifts in the importance of bottom-up and top-down factors on stream periphyton community structure. Oikos 128, 680–691. https://doi.org/10.1111/oik.05844 (2018).Article 
CAS 

Google Scholar 
Hogsden, K. L. & Harding, J. S. Consequences of acid mine drainage for the structure and function of benthic stream communities: A review. Freshw. Sci. 31, 108–120. https://doi.org/10.1899/11-091.1 (2012).Article 

Google Scholar 
Sofi, M. S., Bhat, S. U., Rashid, I. & Kuniyal, J. C. The natural flow regime: A master variable for maintaining river ecosystem health. Ecohydrology 13(8), e2247. https://doi.org/10.1002/eco.2247 (2020).Article 

Google Scholar 
Biggs, B. J. F. Eutrophication of streams and rivers: Dissolved nutrient-chlorophyllrelationship for benthic algae. J. N. Am. Benthol. Soc. 19, 17–31 (2000).
Google Scholar 
Ormerod, S. J., Dobson, M., Hildrew, A. G. & Townsend, C. Multiple stressors in freshwater ecosystems. Freshw. Biol. 55, 1–4 (2010).
Google Scholar 
Poff, et al. The natural flow regime: A paradigm for river conservation and restoration. Bioscience 47, 769–784 (1997).
Google Scholar 
Naiman, R. J., Décamps, H., & McClain, M. E. Riparia: Ecology, Conservation and Management of Streamside Communities, (Elsevier/Academic Press, 2005).Gleick, P. H. Water use. Annu. Rev. Environ. Resour. 28, 275–314 (2003).
Google Scholar 
Jenkins, K. M. & Boulton, A. J. Connectivity in a dryland river: Short-term aquatic macroinvertebrate recruitment following floodplain inundation. Ecology 84(10), 2708–2723 (2003).
Google Scholar 
Biggs, B. J. F. Patterns in benthic algae of streams. In Algal Ecology in Freshwater Benthic Ecosystems (eds. Stevenson, R. J., Bothwell, M. L., & Lowe, R. L.) 31–56 (Academic Press, 1996).Smolar-Žvanut, N. & Mikoš, M. The impact of flow regulation by hydropower dams on the periphyton community in the Soča River, Slovenia. Hydrol. Sci. J. 59(5), 1032–1045. https://doi.org/10.1080/02626667.2013.834339 (2014).Article 
CAS 

Google Scholar 
Curry, C. J. & Baird, D. J. Habitat type and dispersal ability influence spatial structuring of larval Odonata and Trichoptera assemblages. Freshw. Biol. 60, 2142–2152 (2015).
Google Scholar 
Wu, N., Cai, Q. & Fohrer, N. Contribution of microspatial factors to benthic diatom communities. Hydrobiologia 732, 49–60. https://doi.org/10.1007/s10750-014-1843-3 (2014).Article 
CAS 

Google Scholar 
Mueller, M., Pander, J. & Geist, J. The effects of weirs on structural stream habitat and biological communities. J. Appl. Ecol 48(6), 1450–1461. https://doi.org/10.1111/j.1365-2664.2011.02035.x (2011).Article 

Google Scholar 
Davies, P. M. et al. Flow–ecology relationships: closing the loop on effective environmental flows. Mar. Freshw. Res. 65(2), 133–141 (2013).
Google Scholar 
Jun, Y. C. et al. Spatial distribution of benthic macroinvertebrate assemblages in relation to environmental variables in Korean nationwide streams. Water 8(1), 27. https://doi.org/10.3390/w8010027 (2016).Article 

Google Scholar 
Biggs, B. J. F. & Close, M. E. Periphyton biomass dynamics in gravel bed rivers: The relative effects of flows and nutrients. Freshw. Biol. 22, 209–231 (1989).CAS 

Google Scholar 
Jowett, I. & Biggs, B. J. F. Flood and velocity effects on periphyton and silt accumulation in two New Zealand rivers. N. Zeal. J. Mar. Freshw. Res. 31, 287–300 (1997).
Google Scholar 
Biggs, B. J. F., Goring, D. G. & Nikora, V. I. Subsidy and stress responses of stream periphyton to gradients in water velocity as a function of community growth form. J. Phycol. 34, 598–607 (1998).
Google Scholar 
Malmqvist, B. & Englund, G. Effects of hydropower-induced flow perturbations on mayfly (Ephemeroptera) richness and abundance in north Swedish river rapids. Hydrobiologia 341(2), 145–158 (1996).
Google Scholar 
Poff, N. L. & Ward, J. V. Herbivory under different flow regimes: A field experiment and test of a model with a benthic stream insect. Oikos 72, 179–188 (1995).
Google Scholar 
Poff, L. N., Wellnitz, T. & Monroe, J. B. Redundancy among three herbivorous insects across an experimental current velocity gradient. Oecologia 134, 262–269. https://doi.org/10.1007/s00442-002-1086-2 (2003).Article 

Google Scholar 
Vaughn, C. C. The role of periphyton abundance and quality in the microdistribution of a stream grazer, Helicopsyche borealis (Trichoptera: Helicopsychidae). Freshw. Biol. 16, 485–493 (1986).
Google Scholar 
Francoeur, S. N. Meta-analysis of lotic nutrient amendment experiments: Detecting and quantifying subtle responses. J. N. Am. Benthol. Soc. 20, 358–368 (2001).
Google Scholar 
Elser, J. J. et al. Global analysis of nitrogen and phosphorus limitation of primary producers in freshwater, marine and terrestrial ecosystems. Ecol. Lett. 10, 1135–1142 (2007).
Google Scholar 
Hillebrand, H. Meta-analysis of grazer control of periphyton biomass across aquatic ecosystems. J. Phycol. 45, 798–806 (2009).
Google Scholar 
Lamberti, G. A. The role of periphyton in benthic food webs. In Algal Ecology—Freshwater Benthic Ecosystems, 533–572 (eds. Stevenson, R. J., Bothwell, M. L. & Lowe, R. L.) (Academic Press, 1996).Lamberti, G. A. et al. Influence of grazer type and abundance on plant–herbivore interactions in streams. Hydrobiologia 306, 179–188 (1995).
Google Scholar 
Gregory, S. V. Plant–herbivore interactions in stream systems. In Stream Ecology (eds. Barnes, J. R. & Minshall, G. W.) 157–189 (Plenum, 1983).Lamberti, G. A. & Moore, J. W. Aquatic insects as primary consumers. In The Ecology of Aquatic Insects (eds Resh, V. H. & Rosenberg, D. M.) 164–195 (Praeger, 1984).
Google Scholar 
Sterner, R. W., Elser, J. J. & Hessen, D. O. Stoichiometric relationships among producers, consumers and nutrient cycling in pelagic ecosystems. Biogeochemistry 17, 49–67 (1992).CAS 

Google Scholar 
Kahlert, M. & Baunsgaard, M. T. Nutrient recycling—A strategy of a grazer community to overcome nutrient limitation. J. N. Am. Benthol. Soc. 18, 363–369 (1999).
Google Scholar 
Burkholder, J. M., Wetzel, R. G. & Klomparens, K. L. Direct comparison of phosphate uptake by adnate and loosely attached microalgae within and intact biofilm matrix. Appl. Environ. Microbiol. 56, 2882–2890 (1990).CAS 

Google Scholar 
Steinman, A. D. Effects of grazers on freshwater benthic algae. In Algal Ecology: Freshwater Benthic Ecosystems (eds. Stevenson, R. J., Bothwell & Lowe, R. L.) 341–366 (Academic Press, 1996).Smucker, N. J. & Vis, M. L. Spatial factors contribute to benthic diatom structure in streams across spatial scales: Considerations for biomonitoring. Ecol. Indic. 11, 1191–1203 (2011).
Google Scholar 
Myers, et al. Biodiversity hotspots for conservation priorities. Nature 403, 853–858 (2000).Article 
CAS 

Google Scholar 
Wang, J., Pan, F., Soininen, J., Heino, J. & Shen, J. Nutrient enrichment modifies temperature-biodiversity relationship in large scale field experiments. Nat. Commun. 7, 13 (2016).
Google Scholar 
Wu, et al. Flow regimes filter species traits of benthic diatom communities and modify the functional features of lowland streams-a nationwide scale study. Sci. Total Environ. 651, 357–366 (2019).CAS 

Google Scholar 
Nisar, M. A. Geospatial approach to study environmental characterization of a Kashmir wetland (Anchar) catchment with special reference to land use/land cover and changing climate. Ph.D Thesis, Sher-e-Kashmir University of Agricultural Sciences and Technology, Kashmir. Weblink. http://krishikosh.egranth.ac.in/handle/1/91309 (2012).Bhat, S. U., Sofi, A. H., Yaseen, T., Pandit, A. K. & Yousuf, A. R. Macro invertebrate community from Sonamarg streams of Kashmir Himalaya. Pak. J. Biol. Sci. 14(3), 182–194. https://doi.org/10.3923/pjbs.2011.182.194 (2011).Article 
CAS 

Google Scholar 
Baba, A. I., Sofi, A. H., Bhat, S. U., & Pandit, A. K. Periphytic algae of river Sindh in the Sonamarg area of Kashmir valley. J. Phytol. 3(6) (2011).Sofi, M. S., Rautela, K. S., Bhat, S. U., Rashid, I. & Kuniyal, J. C. Application of geomorphometric approach for the estimation of hydro-sedimentological flows and cation weathering rate: Towards understanding the sustainable land use policy for the Sindh Basin, Kashmir Himalaya. Water Air Soil Pollut. 232(7), 1–11. https://doi.org/10.1007/s11270-021-05217-w (2021).Article 
CAS 

Google Scholar 
Romshoo, S. A., & Fayaz, M. Use of high resolution remote sensing for improving environmental Friendly tourism master planning in the Alpine Himalaya: A case study of Sonamarg tourist resort, Kashmir. J. Himalayan Ecol. Sustain. Dev. 14 (2019).Biggs, B. J. F. & Kilroy, C. Stream periphyton monitoring manual. Published by NIWA for Ministry for the Environment, 226 Christchurch, New Zealand: NIWA (2000).APHA. Standard Methods for Examination of Water and Wastewater, 22nd edn. (American Public Health Association, 2012).Cox, E. J. Identification of Freshwater Diatoms from Live Material. (Chapman and Hall, 1996). https://doi.org/10.1017/S0025315400041023.Krammer, K., & Lange-Bertalot, H. Bacillariophyceae, Part 5. English and French Translation of the Keys. (VEB Gustav Fisher Verlag, 2000).Reichardt, E. A remarkable association of diatoms in a spring habitat in the Grazer Bergland, Austria. In Iconographia Diatomologica (ed. Lange-Bertalot, H.) 419–479 (2004).Żelazna-Wieczorek, J. Diatom flora in springs of Lódz Hills (Central Poland). Biodiversity, taxonomy and temporal changes of epipsammic diatom assemblages in springs affected by human impact, 419. Volume 13 of Diatom monographs. Gantner. https://books.google.co.in/books?id=bdxeewAACAAJ (2011).Stark, J. D., Boothroyd, I. K. G., Harding, J. S., Maxted, J. R. & Scarsbrook, M. R. Protocols for sampling macroinvertebrates in wadeable streams. In New Zealand Macroinvertebrate Working Group Report no. 1. Prepared for the Ministry for the Environment. Sustainable Management Fund Project, 5103 (2001).Winterbourn, M. J. Sampling stream invertebrates. In Biological Monitoring of Freshwaters. Proceedings of the Seminar. Water and Soil Miscellaneous Publication No. 83 (eds. Pridmore, R. D., Cooper, A. B.) 241–258. (National Water and Soil Conservation Authority, 1985).Barbour, M. T., Gerritsen, J., Snyder, B. D., Stribling, J. B. Rapid bioassessment protocols for use in streams and wadeable rivers: periphyton, benthic macroinvertebrates and fish, 339. (United States Environmental Protection Agency, Office of Water, 1999).Malmqvist, B. & Hoffsten, P. O. Macroinvertebrate taxonomic richness, community structure and nestedness in Swedish streams. Fundam. Appl. Limnol. 150(1), 29–54. https://doi.org/10.1127/archiv-hydrobiol/150/2000/29 (2000).Article 

Google Scholar 
Ilmonen, J. & Paasivirta, L. Benthic macrocrustacean and insect assemblages in relation to spring habitat characteristics: Patterns in abundance and diversity. Hydrobiologia 533(1–3), 99–113. https://doi.org/10.1007/s10750-004-2399-4 (2005).Article 

Google Scholar 
Munasinghe, D. S. N., Najim, M. M. M., Quadroni, S. & Musthafa, M. M. Impacts of streamflow alteration on benthic macroinvertebrates by mini-hydro diversion in Sri Lanka. Sci. Rep. 11(1), 546. https://doi.org/10.1038/s41598-020-79576-5 (2021).Article 
CAS 

Google Scholar 
Edmondson, W. T. Fresh-Water Biology, 2nd ed. 1050–1056 (Wiley, 1959).Pennak, R. W. Freshwater Invertebrates of United States. (Wiley, 1978).McCafferty, W. P., Provonsha, A. V. Aquatic entomology: The fishermen’s and ecologists’ Illustrated Guide to Insects and their Relatives. (Jones and Bartlett Publishers, 1983).Borror, D., Triplehorn, C., Johnson, N. An Introduction to the Study of Insects, 6th ed. (Saunders College Publishing, 1989).Ward, J. V. Aquatic Insect Ecology, Biology and Habitat. (Wiley, 1992).Engblom, E. & Lingdell, P.E. Analyses of Benthic Invertebrates (ed. Nyman, L.) (1999).Bouchard, R. W. Guide to Aquatic Invertebrates of the Upper Midwest: Identification Manual for Students (University of Minnesota, 2004).
Google Scholar 
Subramanian, K. A. & Sivaramakrishnan, K. G. Aquatic Insects for Biomonitoring Freshwater Ecosystems—A Methodology Manual. (Ashoka Trust for Ecology and Environment (ATREE), 2007).Thorp, J. H., & Covich, A. P. (eds.) Ecology and Classification of North American Freshwater Invertebrates. (Academic Press, 2009).Allan, J. D. & Castillo, M.M. An introduction to fluvial ecosystems. In Stream Ecology: Structure and Function of Running Waters, 1–12 (2007).Oksanen, et al. Vegan: Community ecology package. In: R package version 2.4-3.McCune, B. & Grace, B. Analysis of Ecological Communities (MjM Software Design, 2016).Clarke, K. R. & Gorley, R. N. Primer v6 Permanova+ (Primer-E Ltd., 2006).
Google Scholar 
Salazar, G. EcolUtils: Utilities for Community Ecology Analysis. R package version 0.1 software (2018).Anderson, M. J., Ellingsen, K. E. & McArdle, B. H. Multivariate dispersion as a measure of beta diversity. Ecol. Lett. 9(6), 683–693 (2006).
Google Scholar 
Gardener, M. Community Ecology: Analytical Methods in Using R and Excel. (Pelagic Publishing, 2014).Chao, A. & Bunge, J. Estimating the number of species in a stochastic abundance model. Biometrics 58, 531–539. https://doi.org/10.1111/j.0006-341X.2002.00531.x (2002).Article 
MATH 

Google Scholar 
Peres-Neto, P. R., Legendre, P., Dray, S. & Borcard, D. Variation partitioning of species data matrices: estimation and comparison of fractions. Ecology 87, 2614–2625 (2006).
Google Scholar 
Meng, X. L. et al. Responses of macroinvertebrates and local environment to short-term commercial sand dredging practices in a flood-plain lake. Sci. Total Environ. 631, 1350–1359 (2018).
Google Scholar 
Core Team, R. R: A Language and Environmental for Statistical Computing. (R Foundation for Statistical Computing, 2017).Wood, P. J. & Armitage, P. D. Biological effects of fine sediment in the lotic environment. Environ. Manag. 21, 203–217 (1997).CAS 

Google Scholar 
Marchant, R. Changes in the benthic invertebrate communities of the Thomson River, southeastern Australia, after dam construction. Regul. Rivers Res. Manag. 4, 71–89 (1989).
Google Scholar 
Gray, L. J. & Ward, J. V. Effects of sediment releases from a reservoir on stream macroinvertebrates. Hydrobiologia 96, 177–184 (1982).
Google Scholar 
Sand-Jensen, K., Moller, J. & Olesen, B. H. Biomass regulation of microbenthic algae in Danish lowland streams. Oikos 53, 332–340 (1988).
Google Scholar 
Lewis, M. A., Weber, D. E., Stanley, R. S. & Moore, J. C. Dredging impact on an urbanized Florida bayou: Effects on benthos and algal-periphyton. Environ. Pollut. 115(2), 161–171 (2001).CAS 

Google Scholar 
Biggs, B. J. Algal ecology in freshwater benthic ecosystems geology and landuse to the habitat template of periphyton in stream ecosystems. Freshw. Biol. 33, 419–438 (1995).
Google Scholar 
Taylor, et al. Can diatom-based pollution indices be used for biomonitoring in South Africa? A case study of the Crocodile West and Marico water management area. Hydrobiologia 592, 455–464 (2007).
Google Scholar 
Porter, et al. Efficacy of algal metrics for assessing nutrient and organic enrichment in flowing waters. Freshw. Biol. 53, 1036–1054 (2008).
Google Scholar 
Wetzel, R. G. & Likens, G. E. Limnological analyses, 3rd ed. In Nitrogen, Phosphorus, and Other Nutrients, 85–113. (Springer, 2000). https://doi.org/10.1007/978-1-4757-3250-4.Wetzel, R. G. Attached algal-substrata interactions: Fact or myth, and when and how? vol. 17. In Periphyton of Freshwater Ecosystems (ed. Wetzel, R.) 207–215 (Springer, 1983). https://doi.org/10.1007/978-94-009-7293-3_28.Krajenbrink, H. J. et al. Diatoms as indicators of the effects of river impoundment at multiple spatial scales. PeerJ 7, e8092. https://doi.org/10.7717/peerj.8092 (2019).Article 

Google Scholar 
Poff, N. L., Voelz, N. J., Ward, J. V. & Lee, R. E. Algal colonization under four experimentally-controlled current regimes in a high mountain stream. J. N. Am. Benthol. Soc. 9, 303–318 (1990).
Google Scholar 
Dodds, W. K. & Marra, J. L. Behaviors of the midge, Cricotopus (Diptera; Chironomidae) related to mutualism with Nostoc parmeloides (Cyanobacteria). Aquat. Insects 11, 201–208 (1989).
Google Scholar 
Tang, T., Niu, S. Q. & Dudgeon, D. Responses of epibenthic algal assemblages to water abstraction in Hong Kong streams. Hydrobiologia 703(1), 225–237. https://doi.org/10.1007/s10750-012-1362-z (2013).Article 
CAS 

Google Scholar 
Maheshwari, K., Vashistha, J., Paulose, P. V. & Agarwal, T. Seasonal changes in phytoplankton community of lake Ramgarh, India. Int. J. Curr. Microbiol. Appl. Sci. 4(11), 318–330 (2015).CAS 

Google Scholar 
Luttenton, M. R., & Baisden, C. The relationships among disturbance, substratum size and periphyton community structure. In Advances in Algal Biology: A Commemoration of the Work of Rex Lowe 111–117. (Springer, 2006).Uehlinger, U. Spatial and temporal variability of periphyton biomass in a prealpine river (Necker, Switzerland). Arch. Fur. Hydrobiol. 123, 219–237 (1991).
Google Scholar 
Hill, W. R. Effects of light. In Algal Ecology in Freshwater Benthic Ecosystems. 121–148 (eds. Stevenson, R. J., Bothwell, M. L., Lowe, R. L.) (Academic Press, 1996).DeNichola, D. M. Periphyton responses to temperature at different ecological levels. In Algal Ecology in Freshwater Benthic Ecosystems. (eds. Stevenson, R. J., Bothwell, M. L., Lowe, R. L.) 149–181 (Academic Press, 1996).O’Reilly, C. M. Seasonal dynamics of periphyton in a large tropical lake. Hydrobiologia 553, 293–301. https://doi.org/10.1007/s10750-005-0878-x (2006).Article 

Google Scholar 
Borduqui, M. & Ferragut, C. Factors determining periphytic algae succession in a tropical hypereutrophic reservoir. Hydrobiologia 683, 109–122. https://doi.org/10.1007/s10750-011-0943-6 (2012).Article 
CAS 

Google Scholar 
De Souza, M. L., Pellegrini, B. G. & Ferragut, C. Periphytic algal community structure in relation to seasonal variation and macrophyte richness in a shallow tropical reservoir. Hydrobiologia 755, 183–196. https://doi.org/10.1007/s10750-015-2232-2 (2015).Article 

Google Scholar 
Prowse, T. D. River-ice hydrology. In Encyclopedia of Hydrological Sciences, vol. 4 (ed. Anderson, M. G.). (Wiley, 2005).Rusanov, A. G., Stanislavskaya, E. V. & Ács, É. Periphytic algal assemblages along environmental gradients in the rivers of the Lake Ladoga basin, Northwestern Russia: Implication for the water quality assessment. Hydrobiologia 695(1), 305–327 (2012).CAS 

Google Scholar 
Sofi, M. S., Hamid, A., Bhat, S. U., Rashid, I. & Kuniyal, J. C. Impact evaluation of the run-of-river hydropower projects on the water quality dynamics of the Sindh River in the Northwestern Himalayas. Environ. Monit. Assess. 194(9), 1–6 (2022).
Google Scholar 
MCCormick, P. V. Resource competition and species coexistence in freshwater algal assemblages. In Algal ecology—Freshwater Benthic Ecosystems (eds. Stevenson, R. J., Bothwell, M. L., Lowe, R. L.) 229–252 (Academic, 1996).Hillebrand, H., Worm, B. & Lotze, H. K. Marine microbenthic community structure regulated by nitrogen loading and grazing pressure. Mar. Ecol. Prog. Ser. 204, 27–38 (2000).CAS 

Google Scholar  More

Data were analyzed from samples collected, processed, and published previously [21, 25, 29] and have been summarized here. The present analysis, which consisted of sequence data processing, the calculation of taxon-specific isotopic signatures, and subsequent analyses, reflects original work.Sample collection and isotope incubationTo generate experimental data, three replicate soil samples were collected from the top 10 cm of plant-free patches in four ecosystems along the C. Hart Merriam elevation gradient in Northern Arizona. From low to high elevation, these sites are located in the following environments: desert grassland (GL; 1760 m), piñon-pine juniper woodland (PJ; 2020 m), ponderosa pine forest (PP; 2344 m), and mixed conifer forest (MC; 2620 m). Soil samples 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. This produced three distinct but homogenous soil samples from each of the four ecosystems that were subject to experimental treatments. Three treatments were applied to bring soils to 70% water-holding capacity: water alone (control), water with glucose (C treatment; 1000 µg C g−1 dry soil), or water with glucose and a nitrogen source (CN treatment; [NH4]2SO4 at 100 µg N g−1 dry soil). To track growth through isotope assimilation, both 18O-enriched water (97 atom %) and 13C-enriched glucose (99 atom %) were used. In all treatments isotopically heavy samples were paired with matching “light” samples that received water with a natural abundance isotope signatures. For 18O incubations, this design resulted in three soil samples per ecosystem per treatment (across four ecosystems and three treatments, n = 36) while 13C incubations were limited to only C and CN treatments (n = 24). Previous analyses suggest that three replicates is sufficient to detect growth of 10 atom % 18O in microbial DNA with a power of 0.6 and a growth of 5 atom % 18O with a power of 0.3 (12 and 6 atom % respectively for 13C) [30]. All soils were incubated in the dark for one week. Following incubation, soils were frozen at −80 °C for one week prior to DNA extraction.Quantitative stable isotope probingThe procedure of qSIP (quantitative stable isotope probing) is described here but has been applied to these samples as previously published [17, 21, 25]. DNA extraction was performed on soils using a DNeasy PowerSoil HTP 96 Kit (MoBio Laboratories, Carlsbad, CA, USA) and following manufacturer’s protocol. Briefly, 0.25 g of soils from each sample were carefully added to deep, 96-well plates containing zirconium dioxide beads and a cell lysis solution with sodium dodecyl sulfate (SDS) and shaken for 20 min. Following cell lysis, supernatant was collected and centrifuged three times in fresh 96-well plates with reagents separating DNA from non-DNA organic and inorganic materials. Lastly, DNA samples were collected on silica filter plates, rinsed with ethanol and eluted into 100 µL of a 10 mM Tris buffer in clean 96-well plates. To quantify the degree of 18O or 13C isotope incorporation into bacterial DNA (excess atom fraction or EAF), the qSIP protocol [31] was used, though modified slightly as reported previously [21, 24, 32]. Briefly, microbial growth was quantified as the change in DNA buoyant density due to incorporation of the 18O or 13C isotopes through the method of density fractionation by adding 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. Each post-incubation sample was thus converted from a continuous gradient into approximately 20 fractions (150 µL) using a modified fraction recovery system (Beckman Coulter). The density of each fraction was measured with a Reichart AR200 digital refractometer (Reichert Analytical Instruments, Depew, NY, USA). Fractions with densities between 1.640 and 1.735 g cm−3 were retained as densities outside this range generally did not contain DNA. In all retained fractions, DNA was cleaned and purified using isopropanol precipitation and the abundance of bacterial 16S rRNA gene copies was quantified with qPCR using primers specific to bacterial 16S rRNA genes (Eub 515F: AAT GAT ACG GCG ACC ACC GAG TGC CAG CMG CCG CGG TAA, 806R: CAA GCA GAA GAC GGC ATA CGA GGA CTA CVS GGG TAT CTA AT). Triplicate reactions were 8 µL 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. Reactions were performed on a CFX384 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) under the following cycling conditions: 95 °C at 1 min and 44 cycles at 95 °C (30 s), 64.5 °C (30 s), and 72 °C (1 min). Separate from qPCR, retained sample-fractions were subject to a similar amplification step of the 16S rRNA gene V4 region (515F: GTG YCA GCM GCC GCG GTA A, 806R: GGA CTA CNV GGG TWT CTA AT) in preparation for sequencing with the same reaction mix but differing cycle conditions – 95 °C for 2 min followed by 15 cycles at 95 °C (30 s), 55 °C (30 s), and 60 °C (4 min). The resulting 16S rRNA gene V4 amplicons were sequenced on a MiSeq sequencing platform (Illumina, Inc., San Diego, CA, USA). DNA sequence data and sample metadata have been deposited in the NCBI Sequence Read Archive under the project ID PRJNA521534.Sequence processing and qSIP analysisIndependently from previous publications, we processed raw sequence data of forward and reverse reads (FASTQ) within the QIIME2 environment [33] (release 2018.6) and denoised sequences within QIIME2 using the DADA2 pipeline [34]. We clustered the remaining sequences into amplicon sequence variants (ASVs, at 100% sequence identity) against the SILVA 138 database [35] using a pre-trained open-reference Naïve Bayes feature classifier [36]. We removed samples with less than 3000 sequence reads, non-bacterial lineages, and global singletons and doubletons. We converted ASV sequencing abundances in each fraction to the number of 16S rRNA gene copies per gram dry soil based on qPCR abundances and the known amount of dry soil equivalent added to the initial extraction. This allowed us to express absolute population densities, rather than relative abundances. Across all replicates, we identified 114 543 unique bacterial ASVs.We calculated the 18O and 13C excess atom fraction (EAF) for each bacterial ASV using R version 4.0.3 [37] and data.table [38] with custom scripts available at https://www.github.com/bramstone/. Negative enrichment values were corrected using previously published methods [17]. ASVs that appeared in less than two of the three replicates of an ecosystem-treatment combination (n = 3) and less than three density fractions within those two replicates were removed to avoid assigning spurious estimates of isotope enrichment to infrequent taxa. Any ASVs filtered out of one ecosystem-treatment group were allowed to be present in another if they met the frequency threshold. Applying these filtering criteria, we limited our analysis towards 3759 unique bacterial ASVs which accounted for a small proportion of the total diversity but represented 68.0% of all sequence reads, and encompassed most major bacterial groups (Supplementary Fig. 1).Analysis of life history strategies and nutrient responseAll statistical tests were conducted in R version 4.0.3 [37]. We assessed the ability of phylum-level assignment of life history strategy to predict growth in response to C and N addition, as proxied by the incorporation of heavy isotope during DNA replication [39, 40]. Phylum-level assignments (Table 1) were based on the most frequently observed behavior of lineages with a representative phylum (or subphylum) as compiled previously [23]. We averaged 18O EAF values of bacterial taxa for each treatment and ecosystem and then subtracted the values in control soils from values in C-amended soils to determine C response (∆18O EAFC) and from the 18O EAF of bacteria in CN-amended soils to determine C and N response (Δ18O EAFCN). Because an ASV must have a measurable EAF in both the control and treatment for a valid Δ18O EAF to be calculated, we were only able to resolve the nutrient response for 2044 bacterial ASVs – 1906 in response to C addition and 1427 in response to CN addition.We used Gaussian finite mixture modeling, as implemented by the mclust R package [41], to demarcate plausible multi-isotopic signatures for oligotrophs and copiotrophs. For each treatment, we calculated average per-taxon 13C and 18O EAF values. To compare both isotopes directly, we divided 18O EAF values by 0.6 based on the estimate that this value (designated as µ) represents the fraction of oxygen atoms in DNA derived from the 18O-water, rather than from 16O within available C sources [42]. Two mixture components, corresponding to oligotrophic and copiotrophic growth modes, were defined using the Mclust function using ellipsoids of equal volume and shape. We observed several microorganisms with high 18O enrichment but comparatively low 13C enrichment, potentially indicating growth following the depletion of the added glucose, and that were reasonably clustered as oligotrophs in our mixture model.We tested how frequently mixture model clustering of each microorganism’s growth (based on average 18O–13C EAF in a treatment) could predict its growth across replicates (n = 12 in each treatment—although individual). We applied the treatment-level mixture models defined above to the per-taxon isotope values in each replicate, recording when a microorganism’s life history strategy in a replicate agreed with the treatment-level cluster, and when it didn’t. We used exact binomial tests to test whether the number of “successes” (defined as a microorganism being grouped in the same life history category as its treatment-level cluster) was statistically significant. To account for type I error across all individual tests (one per ASV per treatment), we adjusted P values in each treatment using the false-discovery rate (FDR) method [43].To determine the extent that life history categorizations may be appropriately applied at finer levels of taxonomic resolution, we constructed several hierarchical linear models using the lmer function in the nlme package version 3.1-149 [44]. To condense growth information from both isotopes into a single analysis, 18O and 13C EAF values were combined into a single variable using principal components analysis separately for each treatment. Across the C and CN treatments, the first principal component (PC1) was able to explain – respectively – 86% and 91% of joint variation of 18O and 13C EAF values. In all cases, we applied PC1 as the response variable and treated taxonomy and ecosystem as random model terms to limit the potential of pseudo-replication to bias significance values. We used likelihood ratio analysis and Akaike information criterion (AIC) values to compare models where life history strategy was determined based on observed nutrient responses at different taxonomic levels (Eq. 1) against a model with the same random terms but without any life history strategy data (Eq. 2). Separate models were applied to each treatment. To reduce model overfitting, we removed families represented by fewer than three bacterial ASVs as well as phyla represented by only one order. In addition, we removed bacterial ASVs with unknown taxonomic assignments (following Morrissey et al. [21]). This limited our analysis to 1 049 ASVs in the C amendment and 984 in the CN amendment.$${{{{{rm{PC}}}}}}{1}_{{18{{{{{rm{O}}}}}} – 13{{{{{rm{C}}}}}}}}sim {{{{{rm{strategy}}}}}} + 1|{{{{{rm{phylum}}}}}}/{{{{{rm{class}}}}}}/{{{{{rm{order}}}}}}/{{{{{rm{family}}}}}}/{{{{{rm{genus}}}}}}/{{{{{rm{eco}}}}}}$$
(1)
$${{{{{rm{PC}}}}}}{1}_{{18{{{{{rm{O}}}}}} – 13{{{{{rm{C}}}}}}}}sim 1 + 1|{{{{{rm{phylum}}}}}}/{{{{{rm{class}}}}}}/{{{{{rm{order}}}}}}/{{{{{rm{family}}}}}}/{{{{{rm{genus}}}}}}/{{{{{rm{eco}}}}}}$$
(2)
Here, life history strategy was defined at each taxonomic level using the mixture models above and based on the mean 18O and 13C EAF values of each bacterial lineage (Supplemental Fig. 2). We compared these models with the no-strategy model (Eq. 2) directly using likelihood ratio testing. More