Philippa Kaur

Birks, H. J. B. Strengths and weaknesses of quantitative climate reconstructions based on late-quaternary biological proxies. Open Ecol. J. 3, 68–110 (2011).Article 

Google Scholar 
Chevalier, M. et al. Pollen-based climate reconstruction techniques for late Quaternary studies. Earth-Sci. Rev. 210, 103384 (2020).Article 

Google Scholar 
Bartlein, P. J. et al. Pollen-based continental climate reconstructions at 6 and 21 ka: a global synthesis. Clim. Dyn. 37, 775–802 (2011).Article 

Google Scholar 
Brierley, C. M. et al. Large-scale features and evaluation of the PMIP4-CMIP6 midHolocene simulations. Clim. Past 16, 1847–1872 (2020).Kageyama, M. et al. The PMIP4 Last Glacial Maximum experiments: preliminary results and comparison with the PMIP3 simulations. Clim. Past 17, 1065–1089 (2021).Article 

Google Scholar 
Harrison, S. BIOME 6000 DB classified plotfile version 1. https://doi.org/10.17864/1947.99. (2017).Loarie, S. R. et al. The velocity of climate change. Nature 462, 1052–1055 (2009).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Svenning, J. C. & Sandel, B. Disequilibrium vegetation dynamics under future climate change. Am. J. Bot. 100, 1266–1286 (2013).PubMed 
Article 

Google Scholar 
Neilson, R. P. et al. Forecasting regional to global plant migration in response to climate change. BioScience 55 https://academic.oup.com/bioscience/article/55/9/749/285963 (2005).Normand, S. et al. Postglacial migration supplements climate in determining plant species ranges in Europe. Proc. R. Soc. B: Biol. Sci. 278, 3644–3653 (2011).Article 

Google Scholar 
Seltzer, A. M. et al. Widespread six degrees Celsius cooling on land during the Last Glacial Maximum. Nature 593, 228–232 (2021).Tierney, J. E. et al. Glacial cooling and climate sensitivity revisited. Nature 584, 569 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Ray, N. & Adams, J. M. A GIS-based Vegetation Map of the World at the Last Glacial Maximum (25,000-15,000 BP). Internet Archaeol. 11, https://doi.org/10.11141/ia.11.2 (2001).Birks, H. J. B. & Willis, K. J. Alpines, trees, and refugia in Europe. Plant Ecol. Divers. 1, 147–160 (2008).Article 

Google Scholar 
Roberts, D. R. & Hamann, A. Glacial refugia and modern genetic diversity of 22 western North American tree species. Proc. R. Soc. B: Biol. Sci. 282, 20142903 (2015).Clark, J. S. Why trees migrate so fast: Confronting theory with dispersal biology and the paleorecord. Am. Nat. 152, 204–224 (1998).CAS 
PubMed 
Article 

Google Scholar 
Jackson, S. & Overpeck, J. Responses of plant populations and communities to environmental changes of the late Quaternary. Paleobiology 26, 194–220 (2000).Article 

Google Scholar 
Williams, J. W., Ordonez, A. & Svenning, J.-C. A unifying framework for studying and managing climate-driven rates of ecological change. Nat. Ecol. Evol. 5, 17–26 (2021).Harrison, S. P. & Goñi, M. F. S. Global patterns of vegetation response to millennial-scale variability and rapid climate change during the last glacial period. Quat. Sci. Rev. 29, 2957–2980 (2010).ADS 
Article 

Google Scholar 
Williams, J. W., Post, D. M., Cwynar, L. C., Lotter, A. F. & Levesque, A. J. Rapid and widespread vegetation responses to past climate change in the North Atlantic region. Geology 30, 971–974 (2002).ADS 
CAS 
Article 

Google Scholar 
Giesecke, T., Brewer, S., Finsinger, W., Leydet, M. & Bradshaw, R. H. W. Patterns and dynamics of European vegetation change over the last 15,000 years. J. Biogeogr. 44, 1441–1456 (2017).Article 

Google Scholar 
Ordonez, A. & Williams, J. W. Climatic and biotic velocities for woody taxa distributions over the last 16 000 years in eastern North America. Ecol. Lett. 16, 773–781 (2013).PubMed 
Article 

Google Scholar 
Svenning, J.-C. & Skov, F. Limited filling of the potential range in European tree species. Ecol. Lett. 7, 565–573 (2004).Article 

Google Scholar 
Talluto, M. V., Boulangeat, I., Vissault, S., Thuiller, W. & Gravel, D. Extinction debt and colonization credit delay range shifts of eastern North American trees. Nat. Ecol. Evol. 1, 1–6 (2017).Article 

Google Scholar 
Pearson, R. G. & Dawson, T. P. Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful? Glob. Ecol. Biogeogr. 12, 361–371 (2003).Article 

Google Scholar 
Webb, T. Is vegetation in equilibrium with climate? How to interpret late-Quaternary pollen data. Vegetatio 67, 75–91 (1986).Article 

Google Scholar 
Jackson, S. T. & Williams, J. W. Modern analogs in quaternary paleoecology: Here today, gone yesterday, gone tomorrow? Annu. Rev. Earth Planet. Sci. 32, 495–537 (2004).ADS 
CAS 
Article 

Google Scholar 
Cao, X., Tian, F., Dallmeyer, A. & Herzschuh, U. Northern Hemisphere biome changes ( >30°N) since 40 cal ka BP and their driving factors inferred from model-data comparisons. Quat. Sci. Rev. 220, 291–309 (2019).ADS 
Article 

Google Scholar 
He, F. Simulating transient climate evolution of the last deglaciation with CCSM3 Dissertation at the University of Wisconsin – Madison (2011).Osman, M. B. et al. Globally resolved surface temperatures since the Last Glacial Maximum. Nature 599, 239–244 (2021).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Shakun, J. D. et al. Global warming preceded by increasing carbon dioxide concentrations during the last deglaciation. Nature 484, 49–54 (2012).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Alley, R. B. The Younger Dryas cold interval as viewed from central Greenland. in Quaternary Science Reviews vol. 19 213–226 (Pergamon, 2000).He, C. et al. Hydroclimate footprint of pan-Asian monsoon water isotope during the last deglaciation. Sci. Adv. 7, eabe2611 (2021).ADS 
PubMed 
Article 

Google Scholar 
Reick, C. H., Raddatz, T., Brovkin, V. & Gayler, V. Representation of natural and anthropogenic land cover change in MPI-ESM. J. Adv. Modeling Earth Syst. 5, 459–482 (2013).Prentice, I. C., Guiot, J., Huntley, B., Jolly, D. & Cheddadi, R. Reconstructing biomes from palaeoecological data: a general method and its application to European pollen data at 0 and 6 ka. Clim. Dyn. 12, 185–194 (1996).Article 

Google Scholar 
Dallmeyer, A., Claussen, M. & Brovkin, V. Harmonising plant functional type distributions for evaluating Earth system models. Clim 15, 335–366 (2019).
Google Scholar 
Ni, J., Cao, X., Jeltsch, F. & Herzschuh, U. Biome distribution over the last 22,000 yr in China. Palaeogeogr. Palaeoclimatol. Palaeoecol. 409, 33–47 (2014).Article 

Google Scholar 
Williams, J. W. & Jackson, S. T. Novel climates, no-analog communities, and ecological surprises. Front. Ecol. Environ. 5, 475–482 (2007).Article 

Google Scholar 
Sobol, M. K., Scott, L. & Finkelstein, S. A. Reconstructing past biomes states using machine learning and modern pollen assemblages: a case study from Southern Africa. Quat. Sci. Rev. 212, 1–17 (2019).ADS 
Article 

Google Scholar 
Marinova, E. et al. Pollen‐derived biomes in the Eastern Mediterranean–Black Sea–Caspian‐Corridor. J. Biogeogr. 45, 484–499 (2018).Article 

Google Scholar 
Cao, X. et al. Pollen-based quantitative land-cover reconstruction for northern Asia covering the last 40 ka cal BP. Clim. Past 15, 1503–1536 (2019).Article 

Google Scholar 
Geng, R. et al. Modern pollen assemblages from lake sediments and soil in East Siberia and relative pollen productivity estimates for major taxa. Front. Ecol. Evol. 10, 508 (2022).Article 

Google Scholar 
Cao, X. et al. A taxonomically harmonized and temporally standardized fossil pollen dataset from Siberia covering the last 40 kyr. Earth Syst. Sci. Data 12, 119–135 (2020).ADS 
Article 

Google Scholar 
Sugita, S. Theory of quantitative reconstruction of vegetation I: pollen from large sites REVEALS regional vegetation composition. Holocene 17, 229–241 (2007).ADS 
Article 

Google Scholar 
Githumbi, E. et al. European pollen-based REVEALS land-cover reconstructions for the Holocene: Methodology, mapping and potentials. Earth Syst. Sci. Data 14, 1581–1619 (2022).ADS 
Article 

Google Scholar 
Snell, R. S. et al. Using dynamic vegetation models to simulate plant range shifts. Ecography 37, 1184–1197 (2014).Article 

Google Scholar 
Bullock, J. M. et al. Modelling spread of British wind-dispersed plants under future wind speeds in a changing climate. J. Ecol. 100, 104–115 (2012).Article 

Google Scholar 
Svenning, J. C., Normand, S. & Skov, F. Postglacial dispersal limitation of widespread forest plant species in nemoral Europe. Ecography 31, 316–326 (2008).Article 

Google Scholar 
Herzschuh, U. et al. Glacial legacies on interglacial vegetation at the Pliocene-Pleistocene transition in NE Asia. Nat. Commun. 7, 1–11 (2016).Article 

Google Scholar 
Herzschuh, U. Legacy of the Last Glacial on the present‐day distribution of deciduous versus evergreen boreal forests. Glob. Ecol. Biogeogr. 29, 198–206 (2020).Article 

Google Scholar 
Väliranta, M. et al. Plant macrofossil evidence for an early onset of the Holocene summer thermal maximum in northernmost Europe. Nat. Commun. 6, 1–8 (2015).Article 

Google Scholar 
Schulte, L., Li, C., Livsovski, S. & Herzschuh, U. Forest-permafrost feedbacks and glacial refugia help explain the unequal distribution of larch across continents. J. Biogeogr. 9, 0305–0270 (2022).
Google Scholar 
Davis, M. B., Shaw, R. G. & Etterson, J. R. Evolutionary responses to changing climate. Ecology 86, 1704–1714 (2005).Article 

Google Scholar 
Urban, M. C., Tewksbury, J. J. & Sheldon, K. S. On a collision course: competition and dispersal differences create no-analogue communities and cause extinctions during climate change. Proc. R. Soc. B Biol. Sci. 279, 2072–2080 (2012).Article 

Google Scholar 
Pennington, W. Lags in adjustment of vegetation to climate caused by the pace of soil development. Evidence from Britain. Vegetatio 67, 105–118 (1986).Article 

Google Scholar 
MacDonald, G. M., Kremenetski, K. V. & Beilman, D. W. Climate change and the northern Russian treeline zone. Philos. Trans. R. Soc. B: Biol. Sci. 363, 2285–2299 (2008).CAS 
Article 

