Ecology

1.Food and Agriculture Organization of the United Nations. The State of World Fisheries and Aquaculture http://www.fao.org/3/i2727e/i2727e00.htm (2012).2.Isbell, F. et al. Biodiversity increases the resistance of ecosystem productivity to climate extremes. Nature 526, 574–577 (2015).ADS 
CAS 
PubMed 
Article 

Google Scholar 
3.IPBES. Summary for policymakers of the global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services https://doi.org/10.5281/zenodo.3553579 (2019).4.Tittensor, D. P. et al. A mid-term analysis of progress toward international biodiversity targets. Science 346, 241–244 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
5.Dornelas, M. et al. Assemblage time series reveal biodiversity change but not systematic loss. Science 344, 296–299 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
6.Gonzalez, A. et al. Estimating local biodiversity change: a critique of papers claiming no net loss of local diversity. Ecology 97, 1949–1960 (2016).PubMed 
Article 

Google Scholar 
7.Elahi, R. et al. Recent trends in local-scale marine biodiversity reflect community structure and human impacts. Curr. Biol. 25, 1938–1943 (2015).CAS 
PubMed 
Article 

Google Scholar 
8.Blowes, S. A. et al. The geography of biodiversity change in marine and terrestrial assemblages. Science 366, 339–345 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
9.McCann, K. S. The diversity–stability debate. Nature 405, 228–233 (2000).CAS 
PubMed 
Article 

Google Scholar 
10.Loreau, M. Biodiversity and ecosystem functioning: current knowledge and future challenges. Science 294, 804–808 (2001).ADS 
CAS 
PubMed 
Article 

Google Scholar 
11.Covich, A. P. et al. The role of biodiversity in the functioning of freshwater and marine benthic ecosystems. Bioscience 54, 767–775 (2004).Article 

Google Scholar 
12.Hooper, D. U. et al. Effects of biodiversity on ecosystem functioning: a consensus of current knowledge. Ecol. Monogr. 75, 3–35 (2005).Article 

Google Scholar 
13.Widdicombe, C. E., Eloire, D., Harbour, D., Harris, R. P. & Somerfield, P. J. Long-term phytoplankton community dynamics in the Western English Channel. J. Plankton Res. 32, 643–655 (2010).Article 

Google Scholar 
14.Eloire, D. et al. Temporal variability and community composition of zooplankton at station L4 in the Western Channel: 20 years of sampling. J. Plankton Res. 32, 657–679 (2010).Article 

Google Scholar 
15.Hillebrand, H. et al. In Handbook on Marine Environment Protection (eds Salomon, M. & Markus, T.) 21 (Springer, 2018).16.Bindoff, N. L. et al. In IPCC Special Report on the Ocean and Cryosphere in a Changing Climate (eds H.-O. Pörtner, D. C. Roberts, V. Masson-Delmotte, P. Zhai, M. Tignor, E. Poloczanska, K. Mintenbeck, A. Alegría, M. Nicolai, A. Okem, J. Petzold, B. Rama, N. M. W) Cambridge University Press (2019).17.Barton, A. D., Irwin, A. J., Finkel, Z. V. & Stock, C. A. Anthropogenic climate change drives shift and shuffle in North Atlantic phytoplankton communities. Proc. Natl Acad. Sci. USA 113, 2964–2969 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Pecuchet, L. et al. Spatio‐temporal dynamics of multi‐trophic communities reveal ecosystem‐wide functional reorganization. Ecography 43, 197–208 (2020).Article 

Google Scholar 
19.Poloczanska, E. S. et al. Responses of marine organisms to climate change across oceans. Front. Mar. Sci. 3, 62 (2016).20.Pennekamp, F. et al. Biodiversity increases and decreases ecosystem stability. Nature 563, 109–112 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
21.Duffy, J. E., Godwin, C. M. & Cardinale, B. J. Biodiversity effects in the wild are common and as strong as key drivers of productivity. Nature 549, 261–264 (2017).ADS 
CAS 
PubMed 
Article 

Google Scholar 
22.Worm, B. et al. Impacts of biodiversity loss on ocean ecosystem services. Science 314, 787–790 (2006).ADS 
CAS 
PubMed 
Article 

Google Scholar 
23.Blois, J. L., Zarnetske, P. L., Fitzpatrick, M. C. & Finnegan, S. Climate change and the past, present, and future of biotic interactions. Science 341, 499–504 (2013).ADS 
CAS 
PubMed 
Article 

Google Scholar 
24.Dossena, M. et al. Warming alters community size structure and ecosystem functioning. Proc. R. Soc. B Biol. Sci. 279, 3011–3019 (2012).Article 

Google Scholar 
25.Brander, K. & Kiørboe, T. Decreasing phytoplankton size adversely affects ocean food chains. Glob. Chang. Biol. https://doi.org/10.1111/gcb.15216 (2020).26.Mouw, C. B., Barnett, A., McKinley, G. A., Gloege, L. & Pilcher, D. Phytoplankton size impact on export flux in the global ocean. Glob. Biogeochem. Cycles 30, 1542–1562 (2016).ADS 
CAS 
Article 

Google Scholar 
27.Riahi, K. et al. RCP 8.5—a scenario of comparatively high greenhouse gas emissions. Clim. Change 109, 33–57 (2011).ADS 
CAS 
Article 

Google Scholar 
28.Magnan, A. K. et al. Implications of the Paris agreement for the ocean. Nat. Clim. Chang. 6, 732–735 (2016).ADS 
Article 

Google Scholar 
29.Kuhn, A. M. et al. Temporal and spatial scales of correlation in marine phytoplankton communities. J. Geophys. Res. Ocean. 124, 9417–9438 (2019).ADS 
Article 

Google Scholar 
30.Sonnewald, M., Dutkiewicz, S., Hill, C. & Forget, G. Elucidating ecological complexity: unsupervised learning determines global marine eco-provinces. Sci. Adv. 6, eaay4740 (2020).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
31.Dutkiewicz, S., Boyd, P. W. & Riebesell, U. Exploring biogeochemical and ecological redundancy in phytoplankton communities in the global ocean. Glob. Chang. Biol. 27, 1196–1213 (2021).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
32.Flombaum, P., Wang, W.-L., Primeau, F. W. & Martiny, A. C. Global picophytoplankton niche partitioning predicts overall positive response to ocean warming. Nat. Geosci. 13, 116–120 (2020).ADS 
CAS 
Article 

Google Scholar 
33.Righetti, D., Vogt, M., Gruber, N., Psomas, A. & Zimmermann, N. E. Global pattern of phytoplankton diversity driven by temperature and environmental variability. Sci. Adv. 5, eaau6253 (2019).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Ibarbalz, F. M. et al. Global trends in marine plankton diversity across kingdoms of life. Cell 179, 1084–1097.e21 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Bopp, L. et al. Multiple stressors of ocean ecosystems in the 21st century: projections with CMIP5 models. Biogeosciences 10, 6225–6245 (2013).ADS 
Article 

Google Scholar 
36.Kwiatkowski, L. et al. Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model projections. Biogeosciences 17, 3439–3470 (2020).ADS 
CAS 
Article 

Google Scholar 
37.Cabré, A., Marinov, I. & Leung, S. Consistent global responses of marine ecosystems to future climate change across the IPCC AR5 earth system models. Clim. Dyn. 45, 1253–1280 (2015).Article 

Google Scholar 
38.Bopp, L., Aumont, O., Cadule, P., Alvain, S. & Gehlen, M. Response of diatoms distribution to global warming and potential implications: a global model study. Geophys. Res. Lett. 32, n/a−n/a (2005).Article 
CAS 

Google Scholar 
39.Dutkiewicz, S. et al. Dimensions of marine phytoplankton diversity. Biogeosciences 17, 609–634 (2020).ADS 
Article 

Google Scholar 
40.Dutkiewicz, S., Scott, J. R. & Follows, M. J. Winners and losers: ecological and biogeochemical changes in a warming ocean. Glob. Biogeochem. Cycles 27, 463–477 (2013).ADS 
CAS 
Article 

Google Scholar 
41.Marinov, I., Doney, S. C. & Lima, I. D. Response of ocean phytoplankton community structure to climate change over the 21st century: partitioning the effects of nutrients, temperature, and light. Biogeosciences 7, 3941–3959 (2010).ADS 
Article 

