Ecology

1.Walther, G. R. et al. Ecological responses to recent climate change. Nature 416, 389–395 (2002).ADS 
CAS 
PubMed 
Article 

Google Scholar 
2.Mooij, W. M. et al. The impact of climate change on lakes in the Netherlands: A review. Aquat. Ecol. 39, 381–400 (2005).CAS 
Article 

Google Scholar 
3.Walter, B., Peters, J. & van Beusekom, J. E. E. The effect of constant darkness and short light periods on the survival and physiological fitness of two phytoplankton species and their growth potential after re-illumination. Aquat. Ecol. 51, 591–603 (2017).CAS 
Article 

Google Scholar 
4.Woodward, G., Perkins, D. M. & Brown, L. E. Climate change and freshwater ecosystems: Impacts across multiple levels of organization. Philos. Trans. R. Soc. B Biol. Sci. 365, 2093–2106 (2010).Article 

Google Scholar 
5.Wagner, H., Fanesi, A. & Wilhelm, C. Title: Freshwater phytoplankton responses to global warming. J. Plant Physiol. 203, 127–134 (2016).CAS 
PubMed 
Article 

Google Scholar 
6.Gilbert, J. A. Some phytoplankton like it hot. Nat. Clim. Change 3, 954–955 (2013).ADS 
Article 

Google Scholar 
7.Hense, I., Meier, H. E. M. & Sonntag, S. Projected climate change impact on Baltic Sea cyanobacteria: Climate change impact on cyanobacteria. Clim. Change 119, 391–406 (2013).CAS 
Article 

Google Scholar 
8.Trombetta, T. et al. Water temperature drives phytoplankton blooms in coastal waters. PLoS One 14, e0214933 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
9.Jin, P. & Agustí, S. Fast adaptation of tropical diatoms to increased warming with trade-offs. Sci. Rep. 8, 17771 (2018).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
10.Pinceel, T., Buschke, F., Weckx, M., Brendonck, L. & Vanschoenwinkel, B. Climate change jeopardizes the persistence of freshwater zooplankton by reducing both habitat suitability and demographic resilience. BMC Ecol. 18, 1–9 (2018).Article 

Google Scholar 
11.Shin, H. R. & Kneitel, J. M. Warming interacts with inundation timing to influence the species composition of California vernal pool communities. Hydrobiologia 843, 93–105 (2019).Article 

Google Scholar 
12.Montrone, A. et al. Climate change impacts on vernal pool hydrology and vegetation in northern California. J. Hydrol. 574, 1003–1013 (2019).ADS 
Article 

Google Scholar 
13.Williams, D. D. The biology of temporary waters. Biol. Tempor. Waters https://doi.org/10.1093/acprof:oso/9780198528128.001.0001 (2007).Article 

Google Scholar 
14.Waterkeyn, A., Grillas, P., Vanschoenwinkel, B. & Brendonck, L. Invertebrate community patterns in Mediterranean temporary wetlands along hydroperiod and salinity gradients. Freshw. Biol. 53, 1808–1822 (2008).CAS 
Article 

Google Scholar 
15.Lemmens, P. et al. How to maximally support local and regional biodiversity in applied conservation? Insights from pond management. PLoS One 8, e72538 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Lischeid, G. et al. Natural ponds in an agricultural landscape: External drivers, internal processes, and the role of the terrestrial-aquatic interface. Limnologica 68, 5–16 (2018).CAS 
Article 

Google Scholar 
17.Mancinelli, G., Mali, S. & Belmonte, G. Species richness and taxonomic distinctness of zooplankton in ponds and small lakes from Albania and North Macedonia: The role of bioclimatic factors. Water (Switzerland) 11, 2384 (2019).
Google Scholar 
18.Gołdyn, B., Kowalczewska-Madura, K. & Celewicz-Gołdyn, S. Drought and deluge: Influence of environmental factors on water quality of kettle holes in two subsequent years with different precipitation. Limnologica 54, 14–22 (2015).Article 
CAS 

Google Scholar 
19.Salmaso, N. & Tolotti, M. Phytoplankton and anthropogenic changes in pelagic environments. Hydrobiologia https://doi.org/10.1007/s10750-020-04323-w (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
20.Celewicz, S., Czyż, M. J. & Gołdy, B. Feeding patterns in Eubranchipus grubii (Dybowski 1860) (Branchiopoda: Anostraca) and its potential influence on the phytoplankton communities of vernal pools. J. Limnol. 77, 276–284 (2018).Article 

Google Scholar 
21.Rasconi, S., Winter, K. & Kainz, M. J. Temperature increase and fluctuation induce phytoplankton biodiversity loss—Evidence from a multi-seasonal mesocosm experiment. Ecol. Evol. 7, 2936–2946 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
22.Celewicz-Goldyn, S. & Kuczynska-Kippen, N. Ecological value of macrophyte cover in creating habitat for microalgae (diatoms) and zooplankton (rotifers and crustaceans) in small field and forest water bodies. PLoS One 12, e0177317 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
23.Kozak, A., Celewicz-Gołdyn, S. & Kuczyńska-Kippen, N. Cyanobacteria in small water bodies: The effect of habitat and catchment area conditions. Sci. Total Environ. 646, 1578–1587 (2019).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
24.Iacarella, J. C., Barrow, J. L., Giani, A., Beisner, B. E. & Gregory-Eaves, I. Shifts in algal dominance in freshwater experimental ponds across differing levels of macrophytes and nutrients. Ecosphere 9, e02086 (2018).Article 

Google Scholar 
25.Toseland, A. et al. The impact of temperature on marine phytoplankton resource allocation and metabolism. Nat. Clim. Change 3, 979–984 (2013).ADS 
CAS 
Article 

Google Scholar 
26.Richardson, J. et al. Response of cyanobacteria and phytoplankton abundance to warming, extreme rainfall events and nutrient enrichment. Glob. Change Biol. 25, 3365–3380 (2019).ADS 
Article 