Google Scholar 
Prentice, I. C., Bartlein, P. J. & Webb, T. Vegetation and climate change in eastern North America since the last glacial maximum. Ecology 72, 2038–2056 (1991).Article 

Google Scholar 
Cao, X. Y., Herzschuh, U., Telford, R. J. & Ni, J. A modern pollen-climate dataset from China and Mongolia: Assessing its potential for climate reconstruction. Rev. Palaeobot. Palynol. 211, 87–96 (2014).Article 

Google Scholar 
Leroy, S. A. G., Arpe, K., Mikolajewicz, U. & Wu, J. Climate simulations and pollen data reveal the distribution and connectivity of temperate tree populations in eastern Asia during the Last Glacial Maximum. Clim 16, 2039–2054 (2020).
Google Scholar 
Kaufman, D. et al. A global database of Holocene paleotemperature records. Sci. Data 7, 115 (2020).Mottl, O. et al. Global acceleration in rates of vegetation change over the past 18,000 years. Science 372, 860–864 (2021).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Nolan, C. et al. Past and future global transformation of terrestrial ecosystems under climate change. Science 361, 920–923 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Mauritsen, T. et al. Developments in the MPI-M Earth System Model version 1.2 (MPI-ESM1.2) and Its Response to Increasing CO2. J. Adv. Model. Earth Syst. 11, 998–1038 (2019).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Reick, C. et al. JSBACH 3—The land component of the MPI Earth System Model: documentation of version 3.2. Hamburg: MPI für Meteorologie. Berichte zur Erdsystemforsch. (2021).Brovkin, V., Raddatz, T., Reick, C. H., Claussen, M. & Gayler, V. Global biogeophysical interactions between forest and climate. Geophys. Res. Lett. 36, L07405 (2009).ADS 
Article 

Google Scholar 
Berger, A. L. Long-term variations of daily insolation and Quaternary climatic changes. J. Atmos. Sci. 35, 2361–2367 (1978).ADS 
Article 

Google Scholar 
Köhler, P., Nehrbass-Ahles, C., Schmitt, J., Stocker, T. F. & Fischer, H. A 156 kyr smoothed history of the atmospheric greenhouse gases CO2, CH4, and N2O and their radiative forcing. Earth Syst. Sci. Data 9, 363–387 (2017).ADS 
Article 

Google Scholar 
Tarasov, L., Dyke, A. S., Neal, R. M. & Peltier, W. R. A data-calibrated distribution of deglacial chronologies for the North American ice complex from glaciological modeling. Earth Planet. Sci. Lett. 315–316, 30–40 (2012).ADS 
Article 

Google Scholar 
Loana Meccia, V. & Mikolajewicz, U. Interactive ocean bathymetry and coastlines for simulating the last deglaciation with the Max Planck Institute Earth System Model (MPI-ESM-v1.2). Geosci. Model Dev. 11, 4677–4692 (2018).ADS 
Article 

Google Scholar 
Riddick, T., Brovkin, V., Hagemann, S. & Mikolajewicz, U. Dynamic hydrological discharge modelling for coupled climate model simulations of the last glacial cycle: the MPI-DynamicHD model version 3.0. Geosci. Model Dev. 11, 4291–4316 (2018).ADS 
Article 

Google Scholar 
Kapsch, M., Mikolajewicz, U., Ziemen, F. & Schannwell, C. Ocean response in transient simulations of the last deglaciation dominated by underlying ice‐sheet reconstruction and method of meltwater distribution. Geophys. Res. Lett. 49, e2021GL096767 (2022).ADS 
Article 

Google Scholar 
Murton, J. B., Bateman, M. D., Dallimore, S. R., Teller, J. T. & Yang, Z. Identification of Younger Dryas outburst flood path from Lake Agassiz to the Arctic Ocean. Nature 464, 740–743 (2010).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Rehfeld, K., Marwan, N., Heitzig, J. & Kurths, J. Comparison of correlation analysis techniques for irregularly sampled time series. Nonlinear Process. Geophys. 18, 389–404 (2011).ADS 
Article 

Google Scholar 
Braconnot, P. et al. Evaluation of climate models using palaeoclimatic data. Nat. Clim. Change 2, 417–424 (2012).ADS 
Article 

Google Scholar 
Cao, X. Y., Ni, J., Herzschuh, U., Wang, Y. B. & Zhao, Y. A late Quaternary pollen dataset from eastern continental Asia for vegetation and climate reconstructions: set up and evaluation. Rev. Palaeobot. Palynol. 194, 21–37 (2013).Article 

Google Scholar 
Bigelow, N. H. et al. Climate change and Arctic ecosystems: 1. Vegetation changes north of 55°N between the last glacial maximum, mid-Holocene, and present. J. Geophys. Res. Atmos. 108, 8170 (2003).Ramankutty, N. & Foley, J. A. Estimating historical changes in global land cover: Croplands from 1700 to 1992. Glob. Biogeochem. Cycles 13, 997–1027 (1999).ADS 
CAS 
Article 

Google Scholar 
Berger, A. & Loutre, M. F. Insolation values for the climate of the last 10 million years. Quat. Sci. Rev. 10, 297–317 (1991).ADS 
Article 

Google Scholar 
Deplazes, G. et al. Links between tropical rainfall and North Atlantic climate during the last glacial period. Nat. Geosci. 6, 213–217 (2013).ADS 
CAS 
Article 

Google Scholar 
Wessel, P. et al. Generic mapping tools: improved version released. EOS Trans. AGU 94, 409–410 (2013).ADS 
Article 

Google Scholar  More

