Ecology

Alexandratos, N. & Bruinsma, J. World Agriculture Towards 2030/2050. The 2012 Revision (FAO, 2012).Tilman, D., Balzer, C., Hill, J. & Befort, B. L. Global food demand and the sustainable intensification of agriculture. Proc. Natl Acad. Sci. USA 108, 20260–20264 (2011).CAS 
Article 

Google Scholar 
Mueller, N. D. et al. Closing yield gaps through nutrient and water management. Nature 490, 254–257 (2012).CAS 
Article 

Google Scholar 
Chen, X. et al. Producing more grain with lower environmental costs. Nature 514, 486–489 (2014).CAS 
Article 

Google Scholar 
Fan, M. S. et al. Improving crop productivity and resource use efficiency to ensure food security and environmental quality in China. J. Exp. Bot. 63, 13–24 (2012).CAS 
Article 

Google Scholar 
Godfray, H. C. J. et al. Food security: the challenge of feeding 9 billion people. Science 327, 812–818 (2010).CAS 
Article 

Google Scholar 
Foley, J. A. et al. Solutions for a cultivated planet. Nature 478, 337–342 (2011).CAS 
Article 

Google Scholar 
Porter, J. R. et al. Food Security and Food Production Systems (Cambridge Univ. Press, 2014).Ray, D. K. & Foley, J. A. Increasing global crop harvest frequency: recent trends and future directions. Environ. Res. Lett. 8, 044041 (2013).Article 

Google Scholar 
Lal, R. Restoring soil quality to mitigate soil degradation. Sustainability 7, 5875–5895 (2015).Article 

Google Scholar 
Wall, D. & Six, J. Give soils their due. Science 347, 695 (2015).CAS 
Article 

Google Scholar 
Ray, D. K. et al. Climate variation explains a third of global crop yield variability. Nat. Commun. 6, 5989 (2015).CAS 
Article 

Google Scholar 
Battisti, D. S. & Naylor, R. L. Historical warnings of future food insecurity with unprecedented seasonal heat. Science 323, 240–244 (2009).CAS 
Article 

Google Scholar 
Nelson, G. C. et al. Climate Change: Impact on Agriculture and Costs of Adaptation (International Food Policy Research Institute, 2009).Challinor, A. J., Koehler, A. K., Ramirez-Villegas, J., Whitfield, S. & Das, B. Current warming will reduce yields unless maize breeding and seed systems adapt immediately. Nat. Clim. Change 6, 954–958 (2016).Article 

Google Scholar 
Zhao, C. et al. Temperature increase reduces global yields of major crops in four independent estimates. Proc. Natl Acad. Sci. USA 114, 9326–9331 (2017).CAS 
Article 

Google Scholar 
Schlenker, W., Hanemann, M. & Fisher, A. Will US agriculture really benefit from global warming? Accounting for irrigation in the hedonic approach. Am. Econ. Rev. 95, 395–406 (2005).Article 

Google Scholar 
Piao, S. L. et al. The impacts of climate change on water resources and agriculture in China. Nature 467, 43–51 (2010).CAS 
Article 

Google Scholar 
Ray, D. K. et al. Climate change has likely already affected global food production. PLoS ONE 14, e0217148 (2019).CAS 
Article 

Google Scholar 
Ramankutty, N. et al. The global distribution of cultivable lands: current patterns and sensitivity to possible climate change. Glob. Ecol. Biogeogr. 11, 377–392 (2002).Article 

Google Scholar 
Rosenzweig, C. et al. Assessing agricultural risks of climate change in the 21st century in a global gridded crop model intercomparison. Proc. Natl Acad. Sci. USA 111, 3268–3273 (2014).CAS 
Article 

Google Scholar 
Lobell, D. B. & Burke, M. B. On the use of statistical models to predict crop yield responses to climate change. Agr. For. Meteorol. 150, 1443–1452 (2010).Article 

Google Scholar 
Auffhammer, M. & Schlenker, W. Empirical studies on agricultural impacts and adaptation. Energy Econ. 46, 555–561 (2014).Article 

Google Scholar 
Folberth, C. et al. Uncertainty in soil data can outweigh climate impact signals in global crop yield simulations. Nat. Commun. 7, 11872 (2016).CAS 
Article 

Google Scholar 
Asseng, S. et al. Uncertainty in simulating wheat yields under climate change. Nat. Clim. Change 3, 827–832 (2013).CAS 
Article 

Google Scholar 
Basso, B. et al. Soil organic carbon and nitrogen feedbacks on crop yields under climate change. Agr. Environ. Lett. 3, 180026 (2018).Mϋller, C. et al. Implication of climate mitigation for future agricultural production. Environ. Res. Lett. 10, 125004 (2015).Article 

Google Scholar 
IPCC Climate Change 2022: Impacts, Adaptation, and Vulnerability (eds Pörtner, H. O. et al.) (Cambridge Univ. Press, 2022).Zhang, W. et al. Closing yield gaps in China by empowering smallholder farmers. Nature 537, 671–674 (2016).CAS 
Article 

Google Scholar 
Cui, Z. L. et al. Pursuing sustainable productivity with millions of smallholder farmers. Nature 555, 363–368 (2018).CAS 
Article 

Google Scholar 
Knapp, S. & van der Heijden, M. G. A. A global meta-analysis of yield stability in organic and conservation agriculture. Nat. Commun. 9, 3632 (2018).Article 
CAS 

Google Scholar 
Müller, C. et al. Global Gridded Crop Model evaluation: benchmarking, skills, deficiencies and implications. Geosci. Model Dev. 10, 1403–1422 (2017).Jamieson, P. D., Porter, J. R. & Wilson, D. R. A test of the computer simulation model ARC-WHEAT on wheat crops grown in New Zealand. Field Crops Res. 27, 337–350 (1991).Article 

Google Scholar 
Warszawski, L. et al. The inter-sectoral impact model intercomparison project (ISI–MIP): project framework. Proc. Natl Acad. Sci. USA 111, 3228–3232 (2014).CAS 
Article 

Google Scholar 
Xiong, W. et al. The Impacts of Climate Change on Chinese Agriculture—Phase II National Level Study Final Report (AEA Group, 2008).Liu, B. et al. Similar estimates of temperature impacts on global wheat yield by three independent methods. Nat. Clim. Change 6, 1130–1136 (2016).Article 

Google Scholar 
Tao, F. et al. Global warming, rice production, and water use in China: developing a probabilistic assessment. Agr. For. Meteorol. 148, 94–110 (2008).Article 

Google Scholar 
Xiong, W. et al. Different uncertainty distribution between high and low latitudes in modelling warming impacts on wheat. Nat. Food 1, 63–69 (2020).Article 

Google Scholar 
Fernandez-Illescas, C. P., Porporato, A., Laio, F. & Rodriguez-Iturbe, I. The ecohydrological role of soil texture in a water-limited ecosystem. Water Resour. Res. 37, 2863–2872 (2001).Article 

Google Scholar 
Wang, E. L. et al. Capacity of soils to buffer impact of climate variability and value of seasonal forecasts. Agr. For. Meteorol. 149, 38–50 (2009).Article 

Google Scholar 
Vereecken, H. et al. Modeling soil processes: review, key challenges, and new perspectives. Vadose Zone J. 15, 1–57 (2016).Myers, R. J. K. et al. in The Biological Management of Tropical Soil Fertility (eds Woomer, P.I. & Swift, M.J.) Ch. 4 (Wiley, 1994).Smith, P. & Gregory, P. J. Climate change and sustainable food production. P. Nutr. Soc. 72, 21–28 (2013).Article 

Google Scholar 
Khasawneh, F. E., Sample, E. C. & Kamprath, E. J. The Role of Phosphorus in Agriculture (American Society of Agronomy, 1980).FAOSTAT (Statistics Division of the Food and Agriculture Organization of the United Nations, 2006); http://www.fao.org/faostat/en/#homeFan, M. S. et al. Plant-based assessment of inherent soil productivity and contributions to China’s cereal crop yield increase since 1980. PLoS ONE 8, e74617 (2013).CAS 
Article 

Google Scholar 
Liu, X. & Chen, F. Farming System in China (China Agriculture Press, 2005).Chen, X. P. in Fertilization Technology Highlights, (ed. Zhang, F. S) Ch. 6 (Chinese Agricultural Univ. Press, 2006).Zhang, F. et al. Integrated nutrient management for food security and environmental quality in China. Adv. Agron. 116, 1–40 (2012).CAS 
Article 

Google Scholar 
Bünemann, E. K. et al. Soil quality—a critical review. Soil Biol. Biochem. 120, 105–125 (2018).Article 
CAS 

Google Scholar 
National Soil Survey Office. Chinese Soil (China Agriculture Press, 1998) .Jiang, R. F. & Cui, J. Y. in Fertilization Technology Highlights, (ed. Zhang, F. S.) Ch. 5 (China Agricultural Univ. Press, 2006).Cramer, W. P. & Solomon, A. M. Climatic classification and future global redistribution of agricultural land. Clim. Res. 3, 97–110 (1993).Article 

Google Scholar 
Elith, J., Leathwick, J. R. & Hastie, T. A working guide to boosted regression trees. J. Anim. Ecol. 77, 802–813 (2008).CAS 
Article 

Google Scholar 
Friedman, J. H. Stochastic gradient boosting. Comput. Stat. Data 38, 367–378 (2002).Article 

Google Scholar 
Buston, P. M. & Elith, J. Determinants of reproductive success in dominant pairs of clownfish: a boosted regression tree analysis. J. Anim. Ecol. 80, 528–538 (2011).Article 

Google Scholar 
Friedman, J. H. & Meulman, J. J. Multiple additive regression trees with application in epidemiology. Stat. Med. 22, 1365–1381 (2003).Article 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2018).Kuhn, M. & Johnson, K. Applied Predictive Modeling (Springer, 2013).Yang, J. M., Yang, J. Y., Liu, S. & Hoogenboom, G. An evaluation of the statistical methods for testing the performance of crop models with observed data. Agric. Syst. 127, 81–89 (2014).Article 

Google Scholar 
Loague, K. & Green, R. E. Statistical and graphical methods for evaluating solute transport models: overview and application. J. Contamin. Hydro. 7, 51–73 (1991).CAS 
Article 

Google Scholar 
Akinremi, O. O. et al. Evaluation of LEACHMN under Dryland conditions. I. Simulation of water and solute transport. Can. J. Soil Sci. 85, 223–232 (2005).Article 

Google Scholar 
Palosuo, T. et al. Simulation of winter wheat yield and its variability in different climates of Europe: a comparison of eight crop growth models. Eur. J. Agron. 35, 103–114 (2011).Article 

Google Scholar 
Deng, N. et al. Closing yield gaps for rice self-sufficiency in China. Nat. Commun. 10, 1725 (2019).Article 
CAS 

Google Scholar 
Correndo, A. A. et al. Assessing the uncertainty of maize yield without nitrogen fertilization. Field Crops Res. 260, 107985 (2021).Article 

Google Scholar 
Rattalino Edreira, J. I. et al. Spatial frameworks for robust estimation of yield gaps. Nat. Food 2, 773–779 (2021).Article 

Google Scholar 
Tilman, D., Reich, P. B. & Knops, J. M. H. Biodiversity and ecosystem stability in a decade-long grassland experiment. Nature 441, 629–632 (2006).CAS 
Article 

Google Scholar 
Kuhn, M. Building predictive models in R using the caret package. J. Stat. Softw. 28, 1–26 (2008).Article 

Google Scholar 
IPCC Climate Change 2014: Climate Change: Synthesis Report (eds Core Writing Team, Pachauri, R. K. & Meyer L. A.) (IPCC, 2014).van Vuuren, D. P. et al. The representative concentration pathways: an overview. Clim. Change 109, 5–31 (2011).Article 

Google Scholar 
Hempel, S., Frieler, K., Warszawski, L., Schewe, J. & Piontek, F. A trend-preserving bias correction—the ISI-MIP approach. Earth Syst. Dynam. 4, 219–236 (2013).Article 

Google Scholar 
Chen, H., Sun, J., Lin, W. & Xu, H. Comparison of CMIP6 and CMIP5 models in simulating climate extremes. Sci. Bull. 65, 1415–1418 (2020).Article 

Google Scholar 
China Agriculture Yearbook (China Agriculture Press, 2005). More