Google Scholar 
27.De Senerpont Domis, L. N., Mooij, W. M. & Huisman, J. Climate-induced shifts in an experimental phytoplankton community: A mechanistic approach. Hydrobiologia 584, 403–413 (2007).Article 

Google Scholar 
28.Boyce, D. G., Lewis, M. R. & Worm, B. Global phytoplankton decline over the past century. Nature 466, 591–596 (2010).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
29.Hinder, S. L. et al. Changes in marine dinoflagellate and diatom abundance under climate change. Nat. Clim. Change 2, 271–275 (2012).ADS 
Article 

Google Scholar 
30.Winder, M. & Sommer, U. Phytoplankton response to a changing climate. Hydrobiologia 698, 5–16 (2012).Article 

Google Scholar 
31.Machado, K. B., Vieira, L. C. G. & Nabout, J. C. Predicting the dynamics of taxonomic and functional phytoplankton compositions in different global warming scenarios. Hydrobiologia 830, 115–134 (2019).CAS 
Article 

Google Scholar 
32.O’Neil, J. M., Davis, T. W., Burford, M. A. & Gobler, C. J. The rise of harmful cyanobacteria blooms: The potential roles of eutrophication and climate change. Harmful Algae 14, 313–334 (2012).Article 
CAS 

Google Scholar 
33.Rasconi, S., Gall, A., Winter, K. & Kainz, M. J. Increasing water temperature triggers dominance of small freshwater plankton. PLoS One 10, e0140449 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
34.Wirth, C., Limberger, R. & Weisse, T. Temperature × light interaction and tolerance of high water temperature in the planktonic freshwater flagellates Cryptomonas (Cryptophyceae) and Dinobryon (Chrysophyceae). J. Phycol. 55, 404–414 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Wang, H. et al. High antioxidant capability interacts with respiration to mediate two Alexandrium species growth exploitation of photoperiods and light intensities. Harmful Algae 82, 26–34 (2019).CAS 
PubMed 
Article 

Google Scholar 
36.Fakhri, M., Arifin, N. B., Budianto, B., Yuniarti, A. & Hariati, A. M. Effect of salinity and photoperiod on growth of microalgae Nannochloropsis sp. and Tetraselmis sp. Nat. Environ. Pollut. Technol. 14, 563–566 (2015).
Google Scholar 
37.Torzillo, G., Sacchi, A. & Materassi, R. Temperature as an important factor affecting productivity and night biomass loss in Spirulina platensis grown outdoors in tubular photobioreactors. Bioresour. Technol. 38, 95–100 (1991).Article 

Google Scholar 
38.Shatwell, T., Köhler, J. & Nicklisch, A. Temperature and photoperiod interactions with phosphorus-limited growth and competition of two diatoms. PLoS One 9, e102367 (2014).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
39.Li, G., Talmy, D. & Campbell, D. A. Diatom growth responses to photoperiod and light are predictable from diel reductant generation. J. Phycol. 53, 95–107 (2017).CAS 
PubMed 
Article 

Google Scholar 
40.Reynolds, C. S. Vegetation Processes in the Pelagic: A Model for Ecosystem Theory. Excellence in Ecology Vol. 77 (Ecology Institute, 1997).
Google Scholar 
41.Elliott, J. A., Jones, I. D. & Thackeray, S. J. Testing the sensitivity of phytoplankton communities to changes in water temperature and nutrient load, in a temperate lake. Hydrobiologia 559, 401–411 (2006).CAS 
Article 

Google Scholar 
42.Jöhnk, K. D. et al. Summer heatwaves promote blooms of harmful cyanobacteria. Glob. Change Biol. 14, 495–512 (2008).ADS 
Article 

Google Scholar 
43.Elliott, J. A. Is the future blue–green? A review of the current model predictions of how climate change could affect pelagic freshwater cyanobacteria. Water Res. 46, 1364–1371 (2012).CAS 
PubMed 
Article 

Google Scholar 
44.Ullah, H., Nagelkerken, I., Goldenberg, S. U. & Fordham, D. A. Climate change could drive marine food web collapse through altered trophic flows and cyanobacterial proliferation. PLoS Biol. 16, e2003446 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
45.Hansson, L. A. et al. Food-chain length alters community responses to global change in aquatic systems. Nat. Clim. Change 3, 228–233 (2013).ADS 
Article 

Google Scholar 
46.Burgmer, T. & Hillebrand, H. Temperature mean and variance alter phytoplankton biomass and biodiversity in a long-term microcosm experiment. Oikos 120, 922–933 (2011).Article 

Google Scholar 
47.Hillebrand, H., Burgmer, T. & Biermann, E. Running to stand still: Temperature effects on species richness, species turnover, and functional community dynamics. Mar. Biol. 159, 2415–2422 (2012).Article 

Google Scholar 
48.Lewandowska, A. M. et al. Responses of primary productivity to increased temperature and phytoplankton diversity. J. Sea Res. 72, 87–93 (2012).ADS 
Article 

Google Scholar 
49.Lewandowska, A. M., Hillebrand, H., Lengfellner, K. & Sommer, U. Temperature effects on phytoplankton diversity—The zooplankton link. J. Sea Res. 85, 359–364 (2014).ADS 
Article 

Google Scholar 
50.Bergkemper, V., Stadler, P. & Weisse, T. Moderate weather extremes alter phytoplankton diversity—A microcosm study. Freshw. Biol. 63, 1211–1224 (2018).CAS 
Article 

Google Scholar 
51.McMinn, A. & Martin, A. Dark survival in a warming world. Proc. R. Soc. B Biol. Sci. 280, 20122909 (2013).CAS 
Article 

Google Scholar 
52.Waibel, A., Peter, H. & Sommaruga, R. Importance of mixotrophic flagellates during the ice-free season in lakes located along an elevational gradient. Aquat. Sci. 81, 1–10 (2019).CAS 
Article 