Google Scholar 
42.Dutkiewicz, S., Ward, B. A., Scott, J. R. & Follows, M. J. Understanding predicted shifts in diazotroph biogeography using resource competition theory. Biogeosciences 11, 5445–5461 (2014).ADS 
Article 

Google Scholar 
43.Dutkiewicz, S. et al. Impact of ocean acidification on the structure of future phytoplankton communities. Nat. Clim. Chang. 5, 1002–1006 (2015).ADS 
CAS 
Article 

Google Scholar 
44.Kooijman, S. A. L. M. & Troost, T. A. Quantitative steps in the evolution of metabolic organisation as specified by the dynamic energy budget theory. Biol. Rev. 82, 113–142 (2007).CAS 
PubMed 
Article 

Google Scholar 
45.Lévy, M., Jahn, O., Dutkiewicz, S., Follows, M. J. & D’Ovidio, F. The dynamical landscape of marine phytoplankton diversity. J. R. Soc. Interface 12, 20150481 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
46.Beaugrand, G., Edwards, M., Raybaud, V., Goberville, E. & Kirby, R. R. Future vulnerability of marine biodiversity compared with contemporary and past changes. Nat. Clim. Chang. 5, 695–701 (2015).ADS 
Article 

Google Scholar 
47.Thomas, M. K., Kremer, C. T., Klausmeier, C. A. & Litchman, E. A global pattern of thermal adaptation in marine phytoplankton. Science 338, 1085–1088 (2012).ADS 
CAS 
PubMed 
Article 

Google Scholar 
48.Hillebrand, H. et al. Biodiversity change is uncoupled from species richness trends: consequences for conservation and monitoring. J. Appl. Ecol. 55, 169–184 (2018).Article 

Google Scholar 
49.Litchman, E. & Klausmeier, C. A. Trait-based community ecology of phytoplankton. Annu. Rev. Ecol. Evol. Syst. 39, 615–639 (2008).Article 

Google Scholar 
50.Lindeman, R. L. The trophic-dynamic aspect of ecology. Ecology 23, 399–417 (1942).Article 

Google Scholar 
51.Stock, C. A. et al. Reconciling fisheries catch and ocean productivity. Proc. Natl Acad. Sci. USA 114, E1441–E1449 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
52.Armengol, L., Calbet, A., Franchy, G., Rodríguez-Santos, A. & Hernández-León, S. Planktonic food web structure and trophic transfer efficiency along a productivity gradient in the tropical and subtropical Atlantic Ocean. Sci. Rep. 9, 2044 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
53.Cram, J. A. et al. The role of particle size, ballast, temperature, and oxygen in the sinking flux to the deep sea. Glob. Biogeochem. Cycles 32, 858–876 (2018).ADS 
CAS 
Article 

Google Scholar 
54.Kéfi, S., Dakos, V., Scheffer, M., Van Nes, E. H. & Rietkerk, M. Early warning signals also precede non-catastrophic transitions. Oikos 122, 641–648 (2013).Article 

Google Scholar 
55.Doncaster, C. P. et al. Early warning of critical transitions in biodiversity from compositional disorder. Ecology 97, 3079–3090 (2016).PubMed 
Article 

Google Scholar 
56.Gunderson, L. H. Ecological resilience—in theory and application. Annu. Rev. Ecol. Syst. 31, 425–439 (2000).Article 

Google Scholar 
57.Benedetti, F. et al. The seasonal and inter-annual fluctuations of plankton abundance and community structure in a North Atlantic Marine Protected Area. Front. Mar. Sci. 6, 214 (2019).58.Pannard, A., Bormans, M. & Lagadeuc, Y. Short-term variability in physical forcing in temperate reservoirs: effects on phytoplankton dynamics and sedimentary fluxes. Freshw. Biol. 52, 12–27 (2007).CAS 
Article 

Google Scholar 
59.Vidal, T., Calado, A. J., Moita, M. T. & Cunha, M. R. Phytoplankton dynamics in relation to seasonal variability and upwelling and relaxation patterns at the mouth of Ria de Aveiro (West Iberian Margin) over a four-year period. PLoS One 12, e0177237 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
60.Cermeño, P., de Vargas, C., Abrantes, F. & Falkowski, P. G. Phytoplankton biogeography and community stability in the ocean. PLoS One 5, e10037 (2010).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
61.Allen, S. et al. Interannual stability of phytoplankton community composition in the North-East Atlantic. Mar. Ecol. Prog. Ser. 655, 43–57 (2020).ADS 
Article 

Google Scholar 
62.Barton, A. D., Lozier, M. S. & Williams, R. G. Physical controls of variability in North Atlantic phytoplankton communities. Limnol. Oceanogr. 60, 181–197 (2015).ADS 
Article 

Google Scholar 
63.Collins, S., Rost, B. & Rynearson, T. A. Evolutionary potential of marine phytoplankton under ocean acidification. Evol. Appl. 7, 140–155 (2014).CAS 
PubMed 
Article 

Google Scholar 
64.Lohbeck, K. T., Riebesell, U. & Reusch, T. B. H. Adaptive evolution of a key phytoplankton species to ocean acidification. Nat. Geosci. 5, 346–351 (2012).ADS 
CAS 
Article 

Google Scholar 
65.Irwin, A. J., Finkel, Z. V., Müller-Karger, F. E. & Troccoli Ghinaglia, L. Phytoplankton adapt to changing ocean environments. Proc. Natl Acad. Sci. USA 112, 5762–5766 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
66.Cael, B. B. et al. Marine ecosystem changepoints spread under ocean warming in an Earth System Model. Geophys. Res. Lett.67.Cael, B. B., Dutkiewicz, S. & Henson, S. A. Abrupt shifts in 21st-century plankton communities. Sci. Adv.68.Parmesan, C. & Yohe, G. A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37–42 (2003).ADS 
CAS 
Article 

Google Scholar 
69.Burrows, M. T. et al. The pace of shifting climate in marine and terrestrial ecosystems. Science 334, 652–655 (2011).ADS 
CAS 
PubMed 
Article 

Google Scholar 
70.Chivers, W. J., Walne, A. W. & Hays, G. C. Mismatch between marine plankton range movements and the velocity of climate change. Nat. Commun. 8, 14434 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
71.Poloczanska, E. S. et al. Global imprint of climate change on marine life. Nat. Clim. Chang. 3, 919–925 (2013).ADS 
Article 

Google Scholar 
72.Jonkers, L., Hillebrand, H. & Kucera, M. Global change drives modern plankton communities away from the pre-industrial state. Nature 570, 372–375 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
73.Pond, D. W., Tarling, G. A. & Mayor, D. J. Hydrostatic pressure and temperature effects on the membranes of a seasonally migrating marine copepod. PLoS One 9, e111043 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
74.Mayor, D. J., Sommer, U., Cook, K. B. & Viant, M. R. The metabolic response of marine copepods to environmental warming and ocean acidification in the absence of food. Sci. Rep. 5, 13690 (2015).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
75.Richardson, D. M. & Pyšek, P. Elton, C.S. 1958: The ecology of invasions by animals and plants. London: Methuen. Prog. Phys. Geogr. Earth Environ. 31, 659–666 (2007).Article 

Google Scholar 
76.May, R. M. Qualitative stability in model ecosystems. Ecology 54, 638–641 (1973).Article 

Google Scholar 
77.Lotze, H. K. et al. Global ensemble projections reveal trophic amplification of ocean biomass declines with climate change. Proc. Natl Acad. Sci. USA 116, 12907–12912 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
78.Barange, M. et al. Impacts of climate change on marine ecosystem production in societies dependent on fisheries. Nat. Clim. Chang. 4, 211–216 (2014).ADS 
Article 

Google Scholar 
79.Marañón, E. et al. Unimodal size scaling of phytoplankton growth and the size dependence of nutrient uptake and use. Ecol. Lett. 16, 371–379 (2013).PubMed 
Article 

Google Scholar 
80.Sokolov, A. P. et al. The MIT Integrated Global System Model (IGSM) Version 2: Model Description and Baseline Evaluation Joint Program Report Series, pp. 40 https://globalchange.mit.edu/publication/14579 (2005).81.Dutkiewicz, S. et al. Ocean colour signature of climate change. Nat. Commun. 10, 578 (2019).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
82.Monier, E., Scott, J. R., Sokolov, A. P., Forest, C. E. & Schlosser, C. A. An integrated assessment modeling framework for uncertainty studies in global and regional climate change: the MIT IGSM-CAM (version 1.0). Geosci. Model Dev. 6, 2063–2085 (2013).ADS 
Article 