Chloroplast and repetitive nuclear DNA enrichment in the sedaDNA extractsTo the best of our knowledge, we generated the first large-scale target enriched dataset using sedaDNA extracted from sediments of multiple lakes. Sequencing of two datasets produced 325.5 million (M) quality-filtered paired-end DNA sequences. The first target enriched dataset, targeting both the chloroplast and a set of nuclear genes of Larix on 64 sedaDNA extracts and 19 negative controls from seven lake sediment records resulted in 324 M quality-filtered paired-end sequences. The second target enriched dataset, targeting only the set of nuclear genes of Larix on four samples and two negative controls from an additional lake (Lake CH12) resulted in 1.5 M sequences. Quality-filtering of an additional published target enriched dataset29, targeting the Larix chloroplast genome on the same CH12 samples as applied for the second dataset, added another 54 M sequences.For the chloroplast enrichment, 390 thousand (K) sequences (1%) were classified as Larix at the genus or species level. The average coverage of bait regions was 19% at a mean sequence depth of 0.8. Sequencing of 19 library and extraction blank (negative control) samples resulted in 597 K paired-end sequences, of which 58% quality-filtered and deduplicated sequences remained. Of these, 38% were classified, with 0.03% of them (463 sequences) corresponding to the genus Larix. Negative controls from library preparation resulted in no to very few (0 to 5) sequences mapping to the Larix chloroplast reference genome. Negative controls from DNA extractions, which were in several cases pooled to one library, showed a low number of sequences mapped to Larix (0 to 94 sequences, except 237 sequences in one case). Excluding all sequences in negative controls from the sample analysis had no impact on the patterns resulting from the analysis of sample data. Detailed results and evaluation of negative controls are included in the Supplementary Information (Fig. S5) and Supplementary Data 1 and 2. Samples of all lake records with sufficient sequence coverage showed damage patterns typical of ancient DNA (see Supplementary Data 3).These results are comparable to the results obtained by Schulte et al.29, where 36% of quality-filtered sequences were classified as Viridiplantae with 9% assigned to Larix. In contrast to29, we raised the confidence threshold of taxonomic classification (a parameter defining the number of k-mers needed to produce a match against a taxon in the database), which drastically reduced the number of classified sequences, but increased the confidence in the analysis36.To analyze the enrichment obtained by the nuclear gene bait set, taxonomic classification was repeated using a plant genome database including available Pinaceae genomes. The classification resulted in 716 K sequences assigned to Larix, increasing the previous results by 325 K sequences. However, almost no sequences were mapped against the targeting baits (a maximum of five sequences for some samples). A closer inspection of unmapped sequences assigned to Larix revealed a high content of repetitive DNAs. More specifically, taxonomically classified Larix sequences could be assembled to EulaSat1, the most abundant satellite repeat in the nuclear genome of Larix32,37. This short repeat with a 173 bp long motif is arranged in large arrays of tandemly repeated motifs and is exclusively present in larches32. Analysis of modern L. sibirica, and L. gmelinii (western and eastern range) genomes reveals that EulaSat1 occurs in all species, contributing to 0.62% (L. sibirica), 2.52% (western range L. gmelinii), and 2.39% (eastern range L. gmelinii), of the genomes, respectively (Fig. S2). A comparison of the sequence proportions mapping to the repeat motif in the different datasets of Lake CH12 showed a specific enrichment of the repeat motif by the nuclear gene hybridization probe set (Fig. S3).In total, 17 K sequences mapped to the repeat motif of EulaSat1. The abundance of all sequences mapped per sample is in agreement with the abundance of sequences mapped to the chloroplast genome, confirming the general history of forest development (Fig. 2). Analysis of the nucleotide frequencies in the repeat motif showed a high constancy over all samples (Fig. S4). This suggests high conservation of the EulaSat1 motif in Siberian larches over time and space. Although satellite repeats are reported to have a high sequence turnover, for larches it has been shown that repeat profiles between two geographically well-separated species—the European larch (L. decidua) and the Japanese larch (L. kaempferi)—are very similar32. The main satellite in all larches, EulaSat1, is believed to have greatly multiplied after the split of Larix from Pseudotsuga32. Given the ongoing hybridization between the three Siberian larch species, it is not surprising to find a consistent pattern of nucleotide frequencies in all samples.Fig. 2: Comparison of target enrichment with available DNA metabarcoding and pollen datasets.From left to right: Larix-classified sequence counts mapping to (1) the Larix chloroplast and (2) the EulaSat1 satellite repeat motif, (3) percentage of Larix counts in metabarcoding data, (4) percentage of Larix pollen in pollen assemblages. All data from this study, except metabarcoding data from lakes CH1213 and Bolshoye Shchuchye55 and all pollen data except for several samples of Lake Kyutyunda which were produced in this study56,57,71. Pollen data of Lake Lama and the Holocene part of Lake Kyutyunda are based on parallel sediment cores PG1111 and PG2022, respectively. No available data are marked with crosses, asterisk marks a single Larix pollen grain found in the Bolshoye Shchuchye sediments.Full size imageOff-target sequences in target enriched datasets have already been demonstrated to be useful for the analysis of high-copy DNA such as ribosomal DNA or plastomes34,38,39. A recent study on five modern sedges showed that target enriched sequencing data originally targeting a set of gene exons can also be used to study the repetitive sequence fraction and even infer phylogenetic relationships based on repetitive sequence abundance35. Another study showed that also sequence similarities between homologous repeat motifs can be used to reconstruct phylogenetic relationships among closely related taxa40,41. In the case of Larix satellite EuLaSat1 in our study, no change in nucleotide frequencies, neither related to locations nor in time, could be detected. However, our results show that the off-target fraction in target enriched sedaDNA datasets can hold valuable information and that repeat motifs in more diverse taxon groups could even be a target for enrichment. Specifically enriching for repeat motifs in sedaDNA extracts could enable the study of satellite repeat evolution as well as giving additional information on species abundance and phylogeography.In the two target enriched datasets, sequences taxonomically classified to the genus Larix and mapping to the chloroplast and to the repeat sequence, respectively, show similar patterns of abundance (see Fig. 2). Compared with published metabarcoding and pollen data from the same locations, the Larix abundance patterns can be globally reproduced, underpinning the notion that sequence abundances in target enriched data can be used as good estimates of plant abundances. For older parts of the lake records, target enriched data show Larix where metabarcoding data were unable to detect a clear signal (see Fig. 2, lakes Billyakh, Bolshoye Shchuchye, Kyutyunda, and Lama). This shows that target enrichment is superior to metabarcoding when analyzing one taxonomic group in-depth, as it is less prone to errors by DNA degradation, which can impede primer binding if the molecule becomes too short. Also, independent of age, rare taxa mostly need multiple PCR replicates to be detected by metabarcoding42,43. Target enrichment, however, is more sensitive in identifying one focal taxon group, as the total target length can be much larger (e.g., a complete organellar genome) than for metabarcoding, and the DNA damage patterns are put to use to authenticate ancient DNA. Also, it is limited by molecule length only by the applied threshold in the bioinformatic analysis, for which we used 30 base pairs (bp) as opposed to a minimum of 85 bp molecule length for the Larix metabarcoding marker (for the plant-specific trnL g/h marker44). Similarly, compared to traditional pollen analysis, target enrichment is more accurate at tracing a specific target group, as it is not dependent on pollen productivity. Especially in the case of Larix, pollen productivity is low and preservation poor, resulting in rare findings of its pollen in the sediments22,45. This could explain why for Lake Bolshoye Shchuchye, only a single Larix pollen grain was retrieved throughout the core, whereas target enrichment and metabarcoding show a strong signal in the Holocene sediments (last ~12 ka BP). Target-enriched data also records signals in MIS 2 sediments, however, sequence counts are extremely low, and as it is the only record, where both of the other proxies fail to report a signal, it should be interpreted with caution.A wider pre-glacial distribution of L. sibirica
Chloroplast genomes of L. gmelinii and L. sibirica differ at 157 positions, which can be used to differentiate species in target enriched sedaDNA29. Here, we applied this approach to lake sediment records, which are distributed across Siberia (Fig. 1) and have time ranges back to MIS3, and thereby were able to track species composition in space and time for wide parts of the species ranges.In lakes Billyakh and Kyutyunda, ca. 1500 km east of L. sibirica current range (Fig. 1), we found evidence for a wider distribution of L. sibirica around 32 and 34 ka BP in MIS3 (Fig. 3). Billyakh is situated in the western part of the Verkhoyansk Mountains, and Kyutyunda on the Central Siberian Plateau. Both lakes have low counts of Larix DNA sequences in their oldest samples dated to 51 ka BP (Billyakh) and 38 ka BP (Kyutyunda) with variants of L. gmelinii, but there is a sudden rise in variants attributed to L. sibirica at 34 ka BP (Billyakh) and 32 ka BP (Kyutyunda), which persists in the following samples, but strongly decreases in younger samples (Fig. 3). The rise in the L. sibirica DNA sequence variants coincides with a peak in sequence counts for Lake Kyutyunda. These signals suggest a rapid invasion of L. sibirica into the ranges of L. gmelinii in climatically favorable times and a local depletion or extinction of L. sibirica during the following harsher climates. Lake Billyakh pollen data suggest a moister and warmer climate around 50–30 ka BP than in the latter part of the Last Glacial associated with the MIS3 Interstadial in Siberia46.Fig. 3: Percentage and sequence counts at variable positions along Larix chloroplast genome assigned to species.Left: Alignment of Larix-classified DNA sequences against the chloroplast genome at the 157 variable positions between the species. For each position, the percentage of sequences assigned to a single species is displayed. Each row represents one sample named according to the calibrated age before present. Gray background indicates no coverage at the respective position. Right: Total number of sequences assigned to each of the species per sample.Full size imageStrong support for a wider pre-glacial distribution of L. sibirica comes from genetic analyses which show that it is genetically close to L. olgensis, today occurring on the Korean Peninsula and adjacent areas of China and Russia27,47. It is assumed that the L. sibirica-L. olgensis complex used to share a common range, which was disrupted and displaced when the better cold-adapted L. gmelinii expanded south and southwest during the more continental climatic conditions of the Pleistocene47,48. Furthermore, modern and ancient genetic studies suggest that the L. sibirica zone was recently invaded by L. gmelinii from the east in the hybridization zone of the species, as the climate cooled after the mid-Holocene thermal maximum13,23. Today, pure stands of L. sibirica do not form a continuous habitat, but occur in netted islands5 and morphological features of L. sibirica can be found in populations of L. gmelinii located at least a hundred kilometers east of the closest L. sibirica populations49. Macrofossil findings of L. sibirica in Scandinavia dated to the early Holocene, point to the capability of rapid long-distance jump dispersal of this species50. Fossil L. sibirica cones dated to the end of the Pliocene and in the Pleistocene have also been found far east of its current range in several river valleys including Kolyma, Aldan, and Omolon, and even in the basin of the Sea of Okhotsk9. These indicate long-distance seed dispersal by rivers which may also have assisted in successful establishment since the active-layer depth is deeper close to rivers51,52. As mentioned earlier, L. sibirica is sensitive to permafrost and waterlogged soils. A warmer phase with a deeper thawed layer above the permafrost could have enabled L. sibirica to spread and establish in regions that today are part of the geographic range of L. gmelinii, as L. sibirica is reported to have higher growth rates than L. gmelinii13.
Larix gmelinii formed northern LGM refugia across SiberiaThe possible survival of Larix in high latitude glacial refugia during the LGM is still under discussion4,53 although more and more evidence is reported in favor of the existence of such refugia17,20,21. The question of which of the Larix species formed these populations has hitherto been unanswered, as both pollen and established metabarcoding markers are not able to distinguish between species in the genus Larix, and findings of fossilized cones identifiable to species are rare. By enriching sedaDNA extracts for chloroplast genome sequences, we are, to the best of our knowledge, for the first time, able to distinguish between L. sibirica and L. gmelinii in glacial refugial populations.From Lake Lama, located at the western margin of the Putorana Plateau (Taymyr Peninsula), we obtained a continuous record extending from 23 ka BP to today with varying sequence counts with minima around 18–17 ka BP and 13 ka BP. All samples prior to the Holocene show variations predominantly assigned to L. gmelinii (Fig. 3). Our results suggest a local survival of L. gmelinii in the vicinity of Lake Lama throughout the LGM, which is supported by low numbers of Larix pollen detected through this period. Both target enriched sequence data and pollen indicate an increase from ca. 11 ka BP54. Sparse Larix pollen in the bottom part of the record could be an indication of a possible refugial population (Fig. 2; ref. 54).In Bolshoye Shchuchye, the westernmost lake of the study, situated in the Polar Ural Mountains, all Pleistocene samples show similarly a dominance of L. gmelinii sequence variations (Fig. 3). However, sequence counts for some samples are extremely low and samples from 18 and 10 ka BP had so low counts of mapped DNA sequences that none of the variable positions between the species was covered. Although sequences mapped to the satellite repeat of Larix also showed a Pleistocene signal, this was not repeated in pollen or metabarcoding (Fig. 2) which instead indicates a treeless arctic-alpine flora for the late Pleistocene55,56. Especially for the sample of 20.4 ka, Larix sequence counts are extremely low and new investigations would be needed to confirm a local presence of Larix during the LGM.The record of Lake Billyakh situated in the western Verkhoyansk Mountain Range, likewise shows extremely low counts of sequences mapped to the reference for a range of samples with no sequences covering the studied variable sites (45, 42, and 15 ka BP, 11–56 sequences mapped to non-variable sites). However, the pollen record for the same core shows a quasi-continuous record of Larix with a gap only occurring during the early LGM46 (25–22 ka BP, Fig. 2). Considering the known short-distance dispersal ability and poor preservation of Larix pollen, this strongly supports the presumed existence of a local glacial refugium at Lake Billyakh during that time20. Our samples also show a low but steady presence of Larix throughout the rest of the record, thus making glacial survival probable. The sample closest to the LGM (24 ka BP) indicates a clear dominance of L. gmelinii type variations.The only exception to this general pattern is the record from Lake Kyutyunda, which is located on the Central Siberian Plateau west of the Verkhoyansk Mountain Range. In this record, LGM samples have extremely low counts but show variations assigned to L. sibirica and not to L. gmelinii as in the other lakes. In addition, the preceding sample dated to the MIS3 interstadial shows L. sibirica variation. A possible explanation is that relics of L. sibirica survived during the LGM, but were unable to spread after climate warming, possibly due to genetic depletion or later local extinction. The presence of reworked sediment material can also not be excluded, as suggested by reworked pollen in the record57.In conclusion, our data show almost exclusively L. gmelinii variation for samples covering the most severe LGM climate conditions. This is in agreement with the ecological characteristics describing the species as adapted to extreme cold. In contrast to L. sibirica, it can grow in dwarf forms and propagate clonally and potentially survive thousands of years of adverse climatic conditions58.Postglacial colonization history—differences among larch speciesOf great interest in the Larix history is not only the location and extent of possible high latitude glacial refugia but also if and to what extent these refugia contributed to the recolonization of Siberia after the LGM. Northern refugial populations could have functioned as kernels of postglacial population spread and recolonization, or spreading could have been driven by populations that survived in southern refugia. There are only a few studies on modern populations that report evidence for possible recolonization scenarios of Larix23,27,28. Here, we show that patterns differ between L. sibirica and L. gmelinii.In the western part of our study region, two lakes are situated in the current distribution range of L. sibirica (Figs. 1, 4): Lake Bolshoye Shchuchye in the Polar Ural Mountains and Lake Lama on the Taymyr Peninsula. Despite this, both lakes show L. gmelinii for all Pleistocene samples, and a strong signal of L. sibirica variants only in the Holocene samples, with ages of 5.1 ka BP in Lake Bolshoye Shchuchye and 9.7 ka BP in Lake Lama (Fig. 3). The peak in L. sibirica also coincides with a peak of sequence counts in the respective sample, with a Larix pollen peak in Lake Lama sediments54, and metabarcoding for Lake Bolshoye Shchuchye55. This points to a migration of L. sibirica in its current northern area of northern distribution in the course of climate warming during the early Holocene, whereas glacial refugial populations were consisting of L. gmelinii. Although the local survival of L. gmelinii around Lake Bolshoye Shchuchye remains uncertain due to extremely low sequence counts, it is clear that L. sibirica did not form a refugial population at this site.Fig. 4: Percentage of DNA sequences assigned to references displayed on the geographical locations of the lakes investigated.Samples in the same time frame are averaged. Lake names and current species ranges are annotated in the middle plots. Colors indicate current species distribution (adapted from Semerikov and Lascoux72). The base map is done with ggmap73, map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.Full size imageA range-wide genetic study of L. sibirica analyzing chlorotypes and mitotypes of individuals23 found strong indications for rapid colonization of the West Siberian Plains from populations originating from the foothills of the Sayan Mountains in the south, close to the border of Mongolia, with only limited contribution from local populations. According to our results, the local populations could have been L. gmelinii populations, while the rapid invasion could have been L. sibirica.In the eastern range of the study region, in the current range of L. gmelinii, namely at lakes Emanda, Satagay, and Malaya Chabyda, genetic variations throughout the records are less pronounced. Of the three lake records, only that from Lake Emanda reaches back beyond the LGM, but with a sampling gap for the time of the LGM. Therefore, it remains uncertain whether populations survived the LGM locally, or whether they were invaded or replaced by populations coming from the south with Holocene warming. The restricted variations throughout the record, however, hint at stable populations, which is supported by scarce pollen findings (Fig. 2).Our data suggest that postglacial recolonization of L. sibirica was not started from high latitude glacial refugia, but from southern populations. In contrast, northern glacial populations of L. gmelinii could have potentially enhanced rapid dispersal after the LGM in their current area of distribution.Environment likely plays a more important role than historical factorsThe current boundaries of boreal Larix species arranged from west to east suggest a possible strong influence of the historical species distribution on the current distribution, whereas the gradient of increasing continental climate towards the east assumes a strong influence on the environment. By tracking species distribution in the past, spanning the time of the strongly adverse climate of the LGM, we can give hitherto unprecedented insights into species distribution history.Several lines of evidence suggest a strong influence of the environment on species distribution: (1) Signals for L. sibirica appeared in its current area of distribution as late as the Holocene warming, whereas cold Pleistocene samples are dominated by L. gmelinii type variation; (2) in lakes far east of its modern range, signals of variation typical for L. sibirica coincide with peaks in sequence counts (29 ka BP, Lake Billyakh; 32 ka BP Lake Kyutyunda), which point to more forested vegetation around the lakes and consequently a more favorable climate at that time; and (3) samples dated to the severely cold LGM are dominated by variations of the L. gmelinii type.This is in accordance with the different ecological characteristics described for the species. L. sibirica is sensitive to permafrost and only occurs outside of the zone of continuous permafrost5. In addition, L. sibirica achieves substantially higher growth rates and longer growth periods than L. gmelinii9,13 and can also produce more than twice as many seeds5. This potentially gives L. sibirica the ability to quickly react to climate change and outcompete the other species when the climate becomes more favorable.In contrast, L. gmelinii is adapted to extremely low soil and air temperatures and is able to grow on permafrost with very shallow thaw depths. It’s distribution almost completely coincides with continuous permafrost5, and even a restriction to permafrost areas is discussed as it does not grow well in field trials on warmer soils or where there is a small temperature gradient between air and soil9. Due to this ecology, L. gmelinii is more likely to survive in a high latitude refugium, even during the severe continental climate of the LGM, which was most probably connected to continuous permafrost of low active-layer depths.A study combining mitochondrial barcoding on sedaDNA and a modeling approach on Larix distribution in the Taymyr region around Lake CH12 concluded that the distributions of L. gmelinii and L. sibirica are most strongly influenced by stand density and thus by competition between the species, with L. gmelinii outcompeting L. sibirica at high stand densities13. As our study includes sediment cores reaching further back in time, we see a different trend. Instead of L. gmelinii, it was L. sibirica, which dominated samples with high sequence counts, suggesting high stand density and a more favorable climate. A possible explanation for the different outcomes is the use of different organelle genomes. Epp et al.13 used a marker representing the mitochondrial genome, which is known to introgress more rapidly and as a consequence might show a long past species history59,60.Our findings have potentially important implications for the projections of vegetation-climate feedback. A warming climate in conjunction with a greater permafrost thaw depth could enable the replacement of L. gmelinii by L. sibirica. In contrast to L. gmelinii, L. sibirica is not known to stabilize permafrost thus potentially further promoting permafrost thaw and with it the release of greenhouse gases, creating positive feedback on global warming11. On the other hand, the substantially higher growth rates of L. sibirica in comparison to L. gmelinii would increase carbon sequestration, thus mitigating global warming13. This shows the importance of understanding species-specific reactions to climate change, which can result in great shifts in distribution. Target enrichment applied on sedaDNA is able to reveal the impact of past climate change on populations and the increasing availability of modern reference genomes will further enhance its value of information. More