Google Scholar 
53.Chen, B. Patterns of thermal limits of phytoplankton. J. Plankton Res. 37, 285–292 (2015).Article 

Google Scholar 
54.Reeves, S., McMinn, A. & Martin, A. The effect of prolonged darkness on the growth, recovery and survival of Antarctic sea ice diatoms. Polar Biol. 34, 1019–1032 (2011).Article 

Google Scholar 
55.van de Poll, W. H., Abdullah, E., Visser, R. J. W., Fischer, P. & Buma, A. G. J. Taxon-specific dark survival of diatoms and flagellates affects Arctic phytoplankton composition during the polar night and early spring. Limnol. Oceanogr. 65, 903–914 (2020).ADS 
Article 

Google Scholar 
56.Poniewozik, M. & Juráň, J. Extremely high diversity of euglenophytes in a small pond in eastern Poland. Plant Ecol. Evol. 151, 18–34 (2018).Article 

Google Scholar 
57.Shafik, H. M., Herodek, S., Présing, M. & Vörös, L. Factors effecting growth and cell composition of cyanoprokaryote Cylindrospermopsis raciborskii (Wołoszyńska) Seenayya et Subba Raju. Algol. Stud. Hydrobiol. Suppl. 103, 75–93 (2001).
Google Scholar 
58.Tang, E. P. Y. & Vincent, W. F. Effects of daylength and temperature on the growth and photosynthesis of an Arctic cyanobacterium, Schizothrix calcicola (Oscillatoriaceae). Eur. J. Phycol. 35, 263–272 (2000).Article 

Google Scholar 
59.Agasild, H., Zingel, P., Tõnno, I., Haberman, J. & Nõges, T. Contribution of different zooplankton groups in grazing on phytoplankton in shallow eutrophic Lake Võrtsjärv (Estonia). Hydrobiologia 584, 167–177 (2007).Article 

Google Scholar 
60.Gołdyn, R. & Kowalczewska-Madura, K. Interactions between phytoplankton and zooplankton in the hypertrophic Swarzȩdzkie Lake in western Poland. J. Plankton Res. 30, 33–42 (2008).Article 
CAS 

Google Scholar 
61.Tovar-Sanchez, A., Duarte, C. M., Hernández-León, S. & Sañudo-Wilhelmy, S. A. Krill as a central node for iron cycling in the Southern Ocean. Geophys. Res. Lett. 34, L11601 (2007).ADS 
Article 
CAS 

Google Scholar 
62.Hunt, R. J. & Matveev, V. F. The effects of nutrients and zooplankton community structure on phytoplankton growth in a subtropical Australian reservoir: An enclosure study. Limnologica 35, 90–101 (2005).Article 

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

Google Scholar 
64.Gołdyn, B., Chudzińska, M., Barałkiewicz, D. & Celewicz-Gołdyn, S. Heavy metal contents in the sediments of astatic ponds: Influence of geomorphology, hydroperiod, water chemistry and vegetation. Ecotoxicol. Environ. Saf. 118, 103–111 (2015).PubMed 
Article 
CAS 

Google Scholar 
65.IPCC. Climate Change 2007: The Physical Science Basis (Cambridge University Press, 2007).
Google Scholar 
66.Christensen, J. H. & Christensen, O. B. A summary of the PRUDENCE model projections of changes in European climate by the end of this century. Clim. Change 81, 7–30 (2007).ADS 
Article 

Google Scholar 
67.Beniston, M. et al. Future extreme events in European climate: An exploration of regional climate model projections. Clim. Change 81, 71–95 (2007).Article 

Google Scholar 
68.Anderson, M. J. A new method for non-parametric multivariate analysis of variance. Austral Ecol. 26, 32–46 (2001).
Google Scholar 
69.Arbizu, P. M. pairwiseAdonis: Pairwise Multilevel Comparison Using Adonis. R Packag. version 0.0.1. (2017).70.Rink, B. & Raak, C. J. F. Principal response curves: Analysis of time-dependent multivariate responses of biological community to stress. Environ. Toxicol. Chem. 18, 138–148 (1999).Article 

Google Scholar 
71.Lepš, J. & Šmilauer, P. Multivariate Analysis of Ecological Data using CANOCO. Bulletin of the Ecological Society of America Vol. 87 (Cambridge University Press, 2003).MATH 
Book 

Google Scholar 
72.Jongman, R. H. G., Ter Braak, C. J. F. & van Tongeren, O. F. R. Data Analysis in Community and Landscape Ecology. Data Analysis in Community and Landscape Ecology (Cambridge University Press, 1995). https://doi.org/10.1017/cbo9780511525575.Book 

Google Scholar 
73.ter Braak, J. F. C. & Šmilauer, P. Canoco Reference Manual and CanoDraw for Windows User’s Guide (Microcomputer Power, 2002).
Google Scholar 
74.R Development Core Team. R: A Language and Environment for Statistical Computing (2020).75.Oksanen, J. et al. vegan: Community Ecology Package. R Packag. version 2.5-7 (2020). More