Google Scholar 
83.Moss, R. H. et al. The next generation of scenarios for climate change research and assessment. Nature 463, 747–756 (2010).ADS 
CAS 
PubMed 
Article 

Google Scholar 
84.Buitenhuis, E. T. et al. MAREDAT: towards a world atlas of MARine Ecosystem DATa. Earth Syst. Sci. Data 5, 227–239 (2013).ADS 
Article 

Google Scholar 
85.Ward, B. A. Temperature-correlated changes in phytoplankton community structure are restricted to polar waters. PLoS One 10, e0135581 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
86.Dutkiewicz, S. GUD IGSM depth integrated biomass https://doi.org/10.7910/DVN/LWHQNS (2021).87.Dutkiewicz, S. & Jahn, O. GUD IGSM numerical code and inputs https://doi.org/10.7910/DVN/UA8VNU (2021). More

The model is organized as follows. We simulated water striders as an array of the discrete “objects” which interact one with another according to more or less natural rules. The interaction includes strong short-range repulsion between the animals. Mathematically, short range repulsion means that one animal cannot penetrate inside a private area of other animal. Normally such repulsion appears at relatively short distances corresponding to a radius (R^{repuls}) of their private territory. Besides, there is mutual attraction at longer distances (R^{attract} > R^{repuls}), which biologically corresponds to a tendency to aggregation4. The tendency to aggregation often gives some competitive advantage due to possible collective reactions on the external challenges (for example, on the attacks of predators).The simplest way to simulate an interaction with regulated characteristic distance is to use Gaussian effective potential with some characteristic radius (R_{0}), since the Gaussian potential guaranty the stability in dynamics simulations with relatively large (Delta t)4. Repulsion force in this case looks as follows:$$f_{j}^{rep} = A_{0} left( {vec{r} – vec{r}_{j} } right)expleft[ { – left( {frac{{vec{r} – vec{r}_{j} }}{{R_{0} }}} right)^{2} } right],$$
(2)
where factor A0  R^{repuls}) attraction normally causes a minimum of the interaction energy at some intermediate distance (R^{min }). In infinite empty space, being used in the equations of motion, such a combination of the forces leads to an ordering of the objects with the equilibrium distance corresponding to the position (R_{min }) of the energy minimum. But, if the area ({ [0,Lx],[0,Ly]}) is limited, the equations of motion must be supplied by appropriate boundary conditions which do not allow the animals to leave this area.The simplest way to introduce the boundary conditions is to apply mathematically “soft” but extremely high and narrow walls around the system, which repulse the animals back to the internal space with exponentially growing force$$f_{j}^{bound} = B_{j}^{bound} exp left[ { – left| {frac{{overrightarrow {r}_{j} – overrightarrow {r}_{bound} }}{{R_{{}}^{bound} }}} right|} right]$$
(5)
acting in the direction opposite to that in which the animal occasionally crosses any of the boundaries. The words “extremely high” mean that the amplitude of this force should be supplied by the pre-factor (B_{j}^{bound}) which is much bigger than the amplitude of the repulsion force between the animals (regulated by the amplitude (B_{{_{jk} }}^{repuls})) to be able to overpower their mutual repulsion (B_{j}^{bound} > > B_{j}^{repuls}) pushing them out the boundaries. Therefore also the characteristic distance of the repulsion from the wall should be much shorter than the typical distance between the animals in the population (R^{bond} < < R^{min }). In this sense the boundary should be as narrow as possible.One can expect that mean density of the population is not extremely high and does not force the animals to stay exactly on the minimal distance (R^{min }) between them. In this case total combination of the attracting and repulsing interactions combined with the rejecting boundary conditions normally leads to the specific patterns where relatively dense groups of the animals spread on the distances close to the equilibrium ones are accompanied by almost empty voids between them5. We have checked this hypothesis numerically many times. Below, such patterns will be seen in all the particular visualizations of the numerical results.It should be noted that mathematically, if number of the animals increases inside of the same limited and already populated area, the individuals tend to fill all the voids with the equilibrium (sometimes, even higher than equilibrium) density. It happens in real population as well in that cases when the number of the animals grows quickly for a self-regulation of the population and they become simply forced to occupy every empty space inside the area.However, normally when the population grows too quickly the growth should be regulated by a number of factors. In particular, it will be restricted by a competition for the food and space accompanied by death of some of the participants of the process. From mathematical point of view it means that the equations of motion, which will be written below, must be accompanied by natural generation of the new individuals as well as by their (reasonable) disappearance from the system.It is almost impossible to write such a generation in analytical form, but in numerical model it can be formally presented as a set of natural rules. First of all the length of the array is supposed to be a variable integer. Let us suppose that at every step of calculation (Delta t) the length of the array can potentially increase to the new one (N to N + 1). The potential act of the generation becomes real if the numerically generated random number (varsigma) uniformly distributed in interval [0;1] is less than some threshold (varsigma_{thres} < 1). In the numerical procedure the rate of the generation is certainly regulated by (varsigma_{{{text{thres}}}}) and can be extremely low for example, if (varsigma_{thres} < < 1).New member of the array is generated with random coordinates within the prescribed area ({ [0,Lx],[0,Ly]}) with the mass ((m_{N + 1})) randomly distributed around a starting (basically small) one (m_{s}). Qualitatively it means that at every time moment the array can be supplemented with some probability by a new member which is physically placed in some arbitrary place inside the area ({ [0,Lx],[0,Ly]}). At the same time, we have to introduce a process of disappearance of some of the array members. Natural criterion for this will be established below from the checking of the accumulation or loss of the mass (e.g. reserve fat mass) (m_{j}) by every individual. It is supposed that, if an animal gets critical size (which can be both: maximal (m_{max }) or minimal (m_{min })) it disappears from the system and the length of the array decreases (N to N - 1).Both of such events must be adjusted numerically to some rate natural for a particular biological system. It is quite expected from the very beginning, that there should be a kind of balance between the creation and disappearance rates. If new animas statistically appear too often, the overpopulation will take place. In opposite limit the array will quickly shrink (N to 0) and it will cause a distinction of the population.Obviously both these rates should be naturally regulated by the available food resources. For the particular problem under consideration the resources are generated by a random deposition of the potentially available food onto the water surface. It means that we introduce a new array for the food with coordinates (overrightarrow {r}_{n}^{food}). The length of the array (N^{food}) is also variable and index (n) running in the interval (n = 0,..,N^{food}). In principle, it is possible, and very often happened in our simulations, that at some particular moment available food inside the area can completely disappear. In this case the length of the food array reduces to zero (N^{food} = 0).The food is generated by the random deposition of the food portions distributed inside the area of the system ({ [0,Lx],[0,Ly]}). Generally, it is organized in the same manner as the deposition of the new individuals. If the random number (zeta) uniformly distributed in the interval [0;1] is less than some threshold (varsigma_{thres} < 1) the food portion is deposed at an arbitrary place of the area, in principle at any given time step (Delta t).It is obvious that if the threshold is much smaller than unit: (zeta_{thres}^{food} < < 1) the food portions are produced very rarely. Certainly, the rate of the “food production” in the frames of the model must be properly regulated to make its behaviour natural. Of course, it relates to the size of the portions, but also to the intervals between the food depositions. In any case, these intervals are expected to be much longer than discrete time steps (Delta t) used to solve numerically the equations of motion. At the same time, it must be much shorter than other (biologically reasonable) time-scales of the problem. To control the stability and reasonability of the simulations we have varied this rate in very wide intervals. It was found that at all the reasonable rates food balance in the system tends to the stationary scenario, which corresponds to the expected scenarios in nature.When a portion of the food falls onto the surface the animals which occasionally appear relatively close to it are attracted to this food portion and “eat” it with some characteristic rate. In the particular simulation it was simulated by the additional term of the attraction force:$$f_{{_{jn} }}^{food} (overrightarrow {r}_{j} ,overrightarrow {r}_{n} ) = B_{jn}^{food} (overrightarrow {r}_{j} - overrightarrow {r}_{n} )exp left[ { - left( {frac{{overrightarrow {r}_{j} - overrightarrow {r}_{n} }}{{R_{{_{j} }}^{food} }}} right)^{2} } right]$$ (6) As it is seen from the interactions in the model the animals compete for the food, repulsing one another and reacting faster on the force (f_{{_{jk} }}^{food} (overrightarrow {r}_{j} ,overrightarrow {r}_{k} )). It’s why we apply here non-uniform coefficient (B_{{_{jk} }}^{food}), which is different for different index (j) and in fact depends on some power of mass (B_{jk}^{food} sim m_{j}^{alpha }) with some exponent (alpha). Scaling estimation gives the value of the exponent (alpha) = 2/3. Another form of the competition is related to their simultaneous consumption the food with different rate depending on the size of the different individuals. We suppose that the consumption is proportional to already accumulated mass of the individual: (partial m_{j} /partial t = mu m_{j}). Effective dumping (eta_{j}) on the water surface also depends on the mass of animal (eta_{j} sim m_{j}^{beta }), where estimated exponent (beta) = 1/3.Let us remind that the consumption of the food portion with index (n) is possible only when the animal is close enough to this portion. In other words, one more threshold has to be incorporated: (left| {overrightarrow {r}_{j} - overrightarrow {r}_{n} } right| < R_{thres}^{food}). As we see, the model consists of the dynamic equations of motion and a number of the procedures, which work in parallel and essentially affect both: the dynamic behavior and the results.The equations of motion can be formally written accumulating all the above mentioned forces of the problem:$$m_{j} partial^{2} r_{j} /partial t^{2} + eta_{j} partial r_{j} /partial t = sumlimits_{k} {left[ {f_{{_{jk} }}^{repuls} (overrightarrow {r}_{j} ,overrightarrow {r}_{k} ) + f_{{_{jk} }}^{attract} (overrightarrow {r}_{j} ,overrightarrow {r}_{k} )} right]} + sumlimits_{n} {f_{{_{jn} }}^{food} (overrightarrow {r}_{j} ,overrightarrow {r}_{n} )} + sumlimits_{Boundaries} {f_{j}^{bound} }$$ (7) with their solution combined with the procedures described above. All the initial velocities are zero. This combination makes the solution nontrivial because it is performed for the arrays with changeable lengths and with the varied sizes and forces of the participants. Nevertheless, expected scenarios can be generally predicted and verified later by the large set of numerical experiments at different combinations of the parameters.The longer the time an individual spends near the portion of food and the larger the size it has to the current moment the faster it eats up and faster accumulates additional mass. In general, the larger animals consume more food. However, one should remember that large animals may also have disadvantages. For example, due to inertia of motion faster animals can occasionally jump over a good position near the food. We have observed many such local events during the simulations. Such events partially can be observed in the movies presented below. It is also possible that small animals can appear near the randomly deposed food earlier than the big ones and start consuming it.In some cases such possibilities even emerges due to not just probabilistic but quite regular reason. For example, it can happen due to a further described effect of the “wave of fear”. By influence of such a wave (being scared by somebody moving along the shore), the animals almost synchronically start to run away from one of the boundaries. Big, strong and fast animals escape far away, while the small ones still remain near to the shore. As result, the food deposed into this region will be partially (or even completely) consumed by the weak individuals before the strong ones will return. This effect will be reported further on.Basically, dynamic balance of the population is determined by the relation between a set of the time constants of the problem: how often the food appears, how quick the animals consume it and how efficiently they grow. However, the population structure is also important for the total balance. The more an animal consumes the bigger and stronger it becomes, than faster moves to the new portion of the food and more eats. The individuals are born small and more or less equal. They grow initially depending on their “luck” only. However, population quickly forms some non-uniform distribution of the sizes.In terms of distributions one can say that with the time the individual members of the array shift along the distribution to the larger or smaller size depending on the personal balance of accumulation or losing the mass. At every instant moment the distribution actually demonstrates two opposite fluxes: to the limit of the small and large sizes. Strongest animals evolve to the predefined critical large size (m_{max }), at which they leave the population. Weakest ones tend to the critically small weight mmin and also leave the population. The balance is maintained by the fact that at every time moment the population is incorporated by new young members who either grow or die.In this respect, it is interesting to observe the dynamics of the animal’s sizes distribution. Obviously, the initial population consisting of young animals has histogram localized around small sizes. Later the distribution becomes wider and its edge moves to the larger sizes. The asymptotic histogram shape is determined by the two opposite processes described above. This shape will be compared further with a real distribution found from the field observations.Typical behavior of the model evolution for a population inside a limited area is presented in illustrative movie (video_S2). For a convenience, the population is formally divided into 3 subgroups, which are plotted by the circles having different sizes (small, medium, and big) and colors (blue, green, and red, respectively). The randomly deposed portions of food are shown by large black circles.We start from relatively young population, which contains 100 individuals and allow them to move, grow and disappear according to the rules of game. This behavior demonstrates correlative motions, as well as variation in the population composition, which look self-consistent and rather natural. The interactions between the animals with other animals and with food cause quite predictable motions of the individuals and the changes of their sizes (and colors, respectively) when they leave one of the three groups and join to another. The process stabilizes with the time and seems to become stationary. Below we will study this process by the quantitative time-dependencies and statistical histograms.The bigger an animal the faster it moves in average. Therefore, it is convenient to sort the animals according to the masses and velocities. Such representation is reproduced in the third movie “video_S3.mp4”. The separation between subgroups, their evolution from the initial to a stationary state is clearly seen from the movie in dynamics. Moreover, one can observe even how group of the small animals divides by itself into two subgroups with quite well pronounced gap between them. This separation is caused by the mentioned above two fluxes of the sub-populations, which either grow from the initially small sizes to the medium ones or “in unlucky case” decrease and disappear.Actually, the rate of evolution with fast changes of the masses after very few acts of the interaction with deposed food represented in both movies is strongly overestimated in contrast to the reality. Therefore, as a next step we reduced the rates of accumulation and loss of the mass and proportionally prolonged the simulation time. This process was also recorded in two movies “video_S4.mp4” and “video_S5.mp4”. To reduce the lengths of the videos the time intervals between the frames were made 10 times longer. Because of this reduction the movies look almost stroboscopic. Nevertheless, the movies illustrate quite well long-time dynamics of the system at realistic rate of food deposition and consumption.Information about the system accumulated during long runs including typical pattern formed by the moving animals at some intermediate stage and other plots is presented in Figs. 1, 2, 3 and 4. The spatial pattern in Fig. 1 is taken at some intermediate stage of evolution. Blue, green and red circles correspond to the small medium and large sizes of the individuals. Large black circles mark the food portions, which are recently deposed and not eaten yet.Figure 1Typical spatial pattern formed by the population at intermediate stage of evolution. Blue, green and red circles correspond to the small medium and large sizes of the individuals. One can see the groups packed by the individuals with the distances close to the equilibrium Rmin(R^{min }) and voids. Black circles show food portions remaining to the current moment.Full size imageFigure 2Histograms of the distances between nearest neighbors. Blue line shows some instant distribution. Black line presents histogram accumulated during long run. Maximum of the histogram corresponds to the distance Rmin close to the equilibrium. Long tail of the black curve corresponds to the presence of some individuals inside the voids.Full size imageFigure 3Mass depending values: (a) instant configuration of the velocities of different subgroups marked by the same symbols as in Fig. 1; (b) instant and long run averaged histograms of the masses for complete population; (c) instant and accumulated distributions of the velocities monotonously increasing with masses (volumes) of the individuals.Full size imageFigure 4The comparison between the numerically (curves with points) and experimentally (bars) obtained data. The subplot (a) shows the histogram of the distances between nearest neighbors and the subplot (b) presents the histograms of the volumes of the animals (which are supposed to be approximately proportional to their masses).Full size imageOne can well distinguish in Fig. 1 relatively dense domains packed with the distances between the individuals close to the equilibrium radius (R^{min }) determined by the relation between repulsive (f_{{_{jk} }}^{repuls} (overrightarrow {r}_{j} ,overrightarrow {r}_{k} )) and attraction (f_{{_{jk} }}^{attract} (overrightarrow {r}_{j} ,overrightarrow {r}_{k} )) forces. These domains are separated by voids. To characterize such patterns quantitatively the histograms of the distances between nearest neighbors6 were accumulated during a long simulation run. These histograms are shown in Fig. 2. Thin blue line here reproduces some arbitrary instant distribution. Black line with the circles presents the histogram accumulated during a long run. Maximum of the histogram is close to the equilibrium distance (R^{min }) determined by the balance of interactions. Long tail of the black curve, which extends to the few times longer distances than (R^{min }), exists, since some individuals are located inside voids far-away from all the neighbors.Typical size distributions to which the population evolves with the time are reproduced in Fig. 3. It shows three important size depending values. Figure 3a demonstrates an instant configuration of the velocities. Three different subgroups are marked by the same symbols as in Fig. 1. Figure 3b shows instant histogram of the animal’s masses in complete population and the histogram averaged over long run. They are plotted by the blue and black lines respectively. How mean velocities of the animals depend on their mass is shown in Fig. 3c. As above, the instant and accumulated distributions are plotted by the blue and black lines respectively.It can be noticed from the subplot (c), Fig. 3, that the averaged velocities correlate with the mass of the individuals. In principle, it could be expected according to the model rules, namely due to the nonlinear dependence of strength (f_{{_{jk} }}^{repuls} (overrightarrow {r}_{j} ,overrightarrow {r}_{k} )) as well as effective damping (eta_{j}) for the individual animals on their size (mass (m_{j})). However, taken into account multiple interactions in the system the result was not obvious in advance. The velocity increase is especially pronounced when the animals become large and are close to the final stage of its growth.Further, we would like to compare the numerical and experimentally obtained data. This comparison is presented in Fig. 4 by the curves and bar-plots respectively. From the field videos it is difficult to determine the masses (m_{j}) of the individuals, but we can estimate their volumes (V_{j}) which are approximately proportional to their masses (V_{j} sim m_{j}). Therefore, to compare numerical results with the measurements from videos we have plotted the distributions along the volume coordinate. The subplot (a) in Fig. 4 shows the histograms of the distances between nearest neighbors, and the subplot (b) presents the histograms of the volumes of the animals. The histograms obtained from simulations match very well to the histograms calculated using experimental data. Two sample Kolmogorov–Smirnov test does not reject the hypothesis, that the distributions obtained in experiment and in numerical simulations are from the same continues distribution for both the animal volumes (p = 0.59, D55,500 = 0.105) and for the distances between nearest neighbors (p = 0.84, D55,500 = 0.083).It is important to note that due to stability of the system the distributions independent on their initial shape converge to the same quasi-stationary histograms. To check this convergence, we started from two extreme populations. The first one consisted almost from the large individuals only and another one contained only small individuals. Time depending volume histograms were accumulated into the volume distribution over time gray-scale maps. In Fig. 5 the results for the two cases are shown in the subplots (a) and (b) respectively. Lighter gray color corresponds to the higher density. Both distributions attract to the practically the same final one with the time. The shapes of distributions, which have started from the two extreme initial distributions, almost coincide after long runs, while the distributions flanks shift in different directions as marked by the red and blue arrows in Fig. 6 and visualized in the supplementary movies “video_S6.mp4” and “video_S7.mp4” respectively. As in the static figures the instant and time averaged histograms are shown by the blue and black lines respectively. For the supplementary video “video_S6.mp4” the distribution is initially localized mainly near the right side of the interval and after some quick transient period “jumps” to the distribution close to the final histogram. Comparison with second video “video_S7.mp4” shows how both averaged distributions are attracted after a long simulation run to the similar distribution.Figure 5Gray-scale maps of time-volume histograms accumulated for two limit cases: starting from large and small animals, are shown in the subplots (a,b) respectively. Lighter gray color corresponds to the higher density. Dashed line marks common position of distribution maximums to which they converge during extremely long runs.Full size imageFigure 6Final volume distributions after long runs. It is seen that the maximums coincide already. Remaining trends of the histogram alterations for the distribution started from the large (open circles) and small (closed circles) individuals are marked by the red and blue arrows respectively.Full size imageDuring the system evolution some of the animals leave the population and some new appear. One can accumulate the numbers of died and survived animals for each given time moment. Accumulated numbers of died and survived animals are plotted in Fig. 7a by black and red curves respectively. The balance between these numbers depends on the relation of all the model parameters. Here the same set of the parameters was used as for simulations presented in Figs. 1, 2 and 3.Figure 7Time depending integral populations: (a) accumulated numbers of died (N1) and survived (N2) animals (black and red curves); (b) variations of the numbers of animals in subpopulations of small (Nc1), medium (Nc2) and large (Nc3) animals shown by blue, green and red curves respectively.Full size imageThe curves presenting the integral number of animals are relatively smooth, Fig. 7a, yet the instant numbers of animals in small, medium and large subpopulations shown by blue, green and red curves respectively in Fig. 7b vary much stronger, because at every particular time moment, their numbers are relatively small and the fluctuations are strong. From the figure, it is also obvious that during initial transient time interval the number of small individuals do not exceeds the number of two other populations. It is natural and was expected, because all newborn animals are small. Therefore, subpopulation with small size prevails at stationary stage.As we already announced nontrivial transformation of the spatial pattern and modification of all the resulting distributions can appear due to the “waves of fear”. Such a wave can appear when the animals almost synchronously try to escape from the shoreline if they are being scared by somebody moving along. Mathematically such a wave can be provoked by an additional repulsion force acting from one side of the simulation area to the distance much longer than already incorporated short range reflection from the boundary (Eq. 6). This force influences all the members of the array but appears for relatively limited time. It can be added to the model interaction occurring either randomly or periodically. Second variant is more regular and preferable and allows accumulate statistics faster.Analytically this force has the same form as (f_{j}^{bound}) with the same boundary position (overrightarrow {r}_{bound}):$$f_{j}^{wave} = B_{j}^{wave} exp left[ { - left| {frac{{overrightarrow {r}_{j} - overrightarrow {r}_{bound} }}{{R_{{}}^{wave} }}} right|} right],$$ (8) but it has larger amplitude (B_{j}^{wave} > B_{j}^{bound}) and much longer distance of the exponential decay (R_{{}}^{wave} > > R_{{}}^{bound}). So, it influences not only the individuals occasionally attempting to cross the boundary but almost all others too, even if they are currently relatively far from the boundary.The video in supplementary material “video_S8.mp4” illustrates typical behavior of the system at presence of the “waves of fear” occasionally observed in original sequences (video_S1.avi). One can see how bigger, stronger and faster animals escape quicker and further from the dangerous boundary, while some smallest ones do not react so quickly and remain almost near to it. Food deposition is not correlated with the “waves of fear”. As result, the food deposed near the dangerous boundary is consumed mostly by small animals, since that boundary region is practically depopulated by bigger animals. While stronger individuals return to the shore, such food may be already consumed by the weaker individuals.Static image of such “shifted” to one side pattern for a moment when food and weak members of the array are presented in one side of the area while absolute majority of the population is collected in a center or closer to another side of the area is reproduced in Fig. 8.Figure 8Shifted pattern at the presence of the “wave of fear”. Some food was just deposed inside the depopulated area near to the dangerous boundary. Some small individuals are in the vicinity to the food, while absolute majority of the population is localized in a center or near opposite boundary. The symbols have the same meaning as in Fig. 1.Full size imageBeing regularly applied the “waves of fear” in principle leads to a redistribution of the food and animals in the space. It is quite expected that the density of the animals will spend some time in the areas distanced from the dangerous wall. As result, the food will longer remain available in the places close to the dangerous wall. We have accumulated these densities during long runs at different parameters and obtained such shifted distributions. Typical distributions of the food and animals for the set of parameters which were used in simulations presented in the previous figures are shown in Fig. 9.Figure 9Density distributions of the food and animals at presence of regularly applied “waves of fear”. The densities of the food and animals are shown by the solid black and dashed magenta curves respectively. The curves for the densities of large, middle and small individuals are plotted by the same colors as above. The density shift with increasing size is marked by the arrow.Full size imageMean density of the food is shown by the solid black curve. Total density of the animals is presented by the dashed magenta curve. These densities are anti-correlated. Besides, there is a fine structure of the partial densities of the animals. The curves for the densities of large, middle and small individuals are plotted by the same colors as above. Strongest correlation between the size of the animals and their spatial density is observed in interval between two dash-dotted lines. The density shift with increasing size is marked by the arrow. In close proximity to the dangerous wall, where mean density of the food is abnormally high, the density shift dependence on the individuals size change the direction. Number of the small individuals reduces in average, because they quickly move to the category of middle or even large individuals. More