TRY meta-­analysisWe used the TRY plant trait database27 to identify traits that show systematic differences between the tree and liana growth forms, as a way to narrow the scope of the rest of the analysis. We chose traits to represent major trade­offs within the “economic spectrum” framework, which places plants along a spectrum of strategies from acquisitive, fast return on investment to conservative, slow return on investment according to key functional trait values30. We narrowed traits to those that had observations for at least four tree and liana species. We then compiled our dataset using the following steps during November and December 2019. For each trait, we downloaded the dataset for all species available globally and averaged the observations of the trait to the species level to avoid statistical biases introduced in our growth form comparison due to a high density of observations in a few commercially valuable species. We matched the species ID number with the most frequently used growth form identifier using the TRY “growth form” trait and kept the species with the most frequent identifier of “tree,” “liana,” or “woody vine.” We subsetted the data to keep only species with a majority of observations ascribed to the tree and liana growth forms (i.e., no herbaceous species, ferns, etc.), resulting in observations for 44,222 total species. Finally, we filtered the dataset of 44,222 species by hand to remove species misclassified as trees or lianas; species occurring entirely in temperate to boreal biomes; species from all gymnosperm lineages except the order Gnetales; and entries for taxonomic classifications broader than the genus level (e.g., taxonomic families). We found that hydraulic functional traits in the TRY database (i.e., Ks,max and P50) show systematic differences between growth forms (Supplementary Fig. 1; Supplementary Tables 3 and 4), while there is mixed evidence for differences in the acquisitiveness of trees and lianas in terms of stem anatomical traits (Supplementary Fig. 1; Supplementary Tables 3 and 4) and leaf functional traits (Supplementary Fig. 6; Supplementary Tables 3 and 4), and no evidence of differences between tropical trees and lianas with respect to root functional traits (Supplementary Fig. 7; Supplementary Tables 3 and 4).Extended meta­-analysisWe conducted an additional literature search to supplement the hydraulic trait observations from the TRY database. The additional literature search served two purposes: (1) to fill a major gap identified during our TRY analysis in terms of liana trait observations, and (2) to address the methodological inconsistency of measuring Ks,max and P50 on liana branches shorter than the longest vessel, which incorrectly measures Ks,max and P50 without accounting for end wall resistivity59,60.We conducted a literature search using Web of Science and Google Scholar. We searched the following phrases in combination with “liana:” “hydraulic conductivity,” “hydraulic trait,” “hydraulic efficiency,” and “hydraulic K.” Of the literature we found, we kept only the studies that met the following criteria: (1) reported Ks,max measurements for lianas, (2) measured Ks,max instead of computing Ks,max from xylem conduit dimensions, (3) measured Ks,max on sunlit, terminal branches of mature individuals or saplings, and (4) measured Ks,max on a branch longer than the longest vessel. We considered the authors to have used a branch length longer than maximum vessel length if the authors measured or reported maximum vessel length for the species and a longer branch was used. Because the best methodological practice for measuring P50, especially in species with long vessels, is currently a matter of debate, we additionally removed all observations of P50  > ­0.75. This filtering was performed to reduce the probability that falsely high (i.e., less negative) P50 values were retained in our analysis because of improper measurement technique and is consistent with the P50 filtering performed by Trugman et al.61. Improper measurement technique is a particular concern for lianas, whose wide and long vessels require cautious implementation of the traditional measurement techniques developed for trees. We note that retaining all liana P50 observations (i.e., not filtering out observations  > −0.75) results in a significant difference between trees and lianas (Mann­–Whitney test statistic = 1029, ntree = 61, nliana = 46, p  More

Experimental designTo test whether an 8-year treatment of recurrent summer droughts would change biodiversity effects and species interactions of grassland plants when facing a subsequent drought event, we grew ambient- vs. drought-selected plants of 12 species in a glasshouse. The plants were grown from seeds collected from 40 plots (Supplementary Data 2) under 8-year treatments of yearly summer droughts vs. ambient precipitation in a biodiversity field experiment in Jena, Germany11,41.The Jena Experiment was established in 2002 using a common seed pool of 60 grassland species, with 80 (20times 20,{{{{{rm{m}}}}}}) large plots of species richness levels of 1, 2, 4, 8, 16, and 60 species40. Most of the species are perennial and capable of outcrossing (Supplementary Table 1). The Jena Drought Experiment11,41 was initiated in 2008. Two (1times 1,{{{{{rm{m}}}}}}) subplots were set within each large plot, designated as either drought treatment or ambient control. For the drought treatment, rainout shelters were set up to exclude natural rainfall in mid-summer for 6 weeks. The ambient control treatment got the same shelter construction but rain water was reapplied to not confound the results with artifacts from the shelter60. We repeatedly harvested the aboveground biomass per year, once before and once after the summer drought treatment11,41. The design of the Jena (Drought) Experiment did not allow the exclusion of cross pollination or gene flow between subplots or large plots in the field. Such gene flow may have reduced the possibility for genetic differentiation and for the observed effect sizes of the selection treatment23. We collected seeds from drought and control subplots throughout the 2016 growing season (Fig. 1). We obtained seeds of 17 species, but only used 12 of them, because the other five species had either few seeds or low germination rates. Seeds per species per selection treatment were collected from 4 to 23 (interquartile range: [8.50, 17.00]) maternal plants distributed across 2–10 (interquartile range: [4.75, 9.00]) large plots in Jena Experiment, in which the functional group richness ranged from 1 to 4 (Supplementary Data 2). The 12 plant species represented four functional groups (grass, small herb, tall herb, and legume) (Supplementary Table 1). The detailed classifications of the functional grouping can be found in the design of the Jena Experiment40. Eleven of the 12 species were perennial, and one was annual (Trifolium dubium). The average longevity of the perennial species in the Jena Experiment has been estimated at 3–4 years61, so that multiple generations and sexual reproduction cycles could occur during the 8-year drought treatment. Although each subplot was small, population sizes of each species were estimated to range from 100 to 1000 individuals m−2 in ambient and drought subplots at the beginning of the drought treatment in the field62.We germinated the seeds in Petri dishes and transplanted the seedlings into pots in February 2017 in a glasshouse (day temperature range 20–25 °C, night temperature range 15–21 °C, and humidity range 60–80%) at the University of Zurich, Switzerland. Seedlings were planted individually, in monocultures, or in 2-species mixtures in the glasshouse (Fig. 1). In the glasshouse experiment, both monocultures and mixtures contained four plants within a pot. The pots were (11times 11times 11.5) cm in size and filled with soil composed of 50% collected from a sugar-beet field, 25% sand and 25% perlite. We randomly assigned the pots into four blocks in the glasshouse. To test the effects of drought-induced selection on plant traits, we planted individual seedlings of the 12 species in a fifth block. Within the first 2 weeks, dead individuals were replaced, thereafter dead individuals were not replaced anymore. In total, we established 958 pots: 257 pots of mixtures, 217 pots of monocultures, and 484 pots of individual plants (244 pots of individuals in blocks 1–4, and 240 pots of individuals in block 5; Supplementary Methods). For mixtures, there were 21 species pairs (Supplementary Table 1). Species pairs composed of Crepis biennis or Lotus corniculatus had low numbers of replicates (Supplementary Table 1). However, including or excluding these communities produced qualitatively similar results. Thus, we present the results including these two species in this paper. We provide detailed explanations on the choices of species pairs and regarding the biodiversity treatment history in the Jena Experiment in Supplementary Methods.During a first phase of 3 months in the glasshouse (Fig. 1), pots were watered regularly (“before drought”). After 14–16 weeks, when most of the species had reached peak aboveground biomass, we harvested all individuals in each pot by cutting them 3 cm above the ground, allowing regrowth from the left plant bases (first harvest, “before drought”). The time span for the first harvest included both the time for trait measurements (section “Plant traits” below) and for the immediately following biomass harvest. We completed the biomass harvest of each block within 1–2 days. This allowed us to account for the larger time differences between blocks by fitting block effects in the statistical analyses. After the first harvest of each block, plants were watered regularly and allowed to regrow until the 18th week from planting. This was followed by a second phase of 2 weeks without watering. Soil moisture decreased from more than 40% to less than 10% after 10 days since drought initiation. At the end of the second phase, that is after 20 weeks from planting, we made a second aboveground harvest at 3 cm above the ground (second harvest, “during drought”). During a third phase of 7 weeks, pots were watered regularly again for recovery until most plants reached a new aboveground biomass peak again. At the end of the third phase, that is after 27 weeks from planting, we harvested both above- and belowground plant biomass (third harvest, “after drought”). We checked and confirmed that most plants had reached the full-grown state and peak biomass before each harvest by monitoring their flowering. After each harvest, we cleaned and dried the harvested plant material at 70 °C for 48 h to obtain the dry biomass. We used the aboveground biomass as a proxy for productivity. Although clipping may affect plant responses to the experimental drought in the glasshouse, clipping had the advantage that all plants were “standardized” in height before the experimental drought, thus reducing carry-over effects of differential growth before the experimental drought.Additive partitioningWe used the additive partitioning approach (Eq. 1)17 to decompose the net biodiversity effect (NE) on aboveground biomass into the complementarity effect (CE) and the sampling effect (SE):$$triangle Y={Y}_{O}-{Y}_{E}=N,overline{triangle {RY}},{bar{M}}+N,{{{{{{rm{cov}}}}}}}left({{triangle }}{{{{{bf{RY}}}}}},,{{{{{bf{M}}}}}}right),$$
(1)
where (triangle Y) is the NE; ({Y}_{O}) is the observed yield (productivity) in a mixture; ({Y}_{E}) is the expected yield in the mixture, calculated from the observed yield in monocultures and their corresponding species proportions planted in the mixture, here 0.5; the two additive terms at the right side of the equation represent CE and SE, respectively; N is the number of species in the mixture, here 2. The partitioning is based on the observed and expected relative yield (RY) of species in the mixture. The expected RY of species in the mixture is the proportion planted. ∆({{{{{bf{RY}}}}}}) is the difference between observed and expected RY of species in the mixture; (overline{triangle {RY}}) is the average of ∆({{{{{bf{RY}}}}}}). A positive (overline{triangle {RY}}) indicates a positive CE; a positive covariation between monoculture yield (M), and ∆({{{{{bf{RY}}}}}}) indicates a positive SE. More details about the calculation can be found in Loreau and Hector17. We conducted the partitioning separately for each harvest, selection treatment, and block. We did not perform the partitioning for mixtures with zero biomass63. For monocultures with zero biomass in the second or third harvest, we kept the ones which had positive biomass in the previous harvest but excluded the ones which had zero biomass in the previous harvest. For example, when performing the partitioning for the second harvest, we kept the monocultures that had zero biomass in the second harvest but non-zero biomass in the first harvest; we excluded the monocultures that had zero biomass already in the first harvest. This was to assure that communities that died before the drought could not reappear during or after the drought, and communities that had died during the drought could not reappear after the drought.We used mixed-effects models to assess the influences of drought vs. ambient-selection treatments on biodiversity effects (NEs, CEs, and SEs) separately for each harvest (Fig. 2; Table 1). Block and selection treatment were set as fixed-effects terms, while species composition (identity of species pair) and its interaction with selection treatment were set as random-effects terms. This conservative approach was used to allow for generalizations across all possible species compositions, although a more liberal approach with species composition and its interactions as fixed-effects terms could also have been applied (see Schmid et al.64 for a discussion of defining terms as fixed- vs. random-effects terms, including a justification of preference for treating block as a fixed-effects term). We square-root transformed the CEs and SEs with sign reconstruction (({{{{{{rm{sign}}}}}}}(y)sqrt{y})) prior to analysis to improve the normality of residuals17. The mixed-effects model did not converge in the analysis with CE after the drought event. In this case, we used a general linear model, in which we fitted block, species composition, selection treatment, and species composition by selection treatment interaction in this order. Then we tested the significance of selection treatment using its interaction with species composition as an error term. This procedure is an alternative to mixed-effects models that estimate variance components for random-effects terms with maximum likelihood64.To test whether biodiversity effects on productivity differed from zero, we additionally tested the significance of NEs, CEs, and SEs separately for each selection treatment and harvest (Supplementary Table 3). We set block and species composition as fixed- and random-effects terms, respectively. The model corresponding to CE for ambient-selected plants during the drought event did not converge so that we fitted it with a general linear model, in which we tested the significance of the overall mean (intercept) using species composition as an error term. All statistical analyses were conducted in R 3.6.365. The mixed-effects models were conducted with asreml-R package 4.1.0.11066.Finally, we also tested whether the effects of drought selection on biodiversity effects (NEs, CEs, and SEs) in the glasshouse depended on the history of biodiversity treatment in the Jena Experiment. Most plants in the 2-species communities in the glasshouse originated from mixtures in the Jena Experiment (Supplementary Data 2; whether mixtures in the glasshouse composed of plants originating from monoculture field plots did not affect the effects of drought-selection on biodiversity effects on productivity (Supplementary Data 3)). To increase statistical power, we used functional group richness, ranging from 1 to 4, instead of species richness of the field plots as explanatory variable (Supplementary Methods). We fitted functional group richness either in linear (Supplementary Data 4) or log-linear (Supplementary Data 5) form. We did not detect significant effects of field treatment of functional group richness nor significant interactions between field treatment of functional group richness and the drought-selection history. Therefore, we excluded the history of biodiversity treatments in the field from further analyses.Biomass stability to the drought event in the glasshouseTo assess the temporal responses of community aboveground biomass to the drought event, we calculated three indices representing different facets of stability: biomass resistance, recovery, and resilience (see van Moorsel et al.43 for an example). We calculated resistance as the biomass ratio during vs. before the drought, recovery as the ratio after vs. during the drought and resilience as the ratio after vs. before the drought (see also Isbell et al.9). We log-transformed the indices (plus a half of the minimum positive value to allow taking logs of indices that were originally zero) prior to statistical analyses to improve the normality of residuals. Excluding index values that were originally zero produced qualitatively similar results.To assess the effects of drought-selection on biomass stability, we fitted mixed-effects models with block and selection treatment as fixed-effects terms, and species composition and its interaction with selection treatment as random-effects terms (Supplementary Fig. 3; Supplementary Table 4). We fitted the models separately for mixtures and monocultures. We included the log-transformed biomass at the first harvest as a covariate because biomass stability in response to droughts often depends on plant performance under ambient conditions.In the same way as net biodiversity effects on productivity were calculated for additive partitioning, we calculated biodiversity effects on biomass stability as the difference between each mixture and its corresponding monocultures. Then, we tested the influence of selection treatment on the biodiversity effects on biomass stability. Block and selection treatment were set as fixed-effects terms; species composition and its interaction with selection treatment were set as random-effects terms (Fig. 3; Supplementary Table 5). The log-transformed biomass at the first harvest was also included as a covariate43. To assess the significance of biodiversity effects on biomass stability for each selection treatment, we fitted another set of simplified models, with block and log-transformed biomass as fixed-effects terms, and species composition as random-effects term (Fig. 3).Neighbor interactionsWe assessed interactions between neighboring plants within pots using the metrics of neighbor interaction intensity with multiplicative symmetry (NIntM)44:$${NIn}{t}_{M}=2frac{triangle P}{{P}_{-N}+{P}_{+N}+left|triangle Pright|},$$
(2)
where ({P}_{-N}) and ({P}_{+N}) are the productivities without (individual plant) and with neighbors (monocultures or mixtures), respectively; (triangle P={P}_{+N}-{P}_{-N}). Negative values of NIntM indicate competition and positive values indicate facilitation. NIntM is bounded between –1 (competitive exclusion) and 1 (“obligate” facilitation). For monocultures, we first calculated the per-plant biomass as the ratio between total biomass and planting density, and then used the per-plant value to compare with the corresponding individuals (without neighbor) of the same species with the same selection treatment in the same block. Note that under the reciprocal yield law45, an individual grown alone in a pot should be four times larger than an individual grown with three others in a pot, resulting in a NIntM of –0.75. For 2-species mixtures, we calculated the per-plant biomass separately for each species and took the average NIntM of the two species to measure the interaction intensity of the mixture. We set zero biomass for dead plants in the calculation. Again, if mixtures would also follow the reciprocal yield law independent of species identity, then NIntM = –0.75 would be expected. Values greater than –0.75 indicate some sort of overyielding due to higher density or higher density and higher diversity.To assess how selection treatment modified interactions between plants, we tested the effects of selection treatment on neighbor interaction intensity separately for monocultures and mixtures. We included block and selection treatment as fixed-effects terms, species composition and its interaction with selection treatment as random-effects terms (Supplementary Fig. 4; Supplementary Table 6).We calculated the difference between the heterospecific interaction in a mixture and the conspecific interactions in its two corresponding monocultures. A positive value of this difference indicates a weaker heterospecific than conspecific competition (i.e., niche differentiation) or stronger heterospecific than conspecific facilitation, which may lead to a positive complementarity effect. We tested the effects of selection treatment on interaction difference for each harvest by fitting block and selection treatment as fixed-effects terms, and species composition and its interaction with selection treatment as random-effects terms (Fig. 4; Supplementary Table 8). We also tested the significance of the interaction difference for each selection treatment by fitting block and species composition as fixed- and random-effects term, respectively (Fig. 4; Supplementary Table 7).Plant traitsTo assess whether drought selection would change plant traits, we measured six traits (Supplementary Table 9) closely related to plant usages of water or carbon on plants in pots with one individual from blocks 1–5. We focused on the traits on individual plants without neighbor to evaluate the influence of selection treatment on traits without the impacts of plasticity induced by plant interactions. We measured leaf relative chlorophyll content, leaf area (LA), leaf mass per area (LMA) and leaf osmometric pressure before the drought; leaf stomatal conductance both before and during the drought; and dry biomass ratio between root and shoot after the drought (in the third harvest). Leaf relative chlorophyll content was measured for three mature, fully expanded leaves per plant by using a SPAD-502 Plus chlorophyll meter from Konica Minolta. LA was obtained by scanning 3–4 mature, fully expanded leaves per plant with a LI-3100C Area Meter from LI-COR. LMA was calculated as the ratio between leaf dry mass (oven-dried at 70 °C for 48 h, using the same leaves that for LA) and LA. Leaf osmotic potential at full hydration was considered as an important trait associated with plant tolerance to drought30. We measured leaf osmotic potential with freeze-thaw leaf pieces cut from 1 to 2 mature, fully expanded leaves per plant by using a Wescor vapor pressure osmometer VAPRO (Model 5520) according to the method by Bartlett, et al.30. Plants were fully hydrated 1 day before the leaf sampling for osmotic potential measurement. Leaf stomatal conductance is a measure of exchange rate of carbon dioxide and water vapor through the stomata67. It was measured for 3–5 healthy mature leaves per plant by using a SC-1 Leaf Porometer from Decagon Devices. For grass species, 3 blades were placed adjacent to each other to have a large enough area for the measurement of stomatal conductance. For stomatal conductance during the drought event, we measured the individual plants from block 5 only due to limited time during the drought phase. We harvested aboveground and belowground plant biomass separately for alive individuals at the end of the experiment (after the complete recovery from the drought). The oven-dried (70 °C for 48 h) aboveground and belowground biomass were used to calculate the biomass ratio between root and shoot. We took the average value of each trait of each plant for statistical analyses. Each trait was measured for each block in turn.We used linear mixed-effects models to assess the influence (generalized across species) of selection treatment on trait values (red lines in Supplementary Figs. 5–7). Block and selection treatment were set as fixed-effects terms; species and its interaction with selection treatment were set as random-effects terms. Alternatively, we set species, selection treatment and their interaction as fixed-effects terms to assess whether species responded differently to the selection treatment (Supplementary Table 9). To test whether effects of selection treatment on traits differed between the five trait groups (leaf relative chlorophyll content, leaf area, leaf mass per area, leaf osmometric pressure, and leaf stomatal conductance) measured before the drought event in the glasshouse, we conducted two alternative analyses. First, we performed a principal component analysis with all traits and retained the first two principal axes (PC1 and PC2), which accounted for 39.06% and 22.3% of the total variation, respectively. Then we used PC1 and PC2 as response variables in mixed-effect models, separately. We fitted the models with the same fixed- and random-effects terms as those using each separate trait as the response variable. Effects of selection treatment on PC1 or PC2 were not significant. Second, we pooled the five traits as a single response variable in a mixed-effect model (corresponding to multivariate analysis of variance). Block, trait group (a factor with five levels), selection treatment, and the interaction between trait group and selection treatment were set as fixed-effects terms; species and its interactions with trait group and selection treatment and their three-way interaction were set as random-effects terms. We did not detect significant effects of selection treatment nor its interaction with trait group. Therefore, we did not present the results associated with these multivariate analyses in this paper. LMA, LA, leaf osmotic potential, leaf stomatal conductance, and root-shoot biomass ratio were log-transformed to improve normality of residuals.We also measured leaf relative chlorophyll content, LA and LMA in mixtures before the drought event (Supplementary Table 10) to evaluate the influence of selection treatment on trait dissimilarity between interacting species within communities. We calculated the absolute trait distance between two species in each mixture both separately for each trait and jointly with the three traits. For multi-trait-based dissimilarity, we standardized each trait to mean zero and unit standard deviation and calculated the Euclidean trait distance in standardized three-dimensional trait space.We used linear mixed-effects models to assess the effects of selection treatment on trait dissimilarity in mixtures (Supplementary Table 10). Block and selection treatment were set as fixed-effects terms; species composition and its interaction with selection treatment were set as random-effects terms. The model for LA dissimilarity did not converge so that we fit it with a general linear model, in which we tested the significance of selection treatment using its interaction with species composition as an error term. For the models with LA, LMA, and the joint three traits as dependent variables, we removed one pot (B1P674) because the LA value of Alopecurus pratensis in this pot was extremely small (about 1/3 of the second minimum value of the same species in mixtures). However, including or excluding this pot produced qualitatively similar results.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Trait modelsOur analysis included 491,001 unique trait measurements across 18 traits, encompassing 13,189 tree species from 2313 genera, reflecting ~21% of all known tree species33 (Fig. 1). Traits were measured at 8683 locations across the globe and 373 distinct eco-regions (Supplementary Tables 1, 2), with georeferenced measurements capturing 15% of known tree species in Eurasia, 13% in South America, 9% in Oceania, and 6% in North America and Africa33. The raw data covered 22% of all trait-by-species combinations (Fig. 1b, Supplementary Fig. 2), nearly identical to other large-scale trait analyses across the entire plant kingdom5,17,30. Yet there was considerable variation in coverage across traits, with traits such as specific leaf area and leaf nitrogen measured on more than 60% of all species, versus traits such as crown diameter and conduit diameter, which captured fewer than 5% of species (Fig. 1b, Supplementary Fig. 2). Across all species, 423 had more than 10 unique traits measured, and two species (Picea abies and Pinus sylvestris) had measurements for all 18 traits. In general, there was highly consistent coverage across taxonomic orders and traits (Supplementary Fig. 1), with gymnosperms being slightly overrepresented (comprising 3.1 ± 6.8% of measurements in the database versus ~1% of all known tree species34,35, Fig. 1a), in part reflecting the wider geographic range of many gymnosperms relative to angiosperms36.To explore relationships in functional traits at the individual level, we used random-forest machine-learning models to estimate missing trait values for each individual tree as a function of its environment and phylogenetic history. We also conducted a second set of analyses where trait expression was estimated using phylogenetic information only, which allowed us to include additional non-georeferenced data (Fig. 1), while also quantifying the relative contribution of environmental information on trait expression (Supplementary Fig. 6). Following standard approaches5,15,29,30, all traits were log-transformed and standardized to allow for statistically robust comparisons. Environmental predictors included ten variables encompassing climate37,38,39,40, soil41, topographic42, and geological43 features. Phylogenetic history was incorporated via the first ten phylogenetic eigenvectors44,45 (see Methods). By including environmental information alongside phylogenetic information, this approach not only allowed us to impute species-level traits which have strong phylogenetic signals and weak environmental signals, as is traditionally done17,30 but also to robustly estimate traits which have a weak phylogenetic signal and are instead strongly sensitive to environmental conditions. Moreover, being a non-parametric approach, the random forest makes no a priori assumptions about how trait expression varies across phylogenetic groups or environments.Across all 18 traits, the best-fitting models explained 54 ± 14% of out-of-fit trait variation (VEcv, see Methods), ranging from 26% for stem diameter to 76% of the variation in leaf area (Supplementary Figs. 6, 7). This accuracy was quantified using buffered leave-one-out cross-validation to account for spatial and phylogenetic autocorrelation46, and thus serves as a conservative lower bound for species which are phylogenetically and environmentally distinct from the observations47. There was no significant relationship between out-of-fit cross-validation accuracy and sample size (R2 = 0.06, p = 0.33), highlighting the relatively broad taxonomic coverage for each trait (Fig. 1, Supplementary Fig. 1).Environmental variables and phylogenetic information had approximately equal explanatory power (relative importance of 0.51 vs 0.49 for environment vs. phylogeny), albeit with substantial variation across traits (Supplementary Fig. 9). The inclusion of environmental variables increased the explanatory power of the models by 35%, on average (Supplementary Fig. 6), with crown diameter, crown height, leaf density, and stem diameter exhibiting the largest relative increases (54%, 45%, 73%, and 26%, respectively), mirroring the fact that these traits have comparatively low phylogenetic signal relative to other traits (assessed via Pagel’s λ on the raw data, Fig. 4c). Seed dry mass was the only trait with a substantial increase in accuracy using the phylogeny-only model (25% improvement; Supplementary Fig. 6), reflecting the fact that seed dry mass had the strongest phylogenetic signal of all traits (Fig. 4c), and also because this trait has a substantial amount of additional non-georeferenced data that was included in the phylogeny-only models (Fig. 1b). Wood density was the only trait with nearly identical predictive power whether or not environmental information was included, whereas all other traits exhibited significantly reduced accuracy when environmental information was excluded (Supplementary Fig. 6).Relationships in tree trait expressionUsing the resulting trait models, we imputed missing trait values for every tree with at least one georeferenced trait measurement. For all traits except seed dry mass, we used the random-forest models accounting for environmental and phylogenetic information; for seed dry mass, we used the phylogeny-only model to estimate expression due to its substantially higher data availability and out-of-fit accuracy. For tree height, stem diameter, crown height, crown width, and root depth, we used quantile random forest48 to estimate the upper 90th percentile value for each species in its given location, thereby minimizing ontogenetic variation across a tree’s lifetime (see Methods). We used the resulting trait data to explore the dominant drivers of trait variation using species-weighted principal component analysis, accounting for an unequal number of observations across species.When considering all traits simultaneously, the first two axes of the resulting principal components (PC) capture 41% of the variation in overall trait expression (Fig. 2a; Supplementary Fig. 10; Supplementary Table 5). The first trait axis correlates most strongly with leaf thickness, specific leaf area, and leaf nitrogen (PC loadings of L = 0.77, 0.74, and 0.73, respectively). By capturing key aspects of the leaf-economic spectrum14, these traits reflect various physiological controls on leaf-level resource processing, tissue turnover and photosynthetic rates49. Thick leaves with low specific leaf area (SLA) can help minimize desiccation, frost damage, and nutrient limitation, but at the cost of reduced photosynthetic potential due to primary investment in structural resistance50. Accordingly, leaf nitrogen—a crucial component of Rubisco for photosynthesis51—trades off strongly with leaf thickness. This first axis thus captures the core distinction between “acquisitive” (fast) and “conservative” (slow) life-history strategies across the plant kingdom7,52, reflecting an organismal-level trade-off between the high photosynthetic potential in optimal conditions versus abiotic tolerance in suboptimal conditions. Nevertheless, leaf density—which is related to SLA and is a key feature of the leaf-economic spectrum—loads relatively weakly on this first trait axis compared to other leaf traits (L = −0.28 for axis 1, vs 0.20 for axis 2; Supplementary Table 5), highlighting important aspects of leaf structure that are not captured by this dominant trait axis53.Fig. 2: The dominant trait axes and relationships.Shown are the first two principal component axes capturing trait relationships across the 18 functional traits. a All tree species (n = 30,146 observations), b angiosperms only (n = 24,658), and c gymnosperms only (n = 5498). In a the three variables that load most strongly on each axis are shown in dark black lines, with the remaining variables shown in light grey. These same six variables are highlighted in b and c illustrating how the same relationships extend to angiosperms and gymnosperms (see Supplementary Figs. 10–12 for the full PCAs with all traits visible, and Supplementary Table 5 for the PC loadings).Full size imageThe second trait axis correlates most strongly with maximum tree height (PC loading of L = 0.77), crown height, (L = 0.75), and crown diameter (L = 0.88), highlighting the overarching importance of competition for light and canopy position in forests7 (Fig. 2a; Supplementary Fig. 10; Supplementary Table 5). Large trees and large crowns are critical for light access and for maximizing light interception down through the canopy54. Nevertheless, tall trees with deep crowns also experience greater susceptibility to disturbance and mechanical damage, primarily due to wind and weight25. Because of the massive carbon and nutrient costs required to create large woody structures55,56, larger trees are less viable in nutrient-limited or colder climates57, and in exposed areas with high winds or extreme weather events58. This second axis thus reflects a fundamental biotic/abiotic trade-off related to overall tree size, which is largely orthogonal to leaf-level nutrient-use and photosynthetic capacity.Despite substantial differences in wood and leaf structures between angiosperms and gymnosperms (e.g. vessels vs. tracheids), the two main relationships hold within, as well as across, angiosperms and gymnosperms (Fig. 2b, c; Supplementary Figs. 11, 12). Indeed, angiosperms and gymnosperms are subject to the same physical, mechanical, and chemical processes that determine the ability to withstand various biotic and abiotic pressures59.Collectively, these two primary trait axes capture two dominant ecological trade-offs that underpin tree survival in any given environment: (1) the ability to maximize leaf photosynthetic activity, at the cost of increased risk of leaf desiccation, and (2) the ability to compete for space and maximize light interception, at the cost of increased susceptibility to mechanical damage. By capturing two aspects of conservative-acquisitive life-history strategies, these two relationships closely mirror those seen when considering herbaceous species alongside woody species5,17. However, in line with our expectations, these two axes capture only ~40% of the variation in trait space, versus nearly ~75% of variation when considering only six traits across the entire plant kingdom5. Here, the first seven PC axes are needed to account for 75% of the variation across all 18 traits (Supplementary Table 5). Thus, while this analysis supports the universality of these two primary PC axes, it also demonstrates that the majority of trait variation in trees is unexplained by these two dimensions. As such, quantifying the full dimensionality of trait space by exploring multidimensional trait clusters is needed to better capture the wide breadth of tree form and function.Environmental predictors of trait relationshipsTo examine how environmental variation shapes trait expression across the globe, we next quantified the relationships between environmental conditions and the dominant trait axes. Using Shapley values60, we partitioned the relative influence of each environmental variable on the PC trait axes, controlling for all other variables in the model (see Methods).In line with previous analysis across the plant kingdom61, temperature variables were the strongest drivers of trait relationships (Fig. 3, Supplementary Figs. 17, 18), with annual temperature having the strongest influence both on leaf-economic traits (PC axis 1, Fig. 3c) and on tree-size traits (PC axis 2, Fig. 3d). Leaves face increased frost risk and reduced photosynthetic potential in colder conditions, such that ecological selection should favour thick leaves with low SLA over thin leaves with high SLA and high nutrient-use49. Trees in warm environments are more likely to experience strong biotic interactions, which should increase evolutionary and ecological selection pressures over time62,63, favouring tall species with large crowns that have high competitive ability and efficient light acquisition strategies. Annual temperature thus predominantly reflects the transition from gymnosperm- to angiosperm-dominated ecosystems, with this inflection point occurring at ~15 °C for both axes, demonstrating strong environmental convergence between the dominant axes of trait variation.Fig. 3: The relationship between environmental variables and trait axes.a, b The relative influence of the environmental variables on the two dominant PC axes. The ten variables are sorted by overall variable importance in the models (see Methods). Yellow points are observations which have high values of that environmental variable; blue values are the lowest. Points to the right of zero indicate a positive influence on the PC axis; points to the left indicate a negative influence (see also Supplementary Figs. 17, 18). c–h The relationships between environmental variables and PC axis values for the three variables in a with the strongest influence. Values above zero show a positive influence on PC axis values; values less than zero indicate a negative influence.Full size imageBeyond annual temperature, each trait axis demonstrated different relationships with climate, soil, and topographic variables (Fig. 3a, b, Supplementary Figs. 17, 18). Percent sand content had the second-highest influence on the first trait axis (Fig. 3e), supporting patterns seen across the entire plant kingdom17. Sand content is a strong proxy for soil moisture and soil-available nutrients such as phosphorous, and is therefore closely tied to leaf photosynthetic rates64. In contrast to previous work, however, we find that soil characteristics have correspondingly little effect on the second axis of trait variation (Fig. 3b; Supplementary Fig. 18). Instead, precipitation was the second strongest driver of tree height and crown size (Fig. 3f), with large trees with large crowns becoming consistently more frequent with increasing precipitation. These results highlight that, despite the primary importance of temperature, the main climate stressors to trees (e.g. xylem cavitation and embolism, fire regimes, and leaf desiccation) typically arise via interactions between temperature, soil nutrients, and water availability.For both axes, elevation was the third strongest driver of trait values (Fig. 3g, h), highlighting a critical component of tree functional biogeography that extends beyond climate and soil. Yet the effects of elevation on trait expression differed somewhat across the two axes. For the first axis related to leaf-economic traits, there is little influence at low elevations, followed by a sharp transition at ~2000 m towards gymnosperm-dominated species with thick leaves, low SLA, and low leaf N. For the second trait axis related to tree size, elevation instead has a strong positive influence on tree height and crown size at low elevations, which becomes increasingly less influential past ~500 m. Such results partly reflect the transition from angiosperm to gymnosperm-dominated stands at higher elevations (blue vs. red points, Fig. 3g, h), and potentially the role of environmentally mediated intraspecific variation in traits such as tree height65,66.These results demonstrate close alignment of the dominant trait PC axes across biogeographic regions. Despite the orthogonality of these axes in trait species, environmental conditions place similar constraints on both trait axes, particularly at the environmental extremes (e.g. warm, moist, low elevation vs. cold, dry, high elevation), leading to convergence of the dominant trait axes across environmental gradients.Trait clusters at the global scaleTo better explore the multidimensional nature of trait relationships that are not fully covered by the dominant two axes, we subsequently identified groups of traits that form tightly coupled clusters and which reflect distinct aspects of tree form and function.Our results show that these 18 traits can be grouped into eight trait clusters, each of which reflects a unique aspect of morphology, physiology, or ecology (Fig. 4a, Supplementary Fig. 23). The largest trait cluster (Fig. 4a, pink cluster) demonstrates wood/leaf integration of moisture regulation and photosynthetic activity via the inclusion of leaf area, stem conduit diameter, stomatal conductance, and leaf Vcmax (the maximum rate of carboxylation). Distinct from this cluster are the three traits loading most strongly on PC axis 1 (SLA, leaf thickness and leaf N; Fig. 4a, yellow), highlighting complementary aspects of the leaf-economic spectrum indicative of acquisitive vs. conservative resource use15. The role of leaf K and P in leaf nutrient economies are well established7,67, and yet these traits form a distinct cluster from the other leaf-economic traits (Fig. 4a, light blue) due to their relatively high correlation with tree height and crown size, particularly for leaf K, which loads almost equally on both trait axes (Fig. 4b, Supplementary Table 5).Fig. 4: Trait correlations and functional clusters.a Trait clusters with high average intra-group correlation. The upper triangle gives the species-weighted correlations incorporating intraspecific variation. The lower triangle gives the corresponding correlations among phylogenetic independent contrasts, which adjusts for pseudo-replication due to the non-independence of closely related species. The size of the circle denotes the relative strength of the correlation, with solid circles denoting positive correlations and open circles denoting negative correlations (see Supplementary Fig. 19 for the numeric values). b PC loadings for each trait and each of the first two principal component axes, illustrating which functional trait clusters align most strongly with the dominant axes of trait variation (see Supplementary Table 5 for the full set of PC loadings). c The species-level phylogenetic signal of each trait (Pagel’s λ), calculated using only the raw trait values.Full size imageTree height and crown size form their own distinct cluster (Fig. 4a, dark green), further supporting the inference that these traits reflect key aspects of tree form and function independent of the leaf-economic spectrum. Yet leaf area, despite being part of the cluster reflecting moisture regulation and photosynthetic activity, loads almost equally on PC axes 1 and 2 (Fig. 4b, Supplementary Table 5), highlighting that it serves as an intermediary between the two key aspects of tree size and leaf economics. It is a critical driver of moisture regulation and photosynthetic capacity, while also playing an important role in the light acquisition, leaf-turnover time, and competitive ability54,68.There are two additional two-trait clusters, both of which load relatively poorly on the two primary PC axes: (1) stem diameter and bark thickness (Fig. 4, dark blue), and (2) wood and leaf density (Fig. 4, light green). Bark thickness increases with tree size not only as a result of bark accumulation as trees age, but also due to the functional/metabolic needs of the plant69,70. From an ecological perspective, thick bark can be critical for defense against fire and pest damage (mainly a thick outer bark region), for storage and photosynthate transportation needs (mainly a thick inner bark region)71,72. Yet such relationships are strongly ecosystem-dependent, with tree size emerging as the dominant driver at the global scale70. In contrast, wood density and leaf density are strongly linked to slow/fast life-history strategies, where denser plant parts reduce growth rate and water transport6,15 but protect against pest damage, desiccation, and mechanical breakage6,50,56. As such, leaf density captures fundamentally unique aspects of leaf form and function relative to other leaf traits such as SLA53 (Fig. 4b, Supplementary Table 5), and our results support the inference that these translate into fundamentally different ecological strategies73. Collectively, these two-trait clusters each demonstrate unique and complementary mechanisms that insulate trees against various disturbances and extreme weather events, but at the cost of reduced growth, competitive ability, and productivity under optimal conditions (see Supplementary Notes).Lastly, two traits each comprise their own unique cluster: root depth and seed dry mass (Fig. 4a, purple and orange, respectively). Root growth is subject to a range of belowground processes (e.g. root herbivory, depth to bedrock), and our results confirm previous work demonstrating a clear disconnect between aboveground and belowground traits23,74,75. Root depth accordingly has a relatively weak phylogenetic signal (λ = 0.44, Fig. 4c) but a strong environmental signal (Supplementary Figs. 6, 9), reflecting distinct belowground constraints on trait expression23. In contrast, seed dry mass exhibits the strongest phylogenetic signal (λ = 0.98, Fig. 4c) and weakest environmental signal of any trait (Supplementary Figs. 6, 9), and it accordingly was the only trait where the phylogeny-only model performed substantially better (Supplementary Fig. 6). In line with previous work, seed dry mass has moderate correlations with various other traits underpinning leaf economics and tree size5,28 (e.g. ρ = 0.28, −0.22, and 0.22 for tree height, leaf K, and leaf density, using the raw data), yet it exhibits relatively weak correlation with most other traits, placing it in a distinct functional cluster. Reproductive traits are subject to unique evolutionary pressures26, indicative of different seed dispersal vectors (wind, water, animals) and various ecological stressors that uniquely affect seed viability and germination26. The emergence of root depth and seed dry mass as solo functional clusters thus supports the previous inference that belowground traits74 and reproductive traits26 reflect distinct aspects of tree form and function not fully captured by leaf or wood trait spectrums. More

Benhamou, S. Of scales and stationarity in animal movements. Ecol. Lett. 17, 261–272 (2014).PubMed 
Article 

Google Scholar 
Owen-Smith, N. Effects of temporal variability in resources on foraging behaviour. In Resource Ecology (eds. Prins, H. H. T. & Van Langevelde, F.) 159–181 (Springer Netherlands, 2008).Hutchinson, G. E. The multivariate niche. Cold Spring Harb. Symp. Quant. Biol. 22, 415–421 (1957).Article 

Google Scholar 
Kearney, M. Habitat, environment and niche: what are we modelling? Oikos 115, 186–191 (2006).Article 

Google Scholar 
Tauber, E., Last, K. S., Olive, P. J. W. & Kyriacou, C. P. Clock gene evolution and functional divergence. J. Biol. Rhythms 19, 445–458 (2004).CAS 
PubMed 
Article 

Google Scholar 
Pilorz, V., Helfrich-Förster, C. & Oster, H. The role of the circadian clock system in physiology. Pflug. Arch. – Eur. J. Physiol. 470, 227–239 (2018).CAS 
Article 

Google Scholar 
Levy, O., Dayan, T., Porter, W. P. & Kronfeld-Schor, N. Time and ecological resilience: can diurnal animals compensate for climate change by shifting to nocturnal activity? Ecol. Monogr. 89, e01334 (2019).Article 

Google Scholar 
Gaynor, K. M., Hojnowski, C. E., Carter, N. H. & Brashares, J. S. The influence of human disturbance on wildlife nocturnality. Science 360, 1232–1235 (2018).CAS 
PubMed 
Article 

Google Scholar 
Cox, D. T. C., Gardner, A. S. & Gaston, K. J. Diel niche variation in mammals associated with expanded trait space. Nat. Commun. 12, 1753 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Kronfeld-Schor, N. & Dayan, T. Partitioning of time as an ecological resource. Annu. Rev. Ecol. Evol. Syst. 34, 153–181 (2003).Article 

Google Scholar 
Kronfeld‐Schor, N. et al. On the use of the time axis for ecological separation: diel rhythms as an evolutionary constraint. Am. Nat. 158, 451–457 (2001).PubMed 
Article 

Google Scholar 
Austin, R. W. & Petzold, T. J. Spectral dependence of the diffuse attenuation coefficient of light in ocean waters. OE OPEGAR 25, 253471 (1986).Article 

Google Scholar 
Bandara, K., Varpe, Ø., Wijewardene, L., Tverberg, V. & Eiane, K. Two hundred years of zooplankton vertical migration research. Biol. Rev. 96, 1547–1589 (2021).PubMed 
Article 

Google Scholar 
Brierley, A. S. Diel vertical migration. Curr. Biol. 24, R1074–R1076 (2014).CAS 
PubMed 
Article 

Google Scholar 
Hays, G. C. A review of the adaptive significance and ecosystem consequences of zooplankton diel vertical migrations. In Migrations and Dispersal of Marine Organisms 163–170 (Springer, 2003).Aumont, O., Maury, O., Lefort, S. & Bopp, L. Evaluating the potential impacts of the diurnal vertical migration by marine organisms on marine biogeochemistry. Glob. Biogeochem. Cycles https://doi.org/10.1029/2018GB005886 (2018).Article 

Google Scholar 
Tarrant, A. M., McNamara-Bordewick, N., Blanco-Bercial, L., Miccoli, A. & Maas, A. E. Diel metabolic patterns in a migratory oceanic copepod. J. Exp. Mar. Biol. Ecol. 545, 151643 (2021).Article 

Google Scholar 
Cohen, J. H. & Forward, Jr. R. B. Zooplankton diel vertical migration—a review of proximate control. In Oceanography and Marine Biology (eds Gibson, R. N., Atkinson, R. J. A. & Gordon, J. D. M.) 89–122 (CRC Press, 2009).Benoit-Bird, K. J., Au, W. W. L. & Wisdoma, D. W. Nocturnal light and lunar cycle effects on diel migration of micronekton. Limnol. Oceanogr. 54, 1789–1800 (2009).Article 

Google Scholar 
Last, K. S., Hobbs, L., Berge, J., Brierley, A. S. & Cottier, F. Moonlight drives ocean-scale mass vertical migration of zooplankton during the Arctic Winter. Curr. Biol. 26, 244–251 (2016).CAS 
PubMed 
Article 

Google Scholar 
Omand, M. M., Steinberg, D. K. & Stamieszkin, K. Cloud shadows drive vertical migrations of deep-dwelling marine life. PNAS 118, e2022977118 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Strömberg, J.-O., Spicer, J. I., Liljebladh, B. & Thomasson, M. A. Northern krill, Meganyctiphanes norvegica, come up to see the last eclipse of the millennium? J. Mar. Biol. Assoc. UK 82, 919–920 (2002).Article 

Google Scholar 
Ludvigsen, M. et al. Use of an Autonomous Surface Vehicle reveals small-scale diel vertical migrations of zooplankton and susceptibility to light pollution under low solar irradiance. Sci. Adv. 4, eaap9887 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Häfker, N. S. et al. Circadian clock involvement in zooplankton diel vertical migration. Curr. Biol. 27, 2194–2201 (2017).PubMed 
Article 
CAS 

Google Scholar 
Chen, C. et al. Drosophila Ionotropic Receptor 25a mediates circadian clock resetting by temperature. Nature 527, 516–520 (2015).CAS 
PubMed 
Article 

Google Scholar 
Epifanio, C. E. & Cohen, J. H. Behavioral adaptations in larvae of brachyuran crabs: a review. J. Exp. Mar. Biol. Ecol. 482, 85–105 (2016).Article 

Google Scholar 
Sorek, M. et al. Setting the pace: host rhythmic behaviour and gene expression patterns in the facultatively symbiotic cnidarian Aiptasia are determined largely by Symbiodinium. Microbiome 6, 83 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Hobbs, L., Banas, N. S., Cottier, F. R., Berge, J. & Daase, M. Eat or sleep: availability of winter prey explains mid-winter and spring activity in an Arctic Calanus population. Front. Mar. Sci. 7, 541564 (2020).Article 

Google Scholar 
Urmy, S. S., Horne, J. K. & Barbee, D. H. Measuring the vertical distributional variability of pelagic fauna in Monterey Bay. ICES J. Mar. Sci. 69, 184–196 (2012).Article 

Google Scholar 
Berge, J. et al. Arctic complexity: a case study on diel vertical migration of zooplankton. J. Plankton Res. 36, 1279–1297 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Berge, J. et al. In the dark: A review of ecosystem processes during the Arctic polar night. Prog. Oceanogr. https://doi.org/10.1016/j.pocean.2015.08.005 (2015).Article 

Google Scholar 
Pavlov, A. K. et al. The underwater light climate in Kongsfjorden and Its ecological implications. In The Ecosystem of Kongsfjorden, Svalbard (eds Hop, H. & Wiencke, C.) 137–170 (Springer International Publishing, 2019).Cohen, J. H. et al. Is ambient light during the high arctic polar night sufficient to act as a visual cue for Zooplankton? PLoS ONE 10, e0126247 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Veedin Rajan, V. B. et al. Seasonal variation in UVA light drives hormonal and behavioural changes in a marine annelid via a ciliary opsin. Nat. Ecol. Evol. 5, 204–218 (2021).PubMed 
Article 

Google Scholar 
Vinayak, P. et al. Exquisite light sensitivity of Drosophila melanogaster cryptochrome. PLoS Genet. 9, e1003615 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Verasztó, C. et al. Ciliary and rhabdomeric photoreceptor-cell circuits form a spectral depth gauge in marine zooplankton. eLife 7, e36440 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Hobbs, L. et al. A marine zooplankton community vertically structured by light across diel to interannual timescales. Biol. Lett. 17, 20200810 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Daase, M., Eiane, K., Aksnes, D. L. & Vogedes, D. Vertical distribution of Calanus spp. and Metridia longa at four Arctic locations. Mar. Biol. Res. 4, 193–207 (2008).Article 

Google Scholar 
Irigoien, X., Conway, D. V. P. & Harris, R. P. Flexible diel vertical migration behaviour of zooplankton in the Irish Sea. Mar. Ecol. Prog. Ser. 267, 85–97 (2004).Article 

Google Scholar 
Frost, B. W. & Bollens, S. M. Variability of diel vertical migration in the marine planktonic copepod Pseudocalanus newmani in relation to its predators. Can. J. Fish. Aquat. Sci. 49, 1137–1141 (1992).Article 

Google Scholar 
Tarling, G. A., Jarvis, T., Emsley, S. M. & Matthews, J. B. L. Midnight sinking behaviour in Calanus finmarchicus: a response to satiation or krill predation? Mar. Ecol. Prog. Ser. 240, 183–194 (2002).Article 

Google Scholar 
Hays, G. C., Proctor, C. A., John, A. W. G. & Warner, A. J. Interspecific differences in the diel vertical migration of marine copepods: the implications of size, color, and morphology. Limnol. Oceanogr. 39, 1621–1629 (1994).Article 

Google Scholar 
Gastauer, S., Nickels, C. F. & Ohman, M. D. Body size- and season-dependent diel vertical migration of mesozooplankton resolved acoustically in the San Diego Trough. Limnol. Oceanogr. 67, 300–313 (2021).Article 

Google Scholar 
Hardy, A. C. & Bainbridge, R. Experimental observations on the vertical migrations of plankton animals. J. Mar. Biol. Assoc. UK 33, 409–448 (1954).Article 

Google Scholar 
Musilova, Z. et al. Vision using multiple distinct rod opsins in deep-sea fishes. Science 364, 588–592 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Gornik, S. G. et al. Photoreceptor diversification accompanies the evolution of Anthozoa. Mol. Biol. Evol. 38, 1744–1760 (2021).CAS 
PubMed 
Article 

Google Scholar 
Cohen, J. H. et al. Photophysiological cycles in Arctic krill are entrained by weak midday twilight during the Polar Night. PLoS Biol. 19, e3001413 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Kopperud, K. L. & Grace, M. S. Circadian rhythms of retinal sensitivity in the Atlantic tarpon, Megalops atlanticus. Bull. Mar. Sci. https://doi.org/10.5343/bms.2016.1045 (2017).Article 

Google Scholar 
Ohguro, C., Moriyama, Y. & Tomioka, K. The compound eye possesses a self-sustaining Circadian Oscillator in the Cricket Gryllus bimaculatus. Zool. Sci. 38, 82–89 (2020).Article 

Google Scholar 
Brodrick, E. A., How, M. J. & Hemmi, J. M. Fiddler crab electroretinograms reveal vast circadian shifts in visual sensitivity and temporal summation in dim light. J. Exp. Biol. jeb.243693, https://doi.org/10.1242/jeb.243693 (2022).Kaartvedt, S., Røstad, A., Christiansen, S. & Klevjer, T. A. Diel vertical migration and individual behavior of nekton beyond the ocean’s twilight zone. Deep Sea Res. Part I: Oceanogr. Res. Pap. 103280, https://doi.org/10.1016/j.dsr.2020.103280 (2020).Flôres, D. E. F. L., Jannetti, M. G., Valentinuzzi, V. S. & Oda, G. A. Entrainment of circadian rhythms to irregular light/dark cycles: a subterranean perspective. Sci. Rep. 6, 34264 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Hays, G. C., Kennedy, H. & Frost, B. W. Individual variability in diel vertical migration of a marine copepod: why some individuals remain at depth when others migrate. Limnol. Oceanogr. 46, 2050–2054 (2001).Article 

Google Scholar 
Cohen, J. H. & Forward, R. B. Jr. Photobehavior as an inducible defense in the marine copepod Calanopia americana. Limnol. Oceanogr. 50, 1269–1277 (2005).Article 

Google Scholar 
Kvile, K. Ø., Altin, D., Thommesen, L. & Titelman, J. Predation risk alters life history strategies in an oceanic copepod. Ecology 102, e03214 (2021).PubMed 
Article 

Google Scholar 
Spaak, P. & Ringelberg, J. Differential behaviour and shifts in genotype composition during the beginning of a seasonal period of diel vertical migration. Hydrobiologia 360, 177–185 (1997).Article 

Google Scholar 
Buskey, E. J. & Swift, E. Behavioral responses of oceanic zooplankton to simulated bioluminescence. Biol. Bull. 168, 263–275 (1985).Article 

Google Scholar 
Berndt, A. et al. A novel photoreaction mechanism for the circadian blue light photoreceptor Drosophila cryptochrome. J. Biol. Chem. 282, 13011–13021 (2007).CAS 
PubMed 
Article 

Google Scholar 
Franz-Badur, S. et al. Structural changes within the bifunctional cryptochrome/photolyase CraCRY upon blue light excitation. Sci. Rep. 9, 9896 (2019).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Biscontin, A. et al. Functional characterization of the circadian clock in the Antarctic krill, Euphausia superba. Sci. Rep. 7, 17742 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Piccolin, F. et al. Photoperiodic modulation of circadian functions in Antarctic krill Euphausia superba Dana, 1850 (Euphausiacea). J. Crustacean Biol. 38, 707–715 (2018).
Google Scholar 
Piccolin, F., Pitzschler, L., Biscontin, A., Kawaguchi, S. & Meyer, B. Circadian regulation of diel vertical migration (DVM) and metabolism in Antarctic krill Euphausia superba. Sci. Rep. 10, 16796 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Häfker, N. S., Teschke, M., Hüppe, L. & Meyer, B. Calanus finmarchicus diel and seasonal rhythmicity in relation to endogenous timing under extreme polar photoperiods. Mar. Ecol. Prog. Ser. 603, 79–92 (2018).Article 
CAS 

Google Scholar 
Häfker, N. S. et al. Calanus finmarchicus seasonal cycle and diapause in relation to gene expression, physiology, and endogenous clocks. Limnol. Oceanogr. 63, 2815–2838 (2018).Article 

Google Scholar 
Hüppe, L. et al. Evidence for oscillating circadian clock genes in the copepod Calanus finmarchicus during the summer solstice in the high Arctic. Biol. Lett. 16, 20200257 (2020).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Dmitrenko, I. A. et al. Sea-ice and water dynamics and moonlight impact the acoustic backscatter diurnal signal over the eastern Beaufort Sea continental slope. Ocean Sci. 16, 1261–1283 (2020).CAS 
Article 

Google Scholar 
Hobbs, L., Cottier, F. R., Last, K. S. & Berge, J. Pan-Arctic diel vertical migration during the polar night. Mar. Ecol. Prog. Ser. 605, 61–72 (2018).Article 

Google Scholar 
Chittka, L., Stelzer, R. J. & Stanewsky, R. Daily changes in ultraviolet light levels can synchronize the circadian clock of Bumblebees (Bombus terrestris). Chronobiol. Int. 30, 434–442 (2013).PubMed 
Article 

Google Scholar 
Pauers, M. J., Kuchenbecker, J. A., Neitz, M. & Neitz, J. Changes in the colour of light cue circadian activity. Anim. Behav. 83, 1143–1151 (2012).PubMed 
PubMed Central 
Article 

Google Scholar 
Mouland, J. W., Martial, F., Watson, A., Lucas, R. J. & Brown, T. M. Cones support alignment to an inconsistent world by suppressing mouse circadian responses to the blue colors associated with twilight. Curr. Biol. 29, 4260–4267.e4 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Walmsley, L. et al. Colour as a signal for entraining the mammalian circadian clock. PLoS Biol. 13, e1002127 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Ashley, N. T., Schwabl, I., Goymann, W. & Buck, C. L. Keeping time under the midnight sun: behavioral and plasma melatonin profiles of free-living lapland longspurs (Calcarius lapponicus) during the Arctic Summer. J. Exp. Zool. Part A: Ecol. Genet. Physiol. 319, 10–22 (2013).CAS 
Article 

Google Scholar 
Nordtug, T. & Melø, T. B. Diurnal variations in natural light conditions at summer time in arctic and subarctic areas in relation to light detection in insects. Ecography 11, 202–209 (1988).Article 

Google Scholar 
Cohen, J. H. & Forward, R. B. Jr Diel vertical migration of the marine copepod Calanopia americana. II. Proximate role of exogenous light cues and endogenous rhythms. Mar. Biol. 147, 399–410 (2005).Article 

Google Scholar 
Maas, A. E., Blanco-Bercial, L., Lo, A., Tarrant, A. M. & Timmins-Schiffman, E. Variations in copepod proteome and respiration rate in association with diel vertical migration and circadian cycle. Biol. Bull. 000–000, https://doi.org/10.1086/699219 (2018).Berge, J. et al. Diel vertical migration of Arctic zooplankton during the polar night. Biol. Lett. 5, 69–72 (2009).PubMed 
Article 

Google Scholar 
Dale, T. & Kaartvedt, S. Diel patterns in stage-specific vertical migration of Calanus finmarchicus in habitats with midnight sun. ICES J. Mar. Sci. 57, 1800–1818 (2000).Article 

Google Scholar 
Hut, R. A., van Oort, B. E. H. & Daan, S. Natural entrainment without dawn and dusk: the case of the European Ground Squirrel (Spermophilus citellus). J. Biol. Rhythms 14, 290–299 (1999).CAS 
PubMed 
Article 

Google Scholar 
Williams, C. T., Barnes, B. M., Yan, L. & Buck, C. L. Entraining to the polar day: circadian rhythms in arctic ground squirrels. J. Exp. Biol. 220, 3095–3102 (2017).PubMed 
Article 

Google Scholar 
Daan, S. et al. Lab mice in the field: unorthodox daily activity and effects of a dysfunctional circadian clock allele. J. Biol. Rhythms 26, 118–129 (2011).PubMed 
Article 

Google Scholar 
Gattermann, R. et al. Golden hamsters are nocturnal in captivity but diurnal in nature. Biol. Lett. 4, 253–255 (2008).PubMed 
PubMed Central 
Article 

Google Scholar 
Green, E. W. et al. Drosophila circadian rhythms in seminatural environments: Summer afternoon component is not an artifact and requires TrpA1 channels. PNAS 112, 8702–8707 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Nagy, D. et al. A semi-natural approach for studying seasonal diapause in Drosophila melanogaster reveals robust photoperiodicity. J. Biol. Rhythms 33, 117–125 (2018).CAS 
PubMed 
Article 

Google Scholar 
Prabhakaran, P. M., De, J. & Sheeba, V. Natural conditions override differences in emergence rhythm among closely related Drosophilids. PLoS ONE 8, e83048 (2013).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Ruf, F. et al. Natural Zeitgebers under temperate conditions cannot compensate for the loss of a functional circadian clock in timing of a vital behavior in Drosophila. J. Biol. Rhythms 0748730421998112, https://doi.org/10.1177/0748730421998112 (2021).Dollish, H. K., Kaladchibachi, S., Negelspach, D. C. & Fernandez, F.-X. The Drosophila circadian phase response curve to light: conservation across seasonally relevant photoperiods and anchorage to sunset. Physiol. Behav. 245, 113691 (2022).CAS 
PubMed 
Article 

Google Scholar 
Shaw, B., Fountain, M. & Wijnen, H. Control of daily locomotor activity patterns in Drosophila suzukii by the circadian clock, light, temperature and social interactions. J. Biol. Rhythms 34, 463–481 (2019).CAS 
PubMed 
Article 

Google Scholar 
Chiesa, J. J., Aguzzi, J., García, J. A., Sardà, F. & de la Iglesia, H. O. Light intensity determines temporal niche switching of behavioral activity in deep-water Nephrops norvegicus (Crustacea: Decapoda). J. Biol. Rhythms 25, 277–287 (2010).PubMed 
Article 

Google Scholar 
DeCoursey, P. J. Light-sampling behavior in photoentrainment of a rodent circadian rhythm. J. Comp. Physiol. 159, 161–169 (1986).CAS 
Article 

Google Scholar 
Heard, E. Molecular biologists: let’s reconnect with nature. Nature 601, 9 (2022).CAS 
PubMed 
Article 

Google Scholar 
Deines, K. L. Backscatter estimation using Broadband acoustic Doppler current profilers. In Proc. IEEE Sixth Working Conference on Current Measurement 249–253 (1999).Darnis, G. et al. From polar night to midnight sun: diel vertical migration, metabolism and biogeochemical role of zooplankton in a high Arctic fjord (Kongsfjorden, Svalbard). Limnol. Oceanogr. 62, 1586–1605 (2017).CAS 
Article 

Google Scholar 
Cottier, F. R., Tarling, G. A., Wold, A. & Falk-Petersen, S. Unsynchronized and synchronized vertical migration of zooplankton in a high arctic fjord. Limnol. Oceanogr. 51, 2586–2599 (2006).Article 

Google Scholar 
Johnsen, G. et al. All-sky camera system providing high temporal resolution annual time series of irradiance in the Arctic. Appl. Opt. 60, 6456–6468 (2021).PubMed 
Article 

Google Scholar 
Pan, X. & Zimmerman, R. C. Modeling the vertical distributions of downwelling plane irradiance and diffuse attenuation coefficient in optically deep waters. J. Geophys. Res.: Oceans 115, C08016 (2010).
Google Scholar 
Hersbach, H. et al. The ERA5 global reanalysis. Q. J. R. Meteorol. Soc. 146, 1999–2049 (2020).Article 

Google Scholar 
Tidau, S. et al. Marine artificial light at night: An empirical and technical guide. Methods Ecol. Evol. 12, 1588–1601 (2021).Article 

Google Scholar 
Mobley, C. D. Light and Water: Radiative Transfer in Natural Waters (Academic Press Inc, 1994).Kostakis, I. et al. Development of a bio-optical model for the Barents Sea to quantitatively link glider and satellite observations. Philos. Trans. R. Soc. A: Math., Phys. Eng. Sci. 378, 20190367 (2020).CAS 
Article 

Google Scholar 
Buskey, E. J., Baker, K. S., Smith, R. C. & Swift, E. Photosensitivity of the oceanic copepods Pleuromamma gracilis and Pleuromamma xiphias and its relationship to light penetration and daytime depth distribution. Mar. Ecol. Prog. Ser. 55, 207–216 (1989).Article 

Google Scholar  More

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain
the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in
Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles
and JavaScript. More

Model structureWe used the PoPS (Pest or Pathogen Spread) Forecasting System11 version 2.0.0 to simulate the spread of SLF and calibrated the model (Fig. 6) using Approximate Bayesian Computation (ABC) with sequential Markov chain and a multivariate normal perturbation kernel18,19. We simulated the reproduction and dispersal of SLF groups (at the grid cell level) rather than individuals, as exact measures of SLF populations are not the goal of surveys conducted by USDA and state departments of agriculture. Reproduction was simulated as a Poisson process with mean β that is modified by local conditions. For example, if we have 5 SLF groups in a cell, a β value of 2.2, and a temperature coefficient of 0.7, our modified β value becomes 1.54 and we draw five numbers from a Poisson distribution with a λ value of 1.54. β and dispersal parameters were calibrated to fit the observed patterns of spread. For this application of PoPS, we replaced the long-distance kernel (α2) with a network dispersal kernel based on railroads, along which SLF and tree of heaven are commonly observed7. For each SLF group dispersing, if a railroad is in the grid cell with SLF, we used a Bernoulli distribution with mean of γ (probability of natural dispersal) to determine if an SLF group dispersed via the natural Cauchy kernel with scale (α) or along the rail network. This network dispersal kernel accounts for dispersal along railways if SLF is present in a cell containing a rail line. The network dispersal kernel added three new parameters to the PoPS model: a network file that contained the nodes and edges, minimum distance that each railcar travels, and the maximum distance that each railcar travels. Unlike typical network models, which simulate transport simply between nodes, our approach allows for SLF to disembark a railcar at any point along an edge, more closely mimicking their actual behavior. This network therefore captures the main pathway of SLF long-distance dispersal, i.e., along railways.Fig. 6: Model structure for spotted lanternfly (SLF, Lycorma delicatula).Unused modules in the PoPS model are gray in the equation. a The number of pests that disperse from a single host under optimal environmental conditions (β) is modified by the number of currently infested hosts (I) and environmental conditions in a location (i) at a particular time (t); environmental conditions include seasonality (X) and temperature (T) (see supplementary Fig. 3 for details on temperature). Dispersal is a function of gamma (γ), which is the probability of short-distance dispersal (alpha-1, α1) or long-distance via the rail network (N (dmin, dmax)). For the natural-distance Cauchy kernel, the direction is selected using 0-359 with 0 representing North. For the network kernel, the direction along the rail is selected randomly, and then travel continues in that direction until the drawn distance is reached. Once SLF has landed in a new location, its establishment depends on environmental conditions (X, T) and the availability of suitable hosts (number of susceptible hosts [S] divided by total number of potential hosts [N]). b We used a custom host map for tree of heaven (Ailanthus altissima) to determine the locations of susceptible hosts. The number of newly infested hosts (ψ) is predicted for each cell across the contiguous US.Full size imageSpotted lanternfly model calibrationWe used 2015–2019 data (over 300,000 total observations including both positive and negative surveys) provided by the USDA APHIS and the state Departments of Agriculture of Pennsylvania, New Jersey, Delaware, Maryland, Virginia, and West Virginia to calibrate model parameters (β, α1, γ, dmin, dmax). The calibration process starts by drawing a set of parameters from a uniform distribution. Simulated results for each model run are then compared to observed data within the year they were collected, and accuracy, precision, recall, and specificity are calculated for the simulation period. If each of these statistics is above 65% the parameter set is kept. This process repeats until 10,000 parameter sets are kept; then, the next generation of the ABC process begins: the mean of each accuracy statistic becomes the new accuracy threshold, and parameters are drawn from a multivariate normal distribution based on the means and covariance matrix of the first 10,000 kept parameters. This process repeats for a total of seven generations. Compared to the 2020 and 2021 observation data (over 100,000 total observations including both positive and negative surveys), the model performed well, with an accuracy of 84.4%, precision of 79.7%, recall of 91.55%, and specificity of 77.6%. In contrast, a model run using PoPS’ previous long-distance kernel (α2) instead of the network dispersal kernel had an accuracy of 76.5%, precision of 68.1%, recall of 92.68%, and specificity of 57.2%.We applied the calibrated parameters and their uncertainties (Fig. 7) to forecast the future spread of SLF, using the status of the infestation as of January 1, 2020 as a starting point and data for temperature and the distribution of SLF’s presumed primary host (tree of heaven, Ailanthus altissima) for the contiguous US at a spatial resolution of 5 km.Fig. 7: Parameter distributions.a Reproductive rate (β), b natural dispersal distance (α1), c percent natural dispersal (γ), d minimum distance (dmin), e maximum distance (dmax).Full size imageWeather dataOverwinter survival of SLF egg masses, and therefore spread, is sensitive to temperature (see ref. 2). To run a spread model in PoPS, all raw temperature values are first converted to indices ranging 0–1 to describe their impact on a species’ ability to survive and reproduce. We converted daily Daymet20 temperature into a monthly coefficient ranging 0–1 (Supplementary Fig. 1) and then rescaled from 1 to 5 km by averaging 1-km pixel values. We used weather data 1980–2019 and randomly drew from those historical data to simulate future weather conditions in our simulations, to account for uncertainty in future weather conditions.Tree of heaven distribution mappingSLF is known to feed on >70 species of mainly woody plants7, but tree of heaven is commonly viewed as necessary, or at least highly important, for SLF spread. Young nymphs are host generalists, but older nymphs and adults strongly prefer tree of heaven (in Korea21; in Pennsylvania, US22), and experiments in captivity23 and in situ9 have shown that adult survivorship is higher on the tree of heaven and grapevine than other host plants, likely due to the presence and proportion of sugar compounds important for SLF survival23. Secondary compounds found in tree of heaven also make adult SLF more unpalatable to avian predators24, and researchers have hypothesized that these protective compounds may be passed on to eggs21. For these reasons, tree of heaven is widely considered the primary host for SLF and linked to SLF spread1,25.We, therefore, used tree of heaven as the host in our spread forecast. We estimated the geographic range of tree of heaven using the Maximum Entropy (MaxEnt) model26,27. We chose to use niche modeling because tree of heaven has been in the US for over 200 years and is well past the early stage of invasion at which niche models perform poorly; instead, tree of heaven is well into the intermediate to equilibrium stage of invasion, when niche models perform well28. We obtained 19,282 presences for tree of heaven in the US from BIEN29,30 and EDDmaps31 and selected the most important variables from an initial MaxEnt model of all 19 WorldClim bioclimatic variables32. Our final climate variables were mean annual temperature, precipitation of the coldest quarter, and precipitation of the driest quarter. Given that tree of heaven is non-native and invasive in the US, prefers open and disturbed habitat, and is commonly found along roadsides and in urban landscapes33, we also included distance to major roads and railroads as an additional variable in our model, to account for the presence of disturbed habitat as well as approximate urbanization and anthropogenic degradation. For each 1-km cell in the extent, we calculated distance to the nearest road and nearest railroad using the US Census Bureau’s TIGER data set of primary roads and railroads34. We used our final MaxEnt model to generate the probability of the presence of tree of heaven for each 1-km cell, then reset all cells with a probability ≤0.2 to a value of 0 to minimize overprediction of the tree of heaven locations (because cells ≤0.2 contained less than 1% of the presences used to build the model). We rescaled the remaining probability values 0–1. We used 10% of the tree of heaven presence data to validate the model, which performed well: 95% of the validation data set locations had a probability of presence greater than 65%. We then rescaled the 1-km MaxEnt output to 5 km using the mean value of our 1-km cells, in order to reduce computational time.Forecasting spotted lanternflyWe used the Daymet temperature data and distribution of tree of heaven to simulate SLF spread with PoPS, assuming no further efforts to contain or eradicate either tree of heaven or SLF. We ran the spread simulation 10,000 times from 2020 to 2050 for the contiguous US. After running all 10,000 iterations, we created a probability of occurrence for each cell for each year by dividing the number of simulations in which a cell was simulated as being infested in that year by 10,000 (the total number of simulations). This gave us a probability of occurrence per year. We downscaled our probability of occurrence per year from 5 km to 1 km and set the probability to 0 in 1-km pixels with no tree of heaven occurrence.Data for mapping and comparisonWe compared our probability of occurrence map in 2050 to the SLF suitability map created by Wakie et al.1 using niche modeling to see how well the two modeling approaches would agree if SLF were allowed to spread unmanaged (Fig. 5). Wakie et al.1 categorized pixels below 8.359% as unsuitable, between 8.359% and 26.89% as low risk, between 26.89% and 51.99% as medium risk, and above 51.99% as high risk. To facilitate comparison, we used this same schema to categorize pixels as low, medium, or high probability of spread.We converted the yearly raster probability maps to county-level probabilities in order to examine the yearly risk to crops in counties. We performed this conversion using two methods: (1) the highest probability of occurrence in the county (Supplementary Movie 2) and (2) the mean probability of occurrence in the county (Fig. 1 and Supplementary Movie 1). The first method provides a simple, non-statistical estimate of the probability of SLF presence by assigning the county the value of the highest cell-level probability; the second accounts for all of the probabilities of the cells in the county and typically results in a higher county-level probability. We used USDA county-level production data10 for grapes, almonds, apples, walnuts, cherries, hops, peaches, plums, and apricots to determine the amount of production at risk each year (Fig. 2).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More