Insect preparationWe collected 200–300 larvae and pupae of the diamondback moth from cabbage and cauliflower fields in Wuhan, Beijing, and Shenyang in late September from 2008 to 2012. We mixed all collections into one stock colony because of no geographic differentiation in this species from these sites59. We reared all individuals on cabbage leaves spread evenly across five screen cages (35 × 35 × 15 cm) in growth chambers at constant temperature (25 ± 1 °C) with 15-h light:9-h dark photoperiod, and relative humidity set at 60 ± 10%. We moved any new pupae to new screen cages for adult emergence. Emerged adults were fed with 10% honey solution in cotton balls. To collect eggs, we dipped four small pieces of laboratory film (7 × 5 cm) in fresh cabbage juice for 3–5 s and hung the treated film pieces on the top of each screen cage. To further enlarge the population for our experiments, we reared these insects in the artificial diet in plastic boxes at 25 ± 1 °C. We transferred 200 eggs to the surface of 120 g artificial diet (Southland Products Incorporated, USA) in each plastic box (10 × 10 × 9 cm). The hatched larvae dropped to the surface of the artificial diet and fed on it. Once individuals developed into 3rd or 4th instar larvae, pupae or adults, they were exposed to 10 °C for 24 h (to simulate gradually reduced temperatures in late autumn and allow a thermal acclimation) just before they were placed to the low-temperature regimes for overwintering tests. Overall, we obtained >7000 larvae and >8000 pupae for the laboratory experiment, and >8000 larvae, >8000 pupae, >4000 adults for the field experiment. We have compared the life history traits of insects reared on the artificial diet with natural host plants (cabbage leaves), they performed similarly (Peng and Li, unpublished data).Laboratory experiment of winter survivalSite selectionTo identify what factor determines winter survival under different winter thermal conditions, we conducted a laboratory experiment that simulated temperature regimes of 10 selected sites across a latitudinal gradient in China (Fig. 1, Supplementary Table 1) at which this species is known to occur and damage cruciferous crops during the growing season.Temperature treatmentTo simulate the winter temperatures in the 10 geographically distinct sites (Fig. 1a, Supplementary Table 1), we collected daily mean temperatures during winter (November to next April of 1966–2010) at each site from China Meteorological Data Service Centre (http://data.cma.cn/en). Then, we fitted a polynomial model to the temporal changes of winter daily mean temperatures for each site (Fig. 1b). To simplify the logistics of temperature control procedures, we set all temperature regimes in combinations of linear decline, horizontal maintenance and linear increase to mimic the polynomial changes of winter temperatures in the 10 sites, and adjusted temperature every 10 days as needed (Fig. 1c). We controlled the winter temperature changes of the 10 sites with climate chambers (RXZ-280B, Jiangnan Ltd., Ningbo, China) and refrigerators (Royalstar BCD-246GER) according to curves in Fig. 1c.Experimental protocolsWe conducted a winter survival experiment with 10 low-temperature regimes (Fig. 1c). We exposed 6050 larvae and 7150 pupae to 10 low-temperature regimes according to the experiment design (Fig. 1c). Then we sampled 55 larvae and 65 pupae every 10 days from each temperature regime resulting in 11 sampling points. Sampled larvae were placed at 25 °C for 1.5 days to observe the survival based on if their body kept fresh green60 and appendage moved after touching with a brush34. The pupae were placed at 25 °C, RH 70–80% and photoperiod of 16 L:8D for emergence to determine the survival (emergence rate). These samples were not returned to the temperature treatments. Thus, no individual was measured more than once and each sample interval represents an independent observation.Field survival experiments across 12 geographic sitesTo verify the cold survivals from the laboratory simulation and identify the best predictor under natural conditions, we conducted field experiments to explore the winter survival for multiple years at various geographic sites in China (Fig. 1a, Supplementary Table 1). The diamondback moth overwinters either in remaining cabbage plants or in fallen leaves (post-harvest conditions) in regions without standing cabbage crops in the winter. We tested the winter survival of larvae, pupae and adults in the caged cabbage plants or post-harvest conditions in fallen leaves on the soil surface at each site for 3–4 months. We transferred 30 larvae, 30 pupae or 30 adults from our stock rearing to a cabbage plant in the field. Then each plant was covered with a screen cage to avoid disturbance and contain focal individuals (see Supplementary Fig. 2). We set 6–8 cages for larvae, pupae and adults, respectively, in a field in November or early December. After an exposure of 1, 2, 3 and 4 months, we collected 2 cages of larvae, 2 cages of pupae and 2 cages of adults from the field at each sampling point and kept individuals in the laboratory (25 ± 1 °C, RH 65–75%, L:D = 16:8 h) for two days. We checked the survival status of the larvae based on the change in body coloration (i.e. if the larval body kept fresh green colour)60, pupal survival based on whether adults could emerge from the pupae, and adult survival based on if their appendage moved after touching with a brush.To simulate the field microenvironment of post-harvest conditions in winter, we filled half of a glass jar (diameter = 5.5 cm, height = 14 cm) with moist soil. Then, we transferred 30 larvae or 30 pupae to the soil surface, covered the insects with leaves, and then covered the glass jar with a nylon net (see Supplementary Fig. 2). We buried 6–8 jars for larvae and pupae, respectively, and kept the top of the jar at ground surface level at each site in November and early December. Because almost all adults died in few days within the jar, we did not test the adult survival in post-harvest conditions. After an exposure of 1, 2, 3 and 4 months, we took 2 jars of larvae and 2 jars of pupae per sampling period from the field and placed them in the laboratory with 25 ± 1 °C, RH 65–75%, L:D = 16:8 h for 2 days. The survival status of the larvae and pupae was checked with the same procedures as the overwintering tests on caged cabbage plants. Note that as in the standing plant experiment, no individual was tested more than once assuming that each observation is independent at the replicate level.Modelling and predicting winter survivalModel developmentOur goal was to identify key metrics that best predict the winter survival of the diamondback moth across a climatic gradient. To achieve this goal, we took several steps. First, we fit a set of predictive models to the laboratory experiments to identify which metric and model best describes survival under controlled conditions. We focused on three alternative predictors: the lowest daily mean temperature (MinDTmean), mean temperature (DTmean) combined with exposure days, and low-temperature degree-days (LTDD). The MinDTmean model assumes that survival can simply be predicted as a function of the lowest temperature an individual experienced during its exposure time. The DTmean model assumes that survival depends on both the average temperature individuals experience below the cold threshold for survival (11 °C)32 and exposure duration (note that exposure time varied systematically in 10-day increments). Finally, the LTDD model predicts survival depending on coldness below the cold threshold. We calculated LTDD by summing up negative deviations of daily mean temperatures from the cold threshold (11.0 °C) during each exposure period for each simulated geographic site (Fig. 1c). To detect potential relationships, we fit each model using three different functions, i.e. linear, exponential and sigmoid models to describe the survival probability (Supplementary Table 2, Fig. 1). We estimated parameters of models in SigmaStat 3.5 and compared model fit using R2 and AIC values (see detailed models in Supplementary Table 2).Field validation of survival modelsTo validate winter survival models derived from the laboratory (see models in Supplementary Table 2) for complex and variable field conditions, we compared model predictions to observed survivals in field experiments across 12 different geographic sites over 5 years (Fig. 1a, Supplementary Table 1). To make the connection, we first collected daily mean temperatures recorded at the nearest weather stations to our field sites from China Meteorological Data Service Centre. We then calculated MinDTmean, DTmean and exposure days, and LTDD for each site for each treated period and input these values into these laboratory models to predict winter survival. Note, that because the coefficients were calculated from the laboratory experiment, predictions are completely independent of survival observed under field conditions. During the model validation, we excluded the field data of south China, e.g. Guangzhou, Changsha and Wuhan where the warmer temperatures allowed moths to continue their regular life cycle during the whole winter, resulting in unrealistic winter survival. We also excluded replicates in which glass jars were filled with water and destroyed the tested insects. We used linear regression to compare predicted survival with field observations. The validity of each model was evaluated based on the variance explained, slopes of linear regressions and prediction bias (i.e. deviation from unity slope). Finally, we selected the exponential model driven by LTDD as the model to predict the global distribution of winter survival due to its lowest AIC value (Supplementary Table 2) and the least bias (Supplementary Table 3, Fig. 2) among all models.Global prediction of overwintering range shiftTo extrapolate our winter survival predictions to a global scale under present and future climate conditions, we downloaded global historical daily mean temperature data for 50 years (1967–2016) from Berkeley Earth (1° × 1° grid, http://berkeleyearth.org/data/). We added 1, 2, 3, 4, 5 and 6 °C to mean temperatures of 2012–2016, respectively, to represent the different future warming scenarios37. Then, we calculated the annual LTDD in the northern hemisphere with Eq. (1) and in the southern hemisphere with Eq. (2). For xi,j  x0, we excluded the xi,j for the calculation LTDD. For Eq. (1), we started the calculation of LTDD from July 1st (Julian date 182), ended on June 30th of next year (Julian date 181) to cover the whole low-temperature season in the northern hemisphere cross the calendar year. We used LTDD for every year during past conditions to our validated survival model (LTDD-dependent exponential model) and further calculated the expected corresponding yearly winter survival and 5-year mean survival. Since the diamondback moth only feeds on Brassicaceae plants61, we incorporated host availability to refine the pest distributions. We retrieved Brassicaceae occurrence data during 1967–2016 (3,720,971 records) from the Global Biodiversity Information Facility (GBIF) database (www.gbif.org), and excluded unknown and duplicate records; 919,808 records were retained to model the global distribution of host plants. We used a dataset of eight selected bioclimatic variables as described in a previous Brassicaceae biogeographic study62, including isothermality (bio3), temperature seasonality (bio4), min temperature of coldest month (bio6), mean temperature of wettest quarter (bio8), mean temperature of driest quarter (bio9), precipitation seasonality (bio15), precipitation of warmest quarter (bio18), precipitation of coldest quarter (bio19) from Worldclim dataset63 (http://worldclim.org). We ran the species distribution model using the Maxent algorithm in R package dismo64. Model outputs were presented in grid ranks of host plant presence probability from 0 (unsuitable) to 1 (most suitable). Based on the known distribution of Brassicaceae, we only included grid cells with Brassicaceae presence probability ≥0.3 for our final survival and distribution analysis to ensure the presence of the host plant and mapped them with Arcmap 10.2 (Environmental Systems Research Institute) (see Fig. 3a, e). To show spatial-temporal changes in the geographic distribution of winter survival, we quantified the historical change (expansions or contractions) in the overwintering range based on the total numbers of grids for each year between 1967 and 2016 relative to the baseline area in 1967 (see Fig. 3f) and further calculated average changes of every 5 years (see Fig. 3b). We selected sites in the overwintering marginal belt (with winter survival between 1 and 5%) in the baseline year (1967), calculated the annual LTDD of these sites from 1967 to 2016, and built the linear trend of annual LTDD for years 1967–2016 (see Fig. 3g). We predicted distribution changes for future scenarios (added 1, 2, 3, 4, 5 and 6 °C to the current mean temperatures of 2012–2016) relative to the baseline area of 1967–1971 (see Fig. 3c, d).Meta-analysis linking pesticide resistance to overwintering typeData preparation: literature search and selection criteriaWe performed a comprehensive literature survey to collect data on pesticide resistance of the diamondback moth worldwide. We searched for publications in databases of ISI Web of Science, Scopus and China National Knowledge Infrastructure (CNKI) using keywords “pesticide resistance” in combination with “diamondback moth” or “Plutella xylostella” and expanded references in the selected papers. We reviewed titles, abstracts and in many cases the full articles for relevance and agreement with our inclusion criteria. Studies were included if they (1) monitored the pesticide resistance of field populations, (2) used the leaf dip bioassay method to test pesticide resistance which is the most commonly used method recommended by Insecticide Resistance Action Committee (IRAC, http://www.irac-online.org); (3) provided resistance ratio of field populations. Resistance ratio (abbreviated as RR) is the magnitude of pesticide resistance and is commonly calculated by dividing the median lethal concentration (LC50) of a tested field population by LC50 of the susceptible population (without exposure to pesticide). The LC50 is commonly estimated from a concentration-mortality curve of a given pesticide. The preliminary literature search resulted in 2151 studies out of which 62 matched these criteria. A PRISMA diagram describing details of our literature search is available in Supplementary Fig. 4.Data preparation: data extractionWe extracted data from each selected publication, the names of pesticides, sampling locations and years of field populations, number of tested individuals in a bioassay, resistance ratio of field populations (RR), LC50 of field populations (LC50field) and susceptible populations (LC50susceptible), and 95% confidence intervals (CIs) of LC50field and LC50susceptible. Some studies generated results from multiple types of pesticides with the same field population, each of which was considered as a different entry. Finally, we gathered 1806 entries for pesticide resistance of field populations of the diamondback moth.Data preparation: calculation of the weighted effect sizeWe conducted a meta-analysis to test if pesticide resistance levels vary across different types of overwintering sites. To account for differences in sample sizes and variances in resistance ratios across studies, we calculated the corrected (weighted) resistance ratio for each study following the method in Hedges et al.65. We calculated the logarithm of resistance ratio (logRR) to present the effect size for each entry and further calculated the weighted effect size (wlogRR) by$${{{{{rm{wlogRR}}}}}}={{{{{rm{logRR}}}}}}times w$$
(3)
where w is the weighting factor of each entry, with w = 1/sqrt(VlogRR)66. To consider the contribution from both field and susceptible population, the pooled variance VlogRR was calculated as follows65:$${{{{{{rm{V}}}}}}}_{{{{{{rm{logRR}}}}}}}=frac{{{{{{{{rm{SE}}}}}}}_{{{{{{rm{field}}}}}}}}^{2}}{{n}_{{{{{{rm{field}}}}}}}times {{{{{{{rm{LC50}}}}}}}_{{{{{{rm{field}}}}}}}}^{2}}+frac{{{{{{{{rm{SE}}}}}}}_{{{{{{rm{susceptible}}}}}}}}^{2}}{{n}_{{{{{{rm{susceptible}}}}}}}times {{{{{{{rm{LC50}}}}}}}_{{{{{{rm{susceptible}}}}}}}}^{2}}$$
(4)
where LC50field and LC50susceptible, SEfield and SEsusceptible, nfield and nsusceptible, are LC50, the standard error of LC50 and sample size for field population and susceptible population, respectively. SEfield and SEsusceptible can be calculated from their own confidence intervals (95% CI)67:$${{{{{rm{SE}}}}}}=frac{{{{{{{rm{CI}}}}}}}_{{{{{{rm{upper}}}}}}{{{{{rm{limit}}}}}}}-{{{{{{rm{CI}}}}}}}_{{{{{{rm{lower}}}}}}{{{{{rm{limit}}}}}}}}{2times 1.96}$$
(5)
where CIupper limit is the upper limit and CIlower limit is the lower limit of the 95% CI for LC50.We used the prognostic method68 to estimate VlogRR for entries that miss either 95% CI or LC50 based on the average VlogRR of the other complete entries.Data preparation: potential moderator variablesSeveral factors could influence pesticide resistance besides overwintering temperatures. The effective temperature degree-days (ETDD) may change the annual number of generations, the intensity of pesticide application, and thus the selection stress47, e.g. between 7.4 and 33 °C for the diamondback moth32. In addition, the variety of pesticides used in a study may also affect the resistance levels through their mode of actions (the lethal mechanism) and cross-resistance69,70. To account for these potentially confounding factors, we collected the mode of action for each variety of pesticides from IRAC, and calculated LTDD, ETDD and overwintering type for each of the 1806 original records. We collected data for daily mean temperatures for each site from Berkeley Earth. For each location, we calculated the mean annual LTDD, ETDD and winter survival average across the 5 years before the sample. We split the sampling sites into three types based on predicted winter survival of the diamondback moth: (1) the permanent (overwintering) sites: locations with the mean winter survivals ≥5%, (2) marginal sites: locations with the mean winter survivals 1–5%, (3) transient (non-overwintering) sites: locations with the mean winter survivals More

1.Elmqvist, T. et al. Benefits of restoring ecosystem services in urban areas. Curr. Opin. Environ. Sustain. 14, 101–108 (2015).Article 

Google Scholar 
2.Jones, H. P. et al. Restoration and repair of Earth’s damaged ecosystems. Proc. R. Soc. B Biol. Sci. https://doi.org/10.1098/rspb.2017.2577 (2018).Article 

Google Scholar 
3.Rey Benayas, J. M., Newton, A. C., Diaz, A. & Bullock, J. M. Enhancement of biodiversity and ecosystem services by ecological restoration: A meta-analysis. Science 325, 1121–1124. https://doi.org/10.1126/science.1172460 (2009).ADS 
CAS 
Article 
PubMed 

Google Scholar 
4.Brudvig, L. A. et al. Interpreting variation to advance predictive restoration science. J. Appl. Ecol. 54, 1018–1027. https://doi.org/10.1111/1365-2664.12938 (2017).Article 

Google Scholar 
5.Suding, K. N. Toward an era of restoration in ecology: Successes, failures, and opportunities ahead. Annu. Rev. Ecol. Evol. Syst. 42, 465–487. https://doi.org/10.1146/annurev-ecolsys-102710-145115 (2011).Article 

Google Scholar 
6.Reynolds, H. L., Packer, A., Bever, J. D. & Clay, K. Grassroots ecology: Plant-microbe-soil interactions as drivers of plant community structure and dynamics. Ecology 84, 2281–2291 (2003).Article 

Google Scholar 
7.Van Der Heijden, M. G. A., Bardgett, R. D. & Van Straalen, N. M. The unseen majority: Soil microbes as drivers of plant diversity and productivity in terrestrial ecosystems. Ecol. Lett. 11, 296–310 (2008).Article 

Google Scholar 
8.Hoeksema, J. D. et al. A meta-analysis of context-dependency in plant response to inoculation with mycorrhizal fungi. Ecol. Lett. 13, 394–407. https://doi.org/10.1111/j.1461-0248.2009.01430.x (2010).Article 
PubMed 

Google Scholar 
9.Schultz, P. A. et al. Evidence of a mycorrhizal mechanism for the adaptation of Andropogon gerardii (Poaceae) to high- and low-nutrient prairies. Am. J. Bot. 88, 1650–1656. https://doi.org/10.2307/3558410 (2001).CAS 
Article 
PubMed 

Google Scholar 
10.Koske, R. E., & Gemma, J. N. Mycorrhizae and succession in plantings of beachgrass in sand dunes. Am. J. Bot. 84(1), 118–130 (1997).Article 

Google Scholar 
11.Smith, M. E., Facelli, J. M. & Cavagnaro, T. R. Interactions between soil properties, soil microbes and plants in remnant-grassland and old-field areas: a reciprocal transplant approach. Plant Soil 433, 127–145. https://doi.org/10.1007/s11104-018-3823-2 (2018).CAS 
Article 

Google Scholar 
12.Tipton, A. G., Middleton, E. L., Spollen, W. G. & Galen, C. Anthropogenic and soil environmental drivers of arbuscular mycorrhizal community composition differ between grassland ecosystems. Botany 97, 85–99. https://doi.org/10.1139/cjb-2018-0072 (2019).Article 

Google Scholar 
13.Hamman, S. T. & Hawkes, C. V. Biogeochemical and microbial legacies of non-native grasses can affect restoration success. Restor. Ecol. 21, 58–66. https://doi.org/10.1111/j.1526-100X.2011.00856.x (2013).Article 

Google Scholar 
14.Emery, S. M. & Rudgers, J. A. Beach restoration efforts influenced by plant variety, soil inoculum, and site effects. J. Coast. Res. 27, 636. https://doi.org/10.2112/jcoastres-d-10-00120.1 (2010).Article 

Google Scholar 
15.Sylvia, D. M., Jarstfer, A. G. & Vosátka, M. Comparisons of vesicular-arbuscular mycorrhizal species and inocula formulations in a commercial nursery and on diverse Florida beaches. Biol. Fertil. Soils 16, 139–144. https://doi.org/10.1007/BF00369416 (1993).Article 

Google Scholar 
16.Sylvia, D. M. & Will, M. E. Establishment of vesicular-arbuscular mycorrhizal fungi and other microorganisms on a beach replenishment site in Florida. Appl. Environ. Microbiol. 54, 348–352 (1988).ADS 
CAS 
Article 

Google Scholar 
17.Wubs, E. R. J., van der Putten, W. H., Bosch, M. & Bezemer, T. M. Soil inoculation steers restoration of terrestrial ecosystems. Nat. Plants 2, 16107. https://doi.org/10.1038/nplants.2016.107 (2016).Article 
PubMed 

Google Scholar 
18.Bothe, H., Turnau, K. & Regvar, M. The potential role of arbuscular mycorrhizal fungi in protecting endangered plants and habitats. Mycorrhiza 20, 445–457. https://doi.org/10.5586/asbp.2008.019 (2010).Article 
PubMed 

Google Scholar 
19.Middleton, E. L. & Bever, J. D. Inoculation with a native soil community advances succession in a grassland restoration. Restor. Ecol. 20, 218–226. https://doi.org/10.1111/j.1526-100X.2010.00752.x (2012).Article 

Google Scholar 
20.Crawford, K. M., Busch, M. H., Locke, H. & Luecke, N. C. Native soil microbial amendments generate trade-offs in plant productivity, diversity, and soil stability in coastal dune restorations. Restor. Ecol. https://doi.org/10.1111/rec.13073 (2019).Article 

Google Scholar 
21.Eom, A. H., Hartnett, D. C. & Wilson, G. W. T. Host plant species effects on arbuscular mycorrhizal fungal communities in tallgrass prairie. Oecologia 122, 435–444. https://doi.org/10.1007/s004420050050 (2000).ADS 
Article 
PubMed 

Google Scholar 
22.Brundrett, M. C. & Tedersoo, L. Evolutionary history of mycorrhizal symbioses and global host plant diversity. New Phytol. 220, 1108–1115 (2018).Article 

Google Scholar 
23.Bever, J. D., Mangan, S. A. & Alexander, H. M. Maintenance of plant species diversity by pathogens. Annu. Rev. Ecol. Evol. Syst. 46, 305–325. https://doi.org/10.1146/annurev-ecolsys-112414-054306 (2015).Article 

Google Scholar 
24.Crawford, K. M. et al. When and where plant-soil feedback may promote plant coexistence: a meta-analysis. Ecol. Lett. 22, 13278. https://doi.org/10.1111/ele.13278 (2019).Article 

Google Scholar 
25.Mills, K. E. & Bever, J. D. Maintenance of diversity within plant communities: Soil pathogens as agents of negative feedback. Ecology 79, 1595–1601. https://doi.org/10.1890/0012-9658(1998)079[1595:MODWPC]2.0.CO;2 (1998).Article 

Google Scholar 
26.Koziol, L. et al. The plant microbiome and native plant restoration: The example of native mycorrhizal fungi. Bioscience 68, 996–1006 (2018).Article 

Google Scholar 
27.Maltz, M. R. & Treseder, K. K. Sources of inocula influence mycorrhizal colonization of plants in restoration projects: A meta-analysis. Restor. Ecol. 23, 625–634. https://doi.org/10.1111/rec.12231 (2015).Article 

Google Scholar 
28.Koziol, L. & Bever, J. D. AMF, phylogeny, and succession: Specificity of response to mycorrhizal fungi increases for late-successional plants. Ecosphere https://doi.org/10.1002/ecs2.1555 (2016).Article 

Google Scholar 
29.Middleton, E. L. et al. Locally adapted arbuscular mycorrhizal fungi improve vigor and resistance to herbivory of native prairie plant species. Ecosphere 6, 276. https://doi.org/10.1890/ES15-00152.1 (2015).Article 

Google Scholar 
30.Solís-Domínguez, F. A., Valentín-Vargas, A., Chorover, J. & Maier, R. M. Effect of arbuscular mycorrhizal fungi on plant biomass and the rhizosphere microbial community structure of mesquite grown in acidic lead/zinc mine tailings. Sci. Total Environ. 409, 1009–1016. https://doi.org/10.1016/j.scitotenv.2010.11.020 (2011).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
31.Vogelsang, K. M., Reynolds, H. L. & Bever, J. D. Mycorrhizal fungal identity and richness determine the diversity and productivity of a tallgrass prairie system. New Phytol. 172, 554–562. https://doi.org/10.1111/j.1469-8137.2006.01854.x (2006).Article 
PubMed 

Google Scholar 
32.Larimer, A. L., Bever, J. D. & Clay, K. Consequences of simultaneous interactions of fungal endophytes and arbuscular mycorrhizal fungi with a shared host grass. Oikos 121, 2090–2096. https://doi.org/10.1111/j.1600-0706.2012.20153.x (2012).Article 

Google Scholar 
33.Sikes, B. A., Cottenie, K. & Klironomos, J. N. Plant and fungal identity determines pathogen protection of plant roots by arbuscular mycorrhizas. J. Ecol. 97, 1274–1280. https://doi.org/10.1111/j.1365-2745.2009.01557.x (2009).Article 

Google Scholar 
34.Defeo, O. et al. Threats to sandy beach ecosystems: A review. Estuar. Coast. Shelf Sci. 81, 1–12 (2009).ADS 
Article 

Google Scholar 
35.Feagin, R. A. et al. Going with the flow or against the grain? The promise of vegetation for protecting beaches, dunes, and barrier islands from erosion. Front. Ecol. Environ. 13, 203–210 (2015).Article 

Google Scholar 
36.Feagin, R. A. et al. The role of beach and sand dune vegetation in mediating wave run up erosion. Estuar Coast Shelf Sci. 219, 97–106. https://doi.org/10.1016/j.ecss.2019.01.018 (2019).ADS 
Article 

Google Scholar 
37.Sigren, J. M., Figlus, J. & Armitage, A. R. Coastal sand dunes and dune vegetation: Restoration, erosion, and storm protection. Shore Beach 82, 5–12 (2014).
Google Scholar 
38.Sigren, J. M. et al. The effects of coastal dune volume and vegetation on storm-induced property damage: Analysis from Hurricane Ike. J. Coast Res. 341, 164–173. https://doi.org/10.2112/jcoastres-d-16-00169.1 (2018).Article 

Google Scholar 
39.Silva, R. et al. Response of vegetated dune-beach systems to storm conditions. Coast. Eng. 109, 53–62. https://doi.org/10.1016/j.coastaleng.2015.12.007 (2016).Article 

Google Scholar 
40.Lane, C., Wright, S. J., Roncal, J. & Maschinski, J. Characterizing environmental gradients and their influence on vegetation zonation in a subtropical coastal sand dune system. J. Coast. Res. 4, 213–224. https://doi.org/10.2112/07-0853.1 (2008).CAS 
Article 

Google Scholar 
41.Miller, T. E., Gornish, E. S. & Buckley, H. L. Climate and coastal dune vegetation: Disturbance, recovery, and succession. Plant Ecol. 206, 97–104. https://doi.org/10.1007/s11258-009-9626-z (2010).Article 

Google Scholar 
42.Hewitt, E. J. & Eden, A. Sand and water culture methods used in the study of plant nutrition. Analyst 78, 329–330 (1953).
Google Scholar 
43.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, Austria). https://www.R-project.org/ (2020).
Google Scholar 
44.Farrer, E. C. & Goldberg, D. E. Litter drives ecosystem and plant community changes in cattail invasion. Ecol. Appl. 19, 398–412. https://doi.org/10.1890/08-0485.1 (2009).Article 
PubMed 