1.Szathmáry E, Smith JM. The major evolutionary transitions. Nature 1995;374:227–32.PubMed 
Article 

Google Scholar 
2.West SA, Fisher RM, Gardner A, Kiers ET. Major evolutionary transitions in individuality. Proc Natl Acad Sci USA. 2015;112:10112–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
3.Moran NA. The coevolution of bacterial endosymbionts and phloem-feeding insects. Ann Mo Bot Gard. 2001;88:35–44.Article 

Google Scholar 
4.Bennett GM, Moran NA. Heritable symbiosis: the advantages and perils of an evolutionary rabbit hole. Proc Natl Acad Sci USA. 2015;112:10169–76.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Gil R, Sabater-Munoz B, Latorre A, Silva FJ, Moya A. Extreme genome reduction in Buchnera spp.: toward the minimal genome needed for symbiotic life. Proc Natl Acad Sci USA. 2002;99:4454–8.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
6.Tamames J, Gil R, Latorre A, Pereto J, Silva FJ, Moya A. The frontier between cell and organelle: genome analysis of Candidatus Carsonella ruddii. BMC Evol Biol. 2007;7:181.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
7.Husnik F, Nikoh N, Koga R, Ross L, Duncan RP, Fujie M, et al. Horizontal gene transfer from diverse bacteria to an insect genome enables a tripartite nested mealybug symbiosis. Cell 2013;153:1567–78.CAS 
PubMed 
Article 

