Philippa Kaur

1.Stork, N. E., McBroom, J., Gely, C. & Hamilton, A. J. New approaches narrow global species estimates for beetles, insects, and terrestrial arthropods. PNAS 112, 7519–7752 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
2.Sømme, L. & Block, W. Adaptations to alpine and polar environments in insects and other terrestrial arthropods. In Insects at Low Temperature (eds Lee, R. E. & Denlinger, D. L.) (Springer, 1991).
Google Scholar 
3.Springer, M. Marine insects. In: Marine Biodiversity of Costa Rica, Central America. (ed. Wehrtmann, I. S. & Cortés, J.) Monographiae Biologicae, vol 86. (Springer, 2009).4.Cebeci, H. H. et al. The wood boring insects (Coleoptera: Cerambycidae and Buprestidae) recorded as the new pests for acer Undulatum Pojark from the Babadag mountain (SW Turkey). Fresenius Environ. Bull. 27, 9325–9328 (2018).CAS 

Google Scholar 
5.Becker, G. Untersuchungen über die Ernährungsphysiologie der Hausbockkäferlarven. Z. Vgl. Physiol. 29, 315–388 (1942).CAS 
Article 

Google Scholar 
6.Nardi, J. B., Mackie, R. I. & Dawson, J. O. Could microbial symbionts of arthropod guts contribute significantly to nitrogen fixation in terrestrial ecosystems?. J. Insect Physiol. 48, 751–763 (2002).CAS 
PubMed 
Article 

Google Scholar 
7.Wermelinger, B., Duelli, P. & Obrist, M. K. Dynamics of saproxylic beetles (Coleoptera) in windthrow areas in alpine spruce forests. For. Snow Landsc. 77, 133–148 (2002).
Google Scholar 
8.Evans, H. F., Moraal, L. G. & Pajares, J. A. Biology, ecology and economic importance of Buprestidae and Cerambycidae. In: Bark and Wood Boring Insects in Living Trees in Europe, a Synthesis (ed. Lieutier, F., Day, K. R., Battisti, A., Grégoire, J. C. & Evans, H. F.) 447–474 (Springer, 2004).9.FAO (Food and Agriculture Organization of the United Nations) Global review of forest pests and diseases. FAO Forestry Paper 156 (2009).10.Offenberg, J., Nielsen, J. S. & Damgaard, V. Wood Ant (Formica polyctena) services and disservices in a danish apple plantation. Sociobiology 66, 247–256 (2019).Article 

Google Scholar 
11.Cline, A., Ivie, M. A., Bellamy, C. L. & Scher, J. A Resource for Wood Boring Beetles of the World: Wood Boring Beetle Families, Lucid v. 3.4. USDA/APHIS/PPQ Center for Plant Health Science and Technology, Montana State University, and California Department of Food and Agriculture. http://www.lucidcentral.org/keys/v3/WBB (2009).12.Rudinski, J. A. Ecology of Scolytidae. Annu. Rev. Entomol. 7, 327–348 (1962).Article 

Google Scholar 
13.Huber, D. P. W., Aukema, B. H., Hodgkinson, R. S. & Lindgren, B. S. Successful colonization, reproduction, and new generation emergence in live interior hybrid spruce, Picea engelmannii x glauca, by mountain pine beetle, Dendroctonus ponderosae. Agric. For. Entomol. 11, 83–89 (2009).Article 

Google Scholar 
14.O’neill, K. M., Fultz, J. E. & Ivie, M. A. Distribution of adult Cerambycidae and Buprestidae (Coleoptera) in a subalpine forest under shelterwood management. Coleopt. Bull. 62(1), 27–36 (2008).Article 

Google Scholar 
15.Leather, S. R., Day, K. R. & Salisbury, A. N. The biology and ecology of the large pine weevil, Hylobius abietis (Coleoptera: Curculionidae): A problem of dispersal?. Bull. Entomol. Res. 89, 3–16 (1999).Article 

Google Scholar 
16.Nordlander, G., Bylund, H., Örlander, G. & Wallertz, K. Pine weevil population density and damage to coniferous seedlings in a regeneration area with and without shelterwood. Scand. J. For. Res. 18, 438–448 (2003).Article 

Google Scholar 
17.Nordlander, G., Hellqvist, C., Johansson, K. & Nordenhem, H. Regeneration of European boreal forests: Effectiveness of measures against seedling mortality caused by the pine weevil Hylobius abietis. For. Ecol. Manag. 262(12), 2354–2363 (2011).Article 

Google Scholar 
18.Bale, J. S. et al. Herbivory in global climate change research: Direct effects of rising temperature on insect herbivores. Glob. Change Biol. 8, 1–16 (2002).ADS 
Article 

Google Scholar 
19.Kahuthia-Gathu, R., Kirubi, D. T. & Gitonga, D. Composition and abundance ofwood-boring beetles of Acacia xanthophloea and their associated natural enemies in Thika, Kenya. J. Asia Pac. Biodivers. 11, 248–254 (2018).Article 

Google Scholar 
20.Holmes, T. & Koch, F. Bark beetle epidemics, life satisfaction, and economic well-being. Forests 10, 696 (2019).Article 

Google Scholar 
21.Adlung, K. G. A critical evaluation of the European research on use of red wood ants (Formica rufa group) for the protection of forests against harmful insects. Z. Angew. Entomol. 57, 167–189 (1966).Article 

Google Scholar 
22.Goldazarena, A., Romón, P. & López, S. Bark beetles control in forests of Northern Spain. Integrated Pest Management and Pest Control—Current and Future Tactics, Dr. Sonia Soloneski (Ed), ISBN: 978-953-51-0050-8, In Tech. http://www.intechopen.com/books/integrated-pest-management-and-pest-control-current-and-future-tactics/bark-beetles-control-in-forests-of-northern-spain (2012).23.IUCN 2020. The IUCN Red List of Threatened Species. Version 2020-1. https://www.iucnredlist.org (2020).24.Abdullah, H., Darvishzadeh, R., Skidmore, A. K., Groen, T. A. & Heurich, M. European spruce bark beetle (Ips typographus, L.) green attack affects foliar reflectance and biochemical properties. Int. J. Appl. Earth Obs. 64, 199–209 (2018).Article 

Google Scholar 
25.Bentz, B. J. & Jӧnsson, A. M. Modeling bark beetle responses to climate change. In Bark Beetles: Biology and Ecology of Native and Invasive Species (eds Vega, F. & Hofstetter, R.) 533–553 (Elsevier Academic Press, 2015).Chapter 

Google Scholar 
26.Linder, M. et al. Climate change impacts, adaptive capacity, and vulnerability of European forest ecosystems. For Ecol. Manag. 259, 698–709 (2010).Article 

Google Scholar 
27.Marini, L. et al. Climate drivers of bark beetle outbreak dynamics in Norway spruce forests. Ecography 40, 1426–1435 (2017).Article 

Google Scholar 
28.Schowalter, T. D. Insect ecology: an ecosystem approach, 3rd ed. (ed. Schowalter, T. D.) (Academic Press, 2011).29.Frizzi, F., Masoni, A., Quilghini, G., Ciampelli, P. & Santini, G. Chronicle of an impact foretold: The fate and effect of the introduced Formica paralugubris ant. Biol. Invasions 20, 3575–3589 (2018).Article 

Google Scholar 
30.Leong, M. et al. The habitats humans provide: Factors affecting the diversity and composition of arthropods in houses. Sci. Rep. 7, 15347 (2017).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
31.Schweiger, O. et al. Quantifying the impact of environmental factors on arthropod communities in agricultural landscapes across organisational levels and spatial scales. J. Appl. Ecol. 42, 1129–1139 (2005).Article 

Google Scholar 
32.Trigos-Peral, G. et al. Three categories of urban green areas and the effect of their different management on the communities of ants, spiders and harvestmen. Urban Ecosyst. 23, 803–818 (2020).Article 

Google Scholar 
33.Bellard, C., Bertelsmeier, C., Leadley, P., Thuiller, W. & Courchamp, F. Impacts of climate change on the future of biodiversity. Ecol. Lett. 15, 365–377 (2012).PubMed 
PubMed Central 
Article 

Google Scholar 
34.Czechowski, W., Radchenko, A., Czechowska, W. & Vepsäläinen, K. The ants of Poland with reference to the myrmecofauna of Europe. Fauna Poloniae 4 (Natura optima dux Foundation, 2012).35.Seidl, R. et al. Forest disturbances under climate change. Nat. Clim. Chang. 7, 395–402 (2017).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Sommerfeld, A. et al. Patterns and drivers of recent disturbances across the temperate forest biome. Nat. Commun. 9, 4355 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
37.Sorvari, J. & Hakkarainen, H. Wood ants are wood ants: Deforestation causes population declines in the polydomous wood ant Formica aquilonia. Ecol. Entomol. 32, 707–711 (2007).Article 

Google Scholar 
38.Roura-Pascual, N. et al. Geographic potential of Argentine ants (Linepithema humile Mayr) in the face of global climate change. Proc. R. Soc. Lond. B. 271(1557), 2527–2535 (2004).Article 

Google Scholar 
39.Kroll, J. C. & Fleet, R. R. Impact of woodpecker predation on over-wintering within-tree populations of the southern pine beetle (Dendroctonus frontalis). In The Role of Insectivorous Birds in Forest Ecosystems (eds Dickson, J. G. et al.) 269–281 (Academic Press, 1979).Chapter 

Google Scholar 
40.Khanday, A. L., Sureshan, P. M., Buhroo, A. A., Ranjith, A. P. & Tselikh, E. Pteromalid wasps (Hymenoptera: Chalcidoidea) associated with bark beetles, with the description of a new species from Kashmir, India. J. Asia Pac. Biodivers. 12, 262–272 (2019).Article 

Google Scholar 
41.Reeve, J. D. Predation and bark beetle dynamics. Oecologia 112, 48–54 (1997).ADS 
PubMed 
Article 

Google Scholar 
42.Curtsdotter, A. et al. Ecosystem function in predator–prey food webs confronting dynamic models with empirical data. J. Anim. Ecol. 88, 196–210 (2018).PubMed 
Article 

Google Scholar 
43.Choate, B. & Drummond, F. A. Ants as biological control agents in agricultural cropping systems. Terr. Arthropod Rev. 4, 157–180 (2011).Article 

Google Scholar 
44.Maňák, V., Björklund, N., Lenoir, L. & Nordlander, G. The effect of red wood ant abundance on feeding damage by the pine weevil Hylobius abietis. Agric. For. Entomol. 17, 57–63 (2015).Article 

Google Scholar 
45.Robinson, E. J. H., Stockan, J. A. & Glenn, R. I. Wood ants and their interaction with other organisms. In Wood Ant Ecology and Conservation (eds Stockan, A. & Robinson, E. J. H.) 177–206 (Cambridge Universsity Press, 2016).
Google Scholar 
46.Sorvari, J. Foraging distances and potentiality in forest pest insect control: An example with two candidate ants (Hymenoptera: Formicidae). Myrmecol. News 12, 211–215 (2009).
Google Scholar 
47.Gößwald, K. Die Waldameise im Ökosystem Wald, ihr Nutzen und ihre Hege. Aula Verlag, Wiesbaden, 510 pp (1990). (shortened editon is 978-3-89104-755-2 (ISBN)).48.Juhász, O. et al. Large- and small-scale environmental factors drive distributions of ant mound size across a latitudinal gradient. Insects 11, 350 (2020).PubMed Central 
Article 
PubMed 

Google Scholar 
49.Juhász, O. et al. Consequences of climate change-induced habitat conversions on red wood ants in a Central European mountain: a case study. Animals 10, 1677 (2020).PubMed Central 
Article 
PubMed 

Google Scholar 
50.Risch, A., Ellis, S. & Wiswell, H. Where and why? Wood ant population ecology. In Wood Ant Ecology and Conservation (eds Stockan, A. & Robinson, E. J. H.) 81–105 (Cambridge Universsity Press, 2016).
Google Scholar 
51.R Core Team. R: A Language and Environment for Statistical Computing. https://www.R-project.org/ (2020).52.Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Minchin, P. R., O’Hara, R. B., Simpson, G. L., Solymos, P., Stevens, M. H. H., Szoecs, E., Wagner, H. Vegan: Community Ecology Package (2019).53.Roberts, D.W. (2010) labdsv: Ordination and Multivariate Analysis for Ecology. http://cran.r-project.org/package=labdsv54.Bates, D., Maechler, M., Bolker, B. & Walke, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67(1), 1–48 (2015).Article 

Google Scholar 
55.Hartig, F. DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models. R package version 0.1. 0. CRAN/GitHub (2016).56.Lindén, A. & Mäntyniemi, S. Using the negative binomial distribution to model overdispersion in ecological count data. Ecology 97, 1414–1421 (2011).Article 

Google Scholar 
57.Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag, 2016).MATH 
Book 

Google Scholar 
58.Schowalter, T. D. Ecology and management of bark beetles (Coleoptera: Curculionidae: Scolytinae) in Southern Pine Forests. J. Integ. Pest Mngmt. 3, 1–7 (2012).Article 

Google Scholar 
59.Hölldobler, B. & Wilson, E. O. The Ants (Harvard University Press, 1990).Book 

Google Scholar 
60.Siemann, E., Haarstad, J. & Tilman, D. Dynamics of plant and arthropod diversity during old field succession. Ecography 22, 406–414 (1999).Article 

Google Scholar 
61.Schlyte, F. & Anderbrandt, O. Competition and niche separation between two bark beetles: Existence and mechanisms. Oikos 68, 437–447 (1993).Article 

Google Scholar 
62.Raffa, K. F. et al. Cross-scale drivers of natural disturbances prone to anthropogenic amplification: The dynamics of bark beetle eruptions. BioSci 58, 501–517 (2008).Article 