Google Scholar 
45.Bauer, J. T., Koziol, L. & Bever, J. D. Local adaptation of mycorrhizae communities changes plant community composition and increases aboveground productivity. Oecologia https://doi.org/10.1007/s00442-020-04598-9 (2020).Article 
PubMed 

Google Scholar 
46.Ohsowski, B. M., Klironomos, J. N., Dunfield, K. E. & Hart, M. M. The potential of soil amendments for restoring severely disturbed grasslands. Appl. Soil. Ecol. 60, 77–83. https://doi.org/10.1016/j.apsoil.2012.02.006 (2012).Article 

Google Scholar 
47.Koziol, L. & Bever, J. D. The missing link in grassland restoration: arbuscular mycorrhizal fungi inoculation increases plant diversity and accelerates succession. J. Appl. Ecol. 54, 1301–1309. https://doi.org/10.1111/1365-2664.12843 (2017).Article 

Google Scholar 
48.Bertness, M. D. & Callaway, R. Positive interactions in communities. Trends Ecol. Evol. 9, 191–193. https://doi.org/10.1016/0169-5347(94)90088-4 (1994).CAS 
Article 
PubMed 

Google Scholar 
49.Heneghan, L. et al. Integrating soil ecological knowledge into restoration management. Restor. Ecol. 16, 608–617. https://doi.org/10.1111/j.1526-100X.2008.00477.x (2008).Article 

Google Scholar 
50.Wubs, E. R. J. et al. Single introductions of soil biota and plants generate long-term legacies in soil and plant community assembly. Ecol. Lett. 22, 1145–1151 (2019).Article 

Google Scholar 
51.Hestrin, R., Hammer, E. C., Mueller, C. W. & Lehmann, J. Synergies between mycorrhizal fungi and soil microbial communities increase plant nitrogen acquisition. Commun. Biol. 2, 233–242. https://doi.org/10.1038/s42003-019-0481-8 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar  More

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