Google Scholar 
8.Wilson ACC, Duncan RP. Signatures of host/symbiont genome coevolution in insect nutritional endosymbioses. Proc Natl Acad Sci USA. 2015;112:10255–61.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
9.von Dohlen CD, Kohler S, Alsop ST, McManus WR. Mealybug β-proteobacterial endosymbionts contain γ-proteobacterial symbionts. Nature 2001;412:433–6.Article 

Google Scholar 
10.McCutcheon JP, McDonald BR, Moran NA. Convergent evolution of metabolic roles in bacterial co-symbionts of insects. Proc Natl Acad Sci USA. 2009;106:15394–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
11.Gatehouse LN, Sutherland P, Forgie SA, Kaji R, Christeller JT. Molecular and histological characterization of primary (Betaproteobacteria) and secondary (Gammaproteobacteria) endosymbionts of three mealybug species. Appl Environ Microbiol. 2012;78:1187–97.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
12.Bennett GM, Moran NA. Small, smaller, smallest: the origins and evolution of ancient dual symbioses in a phloem-feeding insect. Genome Biol Evol. 2013;5:1675–88.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
13.Bressan A, Mulligan KL. Localization and morphological variation of three bacteriome-inhabiting symbionts within a planthopper of the genus Oliarus (Hemiptera: Cixiidae): Bacteriome-inhabiting symbionts in Oliarus filicicola. Environ Microbiol Rep. 2013;5:499–505.PubMed 
Article 
PubMed Central 