Google Scholar 
63.Drees, B. M., Jackman, J. A. & Merchant, M. E. Wood-Boring Insects of Trees and Shrub. Texas Agricultural Extension Service. http://hdl.handle.net/1969.1/160397 (1994).64.Ray, C. et al. Patterns of woodboring beetle activity following fires and bark beetle outbreaks in montane forests of California, USA. Fire Ecol. 15, 21 (2019).Article 

Google Scholar 
65.Risch, S. J. & Carroll, C. R. The ecological role of ants in two Mexican agroecosystems. Oecologia 55, 114–119 (1982).ADS 
PubMed 
Article 

Google Scholar 
66.Maák, I. et al. Behaviours indicating cannibalistic necrophagy in ants are modulated by the perception of pathogen infection level. Sci. Rep. 10, 17906 (2020).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
67.Doležal, P. & Sehnal, F. A simple method for the detection of imaginal diapause in beetles. J. Appl. Entomol. 131(3), 221–223 (2007).Article 

Google Scholar 
68.Horstmann, K. Untersuchungen ueber den Nahrungserwerb der Waldameisen (Formica polyctena Foerster) im Eichenwald II Abhaengigkeit vom Jahresverlauf und vom Nahrungsangebot. Oecologia (Berlin) 8, 371–390 (1972).ADS 
Article 

Google Scholar 
69.Overbeck, M. & Schmidt, M. Modelling infestation risk of Norway spruce by Ips typographus (L.) in the Lower Saxon Harz Mountains (Germany). For. Ecol. Manag. 266, 115–125 (2012).Article 

Google Scholar 
70.Ayre, G. Response to movement by Formica polyctena Forst. Nature 199, 405–406 (1963).ADS 
Article 

Google Scholar 
71.Thom, D., Seidl, R., Steyrer, G., Krehan, H. & Formayer, H. Slow and fast drivers of the natural disturbance regime in Central European forest ecosystems. For. Ecol. Manag. 307, 293–302 (2013).Article 

Google Scholar 
72.Czechowski, W. & Vepsäläinen, K. Territory size of wood ants (Hymenoptera: Formicidae): A search for limits of existence of Formica polyctena Först., an inherently polygynic and polycalic species. Ann. Zool. 59, 179–187 (2009).Article 

Google Scholar  More