The pathosystemPlantago lanceolata L. is a perennial monoecious ribwort plantain that reproduces both clonally via the production side rosettes, and sexually via wind pollination. Seeds drop close to the mother plant and usually form a long-term seed bank47. Podospharea plantaginis (Castagne; U. Braun and S. Takamatsu) (Erysiphales, Ascomycota) is an obligate biotrophic powdery mildew that infects only P. lanceolata and requires living host tissue through its life cycle48. It completes its life cycle as localized lesions on host leaves, only the haustorial feeding roots penetrating the leaf tissue to feed nutrients from its host. Infection causes significant stress for host plant and may increase the host mortality31. The interaction between P. lanceolata and P. plantaginis is strain-specific, whereby the same host genotype may be susceptible to some pathogen genotypes while being resistant to others49. The putative resistance mechanism includes two steps. First, resistance occurs when the host plant first recognizes the attacking pathogen and blocks its growth. When the first step fails and infection takes place, the host may mitigate infection development. Both resistance traits vary among host genotypes49.Approximately 4000 P. lanceolata populations form a network covering an area of 50 × 70 km in the Åland Islands, SW of Finland. Disease incidence (0/1) in these populations has been recorded systematically every year in early September since 2001 by approximately 40 field assistants, who record the occurrence of the fungus P. plantaginis in the local P. lanceolata populations30. At this time, disease symptoms are conspicuous as infected plants are covered by white mycelia and conidia. The coverage (m2) of P. lanceolata in the meadows was recorded between 2001 and 2008 and is used as an estimate of host population size. In the field survey two technicians estimate Plantago population size by visually estimating how much ground/other vegetation P. lanceolata foliage covers (m2) in each meadow. The proportion of P. lanceolata plants in each population suffering from drought is also estimated annually in the survey. Data on average rainfall (mm) in July and August was estimated separately for each population using detailed radar-measured rainfall (obtained by Finnish Meteorological Institute) and it was available for years 2001–2008.Host population connectivity (SH)27 for each local population i was computed with the formula that takes into account the area of host coverage (m²) of all host populations surveyed, denoted with (Aj), and their spatial location compared to other host populations. We assume that the distribution of dispersal distances from a location are described by negative exponential distribution. Under this assumption, the following formula (1) quantifies for a focal population i, the effect of all other host populations, taking into account their population sizes and how strongly they are connected through immigration to it:$${S}_{i}^{H}=mathop{sum}limits_{jne i}{{{{{rm{e}}}}}}^{{-alpha d}_{{ij}}}sqrt{{A}_{j}}.$$
(1)
here, dij is the Euclidian distance between populations i and j and 1/α equals the mean dispersal distance, which was set to be two kilometres based on results from a previous study16.The annual survey data has demonstrated that P. plantaginis infects annually 2–16% of all host populations and persists as a highly dynamic metapopulation through extinctions and re-colonizations of local populations16. The number of host populations has remained relatively stable over the study period49. The first visible symptoms of P. plantaginis infection appear in late June as white-greyish lesions consisting of mycelium supporting the dispersal spores (conidia) that are carried by wind to the same or new host plants. Six to eight clonally produced generations follow one another in rapid succession, often leading to local epidemic with substantial proportion of the infected hosts by late summer within the host local population. Podosphaera plantaginis produces resting structures, chasmothecia, that appear towards the end of growing season in August–September31. Between 20% and 90% of the local pathogen populations go extinct during the winter, and thus the recolonization events play an important role in the persistence of the pathogen regionally16.Inoculation assay: Effect of connectivity and disease history on phenotypic disease resistanceHost and pathogen material for the experimentTo examine whether the diversity and level of resistance vary among host populations depending on their degree of connectivity (SH) and disease history, we selected 20 P. lanceolata populations for an inoculation assay. These populations occur in different locations in the host network, and were selected based on their connectivity values (S H of selected populations was 37–110 in isolated and 237–336 in highly connected category, Fig. 1). We did not include host populations in the intermediate connectivity category that was used in the population dynamic analyses in the inoculation assay due to logistic constraints. Podosphaera plantaginis is an obligate biotrophic pathogen that requires living host tissue throughout its life cycle, and obtaining sufficient inoculum for experiments is extremely time and space consuming. In both isolated and highly connected categories, half of the populations (IDs 193, 260, 311, 313, 337, 507, 1821, 1999, 2818 and 5206) were healthy during the study years 2001–2014, while half of the populations (IDs 271, 294, 309, 321, 490, 609, 1553, 1556, 1676 and 1847) were infected by P. plantaginis for several years during the same period. We collected P. lanceolata seeds from randomly selected ten individual plants around the patch area from each host population in August 2014.To acquire inoculum for the assay, we collected the pathogen strains as infected leaves, one leaf from ten plant individuals from four additional host populations (IDs 3301, 4684, 1784, and 3108) in August 2014. None of the pathogen populations were same as the sampled host populations and hence, the strains used in the assay all represent allopatric combinations. Both host and pathogen populations selected for the study were separated by at least two kilometres. The collected leaves supporting infection were placed in Petri dishes on moist filter paper and stored at room temperature until later use.Seeds from ten mother plants from each population were sown in 2:1 mixture of potting soil and sand, and grown in greenhouse conditions at 20 ± 2 °C (day) and 16 ± 2 °C (night) with 16:8 L:D photoperiod. Due to the low germination rate of collected seeds, population 260 (isolated and healthy population) was excluded from the study. Seedlings of ten different mother plants were randomly selected among the germinated plants for each population (n = 190), and grown in individual pots until the plants were eight weeks old.The pathogen strains were purified through three cycles of single colony inoculations and maintained on live, susceptible leaves on Petri dishes in a growth chamber 20 ± 2 °C with 16:8 L:D photoperiod. Every two weeks, the strains were transferred to fresh P. lanceolata leaves. Purified powdery mildew strains (M1–M4), one representing each allopatric population (3301, 4684, 1784 and 3108), were used for the inoculation assay. To produce enough sporulating fungal material, repeated cycles of inoculations were performed before the assay.Inoculation assay quantifying host resistance phenotypesIn order to study how the phenotypic resistance of hosts varies depending on population connectivity and infection history, we scored the resistance of 190 host genotypes, ten individuals from each study populations (n = 19), in an inoculation assay. Here, one detached leaf from each plant was exposed to a single pathogen strain (M1–M4) by brushing spores gently with a fine paintbrush onto the leaf. Leaves were placed on moist filter paper in Petri dishes and kept in a growth chamber at 20 ± 2 with a 16/8D photoperiod. All the inoculations were repeated on two individual Petri plates, leading to 760 host genotype—pathogen genotype combinations and a total of 1520 inoculations (19 populations * 10 plant genotypes * 4 pathogen strains * 2 replicates). We then observed and scored the pathogen infection on day 12 post inoculation, under dissecting microscope. The resulting plant phenotypic response was scored as 0 = susceptible (infection) when mycelium and conidia were observed on the leaf surface, and as 1 = resistance (no infection), when no developing lesions could be detected under a dissecting microscope. A genotype was defined resistant only if both inoculated replicates showed similar response (1), and susceptible if one or both replicates became infected (0).Statistical analysesBayesian spatio-temporal INLA model of the changes in host population sizeTo study how the pathogen infection influences on host population growth, we analyzed the relative change in host population size (m2) (defined as population size (t) − population size (t−1))/population size (t−1)) between consecutive years utilizing data from 2001 to 2008 in response to pathogen presence-absence status at t−1 (Supplementary Table 2). To assess whether this depends on host population connectivity, we estimated the separate effects of pathogen presence/absence in the previous year for connectivity categories—high-, low, and intermediate—that were based on the 0.2 and 0.8 quantiles of the host-connectivity values (Fig. 1A and Supplementary Figs. 1, 2). This allowed us to directly assess and compare the effect of the pathogen on host population growth in the extreme categories between isolated and highly connected host populations which were represented in the sampling for the inoculation study (Fig. 2).As covariates, we included the proportion (0–100%) of dry host plants measured each year within each local population as well as data on the amount of rainfall at the summer months (June, July, and August) obtained from the satellite images, as these were suggested be relevant for this pathosystem in an earlier analysis16. Observations where the change in host population size, or the host population coverage had absolute values larger than their 0.99 quantiles in the whole data, were regarded as outliers and omitted from the analysis. Before the analyses, all the continuous covariates were scaled and centred, and the categorical variables were transformed into binary variables.The relative changes in local host population size between consecutive years was analyzed by a Bayesian spatio-temporal statistical model that simultaneously considers the effects of a set of biologically meaningful predictors. The linear predictor thus consists of two parts (2,3):$$beta {X}_{t}+{z}_{t}{A}_{t}$$
(2)
where (beta) represents the correlation coefficients corresponding to the effects of environmental covariates, ({z}_{t}) corresponds to the spatiotemporal random effect, and ({X}_{t}) and ({A}_{t}) project these to the observation locations. For ({z}_{t}) we assume that the observations from a location in consecutive time points (t−1) and t are described by 1st order autoregressive process:$${z}_{t}=varphi {z}_{t-1}+{w}_{t}$$
(3)
where ({w}_{t}) corresponds to spatially structured zero-mean random noise, for which a Matern covariance function is assumed. Statistical inference then targets jointly the covariate effects (beta), the temporal autocorrelation (varphi), and the hyperparameters describing the spatial autocorrelation in wt. From these the overall variance, as well as spatial range—a distance after which spatial autocorrelation ceases to be significant—can be inferred (Supplementary Fig. 3). For more detailed description of the structure of the statistical model and how to do efficient inference with it using R-INLA, we refer to refs. 16,50.Identification of resistance phenotypesThe phenotype composition of each study population was defined by individual plant responses to the four pathogen strains, where each response could be “susceptible = 0” or “resistant = 1”. For example, a phenotype “1111” refers to a plant resistant to all four pathogen strains. The diversity of distinct resistance phenotypes within populations was estimated using the Shannon diversity index as implemented in the vegan software package51. The Shannon diversity index for all four study groups was then analyzed using a linear model with class predictors population type (well-connected or isolated), infection history (healthy or infected), and their interaction.Analysis of population connectivity and infection history effects on host resistanceTo test whether host population resistance varied depending on connectivity (SH) and infection history, we analyzed the inoculation responses (0 = susceptible, 1=resistant) of each host-pathogen combination by using a logit mixed-effect model in the lme4 package52. The model included the binomial dependent variable (resistance-susceptible; 1/0), and class predictors population type (well-connected or isolated), infection history (healthy or infected), mildew strain (M1, M2, M3, and M4) and their interactions. Plant individual and population were defined as random effects, with plant genotype (sample) hierarchically nested under population. Model fit was assessed using chi-square tests on the log-likelihood values to compare different models and significant interactions, and the best model was selected based on AIC-values. P-values for regression coefficients were obtained by using the car package53. We ran all the analyses in R software54.The metapopulation modelWe model the ecological and co-evolutionary dynamics of host and pathogen metapopulations to understand key features of the experimental system that impact on the qualitative patterns observed. The structure and parameters in our model are therefore not estimated using experimental data, but rather are chosen to cover a range of possibilities (e.g., low vs high transmission rates, variation in trade-off shapes for fitness costs). We construct the metapopulations in two stages to account for relatively well and poorly connected demes. All demes are identical in quality (i.e., no differences in intrinsic birth or death rates between demes) and only differ in their connectivity. Our metapopulation consists of an outer network of 20 demes, equally spaced around the unit square (0.2 units apart), and a 7×7 inner lattice of demes at a minimum distance of 0.2 units from the outer network (Fig. 3A), giving a total of 69 demes. Demes that are separated by a Euclidean distance of at most 0.2 are then connected to each other. This means that populations near the centre of the metapopulation are highly connected, while those on the boundary of the metapopulation are poorly connected. This also has the effect of making connections between well and poorly connected demes assortative (i.e., well/poorly connected demes tend to be connected to well/poorly connected demes). We relax the assumption of assortativity in a second type of network by randomly reassigning connections between demes, while maintaining the same degree distribution. (i.e., the probability of two demes being connected is proportionate to their degree). While well connected demes still have more connections to other well connected demes than to poorly connected demes, they are not more likely to be connected to a well connected deme than by chance based on the degree distribution. In both types of network structure, we classify a deme as well-connected if it is in the top 20% of the degree distribution and poorly connected if it is in the bottom 20%.We model the genetics using a multilocus gene-for-gene framework with haploid host and pathogen genotypes characterized by (L) biallelic loci, where 0 and 1 represent the presence and absence, respectively, of resistance and infectivity alleles. Host genotype (i) and pathogen genotype (j) are represented by binary strings: ({x}_{i}^{1}{x}_{i}^{2}ldots {x}_{i}^{L}) and ({y}_{j}^{1}{y}_{j}^{2}ldots {y}_{j}^{L}). Resistance acts multiplicatively such that the probability of host (i) being infected when challenged by pathogen (j) is ({Q}_{{ij}}={sigma }^{{d}_{{ij}}}), where (sigma) is the reduction in infectivity per effective resistance allele and ({d}_{{ij}}={sum }_{k=1}^{L}{x}_{i}^{k}big(1-{y}_{j}^{k}big)) is the number of effective resistance alleles (i.e., the number of loci where hosts have a resistance allele but pathogens do not have a corresponding infectivity allele). Hosts and pathogens with more resistance or infectivity alleles are assumed to pay higher fitness costs, ({c}_{H}left(iright)) eq. (4) and ({c}_{P}left(jright)) eq. (5) with:$${c}_{H}left(iright)={c}_{H}^{1}left(frac{1-{{{{{rm{e}}}}}}^{frac{{c}_{H}^{2}}{L}{sum }_{k=1}^{L}{x}_{i}^{k}}}{1-{{{{{rm{e}}}}}}^{{c}_{H}^{2}}}right)$$
(4)
and$${c}_{P}left(jright)={c}_{P}^{1}left(frac{1-{{{{{rm{e}}}}}}^{frac{{c}_{P}^{2}}{L}{sum }_{k=1}^{L}{y}_{j}^{k}}}{1-{{{{{rm{e}}}}}}^{{c}_{P}^{2}}}right)$$
(5)
where (0 , < , {c}_{H}^{1},; {c}_{P}^{1},le, 1) control the overall strength of the costs (i.e., the maximum proportional reduction in reproduction (hosts) or transmission rate (pathogens)) and ({c}_{H}^{2},; {c}_{P}^{2}in {{mathbb{R}}}_{ne 0}) control the shape of the trade-off. When ({c}_{H}^{2},; {c}_{P}^{2}, < , 0) the costs decelerate (increasing returns) and when ({c}_{H}^{2},; {c}_{P}^{2}, > , 0) the costs accelerate the costs accelerate (decreasing returns) (Supplementary Fig. 4). This formulation, therefore, allows for a wide-range of trade-off shapes that may occur in nature.The dynamics of the (finite) host and pathogen populations are modelled stochastically using the tau-leap method with a fixed step size of (tau=1). For population (p), the mean host birth rate at time (t) for host (i) (6) is$${B}_{i}^{p}left(tright)=left(aleft(1-{c}_{H}left(iright)right)-q{N}_{p}left(tright)right){S}_{i}^{p}left(tright)$$
(6)
where (a) is the maximum per-capita birth rate, (q) is the strength of density-dependent competition on births, ({N}_{p}left(tright)={S}_{i}^{p}left(tright)+{I}_{icirc }^{p}left(tright)) is the local host population size, ({S}_{i}^{p}left(tright)) and ({I}_{icirc }^{p}left(tright)={sum }_{j=1}^{n}{I}_{{ij}}^{p}left(tright)) are the local sizes of susceptible and infected individuals of genotype (i), and ({I}_{{ij}}^{p}left(tright)) is the local size of hosts of genotype (i) infected by pathogen (j). Host mutations occur at an average rate of ({mu }_{H}) per loci (limited to at most one mutation per time step), so that the mean number of mutations from host type (i) to ({i}^{{prime} }) is ({mu }_{H}{m}_{i{i}^{{prime} }}{B}_{i}^{p}left(tright)), where ({m}_{i{i}^{{prime} }}=1) if genotypes (i) and ({i}^{{prime} }) differ at exactly one locus, and is 0 otherwise.The mean local mortalities for susceptible and infected individuals are (b{S}_{i}^{p}left(tright)) and (left(b+alpha right){I}_{{ij}}^{p}left(tright)), respectively, where (b) is the natural mortality rate and (alpha) is the disease-associated mortality rate. The average number of infected hosts that recover is (gamma {I}_{{ij}}^{p}left(tright)), where (gamma) is the recovery rate.The mean number of new local infections of susceptible host type (i) by pathogen (j) eq. (7) is:$${INF}_{{ij}}^{p}left(tright)=beta left(1-{c}_{P}left(jright)right){Q}_{{ij}}{S}_{i}^{p}left(tright){Y}_{j}^{p}left(tright)$$
(7)
where (beta) is the baseline transmission rate and ({Y}_{j}^{p}left(tright)) is the local number of pathogen propagules following mutation and dispersal. Pathogen mutations occur in a similar manner to host mutations, with mutations from type (j) to ({j}^{{prime} }) occurring at rate ({mu }_{P}{m}_{j{j}^{{prime} }}{I}_{circ j}^{p}left(tright)) where ({mu }_{P}) is the mutation rate per loci (limited to at most one mutation per timestep) and ({I}_{circ j}^{p}left(tright)={sum }_{i=1}^{n}{I}_{{ij}}^{p}left(tright)) is the local number of pathogen (j.) Following mutation, the local number of pathogens of type (j) eq. (8) is:$${W}_{j}^{p}left(tright)={I}_{circ j}^{p}left(tright)left(1-{mu }_{P}Lright)+{mu }_{P}{m}_{j{j}^{{prime} }}{I}_{circ j}^{p}left(tright)$$
(8)
Pathogen dispersal occurs following mutation at a rate of (rho) between connected demes, given by the adjacency matrix ({G}_{{pr}}), with ({G}_{varSigma p}) the total number of connections for deme (p). The mean local number of pathogen propagules following mutation and dispersal eq. (9) is therefore:$${Y}_{j}^{p}left(tright)={W}_{j}^{p}left(tright)left(1-rho {G}_{varSigma p}right)+rho mathop {sum }limits_{r=1}^{{M}_{varSigma }}{G}_{{pr}}{W}_{j}^{r}left(tright)$$
(9)
We focus our parameter sweep on: (i) the structure of the network (assortative or random connections); (ii) the strength (left({c}_{H}^{1},; {c}_{P}^{1}right)) and shape (left({c}_{H}^{2},; {c}_{P}^{2}right)) of the trade-offs; (iii) the transmission rate (left(beta right)); and (iv) the dispersal rate (left(rho right)), fixing the remaining parameters as described in Supplementary Table 1 (preliminary investigations suggested they had less of an impact on the qualitative outcome) and conducting 100 simulations per parameter set. For each simulation we initially seed all populations with the most susceptible host type and place the least infective pathogen type in one of the well-connected populations to minimize the risk of early extinction. We then solve the dynamics for 10,000 time steps (preliminary investigations indicated this was a sufficient period for the metapopulations to reach a quasi-equilibrium in terms of overall resistance). We calculate the average level of resistance (proportion of loci with a resistance allele) between time steps 4001 and 5000 (transient dynamics) and over the final 1000 time steps (long-term dynamics) for well and poorly connected demes, categorized according to whether the disease is present in (infected) or absent from (uninfected) the local population at a given time point and discarding simulations where the pathogen is driven globally extinct.We compare the mean level of resistance in infected/uninfected poorly/well-connected populations across all simulations to the empirical results. We say that a simulation is a qualitative ‘match’ for the empirical findings if: (i) in poorly connected demes, the infected populations are on average at least 5% more resistant than uninfected populations; and (ii) in well-connected demes, the uninfected populations are on average at least 5% more resistant than infected populations. In other words, if ({R}_{{CS}}) is the mean resistance for a population with connectivity (C) ((C=W) and (C=P) for well and poorly connected demes, respectively) and infection status (S) ((S=U) and (S=I) for uninfected and infected populations, respectively), then a parameter set is a qualitative ‘match’ for the empirical findings if ({R}_{{WU}} > 1.05{R}_{{WI}}) and (1.05{R}_{{PI}}, > , 1.05{R}_{{PU}}). If these criteria are not met, then the parameter set is a qualitative ‘mismatch’ for the empirical findings. The model is not intended to be a replica of an empirical metapopulation, but rather is used to reveal the key factors which lead to qualitatively similar distributions of resistance and disease incidences observed in the study of the Åland islands. Hence, the purpose of the model is to determine which biological factors are likely to be crucial to the patterns observed herein.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Animals, housing, and treatmentsAll procedures involving animals were approved by the University of Florida Institutional Animal Care and Use Committee (Protocol #201709925). All methods were performance in accordance with the relevant guidelines and regulations, and permission and informed consent was obtained from the University of Florida (owners) for the use of the steers in this experiment.The experiment was carried out during July and August of 2017 at the Feed Efficiency Facility of University of Florida, North Florida Research and Education Center, located in Marianna, Florida (30°52′N, 85°11″W, 35 m asl). Both ‘Argentine’ bahiagrass and ‘Florigraze’ rhizoma peanut hays were obtained from commercial producers. The hay bales were stored in enclosed barns throughout the duration of the experiment.Twenty-five Brahman × Angus crossbred steers (Bos sp.) were utilized (average BW = 341 ± 17 kg, approx. 16 months of age). The steers were grazing bermudagrass (Cynodon dactylon) pastures, a C4 grass, prior to the start of the study. The day prior to the start of the experiment (e.g. day-1), steers brought to working facilities, where they remained 16 h off feed and water, in order to obtain shrunk bodyweights. On day 0 of the experiment, steers were weighed, blocked by bodyweight, and allocated to five treatments (5 steers per treatment) and housed in grouped pens. Hay intake was recorded utilizing GrowSafe© systems (GrowSafe Systems Ltd., Calgary, AB, Canada), which utilize radio frequency identification to record feed intake by weight change measured to the nearest gram. Water was available ad libitum. Forage treatments were offered ad libitum by providing sufficient hay to maintain full feed troughs throughout each day of the experiment. Treatments were five proportions of ‘Florigraze’ rhizoma peanut hay in ‘Argentine’ bahiagrass hay: (1) 100% bahiagrass hay (0% RP); (2) 25% rhizoma peanut hay + 75% bahiagrass hay (25% RP); (3) 50% rhizoma peanut hay + 50% bahiagrass hay (50% RP); (4) 75% rhizoma peanut hay + 25% bahiagrass hay (75% RP); (5) 100% rhizoma peanut hay (100% RP). Diet chemical composition is presented in Table 1. All treatment proportions were weighed and mixed on as-fed basis. Mixing of diets was done manually; no hay mixers or choppers were used, to minimize leaf shatter.Sample collectionSteers were housed for 32 days and sampling occurred on 0, 8, 16, 24, and 32 days after initiation of treatment diets; exception was for feces, which were collected on d-1 given steers were fasted on d-0 of the experiment. The hay mixtures offered to the steers were collected (10 samples of each diet) and analyzed for nutritive value (Table 1), at the start of the experiment. All sampling occurred between 0700 and 1000 h on each of the sampling days.Fecal samples were collected directly from the rectum and placed in quart-sized plastic bags to avoid contamination. The feces were frozen at −20 °C. All fecal samples were thawed, dried at 55 °C for 72 h, and ground to pass a 2-mm stainless steel screen using a Wiley Mill (Model 4, Thomas-Wiley Laboratory Mill, Thomas Scientific, Swedesboro, NJ, USA). Samples were then ball milled using a Mixer Mill MM400 (Retsch GmbH, Haan, Germany) at 25 Hz for 9 min.Blood was obtained through jugular venipuncture using 10-mL K2 EDTA vials (Becton Dickinson and Company, Franklin Lakes, NJ, USA), and stored in ice until centrifugation. All tubes were centrifuged at 714 G for 20 min using an Allegra X-22R centrifuge (Beckman Coulter, Brea, CA, USA). A 10-mL sample of plasma was collected and placed in a separate glass vial, the remaining plasma, white blood cell, and platelet fractions were discarded. The remaining RBC pellet was re-suspended with 9 vol. 0.9% NaCl solution and mixed at room temperature for 15 min at 2 Hz orbital shaker. The tubes were then centrifuged at 714 G for 20 min. The saline solution from the centrifuged tubes was discarded after centrifugation. The rinse procedure was repeated two more times for a total of three rinses. After the third rinse procedure, a 500-µL sample was removed, frozen at −20 °C, and subsequently freeze-dried for isotopic analyses.Hair clippings were obtained from an area of 20 × 20 cm on the left hindquarter, utilizing veterinary hair clippers (Sunbeam-Oster Inc., Boca Raton, FL, USA). Hair clippings were collected, placed in nylon bags (Ankom Technology, Macedon, NY, USA), and frozen for subsequent analysis. Clippings were always collected in the same location from each animal in order to ensure new hair growth would be analyzed. All hair samples were cleaned using soapy water and defatted following protocols for other keratin-based tissues 31,34. Each sample was sonicated twice for 30 min in a methanol and chloroform solution (2:1, v/v), rinsed with distilled water, and oven dried overnight at 60 °C. Each hair sample was ball milled using a Mixer Mill MM400 (Retsch GmbH, Haan, Germany) at 25 Hz for 9 min.CalculationsAfter processing, all samples were analyzed for total C and N using a CHNS analyzer through the Dumas dry combustion method (Vario MicroCube, Elementar Americas Inc., Ronkonkoma, NY, USA) coupled to an isotope ratio mass spectrometer (IsoPrime 100, Elementar, Elementar Americas Inc., Ronkonkoma, NY, USA). Certified standards of L-glutamic acid (USGS40, USGS41; United States Geological Survey) were used for isotope ratio mass spectrometer calibration. Isotope ratios were as follows: δ13C of −26.39, + 37.63‰, and δ15N of −4.52, 47.57‰ for USGS40 and USGS41, respectively. Calibration of the IRMS was conducted according to Cook, et al. 35, with an accuracy of ≤ 0.06 ‰ for 15N and ≤ 0.08 ‰ for 13C.The isotope ratio for 13C/12C was calculated as:$$delta^{{{13}}} {text{C}} = , left( {^{{{13}}} {text{C}}/^{{{12}}} {text{C}}_{{{text{sample}}}} {-}^{{{13}}} {text{C}}/^{{{12}}} {text{C}}_{{{text{reference}}}} } right)/ , left( {^{{{13}}} {text{C}}/^{{{12}}} {text{C}}_{{{text{reference}}}} times { 1}000} right)$$
(1)

where δ13C is the C isotope ratio of the sample relative to Pee Dee Belemnite (PDB) standard (‰), 13C/12Csample is the C isotope ratio of the sample, and 13C/12Creference is the C isotope ratio of PDB standard 5. The isotope ratio for 15N/14N was calculated as:$$delta^{{{15}}} {text{N}} = , left( {^{{{15}}} {text{N}}/^{{{14}}} {text{N}}_{{{text{sample}}}} -^{{{15}}} {text{N}}/^{{{14}}} {text{N}}_{{{text{reference}}}} } right)/left( {^{{{15}}} {text{N}}/^{{{14}}} {text{N}}_{{{text{reference}}}} times { 1}000} right)$$
(2)
where δ15N is the N isotope ratio of the sample relative to atmospheric nitrogen (‰), 15N/14Nsample is the N isotope ratio of the sample, and 15N/14Nreference is the N isotope ratio of atmospheric N (standard) 5. The fraction factor (Δ) in this study was considered to be the difference between the diet isotopic composition (δ13C or δ15N) and that of the given sample 5.The dietary proportion of rhizoma peanut hay was back-calculated using δ13C and δ15N of each plant on a DM basis 3. This method is advantageous in that it does not require further tissue processing and facilitates implementation at the field scale. The proportion of rhizoma peanut was estimated using Eq. (3), described by Jones et al. 3:$$%RP=100-left{100 times frac{A-C}{B-C}right}$$
(3)
where %RP is the proportion of RP in the diet, A is the δ13C or δ15N of the sample obtained, B is the δ13C or δ15N of bahiagrass, and C is the δ13C or δ15N of RP.Statistical analysisAll response variables were analyzed using linear mixed model procedures as implemented in SAS PROC GLIMMIX (SAS/STAT 15.1, SAS Institute). Individual animals were considered the experimental unit. Treatment, collection day, and their interaction were considered fixed effects, and block was considered a random effect in the model. The data were analyzed as repeated measures, considering collection day as the repeated measure. The best covariance matrix was selected according to the lowest AICC fit statistic. Least squares treatment means were compared through pairwise t test using the PDIFF option of the LSMEAN statement in the aforementioned procedure. Based on the recommendations by Milliken and Johnson 36 and Saville 37, no adjustment for multiple comparisons was made. Orthogonal polynomial contrasts (linear and quadratic effects) were used to test effects of RP inclusion on isotopic responses. The slice option was used when the treatment × collection day interaction was significant (P ≤ 0.05), using collection day as the factor, to test treatment effects across collection days. Significance was declared at P ≤ 0.05. The contrast statement was used to test for linear or quadratic effects. Regression analyses for the dietary predictions were conducted using PROC REG from SAS.Predictions of dietary proportions based on Eq. (3) were generated for 16 subgroups reflecting combinations of isotope (13C or 15N), day (8 or 32), and sample type (feces, plasm, RBC, or hair). Analyses comparing predicted versus actual diet proportions included several components. First, we computed the concordance correlation coefficient (CCC) following the recommendations from Crawford, et al. 38. The CCC is a measure of agreement that encompasses both precision and accuracy, along with corresponding 95% bias accelerated and corrected (BCa) bootstrap confidence intervals. For comparative purposes we calculated the Pearson correlation coefficient which only reflects precision. Both correlation coefficients range from −1.0 to 1.0 and we interpreted values ≥ 0.80 as indicating strong positive agreement/correlation. Next, we regressed the actual percentages on the predicted percentages using linear regression. Perfect prediction corresponds to the estimated regression line having an intercept of zero and a slope of 1.0. We then partitioned the mean square error (MSE) of the predicted (from Eq. (3), not the above linear regression) versus actual percentages as described in Rice and Cochran 39. This partitioning entails calculating the proportion of MSE attributable to three sources of error: the difference in mean predicted and actual values (mean component, denoted “MC”), the error resulting from the slope of the above linear regression deviating from 1.0 (slope component, denoted “SC”), and random error (random component, denoted “RC”). The results from the above analyses were examined to identify subgroups whose predictions were sufficiently good to warrant hypothesis testing. In this context “good” means that the predicted percentages were strongly correlated with the actual percentages and the magnitudes of the predicted percentages were similar to the actual percentages. The objective of the hypothesis testing was to formally evaluate whether dietary proportions could be directly predicted from Eq. (3) (in contrast to generating predictions using the equation from regressing actual dietary percentages on the predicted percentages from Eq. (3)). Paired two one-sided test (TOST) equivalence tests were conducted for the selected subgroups with α = 0.0540. These tests are formulated such that the null hypothesis is “non-equivalence” and the alternative hypothesis is “equivalence”. An input parameter to the test is the equivalence region, a range for which we consider the mean actual minus predicted difference to be unimportant (“equivalent”) from a practical standpoint. For each equivalence test we also computed the 90% confidence interval for the mean actual minus predicted difference which we refer to as the “minimum equivalence region”. Based on the structure of TOST equivalence tests, to reject the null hypothesis at the 0.05 level, the equivalence region specified for the test must completely contain the minimum equivalence region. For example, if the pre-specified equivalence region is (−15%, 15%) and the computed minimum equivalence region is (−16%, −6%) the null hypothesis would not be rejected. Finally, we re-ran all of the analyses described above for the selected subgroups where, prior to analysis, predicted percentages outside of the valid range were assigned the appropriate boundary value (i.e., predicted percentages  100% were assigned a value of 100%). More

Vega, F. E., Rosenquist, E. & Collins, W. Global project needed to tackle coffee crisis. Nature 425, 343 (2003).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Craparo, A. C. W., Van Asten, P. J. A., Läderach, P., Jassogne, L. T. P. & Grab, S. W. Coffea arabica yields decline in Tanzania due to climate change: global implications. Agric. For. Meteorol. 207, 1–10 (2015).ADS 
Article 

Google Scholar 
Davis, A. P. et al. High extinction risk for wild coffee species and implications for coffee sector sustainability. Sci. Adv. 5, eaav3473 (2019).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Davis, A. P., Gole, T. W., Baena, S. & Moat, J. The impact of climate change on indigenous arabica coffee (Coffea arabica): predicting future trends and identifying priorities. PLoS ONE 7, e47981 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Davis, A. P., Mieulet, D., Moat, J., Sarmu, D. & Haggar, J. Arabica-like flavour in a heat-tolerant wild coffee species. Nat. Plants 7, 413–418 (2021).CAS 
PubMed 
Article 

Google Scholar 
Moat, J., Gole, T. W. & Davis, A. P. Least concern to endangered: applying climate change projections profoundly influences the extinction risk assessment for wild Arabica coffee. Global Change Biol. 25, 390–403 (2019).ADS 
Article 

Google Scholar 
Moat, J. et al. Resilience potential of the Ethiopian coffee sector under climate change. Nat. Plants 3, 17081 (2017).PubMed 
Article 

Google Scholar 
Kath, J. et al. Not so robust: Robusta coffee production is highly sensitive to temperature. Global Change Biol. 26, 3677–3688 (2020).ADS 
Article 

Google Scholar 
Liu, L. et al. Soil moisture dominates dryness stress on ecosystem production globally. Nat. Commun. 11, 1–9 (2020).ADS 
CAS 

Google Scholar 
Grossiord, C. et al. Plant responses to rising vapor pressure deficit. New Phytol. 226, 1550–1566 (2020).PubMed 
Article 

Google Scholar 
IPCC Climate Change 2022: Impacts, Adaptation, and Vulnerability (eds. Pörtner, H.-O. et al.) (Cambridge Univ. Press, 2022).Burke, M. et al. Higher temperatures increase suicide rates in the United States and Mexico. Nat. Clim. Change 8, 723–729 (2018).ADS 
Article 

Google Scholar 
Burke, M., Hsiang, S. M. & Miguel, E. Global non-linear effect of temperature on economic production. Nature 527, 235–239 (2015).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Duffy, K. A. et al. How close are we to the temperature tipping point of the terrestrial biosphere? Sci. Adv. 7, eaay1052 (2021).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Scheffer, M. et al. Early-warning signals for critical transitions. Nature 461, 53–59 (2009).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Schneider, S. H. Abrupt non-linear climate change, irreversibility and surprise. Global Environ. Change 14, 245–258 (2004).Article 

Google Scholar 
Lenton, T. M. Early warning of climate tipping points. Nat. Clim. Change 1, 201–209 (2011).ADS 
Article 

Google Scholar 
Lenton, T. M. et al. Climate tipping points—too risky to bet against. Nature. 575, 592–595 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Lobell, D. B., Bänziger, M., Magorokosho, C. & Vivek, B. Nonlinear heat effects on African maize as evidenced by historical yield trials. Nat. Clim. Change 1, 42–45 (2011).ADS 
Article 

Google Scholar 
Lobell, D. B., Deines, J. M. & Tommaso, S. D. Changes in the drought sensitivity of US maize yields. Nat. Food 1, 729–735 (2020).Article 

Google Scholar 
Lobell, D. B. et al. Greater sensitivity to drought accompanies maize yield increase in the US Midwest. Science 344, 516–519 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Rigden, A., Mueller, N., Holbrook, N., Pillai, N. & Huybers, P. Combined influence of soil moisture and atmospheric evaporative demand is important for accurately predicting US maize yields. Nat. Food 1, 127–133 (2020).Article 

Google Scholar 
Schlenker, W. & Roberts, M. J. Nonlinear temperature effects indicate severe damages to US crop yields under climate change. Proc. Natl Acad. Sci. USA 106, 15594–15598 (2009).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
McDowell, N. G. et al. Mechanisms of woody-plant mortality under rising drought, CO2 and vapour pressure deficit. Nat. Rev. Earth Environ. 3, 294–308 (2022).ADS 
CAS 
Article 

Google Scholar 
Sinclair, T. R. et al. Limited-transpiration response to high vapor pressure deficit in crop species. Plant Sci. 260, 109–118 (2017).CAS 
PubMed 
Article 

Google Scholar 
López, J., Way, D. A. & Sadok, W. Systemic effects of rising atmospheric vapor pressure deficit on plant physiology and productivity. Global Change Biol. 27, 1704–1720 (2021).ADS 
Article 

Google Scholar 
McDowell, N. G. & Allen, C. D. Darcy’s law predicts widespread forest mortality under climate warming. Nat. Clim. Change 5, 669–672 (2015).ADS 
Article 

Google Scholar 
Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A. & Hegewisch, K. C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Sci. Data 5, 170191 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
You, L., Wood, S., Wood-Sichra, U. & Wu, W. Generating global crop distribution maps: from census to grid. Agric. Syst. 127, 53–60 (2014).Article 

Google Scholar 
Fong, Y., Huang, Y., Gilbert, P. B. & Permar, S. R. chngpt: threshold regression model estimation and inference. BMC Bioinformatics 18, 1–7 (2017).Article 

Google Scholar 
Qin, Y. et al. Agricultural risks from changing snowmelt. Nat. Clim. Change 10, 459–465 (2020).ADS 
Article 

Google Scholar 
Forster, P. M., Maycock, A. C., McKenna, C. M. & Smith, C. J. Latest climate models confirm need for urgent mitigation. Nat. Clim. Change 10, 7–10 (2020).ADS 
Article 

Google Scholar 
Forster, P. M. et al. Projections of when temperature change will exceed 2 °C above pre-industrial levels. Nat. Clim. Change 10, 407–412 (2011).
Google Scholar 
Joshi, M., Hawkins, E., Sutton, R., Lowe, J. & Frame, D. Projections of when temperature change will exceed 2 °C above pre-industrial levels. Nat. Clim. Change 1, 407–412 (2011).ADS 
Article 

Google Scholar 
IPCC, 2021: Summary for Policymakers. In Climate Change 2021: The Physical Science Basis (eds Masson-Delmotte, V. et al.) (Cambridge Univ. Press, in press).Lobell, D. B. et al. The critical role of extreme heat for maize production in the United States. Nat. Clim. Change 3, 497–501 (2013).ADS 

Google Scholar 
Sinclair, T. R., Hammer, G. L. & Van Oosterom, E. J. Potential yield and water-use efficiency benefits in sorghum from limited maximum transpiration rate. Funct. Plant Biol. 32, 945–952 (2005).PubMed 
Article 

Google Scholar 
Martins, M. Q. et al. Protective response mechanisms to heat stress in interaction with high [CO2] conditions in Coffea spp. Front. Plant Sci. 7, 947 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Rodrigues, W. P. et al. Long‐term elevated air [CO2] strengthens photosynthetic functioning and mitigates the impact of supra‐optimal temperatures in tropical Coffea arabica and C. canephora species. Global Change Biol. 22, 415–431 (2016).ADS 
Article 

Google Scholar 
Ghini, R. et al. Coffee growth, pest and yield responses to free-air CO2 enrichment. Clim. Change 132, 307–320 (2015).ADS 
Article 

Google Scholar 
Rakocevic, M. et al. The vegetative growth assists to reproductive responses of Arabic coffee trees in a long-term FACE experiment. Plant Growth Regul. 91, 305–316 (2020).CAS 
Article 

Google Scholar 
Hammer, G. L. et al. Designing crops for adaptation to the drought and high‐temperature risks anticipated in future climates. Crop Sci. 60, 605–621 (2020).Article 

Google Scholar 
Gennari, P., Rosero-Moncayo, J. & Tubiello, F. N. The FAO contribution to monitoring SDGs for food and agriculture. Nat. Plants 5, 1196–1197 (2019).PubMed 
Article 

Google Scholar 
Lesk, C., Rowhani, P. & Ramankutty, N. Influence of extreme weather disasters on global crop production. Nature 529, 84–87 (2016).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Ortiz-Bobea, A., Ault, T. R., Carrillo, C. M., Chambers, R. G. & Lobell, D. B. Anthropogenic climate change has slowed global agricultural productivity growth. Nat. Clim. Change 11, 306–312 (2021).ADS 
Article 

Google Scholar 
Davis, A. P. et al. Hot coffee: the identity, climate profiles, agronomy, and beverage characteristics of Coffea racemosa and C. zanguebariae. Front. Sustain. Food Syst. 5, 740137 (2021).Article 

Google Scholar 
Sarmiento-Soler, A. et al. Disentangling effects of altitude and shade cover on coffee fruit dynamics and vegetative growth in smallholder coffee systems. Agric. Ecosyst. Environ. 326, 107786 (2022).CAS 
Article 

Google Scholar 
Wood, S. N. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. J. R. Stat. Soc. B 73, 3–36 (2011).MathSciNet 
MATH 
Article 

Google Scholar 
Barton, K. MuMIn: multi-model inference. R-Forge http://r-forge.r-project.org/projects/mumin/ (2009).R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing https://www.r-project.org/ (2021).Harrison, X. A. et al. A brief introduction to mixed effects modelling and multi-model inference in ecology. PeerJ 6, e4794 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Najafi, E., Devineni, N., Khanbilvardi, R. M. & Kogan, F. Understanding the changes in global crop yields through changes in climate and technology. Earths Future 6, 410–427 (2018).ADS 
CAS 
Article 

Google Scholar 
Ovalle-Rivera, O. et al. Assessing the accuracy and robustness of a process-based model for coffee agroforestry systems in Central America. Agrofor. Syst. 94, 2033–2051 (2020).Article 

Google Scholar 
Varma, S. & Simon, R. Bias in error estimation when using cross-validation for model selection. BMC Bioinformatics 7, 1–8 (2006).Article 

Google Scholar 
Yuan, W. et al. Increased atmospheric vapor pressure deficit reduces global vegetation growth. Sci. Adv. 5, eaax1396 (2019).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Son, H. & Fong, Y. Fast grid search and bootstrap-based inference for continuous two-phase polynomial regression models. Environmetrics 32, e2664 (2021).MathSciNet 
Article 

Google Scholar 
Wintgens, J. N. et al. Coffee: Growing, Processing, Sustainable Production. A Guidebook for Growers, Processors, Traders, and Researchers (Wiley, 2004). More

Description of PTEsThe descriptive statistics of the contents of soil PTEs in the study area were shown in Table 1. From Table 1, the mean contents of As and Ni in the oil-affected soils exceeded their corresponding risk screening values33, which may damage the soil ecological environment and affect crop growth. Compared with the secondary standard of soil environmental quality34, the mean contents of As, Cu and Zn were all lower than their corresponding Grade II standard values, but the mean contents of Cd, Cr, Ni and Pb in the oil-affected soils were 1.07, 7.46, 7.14 and 1.36 times of their standard values. In contrast with the background value of Hubei province35, except Mn, the mean contents of As, Cd, Cr, Cu, Ni, Pb, Zn and Ba in the oil-affected soils all exceeded their background values. Meanwhile, the variation coefficient of Cr (1.41) was greater than 1. In general, the soil Cd concentration in the study area was higher than that around Gudao Town, a typical oil-producing region of the Shengli Oilfield in the Yellow River Delta, China12, and from Yellow River Delta, a traditional oil field in China9, but was lower than that around two crude oil flow stations in the Niger Delta, Nigeria36. The concentrations of other PTEs were higher than the corresponding element concentrations, detected in the soil around Gudao Town, a typical oil-producing region of the Shengli Oilfield in the Yellow River Delta, China12, from Yellow River Delta, a traditional oil field in China9, and around two crude oil flow stations in the Niger Delta, Nigeria36. The above analysis exhibited that PTEs in the oil-affected soils had a certain degree of accumulation and may be affected by human activities.Table 1 Statistical characteristics for potential toxic elements in in the study area (mg·kg−1).Full size tableLevels of PTEs enrichment and pollutionThe EF and PLI of soil PTEs in the study area were calculated to evaluate the pollution degree of soil PTEs. The calculation results of EF and PLI were shown in Fig. 2 and Table S4. From Fig. 2, the mean EF values of PTEs were showed as Pb  > Cr  > Ni  > As  > Cd  > Zn  > Cu  > Ba. The mean EFs of all PTEs were greater than 1. Among them, the average EF of Cu, Zn and Ba was between 1 and 2, which was slightly enriched. And As (2.18) and Cd (2.12) were moderately enriched. In particular, the average EF values of Cr, Ni and Pb were 14.23, 8.69 and 15.45, respectively, reaching a significant enrichment level, and all samples of Cr, Ni and Pb were at moderate or above enrichment, of which 10% of the Cr samples were extreme pollution, 85% of Cr samples, 95% of Ni and 5% of Pb (Table S4) were significantly enriched. These proved that these PTEs were generally enriched in the study area, especially Cr, Ni and Pb.Figure 2The map of enrichment factor and contamination factor of PTEs in the study area.Full size imageExcept Mn, the average CF values of other PTEs were all  > 1 (Fig. 2), indicating that the accumulation of Mn in the study area was relatively light, and there was no obvious Mn pollution. The CF values of all samples of As, Cr, Ni and Pb, 80% of Cd samples, 75% of Cu samples, 30% of Mn samples, 65% of Zn samples and 75% of Ba samples (Table S4) were higher than 1. And the mean CF values of Cr, Ni and Pb were 14.21, 7.58 and 12.73, respectively, certifying that the pollution of Cr, Ni and Pb in the study area was considerably serious. PLI was calculated based on the CF value of PTEs, and the results were shown in Fig. 2. The average value of PLI was 2.62, indicating that the soil PTEs in the study area were seriously polluted.Spatial distribution of soil PTEs in the study areaGeostatistical analysis was utilized to do ordinary Kriging interpolation of the PTEs in the study area, the results were shown in Fig. 3. As shown in Fig. 3, the spatial distribution of As, Cr, Ni, Zn and Ba was relatively consistent, and their hot spots were concentrated in the southeast, northwest, and central and eastern parts of the study area where oil wells were distributed. The spatial distribution of Cr and Ni exhibited that there were large-scale hotspots near the oil wells, and the content of Cr and Ni in these hotspots was much higher than second-level environmental quality standards of China, which proved that the content of soil Cr and Ni was significantly affected by the oil production activities of the oil production plant. There were crude oil leaks in B and C, and the contents of Zn and Ba in the vicinity of these two oil wells were relatively high, indicating that soil Zn and Ba in this area may be affected by the crude oil leakage, resulting in a certain degree of accumulation in the soil. The area with the second highest As content mainly resided in the middle of the study area. According to the survey, the herbicides were sprayed every year around the H oil well in the middle of the study area, indicating that the accumulation of As in the soil was not only related to oil extraction activities, but also to the use of pesticides (contains copper arsenate, sodium arsenate, etc.)10, 14. In addition, the hot spots of spatial distribution of Pb, Cd and Mn were concentrated in the southeast, and Cu was mainly concentrated in the southeast and midwest. As analyzed above, in addition to Mn, the PTEs Pb, Cd and Cu all have a certain degree of accumulation. And the investigation found that there were many petroleum machinery manufacturing plants in the central and eastern part of the study area, therefore, the accumulation of Pb, Cd and Cu in the soil may be related to factors such as petroleum extraction, crude oil leakage and machinery manufacturing. The above analysis indicated that the influence of human activities is evident on the distribution of soil PTEs3, 23.Figure 3spatial distribution map of soil PTEs in the study area.Full size imagePotential ecological risk assessmentThe potential ecological risk assessment model after adjusting the threshold was used to evaluate the PER of the oil production plant. The individual potential ecological risk of PTEs was shown in Table 2. From Table 2, the average ({E}_{r}^{i}) values of PTEs were Cr  > Pb  > Cd  > Ni  > As  > Cu  > Zn  > Mn. The average ({E}_{r}^{i}) values of Cr and Pb were 79.62 and 63.64, respectively, reaching a relatively high level of potential ecological risk; the average ({E}_{r}^{i}) values of Cd and Ni were 55.95 and 37.91, respectively, which were at medium potential ecological risk level; the average ({E}_{r}^{i}) values of other PTEs were all lower than 30, with minor potential ecological risk. Specifically, all samples of Cu, Mn and Zn were at slight potential ecological risk level; 5% of As samples, 80% of Cd, 85% of Cr, 80% of Ni and 100% of Pb (Table S5) were at medium and above potential ecological risk. In particular, the potential ecological risks of 35% of Cd samples, 10% of Cr samples, 5% of Ni samples and 80% of Pb samples (Table S5) were relatively high, 10% Cd samples reached high potential ecological risk level, and 10% Cr samples had extremely high potential ecological risk. In summary, Geostatistical analysis shows that the hotspot distribution of all PTEs in the study area is almost related to the distribution of oil wells. In addition, the hotspot distribution of PTEs may also be related to factors such as agricultural and industrial activities3. The average value of PER in the study area was 265.08, and the proportions of the three risk levels of medium, slightly high and high were 5%, 75% and 20%, respectively (Table S5). It proved that the study area was at a higher potential ecological risk. Among them, the PER values of samples A, B, D, E, F, G, H, I and J (Table 2) were all greater than 280, reaching fairly high ecological risk.Table 2 Single ecological risk index and potential ecological risk of soil PTEs in study area.Full size tableHuman health risk assessmentThe non-carcinogenic risk assessment of As, Cd, Cr, Cu, Mn, Ni, Pb, Zn and Ba in the soils of the study area was carried out, and the assessment results were shown in Table 3. The THI values of children and adults under the three exposure routes of soil PTEs in the study area were 7.31 and 1.03, respectively, and the THI values were all  > 1, which indicated that soil PTEs around the oil production plants posed significant non-carcinogenic health risks to children and adults. The non-carcinogenic hazardous quotient (HQ) of children and adults in Table 3 revealed that the HQ of all PTEs for adults under each exposure route was less than 1, while the HQ of Cr and Pb for children under the oral intake route was greater than 1, which were 4.91 and 1.17, respectively. For HQ with different exposure routes of the same PTE, each soil PTE presented the risk of oral ingestion  > oral and nasal inhalation risk  > skin contact risk. The result was in agreement with the reports14, 37. Therefore, oral intake was the main exposure route of non-carcinogenic risk, and oral intake of Cr and Pb caused serious non-carcinogenic risk to children. Statistical analysis of HI for soil PTEs in the study area showed that the HI values of PTEs for children were significantly higher than those of adults, and the HI values of PTEs in children and adults were all Cr  > Pb  >   > As  > Ni  > Mn  > Ba  > Cu  > Zn  > Cd. Among them, the HI values of all PTEs for adults were less than 1, indicating that the non-carcinogenic risks caused by a single PTE did not have a significant impact on adults; while the HI values of Cr and Pb for children were 4.93 and 1.17 greater than 1, indicating that they have caused serious non-carcinogenic risk to local children. In addition, the HI values of As and Ni for children and the HI values of As, Cr and Pb for adults were all greater than 0.1, which requires attention. In summary, children suffered from significant non-carcinogenic risk, and adults suffered from minor non-carcinogenic risk in the study area; soil Cr and Pb were the most important non-carcinogenic risk factors for children and adults in the study area.Table 3 Non-cancer and cancer risk assessment of adults and children under different exposure routes.Full size tableIn this study, soil As, Cd, Cr, Ni and Pb from the study area were assessed for carcinogenic risk, and the results were shown in Table 3. The TCRI of children and adults under the three exposure routes of these five PTEs were 9.44E−04 and 5.75E−04, respectively, indicating that soil PTEs around the oil production plants have caused serious carcinogenic risk to local children and adults. The CR values of children and adults showed that the CR values of Cr (6.33E−04) and Ni (2.64E−04) for children, and Cr (3.87E−04) and Ni (1.49E−04) for adults were all greater than 10–4. In addition, As, Cr and Cd all presented oral intake risk  > oronasal inhalation risk  > skin contact risk. In conclusion, Cr and Ni caused serious carcinogenic risk for children and adults in the study area, and oral intake was also the primary way of carcinogenic risk. The CRI statistics of adults and children exhibited that the CRI values of all PTEs were lower than those of children. The CRI values of the PTEs in adults and children under the three exposure routes were Cr  > Ni  >   > As  > Pb  >   > Cd. Among them, the CRI values of Cr and Ni in children and adults by oral intake were both greater than 10–4, showing a strong carcinogenic risk. It is noteworthy that the assessment based on total concentrations of PTEs in soil might overestimate potential health risks38. The above analysis revealed that both children and adults in the study area suffered from serious carcinogenic risks, and Cr and Ni were the chiefly carcinogenic risk factors. More

We developed the IUCN Global Ecosystem Typology in the following sequence of steps: design criteria; hierarchical structure and definition of levels; generic ecosystem assembly model; top-down classification of the upper hierarchical levels; iterative circumscription of the units and ecosystem-specific adaptations of the assembly model; full description of the units; and map compilation. Some iteration proved necessary, as the description and review process sometimes revealed a need for circumscribing additional units.Design criteria and other typologiesUnder the auspices of the IUCN Commission on Ecosystem Management, we developed six design principles to guide the development of a typology that would meet the needs for global ecosystem reporting, risk assessment, natural capital accounting and ecosystem management: (1) representation of ecological processes and ecosystem functions; (2) representation of biota; (3) conceptual consistency throughout the biosphere; (4) scalable structure; (5) spatially explicit units; and (6) parsimony and utility (see Supplementary Table 1.1 and Supplementary Information, Appendix 1 for definitions and rationale).We assessed 23 existing ecological classifications with global coverage of terrestrial, freshwater, and/or marine environments against these principles to determine their fitness for IUCN’s purpose (Supplementary Information, Appendix 1). These include general classifications of land, water or bioclimate, as well as classifications of units that conform with the definition of ecosystems adopted in the United Nations Convention on Biological Diversity45 or an equivalent definition in the IUCN Red List of Ecosystems30. We reviewed documentation on methods of derivation, descriptions of classification units and maps to assess each classification against the six design principles (Supplementary Table 1.2 for details).Typology structure and ecosystem assemblyWe developed the structure of the Global Ecosystem Typology and the generic ecosystem assembly model at a workshop attended by 48 terrestrial, freshwater and marine ecosystem experts at Kings College London, UK, in May 2017. Participants agreed that a hierarchical structure would provide an effective framework for integrating ecological processes and functional properties (Supplementary Table 1.1, design principle 1), and biotic composition (principle 2) into the typology, while also meeting the requirement for scalability (principle 4). Although neither function nor composition were intended to take primacy within the typology, we reasoned that a hierarchy representing functional features in the upper levels is likely to support generalizations and predictions by leveraging evolutionary convergence13. By contrast, a typology reflecting compositional similarities in its upperlevels is less likely to be stable owing to dynamism of species assemblages and evolving knowledge on species taxonomy and distributions. Furthermore, representation of compositional relationships at a global scale would require many more units in upper levels, and possibly more hierarchical levels. Therefore, we concluded that a hierarchical structure recognizing compositional variants at lower levels within broad functionally based groupings at upper levels would be more parsimonious and robust (principle 6) than one representing composition at upper levels and functions at lower levels.Workshop participants initially agreed that three hierarchical levels for ecosystem function and three levels for biotic composition could be sufficient to represent global variation across the whole biosphere. Participants developed the concepts of these levels into formal definitions (Supplementary Table 3.1), which were reviewed and refined during the development process.To ensure conceptual consistency of the typology and its units throughout the biosphere (principle 3), we drew from community assembly theory to develop a generic model of ecosystem assembly. The traditional community assembly model incorporates three types of filters (dispersal, the abiotic environment and biotic interactions) that determine which biota from a larger pool of potential colonists can occupy and persist in an area13. We extended this model to ecosystems by: (1) defining three groups of abiotic filters (resources, ambient environment and disturbance regimes) and two groups of biotic filters (biotic interactions and human activity); (2) incorporating evolutionary processes that shape characteristic biotic properties of ecosystems over time; (3) defining the outcomes of filtering and evolution in terms of all ecosystem properties including both ecosystem-level functions and species-level traits, rather than only in terms of species traits and composition; and (4) incorporating interactions and feedbacks among filters and selection agents and ecosystem properties to elucidate hypotheses about processes that influence temporal and spatial variability in the properties of ecosystems and their component biota. In community assembly, only a small number of filters are likely to be important in any given habitat13. In keeping with this proposition, we used the generic model to identify biological and physical features that distinguish functionally different groups of ecosystems from one another by focusing on different ecological drivers that come to the fore in structuring their assembly and shaping their properties.Hierarchical levelsThe top level of classification (Fig. 2 and Extended Data Tables 1–4) defines five core realms of the biosphere based on contrasting media that reflect ecological processes and functional properties: terrestrial; freshwaters and inland saline waters (hereafter freshwater); marine; subterranean; and atmospheric. Biome gradient concepts25 highlight continuous variation in ecosystem properties, which is represented in the typology by transitional realms that mark the interfaces between the five core realms (for example, floodplains (terrestrial–freshwater), estuaries (freshwater–marine), and so on). In Supplementary Information, Appendix 3 (pages 3–16) and Supplementary Table 3.1, we describe the five core realms and review the hypothesized assembly filters and ecosystem properties that distinguish different groups within them. The atmospheric realm is included for comprehensive coverage, but we deferred resolution of its lower levels because its biota is poorly understood, sparse, itinerant and represented mainly by dispersive life stages46.Functional biomes (level 2) are components of the biosphere united by one or more major assembly processes that shape key ecosystem functions and ecological processes, irrespective of taxonomic identity (Supplementary Information, Appendix 3, page 17). Our interpretation aligns broadly with ‘functional biomes’ described elsewhere24,25,47, extended here to reflect dominant assembly filters and processes across all realms, rather than the more restricted basis of climate-vegetation relationships that traditionally underpin biome definition on land. Hence, the 25 functional biomes (Supplementary Information, Appendix 4, pages 52–186 and https://global-ecosystems.org/) include some ‘traditional’ terrestrial biomes47, as well as lentic and lotic freshwater systems, pelagic and benthic marine systems, and anthropogenic functional biomes assembled and usually maintained by human activity48.Level 3 of the typology defines 110 ecosystem functional groups described with illustrated profiles in Supplementary Information, Appendix 4 (pages 52–186) and at https://global-ecosystems.org/. These are key units for generalization and prediction, because they include ecosystem types with convergent ecosystem properties shaped by the dominance of a common set of drivers (Supplementary Information, Appendix 3, pages 17–19). Ecosystem functional groups are differentiated along environmental gradients that define spatial and temporal variation in ecological drivers (Figs. 2 and 3 and Supplementary Figs. 3.2 and 3.4). For example, depth gradients of light and nutrients differentiate functional groups in pelagic ocean waters (Fig. 3c and Extended Data Table 4), influencing assembly directly and indirectly through predation. Resource gradients defined by flow regimes (influenced by catchment precipitation and evapotranspiration) and water chemistry, modulated by environmental gradients in temperature and geomorphology, differentiate functional groups of freshwater ecosystems25 (Fig. 3b and Extended Data Table 3). Terrestrial functional groups are distinguished primarily by gradients in water and nutrient availability and by temperature and seasonality (Fig. 3a and Extended Data Table 1), which mediate uptake of those resources and regulate competitive dominance and productivity of autotrophs. Disturbance regimes, notably fire, are important global drivers in assembly of some terrestrial ecosystem functional groups49.Three lower levels of the typology distinguish functionally similar ecosystems based on biotic composition. Our focus in this paper is on global functional relationships of ecosystems represented in the upper three levels of the typology, but the lower levels (Supplementary Information, Appendix 3, pages 19 and 20) are crucial for representing the biota in the typology, and facilitate the scaling up of information from established local-scale typologies that support decisions where most conservation action takes place. These lower levels are being developed progressively through two contrasting approaches with different trade-offs, strengths and weaknesses. First, level 4 units (regional ecosystem subgroups) are ecoregional expressions of ecosystem functional groups developed from the top-down by subdivisions based on biogeographic boundaries (for example, in ref. 50) that serve as simple and accessible proxies for biodiversity patterns51. Second, level 5 units (global ecosystem types) are also regional expressions of ecosystem functional groups, but unlike level 4 units they are explicitly linked to local information sources by bottom-up aggregation52 and rationalization of level 6 units from established subglobal ecological classifications. Subglobal classifications, such as those for different countries (see examples for Chile and Myanmar in Supplementary Tables 3.3 and 3.4), are often developed independently of one another, and thus may involve inconsistencies in methods and thematic resolution of units (that is, broadly defined or finely split). Aggregation of level 6 units to broader units at level 5 based on compositional resemblance is necessary to address inconsistencies among different subglobal classifications and produce compositionally distinctive units suitable for global or regional synthesis.Integrating local classifications into the global typology, rather than replacing them, exploits considerable efforts and investments to produce existing classifications, already developed with local expertise, accuracy and precision. By placing national and regional ecosystems into a global context, this integration also promotes local ownership of information to support local action and decisions, which are critical to ecosystem conservation and management outcomes (Supplementary Information, Appendix 3, page 20). These benefits of bottom-up approaches come at the cost of inevitable inconsistencies among independently developed classifications from different regions, a limitation avoided in the top-down approach applied to level 4.Circumscribing upper-level unitsWe formed specialist working groups (terrestrial/subterranean, freshwater and marine) to develop descriptions of the units within the upper levels of the hierarchy, subdividing realms into functional biomes, and biomes into ecosystem functional groups. We used definitions of the hierarchical levels (Supplementary Table 3.1) and the conceptual model of ecosystem assembly (Fig. 1) to maintain consistency in defining the units at each level during iterative discussions within and between the working groups.Working groups agreed on preliminary lists of functional biomes and ecosystem functional groups by considering variation in major drivers along ecological gradients (Figs. 2 and 3 and Supplementary Figs. 3.2 and 3.4) based on published literature, direct experience and expertise of working group members, and consultation with colleagues in their respective research networks. After the workshop, working groups sought recent global reviews of the candidate units and recent case studies of exemplars to shape descriptions of the major groups of ecosystem drivers and properties for each unit. Circumscriptions and descriptions of the units were reviewed and revised iteratively to ensure clear distinctions among units, with a total of 206 reviews of descriptive profiles undertaken by 60 specialists, a mean of 2.4 reviews per profile (Supplementary Table 5.1). The working groups concurrently adapted the generic model of ecosystem assembly (Fig. 1) to represent working hypotheses on salient drivers and ecosystem properties for each ecosystem functional group.Incorporating human influenceVery few of the ecological typologies reviewed in Supplementary Information, Appendix 1 integrate anthropogenic ecosystems in their classificatory frameworks. Anthropogenic influences create challenges for ecosystem classification, as they may modify defining features of ecosystems to a degree that varies from negligible to major transformation across different locations and times. We addressed this problem by distinguishing transformative outcomes of human activity at levels 2 and 3 of the typology from lesser human influences that may be represented either at levels 5 and 6, or through measurements of ecosystem integrity or condition that reflect divergence from reference states arising from human activity.Anthropogenic ecosystems grouped within levels 2 and 3 were thus defined as those created and sustained by intensive human activities, or arising from extensive modification of natural ecosystems such that they function very differently. These activities are ultimately driven by socio-economic and cultural-spiritual processes that operate across local to global scales of human organization. In many agricultural and aquacultural systems and some others, cessation of those activities may lead to transformation into ecosystem types with qualitatively different properties and organizational processes (see refs. 53,54 for cropland and urban examples, respectively). Indices such as human appropriation of net primary productivity55, combined with land-use maps56, offer useful insights into the distribution of some anthropogenic ecosystems, but further development of indices is needed to adequately represent others, particularly in marine, and freshwater environments. Beyond land-use classification and mapping approaches (Supplementary Information, Appendix 1, page 6), a more comprehensive elaboration of the intensity of human influence underpinning the diverse range of anthropogenic ecosystems requires a multidimensional framework incorporating land-use inputs, outputs, their interactions, legacies of earlier activity and changes in system properties17.Where less intense human activities occur within non-anthropogenic ecosystem types, we focused descriptions on low-impact reference states. Therefore, human activities are not shown as drivers in the assembly models for non-anthropogenic ecosystem groups, even though they may have important influences on the contemporary ecosystem distribution. This approach enables the degree and nature of human influence to be described and measured against these reference states using assessment methods such as the Red List of Ecosystems protocol30, with appropriate data on ecosystem change.Indicative distribution mapsFinally, to produce spatially explicit representations of the units at level 3 of the typology (principle 5), we sought published global maps (sources in Supplementary Table 4.1) that were congruent with the concepts of respective ecosystem functional groups. Where several candidate maps were available, we selected maps with the closest conceptual alignment, finest spatial resolution, global coverage, most recent data and longest time series. The purpose of maps for our study was to visualize global distributions. Prior to applications of map data to spatial analysis, we recommend critical review of methods and validation outcomes reported in each data source to ensure fitness for purpose (Supplementary Information, Appendix 4).Extensive searches of published literature and data archives identified high-quality datasets for some ecosystem functional groups (for example, T1.3 Tropical–subtropical montane rainforests; MT1.4 Muddy shorelines; M1.5 Sea ice) and datasets that met some of these requirements for a number of other ecosystem functional groups (see Supplementary Table 4.1 for details). Where evaluations by authors or reviewers identified limitations in available maps, we used global environmental data layers and biogeographic regionalizations as masks to adjust source maps and improve their congruence to the concept of the relevant functional group (for example, F1.2 Permanent lowland rivers). For ecosystem functional groups with no specific global mapping, we used ecoregions50,57,58 as biogeographic templates to identify broad areas of occurrence. We consulted ecoregion descriptions, global and regional reviews, national and regional ecosystem maps, and applied in situ knowledge of participating experts to identify ecoregions that contain occurrences of the relevant ecosystem functional group (for example, T4.4 Temperate woodlands) (see Supplementary Table 4.1 for details). We mapped ecosystem functional groups as major occurrences where they dominated a landscape or seascape matrix and minor occurrences where they were present, but not dominant in landscape–seascape mosaics, or where dominance was uncertain. Although these two categories in combination communicate more information about ecosystem distribution than binary maps, simple spatial overlays using minor occurrences are likely to inflate spatial statistics. The maps are progressively upgraded in new versions of the typology as explicit spatial models are developed and new data sources become available (see ref. 27 for a current archive of spatial data).The classification and descriptive profiles, including maps, for each functional biome and ecosystem functional group underwent extensive consultation, and targeted peer review and revision through a series of four phases described in Supplementary Information, Appendix 5 (pages 2–4). The reviewer comments and revisions from targeted peer review are documented in Supplementary Table 5.1. In all, more than 100 ecosystem specialists have contributed to the development of v2.1 of the typology.LimitationsUneven knowledge of Earth’s biosphere has constrained the delimitation and description of units within the typology. There is a considerable research bias across the full range of Earth’s ecosystems, with few formal research studies evaluating the relative influence of different ecosystem drivers in many of the functional groups, and abiotic assembly filters generally receiving more attention than biotic and dispersal filters. This poses challenges for developing standardized models of assembly for each ecosystem functional group. The models therefore represent working hypotheses, for which available evidence varies from large bodies of published empirical evidence to informal knowledge of ecosystem experts and their extensive research networks. Large numbers of empirical studies exist for some forest functional groups, savannas, temperate heathlands in Mediterranean-type climates, coral reefs, rocky shores, kelp forests, trophic webs in pelagic waters, small permanent freshwater lakes, and others (see references in the respective profiles (Supplementary Information, Appendix 4)). For example, Bond49 reviewed empirical and modelling evidence on the assembly and function of tropical savannas that make up three ecosystem functional groups, showing that they have a large global biophysical envelope that overlaps with tropical dry forests, and that their distribution and dynamics within that envelope is strongly influenced by top-down regulation via biotic filters (large herbivores and their predators) and recurrent disturbance regimes (fires). Despite the development of this critical knowledge base, savannas suffer from an awareness disparity that hinders effective conservation and management59. In other ecosystems, our assembly models rely more heavily on inferences and generalizations of experts drawn from related ecosystems, are more sensitive to interpretations of participating experts, and await empirical testing and adjustment as understanding improves. Empirical tests could examine hypothesized variation in ecosystem properties along gradients within and between ecosystem functional groups and should return incremental improvements on group delineation and description of assembly processes.High-quality maps at suitable resolution are not yet available for the full set of ecosystem functional groups, which limits current readiness for global analysis. The maps most fit for global synthesis are based on remote sensing and environmental predictors that align closely to the concept of their ecosystem functional group, incorporate spatially explicit ground observations and have low rates of omission and commission errors, ‘high’ spatial resolution (that is, rasters of 1 km2 (30 arcsec) or better), and time series of changes. Sixty of the maps currently in our archive27 aligned directly or mostly with the concept of their corresponding ecosystem functional group, while the remainder were based on indirect spatial proxies, and most were derived from polygon data or rasters of 30 arcsec or finer (Supplementary Table 4.1). Maps for 81 functional groups were based either on known records, or on spatial data validated by quantitative assessments of accuracy or efficacy. Therefore, we suggest that maps currently available for 60–80 of the 110 functional groups are potentially suitable for global spatial analysis of ecosystem distributions. Although, a significant advance on broad proxies such as ecoregions, the maps currently available for ecosystem functional groups would benefit from expanded application of recent advances in remote sensing, environmental datasets, spatial modelling and cloud computing to redress inequalities in reliability and resolution. The most urgent priorities for this work are those identified in Supplementary Table 4.1 as relying on indirect proxies for alignment to concept, qualitative evaluation by experts and coarse resolution ( >1 km2) spatial data.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Johnson, N. et al. (eds) Global Soil Biodiversity Atlas (EU, 2016).
Google Scholar 
FAO et al. State of Knowledge of Soil Biodiversity — Status, Challenges and Potentialities (FAO, 2020).
Google Scholar 
Cameron, E. K. et al. Nature Ecol. Evol. 2, 1042–1043 (2018).PubMed 
Article 

Google Scholar 
van den Hoogen, J. et al. Nature 572, 194–198 (2019).PubMed 
Article 

Google Scholar 
Phillips, H. R. P. et al. Science 366, 480–485 (2019).PubMed 
Article 

Google Scholar 
Guerra, C. A. et al. Nature https://doi.org/10.1038/s41586-022-05292-x (2022).Article 

Google Scholar 
Moore, J. C. & de Ruiter, P. C. Energetic Food Webs: An Analysis of Real and Model Ecosystems (Oxford Univ. Press, 2012).
Google Scholar 
Wolters V. et al. Bioscience 50, 1089–1098 (2000).Article 

Google Scholar 
Schimel, J. P. & Schaeffer, S. M. Front. Microbiol. 3, 348 (2012).PubMed 
Article 

Google Scholar 
IPCC. In Climate Change 2022: Mitigation of Climate Change. Contribution of Working Group III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change Impacts, Adaptation, and Vulnerability: Summary for Policymakers (eds Shukla, P. R. et al.) 50 (Cambridge Univ. Press, 2022).
Google Scholar 
Chenu, C. et al. Soil Till. Res. 188, 41–52 (2019).Article 

Google Scholar 
Liang, C., Schimel, J. P. & Jastrow, J. D. Nature Microbiol. 2, 17105 (2017).PubMed 
Article 

Google Scholar 
Hannula, S. E. & Morriën, E. Geoderma 413, 115767 (2022).Article 

Google Scholar  More