Google Scholar 
14.Bennett GM, Mao M. Comparative genomics of a quadripartite symbiosis in a planthopper host reveals the origins and rearranged nutritional responsibilities of anciently diverged bacterial lineages. Environ Microbiol. 2018;20:4461–72.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
15.von Dohlen CD, Spaulding U, Patch KB, Weglarz KM, Foottit RG, Havill NP, et al. Dynamic acquisition and loss of dual-obligate symbionts in the plant-sap-feeding Adelgidae (Hemiptera: Sternorrhyncha: Aphidoidea). Front Microbiol. 2017;8:1037.Article 

Google Scholar 
16.Mao M, Yang X, Poff K, Bennett G. Comparative genomics of the dual-obligate symbionts from the treehopper, Entylia carinata (Hemiptera: Membracidae), provide insight into the origins and evolution of an ancient symbiosis. Genome Biol Evol. 2017;9:1803–15.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
17.McCutcheon JP, Moran NA. Functional convergence in reduced genomes of bacterial symbionts spanning 200 my of evolution. Genome Biol Evol. 2010;2:708–18.PubMed 
PubMed Central 
Article 

Google Scholar 
18.McCutcheon JP, von Dohlen CD. An interdependent metabolic patchwork in the nested symbiosis of mealybugs. Curr Biol. 2011;21:1366–72.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
19.Sloan DB, Moran NA. Genome reduction and co-evolution between the primary and secondary bacterial symbionts of psyllids. Mol Biol Evol. 2012;29:3781–92.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
20.Hall AAG, Morrow JL, Fromont C, Steinbauer MJ, Taylor GS, Johnson SN, et al. Codivergence of the primary bacterial endosymbiont of psyllids versus host switches and replacement of their secondary bacterial endosymbionts. Environ Microbiol. 2016;18:2591–603.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
21.Tamas I, Klasson L, Canbäck B, Näslund AK, Eriksson A-S, Wernegreen JJ, et al. 50 million years of genomic stasis in endosymbiotic bacteria. Science 2002;296:2376–9.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
22.Shigenobu S, Watanabe H, Hattori M, Sakaki Y, Ishikawa H. Genome sequence of the endocellular bacterial symbiont of aphids Buchnera sp. APS. Nature 2000;407:81–6.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
23.Moran NA, Tran P, Gerardo NM. Symbiosis and insect diversification: an ancient symbiont of sap-feeding insects from the bacterial phylum Bacteroidetes. Appl Environ Microbiol. 2005;71:8802–10.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
24.Gruwell ME, Hardy NB, Gullan PJ, Dittmar K. Evolutionary relationships among primary endosymbionts of the mealybug subfamily Phenacoccinae (Hemiptera: Coccoidea: Pseudococcidae). Appl Environ Microbiol. 2010;76:7521–5.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
25.Koga R, Moran NA. Swapping symbionts in spittlebugs: evolutionary replacement of a reduced genome symbiont. ISME J. 2014;8:1237–46.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
26.Mao M, Bennett GM. Symbiont replacements reset the co-evolutionary relationship between insects and their heritable bacteria. ISME J. 2020;14:1384–95.PubMed 
PubMed Central 
Article 

Google Scholar 
27.Braendle C, Miura T, Bickel R, Shingleton AW, Kambhampati S, Stern DL. Developmental origin and evolution of bacteriocytes in the aphid–Buchnera symbiosis. PLoS Biol. 2003;1:e21.PubMed 
PubMed Central 
Article 

Google Scholar 
28.Weglarz KM, Havill NP, Burke GR, von Dohlen CD. Partnering with a pest: genomes of hemlock woolly adelgid symbionts reveal atypical nutritional provisioning patterns in dual-obligate bacteria. Genome Biol Evol. 2018;10:1607–21.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
29.Toenshoff ER, Penz T, Narzt T, Collingro A, Schmitz-Esser S, Pfeiffer S, et al. Bacteriocyte-associated gammaproteobacterial symbionts of the Adelges nordmannianae/piceae complex (Hemiptera: Adelgidae). ISME J 2012;6:384–96.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
30.Toenshoff ER, Gruber D, Horn M. Co-evolution and symbiont replacement shaped the symbiosis between adelgids (Hemiptera: Adelgidae) and their bacterial symbionts. Environ Microbiol. 2012;14:1284–95.CAS 
PubMed 
Article 

Google Scholar 
31.Toenshoff ER, Szabó G, Gruber D, Horn M. The pine bark adelgid, Pineus strobi, contains two novel bacteriocyte-associated gammaproteobacterial symbionts. Appl Environ Microbiol. 2014;80:878–85.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
32.von Dohlen CD, Spaulding U, Shields K, Havill NP, Rosa C, Hoover K. Diversity of proteobacterial endosymbionts in hemlock woolly adelgid (Adelges tsugae) (Hemiptera: Adelgidae) from its native and introduced range. Environ Microbiol. 2013;15:2043–62.Article 
CAS 