Field survey and samplingField data were collected from 130 study sites spanning a latitudinal gradient of 35.89−50.70° N and a longitudinal gradient of 76.62−122.41°E and covering five provinces across the temperate region in northern China (Xinjiang Autonomous Region, Qinghai Province, Gansu Province, Ningxia Autonomous Region, and Inner Mongolia Autonomous Region; Fig. 2a). Locations for the field study target natural drylands, delineated as areas with aridity level above 0.35 (ref. 30), and represent a large aridity gradient including dry-subhumid (N = 12), semiarid (N = 42), arid (N = 56), and hyperarid (N = 20) regions (Fig. 2a), which are highly vulnerable to expected increases in aridity with human activity and climate change33,71. The aridity level of each site was calculated as 1 – AI, where AI is the ratio of precipitation to potential evapotranspiration38. We obtained AI from the Global Aridity Index and Potential Evapotranspiration Climate database (https://cgiarcsi.community/). The selection of the field sites aimed to minimize the potential impacts of human activity and other disturbances on soil, vegetation, and geomorphological characteristics based on the following three criteria: (i) sites were at least 1 km away from major roads and >50 km from human habitations; (ii) sites were under pristine or unmanaged conditions without visible signs of domestic animal grazing, grass/wood collection, engineering restoration plantings, and infrastructure construction; and (iii) the soil was dry without experiencing rainfall events for at least 3 days prior to sampling. Collectively, our field survey involved a wide range of the abiotic and biotic features of dryland ecosystems across northern China. These sites encompass the 14 soil types, i.e., arenosols, calcisols, cambisols, chernozems, fluvisols, gleysols, greyzems, gypsisols, kastanozems, leptosols, luvisols, phaeozems, solonchaks, and solonetz, and the four main vegetation types44, i.e., typical grassland (dominated by Stipa spp., Leymus spp., Cleistogenes spp., and Agropyron spp.), desert grassland (dominated by Stipa spp., Cleistogenes spp., Suaeda spp., and Artemisia spp.), alpine grassland (dominated by Stipa spp., Leymus spp., Carex spp., and Festuca spp.), and desert (dominated by Reaumuria spp., Salsola spp., Calligonum spp., and Nitraria spp.). Elevation, mean annual temperature, and mean annual precipitation (1970–2000; https://www.worldclim.org/) of the sites varied from 204 to 3,570 m a.s.l. (mean, 1,294 m a.s.l.), from –4.3 to 12.8 °C (mean, 5.0 °C), and from 21 to 453 mm (mean, 195 mm), respectively (Supplementary Table 1).Field sampling was conducted between June and September from 2015 to 2017 (each site was visited once over this period) following well-established standardized protocols as described in refs. 13,34. In brief, three 30 m × 30 m quadrats were established at each site to represent the local vegetation and soil types that covered an area of no less than 10,000 m2. The cover of perennial vegetation was estimated and all perennial plant species were listed by walking steadily along four 1.5 m × 30 m parallel transects (spaced 8 m apart) located within each quadrat using the belt transect method72. Site-level estimate for perennial plant cover was obtained by averaging the values measured in the 12 transects established. After vegetation survey, we located five 1 m × 1 m (for typical grassland, desert grassland, and alpine grassland) or five 5 m × 5 m (for desert) plots within each quadrat (at each corner and the center of the quadrat) to measure site-level plant aboveground and root biomass (g m−2). In each 1 m × 1 m plot, all grasses and dwarf shrubs were harvested to ground level for measurement of aboveground biomass. Five soil cores (7 cm diameter; 0–40 cm depth) per 1-m2 plot were collected randomly, and the roots were removed using a 1-mm sieve and washed cleanly to measure root biomass. All shoot and root samples were dried to constant weight at 65 °C. In each 5 m × 5 m plot, we recorded the number of individuals per dominant shrub species and canopy cover and height of each individual, thereby estimating aboveground and root biomass according to the allometric models developed in previous studies that were conducted in the same regions as sampled here (see Supplementary Table 9 for details). Based on these measurements, we further estimated BNPP. However, BNPP is typically difficult to observe and measure, especially over large spatial scales and environmental gradients as in this study, because the root system is subject to simultaneous growth and turnover73,74. Across our survey areas, ~77–98% of the precipitation occurs between June and September (during the peak-growing season) corresponding to the period of the highest plant above- and belowground biomass34,35,41,75. Therefore, we argue that BNPP can be estimated approximately at each site by the following equation:$$frac{{{{{{rm{Aboveground}}}}}},{{{{{rm{biomass}}}}}}}{{{{{{rm{Root}}}}}},{{{{{rm{biomass}}}}}}}cong frac{{{{{{rm{Aboveground}}}}}},{{{{{rm{net}}}}}},{{{{{rm{primary}}}}}},{{{{{rm{productivity}}}}}},({{{{{rm{ANPP}}}}}})}{{{{{{rm{BNPP}}}}}}}$$
(1)
where both aboveground and root biomass are site-level measurements (g m−2). We used normalized difference vegetation index (NDVI) as a metric for ANPP as explained in recent studies in drylands14,33,70. NDVI data were obtained from the moderate resolution imaging spectroradiometer aboard NASA’s Terra satellites (https://neo.sci.gsfc.nasa.gov/). We used the average NDVI values during our sampling dates as a proxy for ANPP at the site level as described in ref. 14.Five soil cores (0–20 cm depth) per quadrat were then taken randomly under the canopies of the dominant perennial plant species and in bare areas devoid of perennial vegetation, respectively, and then were mixed as one sample for vegetation areas and the other sample for bare ground. When more than one dominant perennial plant species was observed, another three composite samples were collected under the canopies of co-dominant perennial plant species. All vegetation and soil surveys were carried out during the wet season (June to September) when biological activity and productivity are maximal; as such, we do not expect the different sampling times and years or seasonality to be a major factor influencing our conclusions. Collectively, 6–21 soil samples per site were collected, and in total 864 samples were taken and analyzed for each of the seven individual soil functions (see below) and multifunctionality. All soil functions evaluated in the field study were calculated at site level by using a weighted average of the mean values observed in vegetated areas and bare ground by their respective cover13,14,38. After field sampling, the visible pieces of plant material were removed carefully from the soil, which was sieved and divided into three portions. The first portion was air-dried and used for soil organic C, total N, total P, available P, and pH analyses. The second portion was immediately mixed with 2 M KCl and stored at 4 °C for soil ammonium and nitrate analyses. The third portion was immediately frozen at –80 °C for assessing soil microbial diversity.Microcosm experimentIn addition to the large-scale field study described above, we manipulated soil water availability in a microcosm experiment to evaluate the linkages between moisture content, soil microbial diversity, and multifunctionality. It is important to note that our intention is not to directly compare results between these two different approaches [i.e., in the field, measures of soil functions are related to nutrient pools, which we use to associate soil multifunctionality with both plant and soil microbial diversity, whereas in the microcosm experiment the measures of soil functions are related to process rates such as respiration rate and key enzyme activities (see below), which we use to associate soil multifunctionality with microbial diversity in the absence of plants]. Rather, by using an experimental microcosm approach, we aimed to complement the field study and thus further verify the potential increases in aridity to alter the relationship between soil microbial diversity and multifunctionality in the absence of plants. In parallel with the sampling protocols described above, we collected a greater mass of soil (c. 30 kg) under vegetation canopies from one site [i.e., Jingtai country (37.40°N, 104.26°E; Gansu Province, China)]. Soil type, mean annual temperature, mean annual precipitation, and aridity level (1970–2000; https://www.worldclim.org/) of the site is calcisols, 7.9 °C, 205 mm and 0.81, respectively. Following field sampling, the soil was stored immediately at 4 °C until subsequent processing in the laboratory.In brief, a total of 30 experimental microcosms composed of 10 moisture levels with three replicates were established under sterile conditions in a closed incubation chamber (Supplementary Fig. 1a). Each microcosm was filled with 1 kg of soil. These microcosms were incubated at 18.5 °C [the annual mean land surface temperature (1981–2010) for the sampling site; http://data.cma.cn/en], and moisture contents were adjusted and artificially maintained at the ten levels respectively equivalent to 3, 5, 8, 10, 20, 40, 60, 80, 100, and 120% field capacity (27.6%) during the duration of the experiment for 30 days. The corresponding moisture content (%) measured at the end of the experiment varied from 2.03 ± 0.034 to 33.57 ± 1.94, which matched well with differences in moisture conditions among a subset of field soil samples (N = 521; Supplementary Fig. 1b). After incubation, the soil was removed from each microcosm; a portion of the soil was immediately frozen at –80 °C for molecular analysis, and the other fraction was air-dried, sieved, and stored at –20 °C for assessing multiple soil functions as described below.DNA extraction, PCR amplification, and amplicon sequencingFor both the field and experimental studies, we assessed the diversity of soil archaea, bacteria, and fungi using Illumina-based sequencing. Genomic DNA was extracted from 0.5 g of each defrosted soil sample (N = 864 for the field study and N = 30 for the experimental study) using the PowerSoil® DNA Isolation Kit (MO BIO Laboratories, USA) according to the manufacturer’s instructions. For our field study, extracted DNA was pooled at site level, ultimately resulting in 130 composite DNA samples under canopies of vegetation and in bare ground, respectively. Pooling DNA samples may outperform the commonly used method that extracts genomic DNA from mixed soil samples, which could remove large amounts of information on the diversity of soil microorganisms14,22. Negative controls (deionized H2O in place of soil) underwent identical procedures during the extraction to ensure zero contamination in downstream analyses.The V3−V5 regions of the archaeal 16S rRNA gene were amplified using the primer pair Arch344F and Arch915R. Thermal conditions were composed of an initial denaturation of 3 min at 95 °C, ten cycles of touchdown PCR (95 °C for 30 s, annealing temperatures starting at 60 °C for 30 s then decreasing 0.5 °C per cycles, and 72 °C for 1 min), followed by 25 cycles at 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 1 min, with a final extension at 72 °C for 10 min. The primer pair 338F and 806R was used for amplification of the V3−V4 regions of the bacterial 16S rRNA gene. Thermocycling conditions consisted of 3 min at 95 °C and then subjected to 30 amplification cycles of 30 s denaturation at 95 °C, 30 s annealing at 55 °C, followed by 72 °C for 45 s, and a final extension of 72 °C for 10 min. The fungal internal transcribed spacer (ITS) region 1 was amplified using the primer pair ITS1F and ITS2. The amplification conditions involved denaturation at 95 °C for 3 min, 35 cycles of 94 °C for 1 min, 51 °C for 1 min, and 72 °C for 1 min and a final extension at 72 °C for 10 min. Details of primers for each microbial taxa were given in Supplementary Table 10. These primers contained variable length error-correcting barcodes unique to each sample. All amplification reactions were performed in a total volume of 20 μl containing 4 μl of 5× FastPfu Buffer, 2 μl of 2.5 mM dNTPs, 0.8 μl of both the forward and reverse primers, 10 ng of template DNA, and 0.4 μl of FastPfu DNA Polymerase (TransGen Biotech., China). To mitigate individual PCR reaction biases each sample was amplified in triplicate and pooled together. All PCRs were done with the ABI GeneAmp® 9700 Thermal Cycler (Thermo Fisher Scientific, USA). PCR products were evaluated on 2.0% agarose gel with ethidium bromide staining to ensure correct amplicon length, and were gel-purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, USA). Purified amplicons were combined at equimolar concentrations and paired-end sequenced (2 × 300 bp) on an Illumina MiSeq platform (Illumina, USA) at the Majorbio Bio-pharm Technology Co., Ltd. (Shanghai, China) according to standard protocols.Sequence processingInitial sequence processing was conducted with the QIIME pipeline76. Briefly, reads were quality-trimmed with a threshold of an average quality score higher than 20 over 10 bp moving-window sizes and a minimum length of 50 bp. Paired-end reads with at least 10 bp overlap and 2 indicate that the models are different; Supplementary Table 2]. We further assessed whether soil multifunctionality responded more rapidly to aridity than did any individual soil functions. To this end, we explored the presence of aridity thresholds for those relationships that were better fitted by nonlinear regressions (Fig. 2b–i) using the standard protocols developed in ref. 33. The presence of an aridity threshold means that once an aridity level is reached, a given variable either changes abruptly its value (i.e., discontinuous threshold) or its relationship with aridity (i.e., continuous threshold). Hence, a lower aridity threshold indicates that a given variable is more vulnerable to increasing aridity than are others33. We further fitted step (a linear regression that modifies only intercept at a given aridity level) and stegmented (showing changes both in intercept and slope at a given aridity level) regressions for the determination of discontinuous thresholds, and segmented (exhibiting changes only in slope at a given aridity level) regressions for continuous thresholds. Each of these models yields a change point (i.e., threshold) describing the aridity level that evidences the shift in a given nonlinear relationship evaluated. We also used AIC to choose the best threshold model and the corresponding threshold in each case (Supplementary Table 2).We then employed analysis of variance based on type-I sum of squares in a linear mixed-effects model (Eq. (2); Table 1) to test the relationships between the multiple biotic (BNPP, plant species richness, and the soil microbial diversity index) and abiotic (aridity, soil pH, and soil clay content) factors and soil multifunctionality:$${{{{{rm{Soil}}}}}},{{{{{rm{multifunctionality}}}}}} sim; {{{{{rm{Year}}}}}}+{{{{{rm{Plant}}}}}},{{{{{rm{species}}}}}},{{{{{rm{richness}}}}}}\ quad+,{{{{{rm{Soil}}}}}},{{{{{rm{microbial}}}}}},{{{{{rm{diversity}}}}}},{{{{{rm{index}}}}}}\ quad+{{{{{rm{Plant}}}}}},{{{{{rm{species}}}}}},{{{{{rm{richness}}}}}}times {{{{{rm{Soil}}}}}},{{{{{rm{microbial}}}}}},{{{{{rm{diversity}}}}}},{{{{{rm{index}}}}}}+{{{{{rm{Aridity}}}}}}\ quad+,{{{{{rm{Aridity}}}}}}times {{{{{rm{Plant}}}}}},{{{{{rm{species}}}}}},{{{{{rm{richness}}}}}}\ quad+,{{{{{rm{Aridity}}}}}}times {{{{{rm{Soil}}}}}},{{{{{rm{microbial}}}}}},{{{{{rm{diversity}}}}}},{{{{{rm{index}}}}}}+{{{{{rm{BNPP}}}}}}+{{{{{rm{Soil}}}}}},{{{{{rm{pH}}}}}}+{{{{{rm{Soil}}}}}},{{{{{rm{clay}}}}}},{{{{{rm{content}}}}}}\ quad+,{{{{{rm{Elevation}}}}}}+{{{{{rm{Latitude}}}}}}+{{{{{rm{Longitude}}}}}}+(1|{{{{{rm{Soil}}}}}},{{{{{rm{type}}}}}})+(1|{{{{{rm{Vegetation}}}}}},{{{{{rm{type}}}}}})$$
(2)
where × indicates an interaction term. We obtained information on soil clay content (%) from the SoilGrids system (https://soilgrids.org/), and eliminated variation due to different sampling years by first entering the term “Year” into the statistical model41. The elevation, latitude, and longitude of the study sites were included to account for the spatial structure of our dataset13,70. To account for the similarities of soil and vegetation types among study sites we included “Soil type” and “Vegetation type” as random terms.We further simplified the Eq. (2) to focus only on the relationships between aridity, biodiversity, and soil multifunctionality (Eq. (3); Supplementary Fig. 5). We did so because excluding additional biotic and abiotic factors did not change qualitatively the main results presented here (Table 1 and Supplementary Fig. 5), and therefore we used the simplest model to test our hypotheses more clearly. Our simplified model was:$${{{{{rm{Soil}}}}}},{{{{{rm{multifunctionality}}}}}} sim {{{{{rm{Year}}}}}}+{{{{{rm{Plant}}}}}},{{{{{rm{species}}}}}},{{{{{rm{richness}}}}}}\ quad+,{{{{{rm{Soil}}}}}},{{{{{rm{microbial}}}}}},{{{{{rm{diversity}}}}}},{{{{{rm{index}}}}}}\ quad+,{{{{{rm{Aridity}}}}}}+{{{{{rm{Aridity}}}}}}times {{{{{rm{Plant}}}}}},{{{{{rm{species}}}}}},{{{{{rm{richness}}}}}}\ quad+,{{{{{rm{Aridity}}}}}}times {{{{{rm{Soil}}}}}},{{{{{rm{microbial}}}}}},{{{{{rm{diversity}}}}}},{{{{{rm{index}}}}}}\ quad+,{{{{{rm{Aridity}}}}}}times {{{{{rm{Plant}}}}}},{{{{{rm{species}}}}}},{{{{{rm{richness}}}}}}times {{{{{rm{Soil}}}}}},{{{{{rm{microbial}}}}}},{{{{{rm{diversity}}}}}},{{{{{rm{index}}}}}}\ quad+,(1|{{{{{rm{Soil}}}}}},{{{{{rm{type}}}}}})+(1|{{{{{rm{Vegetation}}}}}},{{{{{rm{type}}}}}})$$
(3)
To evaluate how the biodiversity–multifunctionality relationships varied along aridity gradients, we conducted a moving-window analysis as detailed in ref. 69. Briefly, we performed the linear mixed-effects model described in Eq. (3) for a subset window of 60 study sites with the lowest aridity values (this number of sites provided sufficient statistical power for our model), and repeated the same calculations as many times as sites remained (i.e., 70). We then bootstrapped the standardized coefficients of each fixed term within each subset window, which was matched to the average value of aridity across the 60 sites. We fitted linear and nonlinear regressions to the bootstrapped coefficients of biodiversity and its interaction with aridity along aridity gradients (Fig. 3a, b and Supplementary Table 2), and identified the aridity thresholds for the changes in the coefficients of biodiversity (Fig. 3a and Supplementary Table 2) using the same procedure already described above. To provide further support for the aridity thresholds identified here, we also assessed the significance of the bootstrapped standardized coefficients of biodiversity and its interaction with aridity at 95% confidence intervals for each subset window (Fig. 3e). Before fitting threshold regressions, we evaluated whether the variables followed either a unimodal or bimodal distribution using the fitgmdist function in MATLAB (The MathWorks Inc., USA). Our results showed that all variables used for threshold detection presented unimodal distributions (Supplementary Table 11), suggesting that the three threshold regressions mentioned above (i.e., segmented, step, and stegmented) are appropriate in all cases33. We used the chngpt and gam packages in R (http://cran.r-project.org/) to fit segmented/step/stegmented and GAM regressions, respectively. To further check the validity of the thresholds identified, we bootstrapped linear regressions at both sides of each threshold for each variable. We then used the nonparametric Mann–Whitney U-test to compare the slope and the predicted value evaluated before and after each threshold. In all cases, we found significant differences in both of these two parameters (Fig. 3c, d and Supplementary Figs. 2, 3, 6).Given a clear shift in the relationships between plant or microbial diversity and soil multifunctionality occurring at a threshold around an aridity level of 0.8 (Fig. 3), we further used OLS regressions to clarify the relationships between each component of plant or microbial diversity and soil multifunctionality in less and more arid regions separately, as well as across all sites (Fig. 4). To do so, we split our study sites into two groups: sites with aridity 0.8 (more arid regions; N = 76). Moreover, we fitted the mixed-effects model described in Eq. (2) for less and more arid regions separately to ensure the robustness of these bivariate correlations when accounting for multiple biotic and abiotic factors simultaneously, with the exception of using all components of microbial diversity metrics (i.e., soil archaeal, bacterial, and fungal richness) instead of the soil microbial diversity index in the models (Supplementary Table 4). All linear mixed-effects models were performed using the R package lme4. We used a variance inflation factor (VIF) to evaluate the risk of multicollinearity, and selected variables with VIF  0.05). Finally, we also used SEMs to compare the hypothesized direct and indirect relationships between moisture content, microbial diversity, and soil multifunctionality at low and high moisture levels (see an a priori model in Supplementary Figs. 17b and 31c, d). Test of goodness-of-fit for SEMs were same as described above. All the SEM analyses were conducted using AMOS 21.0 (IBM SPSS Inc., USA). Data and code used to perform above analyses are available in figshare98.Reporting SummaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Collection of the samplesTen (10) wood pellet samples were purchased from a different location in B&H, of known suppliers from the market (supermarkets, garden shops, and gas stations). The samples were accompanied by a declaration describing that nine of them were originated from B&H, and one of them was from Italy. Characteristics of collected wood pellet samples (type of wood, energetic value, declared moisture, declared and determined ash amount) are listed in Table 1. All of the samples were analyzed for moisture and ash content. Additionally, in ash samples of mentioned wood pellets, heavy metal concentration (Cd, Co, Cr, Cu, Fe, Mn, Ni, Pb, and Zn) was determined.Table 1 Characteristic of analyzed samples wood pellets.Full size tableAll pellet samples were originated from B&H, purchased from different cities, often used for house heating, instead of sample S3 which was from Italy.Ash determination of wood biomass samplesThe wood pellet samples were oven-dried at 105 °C for 24 h. The content of ash was determined by gravimetric method according to the procedure published by Pan and Eberhardt18 as follows: pellet samples, 1 g (± 0.1 mg) of each was weighed into a previously annealed ceramic pot (m1) and burned in a muffle furnace (Nabertherm) for one hour at 300 °C, following by increasing the temperature to 400 °C for one hour more and then burning the samples for next six hours at 550 °C. The procedure is repeated until a constant mass (m2) was reached. The ash content is determined by the Eq. (1):$${text{Ash content}}, % = frac{{{text{(m}}_{2} – {text{m}}_{{1}} {)}}}{{{text{m}}_{{{text{sample}}}} }} times {100 }{text{.}}$$
(1)
Preparation of samplesThe chemical determinations of the heavy metals in wood pellet ashes (Table 2) were made by wet digestion by soaking the samples in 25 mL of 65% HNO3 in polytetrafluoroethylene (PTFE) vessels. After evaporation of the nitrogen oxides, the vessels were closed and allowed to react for 14 h at 80 °C, following by cooling to room temperature. Then, the digest was filtered, transferred to a 25 mL volumetric flask, and filled up with redistilled water to the mark. All samples and blank were prepared in three replicates19,20,21.Table 2 Heavy metal concentrations (mg kg−1 d.w.) in the wood pellet ashes.Full size tableHeavy metal analysisMetal analyses in ash samples of mentioned wood pellets were performed using a flame atomic absorption spectrometry (Varian AA240FS) for Mn, Fe, Pb, and Zn and graphite furnace (Varian AA240Z) for Cd, Co, Cr, Cu, and Ni. A blank probe was prepared using the same digestion method to avoid the matrix effect. Standard metal solutions used for the calibration graphs were prepared by diluting 1000 mg L−1 stock single-element atomic absorption standard solutions of Cd, Co, Cr, Cu, Fe, Mn, Ni, Pb, and Zn (Certipur Grade, Merck, Germany). Linear calibration graphs with correlation coefficients  > 0.99 were obtained for all analyzed metals. The accuracy of the method was evaluated using the standard reference materials: Fine Fly Ash (CTA-FFA-1, Institute of Nuclear Chemistry and Technology Poland) and Fly Ash from pulverized coal (BCR-038, Institute of reference materials and measurements-IRMM, Belgium). The obtained results were in the range of the reference materials. The detection limit (LOD) and limit of quantification (LOQ) for the nine analyzed metals were calculated based on Xb + 3 SDb and Xb + 10 SDb, respectively, where Xb is the mean concentration of the blank sample (n = 8) and SDb is the standard deviation of the blank for eight readings22. The values of the LOD were: Cd (0.61 µg L−1), Co (0.49 µg L−1), Cr (0.67 µg L−1), Cu (20.10 µg L−1), Fe (83.85 µg L−1), Mn (6.42 µg L−1), Ni (1.12 µg L−1), Pb (23.77 µg L−1), Zn (58.68 µg L−1), and LOQ values were: Cd (1.25 µg L−1), Co (1.41 µg L−1), Cr (1.42 µg L−1), Cu (47.66 µg L−1), Fe (111.2 µg L−1), Mn (16.14 µg L−1), Ni (2.70 µg L−1), Pb (47.73 µg L−1) and Zn (71.05 µg L−1).Pollution evaluationThe metal pollution index (MPI) as the geometric mean of the concentration of all metals found in ashes of wood samples was calculated by the following Eq. (2)23:$${text{MPI}} = left( {{text{C}}_{1} cdot {text{C}}_{2} cdot cdots {text{C}}_{{text{k}}} } right)^{{1/{text{k}}}} ,$$
(2)
where C1 is the concentration value of the first metal, C2 is the concentration value of the second metal, Ck is the concentration value of the kth metal.Evaluation of the presence and the grade of anthropogenic activity were demonstrated through the calculation of the enrichment factor (EF), widely used in environmental issues24. To understand which elements were relatively enriched in the different wood pellet ash samples, the heavy metal enrichment factor was calculated relative to soil values according to Eq. (3)25.$${text{EF}} = frac{{left( {frac{{{text{C}}_{{text{k}}} }}{{{text{E}}_{{{text{ref}}}} }}} right)_{{{text{ashes}}}} }}{{left( {frac{{{text{C}}_{{text{k}}} }}{{{text{E}}_{{{text{ref}}}} }}} right)_{{{text{soil}}}} }},$$
(3)
where Ck is the concentration of the element in the sample or the soil, Eref the concentration of the reference element used for normalization. A reference element is an element commonly stable in the soil characterized by the absence of vertical mobility and/or degradation phenomena. As in many studies as a reference element were Fe, Al, Mn, Sc, or total organic carbon used26,27. Therefore Fe has been chosen as reference material in this study. Iron is one of the major constituents of soil, as well as the average chemical constituent of the upper continental crust26.Health risk assessmentThe general exposure equations used in this study were adapted according to the US Environmental Protection Agency guidance28,29,30. The daily exposure (D) to heavy metals via wood pellet ash was calculated for the three main routes of exposure: (i) direct ingestion of ash particles (Ding); (ii) inhalation of suspended particles via mouth and nose (Dinh); and (iii) dermal absorption to skin adhered ash particles (Ddermal). Equations (4) to (6) were used to calculate exposure via ingestion, inhalation, and dermal route, respectively22,31.$${text{D}}_{{{text{ing}}}} = {text{ C }} cdot frac{{{text{ IngR }} cdot {text{ EF }} cdot {text{ ED}}}}{{{text{BW }} cdot {text{ AT}}}}{ } cdot {text{CF}}1{, }$$
(4)
$${text{D}}_{{{text{inh}}}} = {text{ C }} cdot frac{{{text{ InhR}} cdot {text{ EF }} cdot {text{ ED}}}}{{{text{PEF }} cdot {text{ BW }} cdot {text{ AT}}}}{, }$$
(5)
$${text{D}}_{{{text{dermal}}}} = {text{ C }} cdot frac{{{text{ SA }} cdot {text{ SL }} cdot {text{ABS }} cdot {text{EF }} cdot {text{ ED}}}}{{{text{BW }} cdot {text{ AT}}}}{ } cdot {text{CF}}1{, }$$
(6)

where c (mg kg−1) is the heavy metals concentrations in ash samples; IngR (mg day−1) is the conservative estimates of dust ingestion rates, 50 for adults, 200 for children30,32; InhR (m3 h−1) is the inhalation rate, 2.15 for adults, 1.68 for children32; EF (h year−1) is the exposure frequency, 1225 for adults and children22; ED (years) is the exposure duration, 70 for adults, 6 for children22; BW (kg) is the body weight, 80 for adults, 18.60 for children32; AT (days) is the averaging time, 25,550 for adults, 2190 for children22; PEF is the particle emission factor (m3 kg−1), 6.80 × 108 for adults and children31; SA (cm3) is the exposed skin area, 6840 for adults, 2550 for children32; SL (mg cm−2) is the skin adherence factor, 0.22 for adults, 0.27 for children32; ABS is the dermal absorption factor, 0.001 for adults and children31; CF1 is the unit conversation factor, 10–6 for adults and children22.The potential non-carcinogenic risk for each metal was estimated using the Hazard coefficient (HQ), as suggested by US EPA33. The HQ under various routes of exposure such as ingestion (HQing), inhalation (HQinh), and dermal (HQdermal) was calculated as a ratio of daily exposure (D) to reference dose of each metal (RfD) according to Eq. (7)32.$${text{HQ}}_{{text{k}}} = frac{{{text{D}}_{{text{k}}} }}{{{text{RfD}}}},$$
(7)

where k is ingestion, inhalation, or dermal route. The total hazard index (HI) of heavy metal for all routes of exposure was calculated as a sum of HQing, HQinh, and HQdermal as given in Eq. (8)34.$${text{HI}} = {text{ HQ}}_{{text{ing }}} + {text{ HQ}}_{{text{inh }}} + {text{ HQ}}_{{text{dermal }}} .$$
(8)
The carcinogenic risk (Risk) for potential carcinogenic metals was calculated by multiplying the doses by the corresponding slope factor (SF), as given in Eq. (9)35. The carcinogenic oral, inhalation, and dermal SF, as well as dermal absorption toxicity values, were provided from the Integrated Risk Information System30. The reference doses for Pb were taken from the Guidelines for Drinking Water Quality published by the World Health Organization36.$${text{Risk}} = { }mathop sum limits_{{{text{k}} = 1}}^{{text{n}}} {text{D}}_{{text{k}}} cdot {text{ SF}}_{{text{k}}} ,$$
(9)
where SF is the cancer slope factor for individually metal and k route of exposure (ingestion, inhalation, or dermal path). The total cancer risk (Risktotal) of potential carcinogens was calculated as the sum of the individual risk values using the following Eq. (10).$${text{Risk}}_{{{text{total}}}} = {text{Risk}}_{{{text{ing}}}} + {text{Risk}}_{{{text{inh}}}} + {text{Risk}}_{{{text{dermal}}}} .$$
(10) More

1.Binford, L. R. Constructing Frames of Reference: An Analytical Method for Archaeological Theory Building Using Hunter-Gatherer and Environmental Data Sets (Univ. of California Press, 2001).2.Kelly, R. L. The Lifeways of Hunter-Gatherers: The Foraging Spectrum (Cambridge Univ. Press, 2013).3.Tallavaara, M., Eronen, J. T. & Luoto, M. Productivity, biodiversity, and pathogens influence the global hunter-gatherer population density. Proc. Natl Acad. Sci. USA 115, 1232–1237 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
4.Eriksson, A. et al. Late Pleistocene climate change and the global expansion of anatomically modern humans. Proc. Natl Acad. Sci. USA 109, 16089–16094 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Gurven, M. D. & Davison, R. J. Periodic catastrophes over human evolutionary history are necessary to explain the forager population paradox. Proc. Natl Acad. Sci. USA https://doi.org/10.1073/pnas.1902406116 (2019).6.Tallavaara, M., Luoto, M., Korhonen, N., Järvinen, H. & Seppä, H. Human population dynamics in Europe over the Last Glacial Maximum. Proc. Natl Acad. Sci. USA 112, 8232–8237 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
7.Bradshaw, C. J. A. et al. Minimum founding populations for the first peopling of Sahul. Nat. Ecol. Evol. 3, 1057–1063 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
8.Kavanagh, P. H. et al. Hindcasting global population densities reveals forces enabling the origin of agriculture. Nat. Hum. Behav. 2, 478–484 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
9.Porter, C. C. & Marlowe, F. W. How marginal are forager habitats? J. Archaeol. Sci. 34, 59–68 (2007).Article 

Google Scholar 
10.Reyes-García, V. & Pyhälä, A. Hunter-Gatherers in a Changing World (Springer International Publishing, 2017).11.Lee, R. B. & Daly, R. The Cambridge Encyclopedia of Hunters and Gatherers (Cambridge Univ. Press, 1999).12.Kitanishi, K. Seasonal changes in the subsistence activities and food intake of the Aka hunter-gatherers in northeastern Congo. Afr. Study Monogr. 16, 73–118 (1995).
Google Scholar 
13.Timmermann, A. & Friedrich, T. Late Pleistocene climate drivers of early human migration. Nature 538, 92–95 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
14.Keeley, L. H. Hunter-gatherer economic complexity and ‘population pressure’: a cross-cultural analysis. J. Anthropol. Archaeol. 7, 373–411 (1988).15.Fisher, J. B., Huntzinger, D. N., Schwalm, C. R. & Sitch, S. Modeling the terrestrial biosphere. Annu. Rev. Environ. Resour. 39, 91–123 (2014).Article 

Google Scholar 
16.Pachzelt, A., Forrest, M., Rammig, A., Higgins, S. I. & Hickler, T. Potential impact of large ungulate grazers on African vegetation, carbon storage and fire regimes. Glob. Ecol. Biogeogr. 24, 991–1002 (2015).Article 

Google Scholar 
17.Zhu, D. et al. The large mean body size of mammalian herbivores explains the productivity paradox during the Last Glacial Maximum. Nat. Ecol. Evol. 2, 640–649 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
18.Dyble, M., Thorley, J., Page, A. E., Smith, D. & Migliano, A. B. Engagement in agricultural work is associated with reduced leisure time among Agta hunter-gatherers. Nat. Hum. Behav. 3, 792–796 (2019).19.Hill, K., Kaplan, H., Hawkes, K. & Hurtado, A. M. Men’s time allocation to subsistence work among the Ache of eastern Paraguay. Hum. Ecol. 13, 29–47 (1985).Article 

Google Scholar 
20.Hill, K., Hawkes, K., Hurtado, M. & Kaplan, H. Seasonal variance in the diet of Ache hunter-gatherers in eastern Paraguay. Hum. Ecol. 12, 101–135 (1984).CAS 
Article 

Google Scholar 
21.Marlowe, F. W. et al. Honey, Hadza, hunter-gatherers, and human evolution. J. Hum. Evol. 71, 119–128 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
22.Cordain, L. et al. Plant–animal subsistence ratios and macronutrient energy estimations in worldwide hunter-gatherer diets. Am. J. Clin. Nutr. 71, 682–692 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
23.Gurven, M. & Kaplan, H. Longevity among hunter-gatherers: a cross-cultural examination. Popul. Dev. Rev. 33, 321–365 (2007).Article 

Google Scholar 
24.Klein Goldewijk, K., Beusen, A. & Janssen, P. Long-term dynamic modeling of global population and built-up area in a spatially explicit way: HYDE 3.1. Holocene 20, 565–573 (2010).Article 

Google Scholar 
25.Marlowe, F. W. Hunter-gatherers and human evolution. Evol. Anthropol. 14, 54–67 (2005).Article 

Google Scholar 
26.Burger, J. R. & Fristoe, T. S. Hunter-gatherer populations inform modern ecology. Proc. Natl Acad. Sci. USA 115, 1137–1139 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
27.Hurtado, A. M. & Hill, K. R. Seasonality in a foraging society: variation in diet, work effort, fertility, and sexual division of labor among the Hiwi of Venezuela. J. Anthropol. Res. 46, 293–346 (1990).Article 

Google Scholar 
28.Wilmsen, E. N. Studies in diet, nutrition, and fertility among a group of Kalahari Bushmen in Botswana. Soc. Sci. Inf. 21, 95–125 (1982).Article 

Google Scholar 
29.Lee, R. B. in Man the Hunter (eds. Lee, R. B. & DeVore, I.) 30–48 (Aldine de Gruyter, 1968).30.Hamilton, M. J., Milne, B. T., Walker, R. S. & Brown, J. H. Nonlinear scaling of space use in human hunter-gatherers. Proc. Natl Acad. Sci. USA 104, 4765–4769 (2007).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
31.Messer, E. Anthropological perspectives on diet. Annu. Rev. Anthropol. 13, 205–249 (1984).Article 

Google Scholar 
32.Testart, A. et al. The significance of food storage among hunter-gatherers: residence patterns, population densities, and social Inequalities [and Comments and Reply]. Curr. Anthropol. 23, 523–537 (1982).Article 

Google Scholar 
33.Winterhalder, B. Diet choice, risk, and food sharing in a stochastic environment. J. Anthropol. Archaeol. 5, 369–392 (1986).Article 

Google Scholar 
34.Kelly, R. L., Pelton, S. R. & Robinson, E. in Towards a Broader View of Hunter-Gatherer Sharing (eds Lavi, N. & Friesem, D. E.) Ch. 10 (McDonald Institute for Archaeological Research, 2019).35.Joannes-Boyau, R. et al. Elemental signatures of Australopithecus africanus teeth reveal seasonal dietary stress. Nature 572, 112–115 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Smits, S. A. et al. Seasonal cycling in the gut microbiome of the Hadza hunter-gatherers of Tanzania. Science 357, 802–806 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Barnosky, A. D. Assessing the causes of Late Pleistocene extinctions on the continents. Science 306, 70–75 (2004).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
38.Henrich, J. Demography and cultural evolution: how adaptive cultural processes can produce maladaptive losses—the Tasmanian case. Am. Antiq. 69, 197–214 (2004).Article 

Google Scholar 
39.Powell, A., Shennan, S. & Thomas, M. G. Late Pleistocene demography and the appearance of modern human behavior. Science 324, 1298–1301 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
40.D’Alpoim Guedes, J. A., Crabtree, S. A., Bocinsky, R. K. & Kohler, T. A. Twenty-first century approaches to ancient problems: climate and society. Proc. Natl Acad. Sci. USA 113, 14483–14491 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
41.Cegielski, W. H. & Rogers, J. D. Rethinking the role of Agent-based modeling in archaeology. J. Anthropol. Archaeol. 41, 283–298 (2016).Article 

Google Scholar 
42.Axtell, R. L. et al. Population growth and collapse in a multiagent model of the Kayenta Anasazi in long house valley. Proc. Natl Acad. Sci. USA 99, 7275–7279 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
43.Hayden, B. Research and development in the stone age: technological transitions among hunter-gatherers. Curr. Anthropol. 22, 519–548 (1981).Article 

Google Scholar 
44.Itkonen, T. I. Suomen Lappalaiset Vuoteen 1945. Ensimmäinen Osa (WSOY, 1848).45.Kirby, K. R. et al. D-PLACE: a global database of cultural, linguistic and environmental diversity. PLoS ONE 11, e0158391 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
46.MODIS NPP (MOD17A3) (NTSG, accessed 12 March 2015); http://files.ntsg.umt.edu/data/NTSG_Products/MOD17/47.Defries, R.S. et al. ISLSCP II Continuous Fields of Vegetation Cover, 1992–1993 (ORNL DAAC, 2009); https://doi.org/10.3334/ORNLDAAC/93148.Bodesheim, P., Jung, M., Gans, F., Mahecha, M. D. & Reichstein, M. Upscaled diurnal cycles of land–atmosphere fluxes: a new global half-hourly data product. Earth Syst. Sci. Data 10, 1327–1365 (2018).Article 

Google Scholar 
49.Šímová, I. & Storch, D. The enigma of terrestrial primary productivity: measurements, models, scales and the diversity–productivity relationship. Ecography 40, 239–252 (2017).Article 

Google Scholar 
50.Bontemps, S. et al. Consistent global land cover maps for climate modelling communities: current achievements of The ESA Land Cover CCI. In Proc. ESA Living Planet Symposium 2013 (ESA, 2013).51.Bliege Bird, R. & Bird, D. W. Why women hunt—risk and contemporary foraging in a western desert aboriginal community. Curr. Anthropol. 49, 655–693 (2008).Article 

Google Scholar 
52.Bliege Bird, R., Codding, B. F. & Bird, D. W. What explains differences in men’s and women’s production? Hum. Nat. 20, 105–129 (2009).PubMed 
Article 
PubMed Central 

Google Scholar 
53.Reyes-García, V., Díaz-Reviriego, I., Duda, R., Fernández-Llamazares, Á. & Gallois, S. ‘Hunting otherwise’—women’s hunting in two contemporary forager-horticulturalist societies. Hum. Nat. 31, 203–221 (2020).PubMed 
Article 
PubMed Central 

Google Scholar 
54.Krinner, G. et al. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. Global Biogeochem. Cycles 19, GB1015 (2005).55.Hamilton, M. J., Lobo, J., Rupley, E., Youn, H. & West, G. B. The ecological and evolutionary energetics of hunter‐gatherer residential mobility. Evol. Anthropol. 25, 124–132 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
56.Abrams, P. A. & Ginzburg, L. R. The nature of predation: prey dependent, ratio dependent or neither? Trends Ecol. Evol. 15, 337–341 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
57.Winterhalder, B., Baillargeon, W., Cappelletto, F., Randolph Daniel, I. & Prescott, C. The population ecology of hunter-gatherers and their prey. J. Anthropol. Archaeol. 7, 289–328 (1988).Article 

Google Scholar 
58.Illius, A. W. & O’Connor, T. G. Resource heterogeneity and ungulate population dynamics. Oikos 89, 283–294 (2000).Article 

Google Scholar 
59.Golley, F. B. Energy values of ecological materials. Ecology 42, 581–584 (1961).Article 

Google Scholar 
60.Herbers, J. M. Time resources and laziness in animals. Oecologia 49, 252–262 (1981).PubMed 
Article 
PubMed Central 

Google Scholar 
61.Raichlen, D. A. et al. Sitting, squatting, and the evolutionary biology of human inactivity. Proc. Natl Acad. Sci. USA https://doi.org/10.1073/pnas.1911868117 (2020).62.Abrams, H., Jr. in Food and Evolution (eds Harris, M. & Ross, E.) 207–223 (Temple Univ. Press, 1987).63.Hanya, G. & Aiba, S. Fruit fall in tropical and temperate forests: implications for frugivore diversity. Ecol. Res. 25, 1081–1090 (2010).Article 

Google Scholar 
64.Gherardi, L. A. & Sala, O. E. Global patterns and climatic controls of belowground net carbon fixation. Proc. Natl Acad. Sci. USA https://doi.org/10.1073/pnas.2006715117 (2020).65.van Zonneveld, M. et al. Human diets drive range expansion of megafauna-dispersed fruit species. Proc. Natl Acad. Sci. USA 115, 3326–3331 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
66.Max R., Hannah R. & Esteban Ortiz-Ospina World Population Growth (GCDL, 2020); https://ourworldindata.org/world-population-growth67.Pontzer, H. et al. Metabolic acceleration and the evolution of human brain size and life history. Nature 533, 390 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
68.Viovy, N. CRUNCEP Version 7—Atmospheric Forcing Data for the Community Land Model (Research Data Archive at the National Center for Atmospheric Research, Computational and Information Systems Laboratory, 2018); https://doi.org/10.5065/PZ8F-F01769.Sobol, I. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math. Comput. Simul. 55, 271–280 (2001).Article 

Google Scholar  More

Two databases including yield, management, and weather data for maize (n = 17,013) and soybean (n = 24,848) involving US crop performance trials conducted in 28 states between 2016 to 2018 for maize and between 2014 to 2018 for soybean, were developed (Fig. 1). Crop yield and management data were obtained from publicly available variety performance trials which are typically performed yearly in several locations across each state (see methods for more information). Final databases were separated in training (80% of database) and testing (20% of database) datasets using stratified sampling by year, use of irrigation, and soil type. For each crop, an extreme gradient boosting (XGBoost, see methods for more information) algorithm to estimate yield based on soil type and weather conditions (E), seed traits (G) and management practices (M) was developed (see variables listed in Tables S1 and S2 for maize and soybean, respectively, and data science workflow in Fig. S1).Figure 1Locations where maize and soybean trials were performed during the examined period. The map was developed in ArcGIS Pro 2.8.0 (https://www.esri.com).Full size imageThe developed algorithms exhibited a high degree of accuracy when estimating yield in independent datasets (test dataset not used for model calibration) (Fig. 2). For maize, the root mean square error (RMSE) and mean absolute error (MAE) was a respective 4.7 and 3.6% of the dataset average yield (13,340 kg/ha). For soybean, the respective RMSE and MAE was 6.4 and 4.9% of the dataset average yield (4153 kg/ha). As is evident in the graphs (Fig. 2), estimated yields exhibited a high degree of correlation with actual yields for both algorithms in the independent datasets. For maize and soybean, 72.3 and 60% of cases in the test dataset deviated less than 5% from actual yields, respectively. Maximum deviation for maize and soybean reached 43 and 70%, respectively. Data points with deviations greater than 15% from actual yield were 1.5% in maize and 3.6% in soybean databases. These results suggest that the developed algorithms can accurately estimate maize and soybean yields utilizing database-generated information involving reported environmental, seed genetic, and crop management variables.Figure 2Actual versus algorithm-derived maize (left) and soybean (right) yield in test datasets. Black solid line indicates y = x, red short-dashed lines, black dashed lines, and red long-dashed lines indicate ± 5, 10, and 15% deviation from the y = x line. RMSE, root mean square error; MAE, mean absolute error; r2, coefficient of determination; n = number of observations. Each observation corresponds to a yield of an individual cropping system in a specific environment (location-year).Full size imageIn contrast to statistical models, ML algorithms can be complex, and the effect of single independent variables may not obvious. However, accumulated local effects (ALE) plots14 can aid the understanding and visualization of important and possibly correlated features in ML algorithms. For both crops, indicatively important variables included sowing date, seeding rate, nitrogen fertilizer (for maize), row spacing (for soybean) and June to September cumulative precipitation (Fig. 3). Across the entire region and for both crops, the algorithm-derived trends suggest that above average yields occur in late April to early May sowing dates, but sharply decrease thereafter. Similar responses have been observed in many regional studies across the US for both, maize15,16,17,18 and soybean19. Similarly, simulated yield curves due to increasing seeding rate are in close agreement with previous maize20,21 and soybean22 studies. The maize algorithm has captured the increasing yield due to increasing N fertilizer rate. The soybean algorithm suggests that narrower row spacing resulted in above average yield compared to wider spacing. Such response has been observed in many regions across the US23. Season cumulative precipitation between 400 and 700 mm resulted in above average yields for both crops.Figure 3Accumulated local effect plots for maize sowing date (A), seeding rate (B), Nitrogen fertilizer rate (C), and cumulative precipitation between June and September (mm) (D), and soybean sowing date (E), seeding rate (F), row spacing (G), and cumulative precipitation between June and September (mm) (H).Full size imageThe responses in the ALE plots (Fig. 3) suggest that these algorithms have captured the general expected average responses for important single features. Nevertheless, our databases include hundreds of locations with diverse environments across the US and site-specific crop responses which may vary due to components of the G × E × M interaction. We argue that, instead of examining a single or low-order management interactions, site-specific evaluation of complex high order interactions (a.k.a. cropping systems) can reveal yield differences that current research approaches cannot fully explore and quantify. For example, sowing date exerts a well-known impact on maize and soybean yield. For each crop separately, by creating a hypothetical cropping system (a single combination of all management and traits in Tables S1 and S2) in a randomly chosen field in south central Wisconsin (latitude = 43.34, longitude = -89.38), and by applying the developed algorithms, we can generate estimates of maize and soybean yield. For that specific field and cropping system (out of the vast number of management combinations a farmer can choose from), maize yield with May 1st sowing was 711 kg/ha greater (6% increase) than June sowing (Fig. 4A). By creating scenarios with 256 background cropping system choices (Table S3), the resultant algorithm-derived yield estimate difference for the same sowing date contrast (averaged across varying cropping systems) was smaller but still positive (3% increase), although the range of possible yield differences was wider (Fig. 4B). However, when comparing, instead of averaging, the estimated yield potential among the simulated cropping systems, a 2903 kg/ha yield difference (25% difference) was observed (Fig. 4C). Interestingly, when focusing on the early sown fields that were expected to exhibit the greatest yield, the same yield difference was observed (Fig. 4D). This result shows that sub-optimal background management can mitigate the beneficial effect of early sowing (Table S4).Figure 4Maize yield difference (in kg/ha and percentage) due to sowing date (May 1st vs. June 1st) for a single identical background cropping system (A), maize yield difference due to sowing date when averaged across 256 (3 years × 256 cropping systems = 768 year-specific yields) (B), maize yield variability in each of the 256 cropping systems (C), and maize yield variability in each of the 128 cropping systems with early sowing (D). Soybean yield difference due to sowing date (May 1st vs June 1st) for a single identical background cropping system (E), soybean yield difference due to sowing date when averaged across 128 (5 years × 128 cropping systems = 640 year-specific yields) (F), soybean yield in each of the 128 cropping systems (G) and soybean yield variability due in each of the 64 cropping systems with early sowing (H). Within each panel, the horizontal red and grey lines indicate the boxplot with maximum and minimum yield, respectively. In the left four panels, boxes delimit first and third quartiles; solid lines inside boxes indicate median and green triangles indicate means. Upper and lower whiskers extend to maximum and minimum yields. Each maize and soybean cropping system is a respective 8-way and a 7-way interaction of management practices in a randomly chosen field in Wisconsin, USA (Table S3 and S5, respectively).Full size imageIn the case of soybean, a May 1st sowing resulted in greater yield (588 kg/ha; a 14% increase) than a June 1st in the single background cropping system (Fig. 4E). The result was consistent when yield differences due to sowing date were averaged across 128 background cropping system choices (Table S5) (Fig. 4F). Similar to what was observed in maize, among all cropping systems, yield varied by 1704 kg/ha (44% difference) (Fig. 4G). When focusing only on the early sown fields, a 1181 kg/ha yield difference (27% yield increase) was observed (Fig. 4H). In agreement with maize, this result highlights the importance of accounting for sub-optimal background management which can mitigate the beneficial effect of early sowing (Table S6).We note here the ability of farmers to change management practices can be limited due to an equipment constraint (e.g., change planter unit row width) or simply impossible (e.g., change the previous year’s crop). Thus, recommended management practices that were evaluated in studies that used specific background management may not be applicable in some instances. The benefits of the foregoing approach, which involves extensive up-to-date agronomic datasets and high-level computational programing, can have important and immediate implications in future agricultural trials. Our approach allows for more precise examination of complex management interactions in specific environments (soil type and growing season weather) across the US (region covered in Fig. 1). The ability to extract single management practice information (even across cropping systems) is also possible by utilizing ALE plots, or by calculation of the frequency at which a given level/rate of a management practice appeared among the highest yielding cropping systems (Tables S4 and S6).Among all available 30-d weather variables, many were strongly correlated in both crop databases (Figs. S2 and S3 for maize and soybean, respectively). Models using all 30-d interval variables with r  More

Sampling and genotypingWe collected hair and cheek swab samples from 77 men from geographically and linguistically distinct groups of the Caucasus: Kartvelian speakers from Georgia and Turkey, Northeast Caucasian speakers and Turkic speakers from the Russian Federation and Armenian speakers from Georgia’s southern province of Javakheti, descendants of the families displaced from Mush and Erzurum provinces of eastern Turkey in the early nineteenth century (Table 1, Fig. 1). To maximize the representativeness of the genetic signature of each population, the samples were collected from locals with no ancestors from outside of the respective ethnic/geographic population over the last three generations. DNA was extracted from follicles of 10–12 male chest hairs and cheek swab samples. Extraction was performed using Qiagen DNeasy Blood and Tissue kit, following the manufacturer’s recommendations (Qiagen, Valencia, CA, USA). The DNA samples were genotyped for 693,719 autosomal and 17,678 X-chromosomal SNPs by Family Tree DNA (FTDNA—Gene By Gene, Ltd, Houston, TX, www.familytreedna.com).Table 1 Modern study populations of the Caucasus. Latitude and longitude georeference population hubs.Full size tableFigure 1The distribution of the study populations: averaged centroids of ancient populations (uniquely colored points in the main map, see Table 2 for details) and hubs of the modern Caucasian populations (identified in the inset map, see Table 1 for details). Glacial human refugia extracted from Gavashelishvili and Tarkhnishvili5 are shaded in purple. The map is generated using QGIS Desktop 3.10.6-A Coruña (https://qgis.org).Full size imageOur dataset of modern Caucasian genotypes was supplemented with published 10 modern Mbuti (Supplementary Table S1) and 122 Upper Paleolithic-Mesolithic human genotypes, retrieved as a part of 1240 K dataset from David Reich’s Lab website, Harvard University (https://reich.hms.harvard.edu/downloadable-genotypes-present-day-and-ancient-dna-data-compiled-published-papers; see Supplementary Table S2 for details). The ancient genotypes were selected such that they either dated from the LGM or fell within the glacial refugia identified by Gavashelishvili and Tarkhnishvili5. We did so in order to maximize the genetic signature of potential refugial populations in our analysis. We divided the ancient genotypes into 2000-year-long intervals, and then grouped each of these intervals into geographic units (hereafter ancient populations, Table 2, Fig. 1). The modern and ancient genotypes were merged using PLINK 1.9 (PLINK 1.9: www.cog-genomics.org/plink/1.9/27.Table 2 Ancient study populations. The ancient genotypes are divided into 2000-year-long intervals, and then each of these intervals is grouped into geographic units (i.e. ancient populations). Age, latitude and longitude are averaged across each ancient population (see Supplementary Table S2 for details).Full size tableEthics statementThe research team members, through their contacts in the studied communities, inquired whether locals would voluntarily participate in genetic research that would help clarify the genetic makeup of the Caucasus. A verbal agreement was made with volunteer donors of DNA samples, according to which the results would be communicated, electronically or in hard copy, with participants individually. Participants were informed that, upon the completion of the lab work, the research would be published without mentioning the names of sample donors. Those who agreed provided us with the envelopes containing their chest hairs or cheek swab samples, with the birthplace of their ancestors (last three generations) written on the envelope or a piece of paper. In accordance with the preferences of the sample donors, the agreement was verbal and not written. The envelopes and papers are stored as evidence of voluntary provision of the samples and the related information. Analysis of data was done anonymously, using only location and ethnic information; only the first and third authors of the manuscript had access to names associated with the samples. Therefore, this study was based on noninvasive and nonintrusive sampling (volunteers provided hair and swab samples they collected themselves), and the information destined for open publication does not contain any personal information. The study methodology and the procedure of verbal consent was discussed in detail with and approved by the members of the Ilia State University Commission for Ethical Issues before the field survey started, and the commission decided that formal ethical approval was not needed for conducting this study. This is confirmed in a letter from the commission chairman, a copy of which has been provided to the journal editor as part of the submission process.Genetic affinity and geographyFirst, we measured genetic affinity between the modern Caucasian populations, and between the modern populations and the ancient populations of hunter-gatherers, and then tested whether the genetic affinity between these populations was determined by geographic features. Data were mapped using QGIS Desktop 3.10.6-A Coruña, whereas graphs were created using the “ggplot2” package28 in R version 3.5.229.To evaluate genetic affinities and structure of the modern populations, we used Wright’s fixation index (Fst), inbreeding coefficient, admixture analysis and the principal component analysis (PCA). For these procedures we filtered the raw SNP genotypes in PLINK 1.9, first removing all SNPs with the minor allele frequency  0.3, calculated in windows of 50 bp size and 10 bp steps (–maf 0.05 –indep-pairwise 50 10 0.3). Since all individuals in our dataset possess a single copy of the X-chromosome, we did not expect any differential ploidy bias, and SNPs on the X were treated similarly to those on the autosomes. Fst pairwise values were calculated using the smartpca program of EIGENSOFT30 with default parameters, inbreed: YES, and fstonly: YES. The relationship between the modern populations based on Fst values was visualized by constructing a neighbor-joining tree using the “ape” package31 in R version 3.5.2. The average and standard deviation of the inbreeding coefficient for each population was calculated using “fhat2” estimate of PLINK 1.9. The LD pruned genotypes were used in ADMIXTURE 1.3.032, performed in unsupervised mode in order to infer the population structure from the modern individuals. The number of clusters (k) was varied from 2 to 7 and the fivefold cross-validation error was calculated for each k33. We conducted principal components analysis in the smartpca program of EIGENSOFT30, using default parameters and the lsqproject: YES and numoutlieriter: 0 options. Eigenvectors of principal components were inferred with the modern populations from the Caucasus, while the ancient populations were then projected onto the PCA plots. We also assessed the relatedness between sampled individuals using kinship coefficients estimated by KING34.To quantify genetic affinities between the modern and ancient populations, we used the programs qp3Pop and qpDstat in the ADMIXTOOLS suite (https://github.com/DReichLab35 for f3- and f4-statistics, respectively. f3-statistics of the form f3(X,Y,Outgroup) measure the amount of shared genetic drift of populations X and Y after their divergence from an outgroup. We used an ancient population and a modern Caucasian population for X, Y and Mbuti as an outgroup. f4-statistics of the form f4(Outgroup,Test;X,Y) show if population Test is equally related to X and Y or shares an excess of alleles with either of the two. In the f4-statistic calculation we used Mbuti for Outgroup, a modern population of the Caucasus for Test, and X and Y for contemporaneous ancient populations. This meant that f4  0 indicated higher genetic affinity between the test population and Y.To quantify geographic features, we derived least-cost paths and measured least-cost distances (LCD) between the modern and ancient populations using the Least Cost Path Plugin for QGIS. The computation of LCD considers a friction grid that is a raster map where each cell indicates the relative difficulty (or cost) of moving through that cell. A least-cost path minimizes the sum of frictions of all cells along the path, and this sum is the least-cost distance (LCD). For impedance to human movement and expansion, we used 15 geographic features (Table 3). All gridded geographic features (i.e. raster layers) were resampled to a resolution of 1 km using the nearest-neighbor assignment technique. All possible subsets of the 15 geographic features, that did not cancel out each other, were used to calculate different variables of LCD. We assumed that most human movements occurred during climate warming events when the earth’s surface was not dramatically different from that of today, and hence used the current data of the geographic features.Table 3 Geographic features used in combinations to calculate least-cost distances (LCD) between ancient populations and modern Caucasians.Full size tableLinking genetic affinity and geographyGeneralized additive models (GAMs) were used to fit the outgroup f3-statistic to time and variously calculated LCD between the modern and ancient populations using the “mgcv” package36 in R version 3.5.2. Time between the modern and ancient populations was measured in BP (years before present, defined by convention as years before 1950 CE). We used GAMs because without any assumptions they are able to find nonlinear and non-monotonic relationships. GAMs were fitted using a Gamma family with a log link function. Penalized thin plate regression splines were used to represent all the smooth terms. The restricted maximum likelihood (REML) estimation method was implemented to estimate the smoothing parameter because it is the most robust of the available GAM methods36.Model and variable selection were performed by exploring LCD, time BP and the interaction term. The predictive power of the models was evaluated through a tenfold cross-validation. The cross-validation of many models was handled through R’s parallelization capabilities37,38. The best model was selected by the mean squared error of the cross-validation. Akaike’s Information Criterion (AIC) is generally used as a means for model selection. However, we preferred cross-validation for model selection because AIC a priori assumes that simpler models with the high goodness of fit are more likely to have the higher predictive power, while cross-validation without any a priori assumptions measures the predictive performance of a model by efficiently running model training and testing on the available data.We additionally validated the effect of different subsets of geographic features by assessing the relationship between statistically significant values of f4-statistic (i.e. |Z| > 3) and each subset. The relationship between f4-statistic of the form of f4(Outgroup,Test;X,Y) and geographic features was determined by measuring the agreement between the negative/positive signs of f4-statistic and the difference in LCD (LCD.D) for each pair of contemporaneous ancient populations X and Y. LCD.D was calculated as (LCD1–LCD2), where LCD1 was least-cost distance between the test population and X, and LCD2 was least-cost distance between the test population and Y. LCD.D  0 indicated less least-cost distance between Test and Y. So, the same sign of f4 and LCD.D values indicated agreement between geographic proximity and genetic affinity. We used Cohen’s kappa39 to measure the agreement.In order to test if geographic features (Table 3) accounted for present-day genetic differentiation in the Caucasus, we measured the relationship between Fst and LCD across the modern populations using the Mantel test in the “vegan” package40 in R version 3.5.2. In addition, we checked whether contribution from ancient samples was related to today’s genetic differentiation. To do so, we calculated median of f3-statistic of ancient populations of each geographic grouping (e.g. the following 6 populations made up one group: Balkans 39,950–41,950 BP, Balkans 37,950–39,950 BP, Balkans 31,950–33,950 BP, Balkans 9950–11,950 BP, Balkans 7950–9950 BP, Balkans 5950–7950 BP). Then we measured the manhattan distance of f3 median values of all combinations of the geographic groupings between the modern populations and compared the results to Fst and LCD using the Mantel test. More

Regular structures and isothermal theoremFor networks where each node represents a single individual, the isothermal theorem of evolutionary graph theory shows that the fixation probability is the same as the fixation probability of a well-mixed population if the temperature distribution is homogeneous across the whole population1. The temperature of a node defined as the sum over all the weights leads to that node. This theorem extends to structured meta-populations for any migration probability (lambda ): If the underlying structure of the meta-population that connects the patches is a regular network and the local population size is identical in each patch, the temperature of all individuals is identical, regardless of the value of the migration probability. Therefore, the fixation probability in a population with such a structure is the same as the fixation probability in a well-mixed population of the same total population size (N=sum _{j=1}^M N_j), given by ( phi _{mathrm{wm}}^N(r)).Small migration regimeIf the migration probability is small enough such that the time between two subsequent migration events (( sim frac{1}{lambda } )) is much longer than the absorption time within any patch, then at the time of each migration event we may suppose that the meta-population is in a homogeneous configuration22,28. In other words, the low migration regime is an approximation in which we neglect the probability that the meta-population is not in a homogeneous configuration at the time of migration events. We define a homogeneous configuration of the meta-population as a configuration in which in all patches either all individuals are mutants, or all are wild-types.Therefore, instead of having (2^N) states, where N is the population size, the system has only (2^M) states, where M is the number of patches. Thus, we can calculate the fixation probability exactly as in the case of a standard evolutionary graph model where each node represents a single individual but with a modified transition probabilities.In a network with homogeneous patches, in order to increase the number of homogeneous mutant-patches one individual mutant needs to migrate to one of its neighbouring homogeneous wild-type-patches and reaches fixation there. For example if node j is occupied by mutants and one of its neighbouring patches, node k, is occupied by wild-types, the probability that one mutant individual from patch j migrates to patch k and reaches fixation there is (frac{lambda }{mathrm{deg} (j)}phi _{mathrm{wm}}^{N_{k}}(r) ), where (mathrm{deg} (j) ) is the degree of node j to take into account that the mutant can move to different patches. This is analogous to the probability that one mutant in node j replaces one wild-type in node k ,(T^{jrightarrow k}), in the network of individuals.Similarly, if node j is occupied by wild-types and one of its neighbouring patches, node j, is occupied by mutants the probability that one wild-type individual from patch j migrates to patch k and reaches fixation there equals to (frac{lambda }{mathrm{deg} (j)}phi _{mathrm{wm}}^{N_{k}}(1/r) ) where (mathrm{deg} (j) ). Overall, we can move from network of individuals to the network of homogeneous patches by replacing the transition probabilities with the product of migration and fixation probabilities.Two-patch meta-populationThe simplest non-trivial case is the fixation probability in a two-patch meta-population with different local size for small migration probability (lambda ). If the migration probability (lambda ) is very small, a new mutant first needs to take over its own patch and only then the first migrant arrives in the second patch. To be more precise, the time between two migration events has to be much higher than the typical time that it takes for the migrant to take over the patch or go extinct again38. In this case, we can divide the dynamics into two phases: A first phase in which a mutant invades one patch and a second phase in which a homogeneous patch of mutants invades the whole meta-population. Assume a new mutation arises in patch 1. Only if this mutant reaches fixation in patch 1, it also has a chance to reach fixation in patch 2. When patch 1 consists of only mutants and patch 2 consists of only wild-types, there are two possibilities for the ultimate fate of the mutant:

(i)

Eventually, the offspring of one mutant selected from patch 1 for reproduction will migrate to patch 2 and reach fixation there. The wild-type goes extinct. This happens with probability ( frac{N_1 r}{N_1 r+N_2} phi _{mathrm{wm}}^{N_2}(r)).

(ii)

Eventually, the offspring of one wild-type selected from patch 2 for reproduction will migrate to patch 1 and the mutant goes extinct. This occurs with probability ( frac{N_2}{N_1r+N_2} phi _{mathrm{wm}}^{N_1}(tfrac{1}{r})).

Therefore, the probability that a single mutant arising in patch 1 reaches fixation in the entire population is $$begin{aligned} phi _{mathrm{wm}}^{N_1}(r) frac{frac{N_1 r}{N_1 r+N_2} phi _{mathrm{wm}}^{N_2}(r)}{frac{N_1 r}{N_1 r+N_2} phi _{mathrm{wm}}^{N_2}(r)+frac{N_2}{N_1r+N_2} phi _{mathrm{wm}}^{N_1}left( tfrac{1}{r}right) }=phi _{mathrm{wm}}^{N_1}(r) phi _{mathrm{wm}}^{N_2}(r) frac{1 }{ phi _{mathrm{wm}}^{N_2}(r) +frac{N_2}{N_1} frac{1}{r}phi _{mathrm{wm}}^{N_1} left( tfrac{1}{r}right) }. end{aligned}$$
(3a)
Similarly the probability that a mutant arising in patch 2 takes over the whole population equals$$begin{aligned} phi _{mathrm{wm}}^{N_2}(r) phi _{mathrm{wm}}^{N_1}(r) frac{1 }{phi _{mathrm{wm}}^{N_1}(r)+frac{N_1}{N_2} frac{1}{r} phi _{mathrm{wm}}^{N_2}left( tfrac{1}{r}right) }. end{aligned}$$
(3b)
If we assume that the mutant arises in a patch with a probability proportional to the patch size, the average fixation probability (phi _{bullet !!-!!bullet }) in a two patch population for small migration probability is the weighted sum of Eqs. (3a) and (3b),$$begin{aligned} phi _{bullet !!-!!bullet }&= phi _{mathrm{wm}}^{N_1}(r) phi _{mathrm{wm}}^{N_2}(r) nonumber \&quad times left( frac{frac{N_1}{N_1+N_2} }{ phi _{mathrm{wm}}^{N_2}(r) +frac{N_2}{N_1} frac{1}{r}phi _{mathrm{wm}}^{N_1}left( tfrac{1}{r}right) } +frac{frac{N_2}{N_1+N_2} }{ phi _{mathrm{wm}}^{N_1}(r) +frac{N_1}{N_2} frac{1}{r} phi _{mathrm{wm}}^{N_2}left( tfrac{1}{r}right) }right) . end{aligned}$$
(4)
In the case of neutrality, (r=1), we recover (phi _{bullet !!-!!bullet } = frac{1}{N_1+N_2})—the fixation probability in a population of the total size of the two patches. For identical patch sizes, ( N_1=N_2 ), Eq. (4) simplifies to$$begin{aligned} phi _{bullet !!-!!bullet } = left( phi _{mathrm{wm}}^{N_1}(r)right) ^2 frac{1}{phi _{mathrm{wm}}^{N_1}(r)+frac{1}{r} phi _{mathrm{wm}}^{N_1}left( tfrac{1}{r}right) } = phi _{mathrm{wm}}^{2 N_1}(r), end{aligned}$$
(5)
where the simplification to the fixation probability within a single population of size (2N_1) reflects the validity of the isothermal theorem.For (N_1 ne N_2), we approximate Eq. (4) for weak and strong selection. Let us first consider highly advantageous mutants, (r gg 1). In this case, we have (phi _{mathrm{wm}}^{N_1}(r) gg phi _{mathrm{wm}}^{N_1}(tfrac{1}{r})) and thus we can neglect the possibility that a wild-type takes over a mutant patch if patch sizes are sufficiently large. The probability (phi _{bullet !!-!!bullet } ) then becomes a weighted average reflecting patch sizes. For identical patch size (N_1=N_2 = N/2), it reduces to (phi _{bullet !!-!!bullet } approx phi _{mathrm{wm}}^{N_1}(r)=phi _{mathrm{wm}}^{N/2}(r)). In other words, taking over the first patch is sufficient to make fixation in the entire population certain. For patches of very different size, (N_1 gg N_2), we have (N approx N_1) and find (phi _{bullet !!-!! bullet } approx phi _{mathrm{wm}}^{N}(r), ) which implies that fixation is driven by the fixation process in the larger patch, regardless of where the mutant arises. Note that there is a difference between the case of identical patch size and very different patch size . The case of highly disadvantageous mutants, (r ll 1), can be handled in a very similar way.Next, we consider weak selection, (r approx 1). We can approximate the fixation probability as (phi _{mathrm{wm}}^{N}(r^{pm 1}) approx frac{1}{N} pm frac{N-1}{2N} (r-1)). With this, we find$$begin{aligned} phi _{bullet !!-!!bullet } approx frac{1}{N_1+N_2} +frac{1}{2} left( 1 – frac{1}{N_1+N_2} -frac{(N_1-N_2)^2}{(N_1^2+N_2^2)^2} N_1 N_2right) (r-1). end{aligned}$$
(6)
For identical patch size (N_1=N_2 = N/2), this reduces to$$begin{aligned} phi _{bullet !!-!!bullet } approx tfrac{1}{N} +tfrac{N-1}{2N} (r-1), end{aligned}$$
(7)
which is the known result for a single population of size (N=N_1+N_2). When patches have very different size, (N_1 gg N_2) such that (N approx N_1), we recover the same result. Thus, the difference between the fixation probability of a two-patch meta-population with identical patch size and the fixation probability of a two-patch meta-population with very different patch size that we found for highly advantageous mutants is no longer observed for weak selection.When migration probabilities become larger, our approximation is no longer valid and we need to rely on numerical approaches. Figure 2 illustrates the difference between the fixation probability of a two-patch structure meta-population and the equivalent well-mixed population of size (N_1+N_2 ) when migration is low using Eq. (4) and comparing with the numerical approach in Ref.39.While the fixation probability of the two-patch meta-population is very close to the fixation probability of the well-mixed population40, a close inspection reveals an interesting property: For low migration probabilities and (N_1 ne N_2), the two patch structure is a suppressor of selection in the original sense of Lieberman et al.1: For advantageous mutations, (r >1), it decreases the fixation probability, whereas for disadvantageous mutations, (r1) and negative for (r1 ) the minimum fixation probability occurs when the two patch sizes are identical, ( N_1=N_2=N/2 ). Similarly, for fitness values ( r More

1.Liu, M. et al. Evaluation of high-resolution satellite rainfall products using rain gauge data over complex terrain in southwest China. Theoret. Appl. Climatol. 119, 203–219. https://doi.org/10.1007/s00704-014-1092-4 (2014).ADS 
Article 

Google Scholar 
2.Piao, S. et al. Leaf onset in the northern hemisphere triggered by daytime temperature. Nat. Commun. 6, 6911. https://doi.org/10.1038/ncomms7911 (2015).ADS 
CAS 
Article 
PubMed 

Google Scholar 
3.Zheng, Z. et al. Continuous but diverse advancement of spring-summer phenology in response to climate warming across the Qinghai-Tibetan Plateau. Agric. For. Meteorol. 223, 194–202. https://doi.org/10.1016/j.agrformet.2016.04.012 (2016).ADS 
Article 

Google Scholar 
4.Shen, M. et al. Can changes in autumn phenology facilitate earlier green-up date of northern vegetation?. Agric. For. Meteorol. https://doi.org/10.1016/j.agrformet.2020.108077 (2020).Article 

Google Scholar 
5.Zhang, G., Zhang, Y., Dong, J. & Xiao, X. Green-up dates in the Tibetan Plateau have continuously advanced from 1982 to 2011. Proc. Natl. Acad. Sci. U. S. A. 110, 4309–4314. https://doi.org/10.1073/pnas.1210423110 (2013).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
6.Zhu, Z. C. et al. Greening of the Earth and its drivers. Nat. Clim. Chang. 6, 791–795. https://doi.org/10.1038/nclimate3004 (2016).ADS 
CAS 
Article 

Google Scholar 
7.Sun, Q., Li, B., Zhou, G., Jiang, Y. & Yuan, Y. Delayed autumn leaf senescence date prolongs the growing season length of herbaceous plants on the Qinghai-Tibetan Plateau. Agric. For. Meteorol. 284, 1. https://doi.org/10.1016/j.agrformet.2019.107896 (2020).Article 

Google Scholar 
8.Gao, Q. et al. Climatic change controls productivity variation in global grasslands. Sci. Rep. 6, 26958. https://doi.org/10.1038/srep26958 (2016).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
9.Zhang, K. et al. Vegetation greening and climate change promote multidecadal rises of global land evapotranspiration. Sci. Rep. 5, 15956. https://doi.org/10.1038/srep15956 (2015).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
10.Li, Z., Chen, Y., Wang, Y. & Fang, G. Dynamic changes in terrestrial net primary production and their effects on evapotranspiration. Hydrol. Earth Syst. Sci. 20, 2169–2178. https://doi.org/10.5194/hess-20-2169-2016 (2016).ADS 
Article 

Google Scholar 
11.Jeong, S.-J., Ho, C.-H., Gim, H.-J. & Brown, M. E. Phenology shifts at start vs. end of growing season in temperate vegetation over the Northern Hemisphere for the period 1982–2008. Global Change Biol. 17, 2385–2399. https://doi.org/10.1111/j.1365-2486.2011.02397.x (2011).ADS 
Article 

Google Scholar 
12.Wang, Y., Gao, Q., Liu, T., Tian, Y. & Yu, M. The greenness of major shrublands in china increased from 2001 to 2013. Remote Sens. https://doi.org/10.3390/rs8020121 (2016).Article 

Google Scholar 
13.Xu, X. et al. Plant community structure regulates responses of prairie soil respiration to decadal experimental warming. Glob. Change Biol. 21, 3846–3853. https://doi.org/10.1111/gcb.12940 (2015).ADS 
Article 

Google Scholar 
14.Gang, C. et al. Drought-induced dynamics of carbon and water use efficiency of global grasslands from 2000 to 2011. Ecol. Ind. 67, 788–797. https://doi.org/10.1016/j.ecolind.2016.03.049 (2016).CAS 
Article 

Google Scholar 
15.Yao, J., Yang, Q., Mao, W., Zhao, Y. & Xu, X. Precipitation trend–Elevation relationship in arid regions of the China. Glob. Planet. Change 143, 1–9. https://doi.org/10.1016/j.gloplacha.2016.05.007 (2016).ADS 
Article 

Google Scholar 
16.Yuan, X. et al. Vegetation changes and land surface feedbacks drive shifts in local temperatures over Central Asia. Sci. Rep. 7, 3287. https://doi.org/10.1038/s41598-017-03432-2 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
17.Liu, Y., Kumar, M., Katul, G. G., Feng, X. & Konings, A. G. Plant hydraulics accentuates the effect of atmospheric moisture stress on transpiration. Nat. Clim. Chang. 10, 691–695. https://doi.org/10.1038/s41558-020-0781-5 (2020).ADS 
CAS 
Article 

Google Scholar 
18.Li, Y. et al. Evaluation and projection of snowfall changes in High Mountain Asia based on NASA’s NEX-GDDP high-resolution daily downscaled dataset. Environ. Res. Lett. https://doi.org/10.1088/1748-9326/aba926 (2020).Article 

Google Scholar 
19.Yao, T. et al. Chained impacts on modern environment of interaction between Westerlies and Indian Monsoon on Tibetan Plateau. Bull. Chin. Acad. Sci. 32, 976–984. https://doi.org/10.16418/j.issn.1000-3045.2017.09.007 (2017).Article 

Google Scholar 
20.Pritchard, H. D. Asia’s shrinking glaciers protect large populations from drought stress. Nature 569, 649–654. https://doi.org/10.1038/s41586-019-1240-1 (2019).ADS 
CAS 
Article 
PubMed 

Google Scholar 
21.He, Z. et al. Assessing temperature sensitivity of subalpine shrub phenology in semi-arid mountain regions of China. Agric. For. Meteorol. 213, 42–52. https://doi.org/10.1016/j.agrformet.2015.06.013 (2015).ADS 
Article 

Google Scholar 
22.Zhou, J. et al. Alpine vegetation phenology dynamic over 16years and its covariation with climate in a semi-arid region of China. Sci. Total Environ. 572, 119–128. https://doi.org/10.1016/j.scitotenv.2016.07.206 (2016).ADS 
CAS 
Article 
PubMed 

Google Scholar 
23.Zhao, J. et al. Increased precipitation offsets the negative effect of warming on plant biomass and ecosystem respiration in a Tibetan alpine steppe. Agric. For. Meteorol. https://doi.org/10.1016/j.agrformet.2019.107761 (2019).Article 
PubMed 

Google Scholar 
24.Deng, H., Pepin, N. C. & Chen, Y. Changes of snowfall under warming in the Tibetan Plateau. J. Geophys. Res. Atmos. 122, 7323–7341. https://doi.org/10.1002/2017jd026524 (2017).ADS 
Article 

Google Scholar 
25.Yao, T. Tackling on environmental changes in Tibetan Plateau with focus on water, ecosystem and adaptation. Sci. Bull. 64, 1. https://doi.org/10.1016/j.scib.2019.03.033 (2019).Article 

Google Scholar 
26.Shen, M. et al. Strong impacts of daily minimum temperature on the green-up date and summer greenness of the Tibetan Plateau. Glob. Change Biol. 22, 3057–3066. https://doi.org/10.1111/gcb.13301 (2016).ADS 
Article 

Google Scholar 
27.Piao, S. et al. Plant phenology and global climate change: Current progresses and challenges. Glob. Change Biol. 25, 1922–1940. https://doi.org/10.1111/gcb.14619 (2019).ADS 
Article 

Google Scholar 
28.Xu, M. & Xue, X. A research on summer vegetation characteristics & short-time responses to experimental warming of alpine meadow in the Qinghai-Tibetan Plateau. Acta Ecol. Sin. 33, 2071–2083. https://doi.org/10.5846/stxb201112201935 (2013).Article 

Google Scholar 
29.Huang, N., He, J. S., Chen, L. & Wang, L. No upward shift of alpine grassland distribution on the Qinghai-Tibetan Plateau despite rapid climate warming from 2000 to 2014. Sci. Total Environ. 625, 1361–1368. https://doi.org/10.1016/j.scitotenv.2018.01.034 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
30.Piao, S. et al. Evidence for a weakening relationship between interannual temperature variability and northern vegetation activity. Nat. Commun. 5, 5018. https://doi.org/10.1038/ncomms6018 (2014).ADS 
CAS 
Article 
PubMed 

Google Scholar 
31.Huang, M. et al. Air temperature optima of vegetation productivity across global biomes. Nat. Ecol. Evol. 3, 1. https://doi.org/10.1038/s41559-019-0838-x (2019).CAS 
Article 

Google Scholar 
32.Liu, H., Zhang, M., Lin, Z. & Xu, X. Spatial heterogeneity of the relationship between vegetation dynamics and climate change and their driving forces at multiple time scales in Southwest China. Agric. For. Meteorol. 256–257, 10–21. https://doi.org/10.1016/j.agrformet.2018.02.015 (2018).ADS 
Article 

Google Scholar 
33.Chen, Z., Wang, W. & Fu, J. Vegetation response to precipitation anomalies under different climatic and biogeographical conditions in China. Sci. Rep. 10, 830. https://doi.org/10.1038/s41598-020-57910-1 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
34.Guo, H. et al. Space-time characterization of drought events and their impacts on vegetation in Central Asia. J. Hydrol. 564, 1165–1178. https://doi.org/10.1016/j.jhydrol.2018.07.081 (2018).ADS 
Article 

Google Scholar 
35.Li, P., Hu, Z. & Liu, Y. Shift in the trend of browning in Southwestern Tibetan Plateau in the past two decades. Agric. For. Meteorol. https://doi.org/10.1016/j.agrformet.2020.107950 (2020).Article 

Google Scholar 
36.Liu, Z., Li, C., Zhou, P. & Chen, X. A probabilistic assessment of the likelihood of vegetation drought under varying climate conditions across China. Sci. Rep. 6, 35105. https://doi.org/10.1038/srep35105 (2016).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
37.Gao, Q.-Z., Li, Y., Xu, H.-M., Wan, Y.-F. & Jiangcun, W.-Z. Adaptation strategies of climate variability impacts on alpine grassland ecosystems in Tibetan Plateau. Mitig. Adapt. Strat. Glob. Change 19, 199–209. https://doi.org/10.1007/s11027-012-9434-y (2012).CAS 
Article 

Google Scholar 
38.Guo, Y. & Wang, C. Trends in precipitation recycling over the Qinghai-Xizang Plateau in last decades. J. Hydrol. 517, 826–835. https://doi.org/10.1016/j.jhydrol.2014.06.006 (2014).ADS 
Article 

Google Scholar 
39.Schlaepfer, D. R. et al. Climate change reduces extent of temperate drylands and intensifies drought in deep soils. Nat. Commun. 8, 14196. https://doi.org/10.1038/ncomms14196 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
40.Yao, J. et al. Climatic and associated atmospheric water cycle changes over the Xinjiang, China. J. Hydrol. 585, 1. https://doi.org/10.1016/j.jhydrol.2020.124823 (2020).Article 

Google Scholar 
41.Sun, A. et al. Quantified hydrological responses to permafrost degradation in the headwaters of the Yellow River (HWYR) in High Asia. Sci. Total Environ. 712, 135632. https://doi.org/10.1016/j.scitotenv.2019.135632 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
42.Brun, F., Berthier, E., Wagnon, P., Kaab, A. & Treichler, D. A spatially resolved estimate of High Mountain Asia glacier mass balances, 2000–2016. Nat. Geosci. 10, 668–673. https://doi.org/10.1038/NGEO2999 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
43.Luo, D., Liu, L., Jin, H., Wang, X. & Chen, F. Characteristics of ground surface temperature at Chalaping in the Source Area of the Yellow River, northeastern Tibetan Plateau. Agric. For. Meteorol. https://doi.org/10.1016/j.agrformet.2019.107819 (2020).Article 

Google Scholar 
44.Che, M. et al. Spatial and temporal variations in the end date of the vegetation growing season throughout the Qinghai-Tibetan Plateau from 1982 to 2011. Agric. For. Meteorol. 189–190, 81–90. https://doi.org/10.1016/j.agrformet.2014.01.004 (2014).ADS 
Article 

Google Scholar 
45.Ji, Z. et al. Investigation of mineral aerosols radiative effects over High Mountain Asia in 1990–2009 using a regional climate model. Atmos. Res. 178–179, 484–496. https://doi.org/10.1016/j.atmosres.2016.05.003 (2016).CAS 
Article 

Google Scholar 
46.Wang, X. et al. Spring temperature change and its implication in the change of vegetation growth in North America from 1982 to 2006. Proc. Natl. Acad. Sci. U. S. A. 108, 1240–1245. https://doi.org/10.1073/pnas.1014425108 (2011).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
47.Piao, S. et al. Responses and feedback of the Tibetan Plateau’s alpine ecosystem to climate change. Chin. Sci. Bull. 64, 2842–2855. https://doi.org/10.1360/TB-2019-0074 (2019).Article 

Google Scholar 
48.Zeng, Z. et al. Climate mitigation from vegetation biophysical feedbacks during the past three decades. Nat. Clim. Change 7, 432–436. https://doi.org/10.1038/nclimate3299 (2017).ADS 
Article 

Google Scholar 
49.Xu, H. J., Wang, X. P. & Yang, T. B. Trend shifts in satellite-derived vegetation growth in Central Eurasia, 1982–2013. Sci. Total Environ. 579, 1658–1674. https://doi.org/10.1016/j.scitotenv.2016.11.182 (2017).ADS 
CAS 
Article 
PubMed 

Google Scholar 
50.Zhang, Y. et al. Satellite-observed global terrestrial vegetation production in response to water availability. Remote Sens. 13, 1. https://doi.org/10.3390/rs13071289 (2021).Article 

Google Scholar 
51.Curio, J. & Scherer, D. Seasonality and spatial variability of dynamic precipitation controls on the Tibetan Plateau. Earth Syst. Dyn. Discus. https://doi.org/10.5194/esd-2016-1,10.5194/esd-2016-1 (2016).Article 

Google Scholar 
52.Li, J., Sun, C. & Jin, F. F. NAO implicated as a predictor of Northern Hemisphere mean temperature multidecadal variability. Geophys. Res. Lett. 40, 5497–5502. https://doi.org/10.1002/2013gl057877 (2013).ADS 
Article 

Google Scholar 
53.Turner, A. G. & Annamalai, H. Climate change and the South Asian summer monsoon. Nat. Clim. Chang. 2, 587–595. https://doi.org/10.1038/nclimate1495 (2012).ADS 
Article 

Google Scholar 
54.Crimmins, T. M., Crimmins, M. A. & DavidBertelsen, C. Complex responses to climate drivers in onset of spring flowering across a semi-arid elevation gradient. J. Ecol. 98, 1042–1051. https://doi.org/10.1111/j.1365-2745.2010.01696.x (2010).Article 

Google Scholar 
55.Du, J. et al. Interacting effects of temperature and precipitation on climatic sensitivity of spring vegetation green-up in arid mountains of China. Agric. For. Meteorol. 269–270, 71–77. https://doi.org/10.1016/j.agrformet.2019.02.008 (2019).ADS 
Article 

Google Scholar 
56.Huang, J. et al. Global semi-arid climate change over last 60 years. Clim. Dyn. 46, 1131–1150. https://doi.org/10.1007/s00382-015-2636-8 (2015).Article 

Google Scholar 
57.Sun, J., Qin, X. & Yang, J. The response of vegetation dynamics of the different alpine grassland types to temperature and precipitation on the Tibetan Plateau. Environ. Monit. Assess. 188, 20. https://doi.org/10.1007/s10661-015-5014-4 (2016).Article 
PubMed 

Google Scholar 
58.Ganjurjav, H. et al. Differential response of alpine steppe and alpine meadow to climate warming in the central Qinghai-Tibetan Plateau. Agric. For. Meteorol. 223, 233–240. https://doi.org/10.1016/j.agrformet.2016.03.017 (2016).ADS 
Article 

Google Scholar 
59.Xu, M. et al. Year-round warming and autumnal clipping lead to downward transport of root biomass, carbon and total nitrogen in soil of an alpine meadow. Environ. Exp. Bot. 109, 54–62. https://doi.org/10.1016/j.envexpbot.2014.07.012 (2015).CAS 
Article 

Google Scholar 
60.Xie, J. et al. Land surface phenology and greenness in Alpine grasslands driven by seasonal snow and meteorological factors. Sci. Total Environ. 725, 138380. https://doi.org/10.1016/j.scitotenv.2020.138380 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
61.Zhang, Y. et al. Vegetation dynamics and its driving forces from climate change and human activities in the Three-River Source Region, China from 1982 to 2012. Sci. Total Environ. 563–564, 210–220. https://doi.org/10.1016/j.scitotenv.2016.03.223 (2016).ADS 
CAS 
Article 
PubMed 

Google Scholar 
62.Liu, L. et al. Elevation-dependent decline in vegetation greening rate driven by increasing dryness based on three satellite NDVI datasets on the Tibetan Plateau. Ecol. Indic. https://doi.org/10.1016/j.ecolind.2019.105569 (2019).Article 

Google Scholar 
63.Piao, S. et al. Altitude and temperature dependence of change in the spring vegetation green-up date from 1982 to 2006 in the Qinghai-Xizang Plateau. Agric. For. Meteorol. 151, 1599–1608. https://doi.org/10.1016/j.agrformet.2011.06.016 (2011).ADS 
Article 

Google Scholar 
64.Zhang, X., Tarpley, D. & Sullivan, J. T. Diverse responses of vegetation phenology to a warming climate. Geophys. Res. Lett. https://doi.org/10.1029/2007gl031447 (2007).Article 

Google Scholar 
65.Gao, Y. et al. Vegetation net primary productivity and its response to climate change during 2001–2008 in the Tibetan Plateau. Sci. Total Environ. 444, 356–362. https://doi.org/10.1016/j.scitotenv.2012.12.014 (2013).ADS 
CAS 
Article 
PubMed 

Google Scholar 
66.Shen, M. et al. Influences of temperature and precipitation before the growing season on spring phenology in grasslands of the central and eastern Qinghai-Tibetan Plateau. Agric. For. Meteorol. 151, 1711–1722. https://doi.org/10.1016/j.agrformet.2011.07.003 (2011).ADS 
Article 

Google Scholar 
67.Chen, N. et al. The compensation effects of post-drought regrowth on earlier drought loss across the tibetan plateau grasslands. Agric. For. Meteorol. https://doi.org/10.1016/j.agrformet.2019.107822 (2020).Article 

Google Scholar 
68.Zhao, W. et al. Contributions of climatic factors to interannual variability of the vegetation index in Northern China Grasslands. J. Clim. 33, 175–183. https://doi.org/10.1175/jcli-d-18-0587.1 (2020).ADS 
Article 

Google Scholar 
69.Liang, J. et al. Where will threatened migratory birds go under climate change? Implications for China’s national nature reserves. Sci. Total Environ. 645, 1040–1047. https://doi.org/10.1016/j.scitotenv.2018.07.196 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
70.Qu, S. et al. What drives the vegetation restoration in Yangtze River basin, China: Climate change or anthropogenic factors?. Ecol. Ind. 90, 438–450. https://doi.org/10.1016/j.ecolind.2018.03.029 (2018).Article 

Google Scholar 
71.Yin, L. et al. What drives the vegetation dynamics in the Hengduan Mountain region, southwest China: Climate change or human activity?. Ecol. Ind. 112, 106013. https://doi.org/10.1016/j.ecolind.2019.106013 (2020).Article 

Google Scholar 
72.Zhou, X., Yamaguchi, Y. & Arjasakusuma, S. Distinguishing the vegetation dynamics induced by anthropogenic factors using vegetation optical depth and AVHRR NDVI: A cross-border study on the Mongolian Plateau. Sci. Total Environ. 616–617, 730–743. https://doi.org/10.1016/j.scitotenv.2017.10.253 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
73.Li, Y. et al. The effects of fencing on carbon stocks in the degraded alpine grasslands of the Qinghai-Tibetan Plateau. J. Environ. Manag. 128, 393–399. https://doi.org/10.1016/j.jenvman.2013.05.058 (2013).CAS 
Article 

Google Scholar 
74.Liu, X. et al. How does grazing exclusion influence plant productivity and community structure in alpine grasslands of the Qinghai-Tibetan Plateau?. Glob. Ecol. Conserv. https://doi.org/10.1016/j.gecco.2020.e01066 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
75.Li, W. et al. Effects of grazing regime on vegetation structure, productivity, soil quality, carbon and nitrogen storage of alpine meadow on the Qinghai-Tibetan Plateau. Ecol. Eng. 98, 123–133. https://doi.org/10.1016/j.ecoleng.2016.10.026 (2017).Article 

Google Scholar 
76.Deng, L. et al. Effects of grazing exclusion on carbon sequestration in China’s grassland. Earth Sci. Rev. 173, 84–95. https://doi.org/10.1016/j.earscirev.2017.08.008 (2017).ADS 
CAS 
Article 

Google Scholar 
77.Yu, L. et al. Effects of grazing exclusion on soil carbon dynamics in alpine grasslands of the Tibetan Plateau. Geoderma 353, 133–143. https://doi.org/10.1016/j.geoderma.2019.06.036 (2019).ADS 
CAS 
Article 

Google Scholar 
78.Shao, Q. et al. Effects of an ecological conservation and restoration project in the Three-River Source Region, China. J. Geograph. Sci. 27, 183–204. https://doi.org/10.1007/s11442-017-1371-y (2016).Article 

Google Scholar 
79.Sun, Q. et al. A systematic review of research studies on the estimation of net primary productivity in the Three-River Headwater Region, China. J. Geograph. Sci. 27, 161–182. https://doi.org/10.1007/s11442-017-1370-z (2016).Article 

Google Scholar 
80.Shen, X. et al. Marshland loss warms local land surface temperature in China. Geophys. Res. Lett. https://doi.org/10.1029/2020GL087648 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
81.Shen, X. et al. Aboveground biomass and its spatial distribution pattern of herbaceous marsh vegetation in China. Sci. China Earth Sci. 64, 1115–1125. https://doi.org/10.1007/s11430-020-9778-7 (2021).ADS 
Article 

Google Scholar 
82.Wang, Y. et al. Spatiotemporal change of aboveground biomass and its response to climate change in marshes of the Tibetan Plateau. Int. J. Appl. Earth Observ. Geoinf. https://doi.org/10.1016/j.jag.2021.102385 (2021).Article 

Google Scholar 
83.Jeong, S.-J., Ho, C.-H. & Jeong, J.-H. Increase in vegetation greenness and decrease in springtime warming over east Asia. Geophys. Res. Lett. https://doi.org/10.1029/2008gl036583 (2009).Article 

Google Scholar 
84.Shen, M. et al. Evaporative cooling over the Tibetan Plateau induced by vegetation growth. Proc. Natl. Acad. Sci. U. S. A. 112, 9299–9304. https://doi.org/10.1073/pnas.1504418112 (2015).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
85.Shen, X. et al. Asymmetric effects of daytime and nighttime warming on spring phenology in the temperate grasslands of China. Agric. For. Meteorol. 259, 240–249. https://doi.org/10.1016/j.agrformet.2018.05.006 (2018).ADS 
Article 

Google Scholar 
86.Shen, X. et al. Spatiotemporal variation in vegetation spring phenology and its response to climate change in freshwater marshes of Northeast China. Sci. Total Environ. 666, 1169–1177. https://doi.org/10.1016/j.scitotenv.2019.02.265 (2019).ADS 
CAS 
Article 
PubMed 

Google Scholar 
87.Niittynen, P. et al. Fine-scale tundra vegetation patterns are strongly related to winter thermal conditions. Nat. Clim. Chang. https://doi.org/10.1038/s41558-020-00916-4 (2020).Article 

Google Scholar 
88.Wu, D. et al. Evaluation of spatiotemporal variations of global fractional vegetation cover based on GIMMS NDVI data from 1982 to 2011. Remote Sens. 6, 4217–4239. https://doi.org/10.3390/rs6054217 (2014).ADS 
Article 

Google Scholar 
89.Zhang, H. et al. Calculation of evapotranspiration in different climatic zones combining the long-term monitoring data with bootstrap method. Environ. Res. 191, 110200. https://doi.org/10.1016/j.envres.2020.110200 (2020).CAS 
Article 
PubMed 

Google Scholar 
90.Kalisa, W. et al. Assessment of climate impact on vegetation dynamics over East Africa from 1982 to 2015. Sci. Rep. 9, 16865. https://doi.org/10.1038/s41598-019-53150-0 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
91.Chen, Y. Geographical data analysis with Matlab 202–220 (Chen, 2012). More