Google Scholar 
33.Havelka J, Danilov J, Rakauskas R. Relationships between aphid species of the family Adelgidae (Hemiptera Adelgoidea) and their endosymbiotic bacteria: a case study in Lithuania. Bull Insectology. 2021;74:1–10.
Google Scholar 
34.Favret C, Havill NP, Miller GL, Sano M, Victor B. Catalog of the adelgids of the world (Hemiptera, Adelgidae). Zookeys 2015;534:35–54.Article 

Google Scholar 
35.Blackman RL, Eastop VF Aphids on the world’s trees: an identification and information guide. 1994. CAB International.36.Havill NP, Foottit RG. Biology and evolution of Adelgidae. Ann Rev Ento. 2007;52:325–49.CAS 
Article 

Google Scholar 
37.Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 2014;30:2114–20.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
38.Zhang J, Kobert K, Flouri T, Stamatakis A. PEAR: a fast and accurate Illumina paired-end read mergeR. Bioinformatics 2014;30:614–20.CAS 
PubMed 
Article 

Google Scholar 
39.Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J Comp Bio. 2012;19:455–77.CAS 
Article 

Google Scholar 
40.Walker BJ, Abeel T, Shea T, Priest M, Abouelliel A, Sakthikumar S, et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS ONE. 2014;9:e112963.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
41.Kolmogorov M, Yuan J, Lin Y, Pevzner PA. Assembly of long, error-prone reads using repeat graphs. Nat Biotechnol. 2019;37:540–6.CAS 
PubMed 
Article 

Google Scholar 
42.Laetsch DR, Blaxter ML. BlobTools: Interrogation of genome assemblies. F1000Research. 2017;6:1287.Article 

Google Scholar 
43.Boetzer M, Henkel CV, Jansen HJ, Butler D, Pirovano W. Scaffolding pre-assembled contigs using SSPACE. Bioinformatics 2011;27:578–9.CAS 
PubMed 
Article 

Google Scholar 
44.Chu C, Li X, Wu Y. GAPPadder: a sensitive approach for closing gaps on draft genomes with short sequence reads. BMC Genomics. 2019;20:426.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
45.Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics 2014;30:2068–9.CAS 
PubMed 
Article 

Google Scholar 
46.Varani AM, Siguier P, Gourbeyre E, Charneau V, Chandler M. ISsaga is an ensemble of web-based methods for high throughput identification and semi-automatic annotation of insertion sequences in prokaryotic genomes. Genome Biol. 2011;12:R30.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
47.Karp PD, Billington R, Caspi R, Fulcher CA, Latendresse M, Kothari A, et al. The BioCyc collection of microbial genomes and metabolic pathways. Brief Bioinform. 2019;20:1085–93.CAS 
PubMed 
Article 

Google Scholar 
48.Karp PD, Ong WK, Paley S, Billington R, Caspi R, Fulcher C, et al. The EcoCyc database. EcoSal Plus. 2018;8:10.1128.Article 

Google Scholar 
49.Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007;35:W182–5.PubMed 
PubMed Central 
Article 

Google Scholar 
50.Huerta-Cepas J, Forslund K, Coelho LP, Szklarczyk D, Jensen LJ, von Mering C, et al. Fast genome-wide functional annotation through orthology assignment by eggNOG-Mapper. Mol Biol Evol. 2017;34:2115–22.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
51.Tatusov RL, Galperin MY, Natale DA, Koonin EV. The COG database: a tool for genome-scale analysis of protein functions and evolution. Nucleic Acids Res. 2000;28:33–36.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
52.Wang Y, Tang H, DeBarry JD, Tan X, Li J, Wang X, et al. MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity. Nucleic Acids Res. 2012;40:e49–e49.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Xu L, Dong Z, Fang L, Luo Y, Wei Z, Guo H, et al. OrthoVenn2: a web server for whole-genome comparison and annotation of orthologous clusters across multiple species. Nucleic Acids Res. 2019;47:W52–W58.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
54.Xu Y, Bi C, Wu G, Wei S, Dai X, Yin T, et al. VGSC: a web-based vector graph toolkit of genome synteny and collinearity. Biomed Res Int. 2016;2016:7823429.PubMed 
PubMed Central 

Google Scholar 
55.Adeolu M, Alnajar S, Naushad S, S Gupta R. Genome-based phylogeny and taxonomy of the ‘Enterobacteriales’: proposal for Enterobacterales ord. nov. divided into the families Enterobacteriaceae, Erwiniaceae fam. nov., Pectobacteriaceae fam. nov., Yersiniaceae fam. nov., Hafniaceae fam. nov., Morganellaceae fam. nov., and Budviciaceae fam. nov. Int J Syst Evol Microbiol. 2016;66:5575–99.CAS 
PubMed 
Article 

Google Scholar 
56.Guy L. phyloSkeleton: taxon selection, data retrieval and marker identification for phylogenomics. Bioinformatics 2017;33:1230–2.CAS 
PubMed 
PubMed Central 

Google Scholar 
57.Eddy SR. Accelerated profile HMM searches. PLoS Comput Biol. 2011;7:e1002195.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
58.Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30:772–80.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
59.Capella-Gutiérrez S, Silla-Martínez JM, Gabaldón T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics 2009;25:1972–3.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
60.Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 2014;30:1312–3.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
61.Lartillot N, Rodrigue N, Stubbs D, Richer J. PhyloBayes MPI: phylogenetic reconstruction with infinite mixtures of profiles in a parallel environment. Syst Biol. 2013;62:611–5.CAS 
PubMed 
Article 

Google Scholar 
62.Husník F, Chrudimský T, Hypša V. Multiple origins of endosymbiosis within the Enterobacteriaceae (γ-Proteobacteria): convergence of complex phylogenetic approaches. BMC Biology. 2011;9:1–17.Article 

Google Scholar 
63.Emms DM, Kelly S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol. 2019;20:238.PubMed 
PubMed Central 
Article 

Google Scholar 
64.Hesse C, Schulz F, Bull CT, Shaffer BT, Yan Q, Shapiro N, et al. Genome-based evolutionary history of Pseudomonas spp. Environ Microbiol. 2018;20:2142–59.CAS 
PubMed 
Article 

Google Scholar 
65.Burke GR, Normark BB, Favret C, Moran NA. Evolution and diversity of facultative symbionts from the aphid subfamily Lachninae. Appl Environ Microbiol. 2009;75:5328–35.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
66.Manzano‐Marín A, Szabó G, Simon J, Horn M, Latorre A. Happens in the best of subfamilies: establishment and repeated replacements of co‐obligate secondary endosymbionts within Lachninae aphids: co-obligate endosymbiont dynamics in the Lachninae. Environ Microbiol. 2017;19:393–408.PubMed 
Article 
CAS 

Google Scholar 
67.Castresana J. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol Biol Evol. 2000;17:540–52.CAS 
PubMed 
Article 

Google Scholar 
68.ggplot2. Create elegant data visualisations using the grammar of graphics. https://ggplot2.tidyverse.org/. Accessed Apr 2021.69.Manzano-Marín A, Oceguera-Figueroa A, Latorre A, Jiménez-García LF, Moya A. Solving a bloody mess: B-vitamin independent metabolic convergence among gammaproteobacterial obligate endosymbionts from blood-feeding arthropods and the leech Haementeria officinalis. Genome Biol Evol. 2015;7:2871–84.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
70.Janda JM, Abbott SL. The genus Hafnia: from soup to nuts. Clin Microbiol Rev. 2006;19:12–18.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
71.Szabó G, Schulz F, Manzano-Marín A, Toenshoff ER, Horn M Evolutionary recent dual obligatory symbiosis among adelgids indicates a transition between fungus and insect associated lifestyles. bioRxiv. 2020; e-pub ahead of print 16 October 2020; https://doi.org/10.1101/2020.10.16.342642.72.Wilson ACC, Ashton PD, Calevro F, Charles H, Colella S, Febvay G, et al. Genomic insight into the amino acid relations of the pea aphid, Acyrthosiphon pisum, with its symbiotic bacterium Buchnera aphidicola. Insect Mol Biol. 2010;19:249–58.CAS 
PubMed 
Article 

Google Scholar 
73.Sloan DB, Nakabachi A, Richards S, Qu J, Murali SC, Gibbs RA, et al. Parallel histories of horizontal gene transfer facilitated extreme reduction of endosymbiont genomes in sap-feeding insects. Mol Biol Evol. 2014;31:857–71.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
74.Hansen AK, Moran NA. The impact of microbial symbionts on host plant utilization by herbivorous insects. Mol Ecol. 2014;23:1473–96.PubMed 
Article 

Google Scholar 
75.Manzano-Marı́n A, Coeur d’acier A, Clamens A-L, Orvain C, Cruaud C, Barbe V, et al. Serial horizontal transfer of vitamin-biosynthetic genes enables the establishment of new nutritional symbionts in aphids’ di-symbiotic systems. ISME J. 2020;14:259–73.Article 
CAS 

Google Scholar 
76.Lo W-S, Huang Y-Y, Kuo C-H. Winding paths to simplicity: genome evolution in facultative insect symbionts. FEMS Microbiol Rev. 2016;40:855–74.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
77.Toh H, Weiss BL, Perkin SAH, Yamashita A, Oshima K, Hattori M, et al. Massive genome erosion and functional adaptations provide insights into the symbiotic lifestyle of Sodalis glossinidius in the tsetse host. Genome Res. 2006;16:149–56.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
78.Cole ST, Eiglmeier K, Parkhill J, James KD, Thomson NR, Wheeler PR, et al. Massive gene decay in the leprosy bacillus. Nature 2001;409:1007–11.CAS 
PubMed 
Article 

Google Scholar 
79.Moran NA, Bennett GM. The tiniest tiny genomes. Annu Rev Microbiol. 2014;68:195–215.CAS 
PubMed 
Article 

Google Scholar 
80.Bennett GM, McCutcheon JP, MacDonald BR, Romanovicz D, Moran NA. Differential genome evolution between companion symbionts in an insect-bacterial symbiosis. mBio 2014;5:e01697–14.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
81.Degnan PH, Ochman H, Moran NA. Sequence conservation and functional constraint on intergenic spacers in reduced genomes of the obligate symbiont Buchnera. PLoS Genet. 2011;7:e1002252.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
82.Van Leuven JT, Meister RC, Simon C, McCutcheon JP. Sympatric speciation in a bacterial endosymbiont results in two genomes with the functionality of one. Cell 2014;158:1270–80.PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
83.Gomez-Valero L. The evolutionary fate of nonfunctional DNA in the bacterial endosymbiont Buchnera aphidicola. Mol Biol Evol. 2004;21:2172–81.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
84.Manzano-Marı́n A, Coeur d’acier A, Clamens A-L, Orvain C, Cruaud C, Barbe V, et al. A freeloader? The highly eroded yet large genome of the Serratia symbiotica symbiont of Cinara strobi. Genome Biol Evol. 2018;10:2178–89.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
85.Santos-Garcia D, Silva FJ, Morin S, Dettner K, Kuechler SM. The all-rounder Sodalis: a new bacteriome-associated endosymbiont of the lygaeoid bug Henestaris halophilus (Heteroptera: Henestarinae) and a critical examination of its evolution. Genome Biol Evol. 2017;9:2893–910.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
86.Havill NP, Foottit RG, von Dohlen CD. Evolution of host specialization in the Adelgidae (Insecta: Hemiptera) inferred from molecular phylogenetics. Mol Phylogenet. 2007;44:357–70.CAS 
Article 

Google Scholar 
87.Manzano-Marı́n A, Latorre A. Snapshots of a shrinking partner: genome reduction in Serratia symbiotica. Sci Rep. 2016;6:32590.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
88.Monnin D, Jackson R, Kiers ET, Bunker M, Ellers J, Henry LM. Parallel evolution in the integration of a co-obligate aphid symbiosis. Curr Biol. 2020;30:1949–57. e6CAS 
PubMed 
Article 

Google Scholar 
89.Husnik F, McCutcheon JP. Repeated replacement of an intrabacterial symbiont in the tripartite nested mealybug symbiosis. Proc Natl Acad Sci USA. 2016;113:e5416–24.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
90.Moran NA, McCutcheon JP, Nakabachi A. Genomics and evolution of heritable bacterial symbionts. Annu Rev Genet. 2008;42:165–90.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
91.Degnan PH, Leonardo TE, Cass BN, Hurwitz B, Stern D, Gibbs RA, et al. Dynamics of genome evolution in facultative symbionts of aphids. Environ Microbiol. 2010;12:2060–9.CAS 
PubMed 
PubMed Central 

Google Scholar 
92.Burke GR, Moran NA. Massive genomic decay in Serratia symbiotica, a recently evolved symbiont of aphids. Genome Biol Evol. 2011;3:195–208.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
93.Munson MA, Baumann P, Clark MA, Baumann L, Moran NA, Voegtlin DJ, et al. Evidence for the establishment of aphid-eubacterium endosymbiosis in an ancestor of four aphid families. J Bacteriol. 1991;173:6321–4.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
94.Moran NA, Munson MA, Baumann P, Ishikawa H. A molecular clock in endosymbiotic bacteria is calibrated using the insect hosts. Proc R Soc B 1993;253:167–71.Article 

Google Scholar 
95.Kuechler SM, Gibbs G, Burckhardt D, Dettner K, Hartung V. Diversity of bacterial endosymbionts and bacteria-host co-evolution in Gondwanan relict moss bugs (Hemiptera: Coleorrhyncha: Peloridiidae). Environ Microbiol. 2013;15:2031–42.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
96.Thao ML, Moran NA, Abbot P, Brennan EB, Burckhardt DH, Baumann P. Cospeciation of psyllids and their primary prokaryotic endosymbionts. Appl Environ Microbiol. 2000;66:2898–905.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
97.Thao ML, Baumann P. Evolutionary relationships of primary prokaryotic endosymbionts of whiteflies and their hosts. Appl Environ Microbiol. 2004;70:3401–6.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
98.Meseguer AS, Manzano-Marín A, Coeur d’Acier A, Clamens AL, Godefroid M, Jousselin E. Buchnera has changed flatmate but the repeated replacement of co-obligate symbionts is not associated with the ecological expansions of their aphid hosts. Mol Ecol. 2017;26:2363–78.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
99.McCutcheon JP, Moran NA. Parallel genomic evolution and metabolic interdependence in an ancient symbiosis. Proc Natl Acad Sci USA. 2007;104:19392–7.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
100.Rao Q, Rollat-Farnier PA, Zhu DT, Santos-Garcia D, Silva FJ, Moya A, et al. Genome reduction and potential metabolic complementation of the dual endosymbionts in the whitefly Bemisia tabaci. BMC Genomics. 2015;16:226.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
101.Rosenblueth M, Sayavedra L, Sámano-Sánchez H, Roth A, Martínez-Romero E. Evolutionary relationships of flavobacterial and enterobacterial endosymbionts with their scale insect hosts (Hemiptera: Coccoidea). J Evol Biol. 2012;25:2357–68.PubMed 
Article 
PubMed Central 

Google Scholar 
102.Michalik K, Szklarzewicz T, Kalandyk-Kołodziejczyk M, Jankowska W, Michalik A. Bacteria belonging to the genus Burkholderia are obligatory symbionts of the eriococcids Acanthococcus aceris Signoret, 1875 and Gossyparia spuria (Modeer, 1778) (Insecta, Hemiptera, Coccoidea). Arthropod Struct Dev. 2016;45:265–72.PubMed 
Article 
PubMed Central 

Google Scholar 
103.Van Ham RC, Kamerbeek J, Palacios C, Rausell C, Abascal F, Bastolla U, et al. Reductive genome evolution in Buchnera aphidicola. Proc Natl Acad Sci USA. 2003;100:581–6.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
104.Vogel KJ, Moran NA. Effect of host genotype on symbiont titer in the aphid-Buchnera symbiosis. Insects 2011;2:423–34.PubMed 
PubMed Central 
Article 

Google Scholar 
105.Bennett GM, McCutcheon JP, McDonald BR, Moran NA. Lineage-specific patterns of genome deterioration in obligate symbionts of sharpshooter leafhoppers. Genome Biol Evol. 2015;8:296–301.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
106.Havill NP, Griffin BP, Andersen JC, Foottit RG, Justesen MJ, Caccone A, et al. Species delimitation and invasion history of the balsam woolly adelgid, Adelges (Dreyfusia) piceae (Hemiptera: Aphidoidea: Adelgidae), species complex. Syst Entomol. 2021;46:186–204.Article 

Google Scholar  More

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Google Scholar  More

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

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

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

(i)

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

(ii)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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