Philippa Kaur

1.Allen, T. et al. Global hotspots and correlates of emerging zoonotic diseases. Nat. Commun. 8, 1124 (2017).2.Gibb, R. et al. Zoonotic host diversity increases in human-dominated ecosystems. Nature 584, 398–402 (2020).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
3.Hassell, J. M., Begon, M., Ward, M. J. & Fèvre, E. M. Urbanization and disease emergence: dynamics at the wildlife–livestock–human interface. Trends Ecol. Evol. 32, 55–67 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
4.Plowright, R. K. et al. Pathways to zoonotic spillover. Nat. Rev. Microbiol. 15, 502–510 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Holmes, E. C., Rambaut, A. & Andersen, K. G. Pandemics: spend on surveillance, not prediction comment. Nature 558, 180–182 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
6.Tabachnick, W. J. Challenges in predicting climate and environmental effects on vector-borne disease episystems in a changing world. J. Exp. Biol. 213, 946–954 (2010).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
7.Franklinos, L. H. V., Jones, K. E., Redding, D. W. & Abubakar, I. The effect of global change on mosquito-borne disease. Lancet Inf. Dis. 19, e302–e312 (2019).Article 

Google Scholar 
8.Seabloom, E. W. et al. The community ecology of pathogens: coinfection, coexistence and community composition. Ecol. Lett. 18, 401–415 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
9.Johnson, P. T. J., De Roode, J. C. & Fenton, A. Why infectious disease research needs community ecology. Science 349, 1259504 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
10.Parker, I. M. et al. Phylogenetic structure and host abundance drive disease pressure in communities. Nature 520, 542–544 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Telfer, S. et al. Species interactions in a parasite community drive infection risk in a wildlife population. Science 330, 243–246 (2010).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
12.Dallas, T. A., Laine, A.-L. L. & Ovaskainen, O. Detecting parasite associations within multi-species host and parasite communities. Proc. R. Soc. B 286, 20191109 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
13.Weinstein, S., Titcomb, G., Agwanda, B., Riginos, C. & Young, H. Parasite responses to large mammal loss in an African savanna. Ecology 98, 1839–1848 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
14.Anderson, R. & May, R. Infectious Diseases of Humans: Dynamics and Control (Oxford Univ. Press, 1992).15.Keeling, M. J. & Rohani, P. Modeling Infectious Diseases in Humans and Animals (Princeton Univ. Press, 2011).16.Buhnerkempe, M. G. et al. Eight challenges in modelling disease ecology in multi-host, multi-agent systems. Epidemics 10, 26–30 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
17.Cross, P. C., Prosser, D. J., Ramey, A. M., Hanks, E. M. & Pepin, K. M. Confronting models with data: the challenges of estimating disease spillover. Philos. Trans. R. Soc. B 374, 20180435 (2019).Article 

Google Scholar 
18.Johnson, E. E., Escobar, L. E. & Zambrana-Torrelio, C. An ecological framework for modeling the geography of disease transmission. Trends Ecol. Evol. 34, 655–668 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
19.Warton, D. I. et al. So many variables: joint modeling in community ecology. Trends Ecol. Evol. 30, 766–779 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
20.Carlson, C. J. et al. The global distribution of Bacillus anthracis and associated anthrax risk to humans, livestock and wildlife. Nat. Microbiol. 4, 1337–1343 (2019).CAS 
PubMed 
Article 

Google Scholar 
21.Sutherland, W. J. Predicting the ecological consequences of environmental change: a review of the methods. J. Appl. Ecol. 43, 599–616 (2006).Article 

Google Scholar 
22.Getz, W. M. et al. Making ecological models adequate. Ecol. Lett. 21, 153–166 (2018).PubMed 
Article 

Google Scholar 
23.Carlson, C. J., Chipperfield, J. D., Benito, B. M., Telford, R. J. & O’Hara, R. B. Species distribution models are inappropriate for COVID-19. Nat. Ecol. Evol. 4, 770–771 (2020).PubMed 
Article 
PubMed Central 

Google Scholar 
24.Evans, M. R. et al. Predictive systems ecology. Proc. R. Soc. B 280, 20131452 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
25.Purves, D. & Pacala, S. Predictive models of forest dynamics. Science 320, 1452–1453 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
26.Harfoot, M. B. J. et al. Emergent global patterns of ecosystem structure and function from a mechanistic general ecosystem model. PLoS Biol. 12, e1001841 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
27.Lafferty, K. D. et al. Parasites in food webs: the ultimate missing links. Ecol. Lett. 11, 533–546 (2008).PubMed 
PubMed Central 
Article 

Google Scholar 
28.Redding, D. W. et al. Impacts of environmental and socio-economic factors on emergence and epidemic potential of Ebola in Africa. Nat. Commun. 10, 4531 (2019).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
29.Redding, D. W., Moses, L. M., Cunningham, A. A., Wood, J. & Jones, K. E. Environmental-mechanistic modelling of the impact of global change on human zoonotic disease emergence: a case study of Lassa fever. Methods Ecol. Evol. 7, 646–655 (2016).Article 

Google Scholar 
30.Carlson, C. J., Dallas, T. A., Alexander, L. W., Phelan, A. L. & Phillips, A. J. What would it take to describe the global diversity of parasites? Proc. R. Soc. B 287, 20201841 (2020).PubMed 
Article 
PubMed Central 

Google Scholar 
31.Rynkiewicz, E. C., Pedersen, A. B. & Fenton, A. An ecosystem approach to understanding and managing within-host parasite community dynamics. Trends Parasitol. 31, 212–221 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
32.Lello, J. & Hussell, T. Functional group/guild modelling of inter-specific pathogen interactions: a potential tool for predicting the consequences of co-infection. Parasitology 135, 825–839 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
33.Purves, D. et al. Time to model all life on Earth. Nature 493, 295–297 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
34.Kalka, M. B., Smith, A. R. & Kalko, E. K. V. Bats limit arthropods and herbivory in a tropical forest. Science 320, 71 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
35.Lafferty, K. D. et al. A general consumer-resource population model. Science 349, 854–857 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
36.Friedman, J., Higgins, L. M. & Gore, J. Community structure follows simple assembly rules in microbial microcosms. Nat. Ecol. Evol. 1, 0109 (2017).Article 

Google Scholar 
37.Goldford, J. E. et al. Emergent simplicity in microbial community assembly. Science 361, 469–474 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
38.Hatton, I. A. et al. The predator-prey power law: biomass scaling across terrestrial and aquatic biomes. Science 349, 6252 (2015).Article 
CAS 

Google Scholar 
39.Hatton, I. A., Dobson, A. P., Storch, D., Galbraith, E. D. & Loreau, M. Linking scaling laws across eukaryotes. Proc. Natl Acad. Sci. USA 116, 21616 (2019).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
40.Locey, K. J. & Lennon, J. T. Scaling laws predict global microbial diversity. Proc. Natl Acad. Sci. USA 113, 5970–5975 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
41.Godon, J. J., Arulazhagan, P., Steyer, J. P. & Hamelin, J. Vertebrate bacterial gut diversity: size also matters. BMC Ecol. 16, 12 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
42.Faust, C. L. et al. Null expectations for disease dynamics in shrinking habitat: dilution or amplification? Philos. Trans. R. Soc. B 372, 20160173 (2017).Article 

Google Scholar 
43.De Leo, G. A. & Dobson, A. P. Allometry and simple epidemic models for microparasites. Nature 379, 720–722 (1996).PubMed 
Article 
PubMed Central 

Google Scholar 
44.Strauss, A. T., Shoemaker, L. G., Seabloom, E. W. & Borer, E. T. Cross‐scale dynamics in community and disease ecology: relative timescales shape the community ecology of pathogens. Ecology 100, e02836 (2019).45.Handel, A. & Rohani, P. Crossing the scale from within-host infection dynamics to between-host transmission fitness: A discussion of current assumptions and knowledge. Philos. Trans. R. Soc. B 370, 20140302 (2015).46.Tibayrenc, M. & Ayala, F. J. Reproductive clonality of pathogens: a perspective on pathogenic viruses, bacteria, fungi, and parasitic protozoa. Proc. Natl Acad. Sci. USA 109, E3305–E3313 (2012).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
47.Gorter, F. A., Manhart, M. & Ackermann, M. Understanding the evolution of interspecies interactions in microbial communities. Philos. Trans. R. Soc. B 375, 1798 (2020).Article 
CAS 

Google Scholar 
48.Zeng, Q., Wu, S., Sukumaran, J. & Rodrigo, A. Models of microbiome evolution incorporating host and microbial selection. Microbiome 5, 127 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
49.Zeng, Q., Sukumaran, J., Wu, S. & Rodrigo, A. Neutral models of microbiome evolution. PLoS Comput. Biol. 11, e1004365 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
50.Liautaud, K., van Nes, E. H., Barbier, M., Scheffer, M. & Loreau, M. Superorganisms or loose collections of species? A unifying theory of community patterns along environmental gradients. Ecol. Lett. 22, 1243–1252 (2019).PubMed 
PubMed Central 

Google Scholar 
51.Coyte, K. Z., Schluter, J. & Foster, K. R. The ecology of the microbiome: Networks, competition, and stability. Science 350, 663–666 (2015).CAS 
Article 
PubMed 

Google Scholar 
52.Wright, E. S. & Vetsigian, K. H. Inhibitory interactions promote frequent bistability among competing bacteria. Nat. Commun. 7, 11274 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Sanchez, A. & Gore, J. Feedback between population and evolutionary dynamics determines the fate of social microbial populations. PLoS Biol. 11, e1001547 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
54.Lopatkin, A. J. & Collins, J. J. Predictive biology: modelling, understanding and harnessing microbial complexity. Nat. Rev. Microbiol. 18, 507–520 (2020).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
55.Wu, F. et al. A unifying framework for interpreting and predicting mutualistic systems. Nat. Commun. 10, 242 (2019).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
56.Restif, O. et al. Model-guided fieldwork: practical guidelines for multidisciplinary research on wildlife ecological and epidemiological dynamics. Ecol. Lett. 15, 1083–1094 (2012).PubMed 
PubMed Central 
Article 

Google Scholar 
57.Herzog, S. A., Blaizot, S. & Hens, N. Mathematical models used to inform study design or surveillance systems in infectious diseases: a systematic review. BMC Infect. Dis. 17, 1–10 (2017).Article 

Google Scholar 
58.Cotterill, G. G. et al. Winter feeding of elk in the Greater Yellowstone Ecosystem and its effects on disease dynamics. Philos. Trans. R. Soc. B 373, 20170093 (2018).Article 

Google Scholar 
59.Cross, P. C. et al. Estimating distemper virus dynamics among wolves and grizzly bears using serology and Bayesian state-space models. Ecol. Evol. 8, 8726–8735 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
60.Hopkins, J. B., Ferguson, J. M., Tyers, D. B. & Kurle, C. M. Selecting the best stable isotope mixing model to estimate grizzly bear diets in the Greater Yellowstone Ecosystem. PLoS ONE 12, e0174903 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
61.Schwartz, C. C. et al. Body and diet composition of sympatric black and grizzly bears in the Greater Yellowstone Ecosystem. J. Wildl. Manag. 78, 68–78 (2014).Article 

Google Scholar 
62.Chester, C. C. Yellowstone to Yukon: transborder conservation across a vast international landscape. Environ. Sci. Policy 49, 75–84 (2015).Article 

Google Scholar 
63.Young, H. S. et al. Interacting effects of land use and climate on rodent-borne pathogens in central Kenya. Philos. Trans. R. Soc. B 372, 20160116 (2017).Article 

Google Scholar 
64.Sitters, J., Kimuyu, D. M., Young, T. P., Claeys, P. & Olde Venterink, H. Negative effects of cattle on soil carbon and nutrient pools reversed by megaherbivores. Nat. Sustain. 3, 360–366 (2020).Article 

Google Scholar 
65.Sethi, S. S., Ewers, R. M., Jones, N. S., Orme, C. D. L. & Picinali, L. Robust, real‐time and autonomous monitoring of ecosystems with an open, low‐cost, networked device. Methods Ecol. Evol. 9, 2383–2387 (2018).Article 

Google Scholar 
66.Alfano, N., Dayaram, A. & Tsangaras, K. Non-invasive surveys of mammalian viruses using environmental DNA. Preprint at bioRxiv https://doi.org/10.1101/2020.03.26.009993 (2020)67.Coyte, K. Z. & Rakoff-Nahoum, S. Understanding competition and cooperation within the mammalian gut microbiome. Curr. Biol. 29, R538–R544 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
68.Murray, M. H. et al. Gut microbiome shifts with urbanization and potentially facilitates a zoonotic pathogen in a wading bird. PLoS ONE 15, e0220926 (2020).69.McCallum, H. I. et al. Does terrestrial epidemiology apply to marine systems? Trends Ecol. Evol. 19, 585–591 (2004).Article 

Google Scholar 
70.Wu, S., Carvalho, P. N., Müller, J. A., Manoj, V. R. & Dong, R. Sanitation in constructed wetlands: a review on the removal of human pathogens and fecal indicators. Sci. Total Environ. 541, 8–22 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
71.Lamb, J. B. et al. Seagrass ecosystems reduce exposure to bacterial pathogens of humans, fishes, and invertebrates. Science 355, 731–733 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
72.Janssen, M. A., Walker, B. H., Langridge, J. & Abel, N. An adaptive agent model for analysing co-evolution of management and policies in a complex rangeland system. Ecol. Model. 131, 249–268 (2000).Article 

Google Scholar 
73.Ngonghala, C. N. et al. General ecological models for human subsistence, health and poverty. Nat. Ecol. Evol. 1, 1153–1159 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
74.Hosseini, P. R. et al. Does the impact of biodiversity differ between emerging and endemic pathogens? The need to separate the concepts of hazard and risk. Philos. Trans. R. Soc. B 372, 20160129 (2017).Article 

Google Scholar 
75.Washburne, A. D. et al. Percolation models of pathogen spillover. Philos. Trans. R. Soc. B 374, 20180331 (2019).Article 

Google Scholar 
76.Olival, K. J. et al. Host and viral traits predict zoonotic spillover from mammals. Nature 546, 646–650 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
77.Dobson, A. et al. Habitat loss, trophic collapse, and the decline of ecosystem services. Ecology 87, 1915–1924 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
78.Faust, C. L. et al. Pathogen spillover during land conversion. Ecol. Lett. 21, 471–483 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
79.Sokolow, S. H. et al. Ecological interventions to prevent and manage zoonotic pathogen spillover. Philos. Trans. R. Soc. B 374, 20180342 (2019).Article 

Google Scholar 
80.Kauffman, M. J. et al. Landscape heterogeneity shapes predation in a newly restored predator-prey system. Ecol. Lett. 10, 690–700 (2007).PubMed 
Article 
PubMed Central 

Google Scholar 
81.Smith, D. W., Peterson, R. O. & Houston, D. B. Yellowstone after wolves. BioScience 53, 330–340 (2003).Article 

Google Scholar 
82.McNaughton, S. J. Ecology of a grazing ecosystem: the Serengeti. Ecol. Monogr. 55, 259–294 (1985).Article 

Google Scholar 
83.Atkins, J. L. et al. Cascading impacts of large-carnivore extirpation in an African ecosystem. Science 364, 173–177 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
84.Cross, P. C., Edwards, W. H., Scurlock, B. M., Maichak, E. J. & Rogerson, J. D. Effects of management and climate on elk brucellosis in the Greater Yellowstone Ecosystem. Ecol. Appl. 17, 957–964 (2007).PubMed 
Article 
PubMed Central 

Google Scholar 
85.Almberg, E. S., Cross, P. C. & Smith, D. W. Persistence of canine distemper virus in the Greater Yellowstone Ecosystem’s carnivore community. Ecol. Appl. 20, 2058–2074 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
86.Holdo, R. M. et al. A disease-mediated trophic cascade in the Serengeti and its implications for ecosystem C. PLoS Biol. 7, e1000210 (2009).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
87.Borer, E. T. in Unsolved Problems in Ecology (eds Dobson, A. P. et al.) 3–15 (Princeton Univ. Press, 2020).88.Kao, R. H. et al. NEON terrestrial field observations: designing continental-scale, standardized sampling. Ecosphere https://doi.org/10.1890/ES12-00196.1 (2012).89.Springer, Y. P. et al. Tick-, mosquito-, and rodent-borne parasite sampling designs for the National Ecological Observatory Network. Ecosphere 7, e01271 (2016).Article 

Google Scholar 
90.Anderson-Teixeira, K. J. et al. CTFS-ForestGEO: a worldwide network monitoring forests in an era of global change. Glob. Chang. Biol. 21, 528–549 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
91.Dobson, A. P. et al. Ecology and economics for pandemic prevention. Science 369, 379–381 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
92.Kress, W. J., Mazet, J. A. K. & Hebert, P. D. N. Opinion: intercepting pandemics through genomics. Proc. Natl Acad. Sci. USA 117, 202009508 (2020).
Google Scholar 
93.Durmuş, S. & Ülgen, K. Comparative interactomics for virus–human protein–protein interactions: DNA viruses versus RNA viruses. FEBS Open Bio 7, 96–107 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
94.Becker, D. J. & Albery, G. F. Expanding host specificity and pathogen sharing beyond viruses. Mol. Ecol. 29, 3170–3172 (2020).PubMed 
Article 
PubMed Central 

Google Scholar 
95.Rittershaus, E. S. C., Baek, S. H. & Sassetti, C. M. The normalcy of dormancy: common themes in microbial quiescence. Cell Host Microbe 13, 643–651 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
96.Pedersen, A. B. & Fenton, A. Emphasizing the ecology in parasite community ecology. Trends Ecol. Evol. 22, 133–139 (2007).PubMed 
Article 
PubMed Central 

Google Scholar 
97.Schmid-Hempel, P. Immune defence, parasite evasion strategies and their relevance for ‘macroscopic phenomena’ such as virulence. Philos. Trans. R. Soc. B 364, 85–98 (2009).Article 

Google Scholar 
98.Hekstra, D. R. & Leibler, S. Contingency and statistical laws in replicate microbial closed ecosystems. Cell 149, 1164–1173 (2012).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
99.Plowright, R. K. et al. Land use-induced spillover: a call to action to safeguard environmental, animal, and human health. Lancet Planet. 5, E237–E245 (2021).Article 

Google Scholar 
100.Barychka, T., Mace, G. & Purves, D. The Madingley General Ecosystem Model predicts bushmeat yields, species extinction rates and ecosystem-level impacts of bushmeat harvesting. Preprint at biorXiv https://doi.org/10.1101/2020.03.02.959718 (2020). More

General framework for eco-coevolutionary transitions from antagonism to mutualismWe develop a general framework in which we model interactions between host species i (density Hi) and its partner species k (density Fk), which are initially purely antagonistic. The model is general, but could be applied broadly to bacterial hosts and parasitic phages or plant hosts and animal or fungal partners, for example. The ecological dynamics of this community (without evolution) are given by:$$frac{d{H}_{i}}{{dt}}={g}_{i}{H}_{i}left(1-mathop{sum}limits_{j}{q}_{{ij}}{H}_{j}right)+{sum }_{k}{f}!_{ik}left[beta left({H}_{i},{F}_{k}right),alpha left({{H}_{i},F}_{k}right)right]$$
(4a)
$$frac{d{F}_{k}}{{dt}}=mathop{sum}limits_{i}{f}_{ki}left[beta left({H}_{i},{F}_{k}right),alpha left({{H}_{i},F}_{k}right)right]-{delta }_{k}{F}_{k}$$
(4b)
The first term of Eq. (4a) describes host population growth in the absence of partner species, where gi is its intrinsic per capita growth rate and qij is the competitive effect of host j on host i for other limiting factors. The general function fik describes the effects of interactions with partner k on host i: β(Hi, Fk) gives the potential mutualism and α(Hi, Fk) describes the antagonism. In Eq. (4b), the general function fki gives the effects of interactions with host i on partner k and δk is the partner’s per capita mortality rate.To derive an explicit eco-coevolutionary model, we apply Equation (4) to model interactions between a single host species and its exclusive partner species (for the sake of simplicity) in terms of host traits xi and partner traits yi (involved in interactions with host i); the ecological dynamics of which are given by:$$frac{1}{{H}_{i}}frac{d{H}_{i}}{{dt}}={g}_{i}left(1-{q}_{i}{H}_{i}right)+frac{bleft[{x}_{i}^{B}right]vleft[{x}_{i}^{V},{y}_{i}^{V}right]{F}_{k}}{{S}_{i}+vleft[{x}_{i}^{V},{y}_{i}^{V}right]{F}_{k}}-hleft[{x}_{i}^{H},{y}_{i}^{H}right]vleft[{x}_{i}^{V},{y}_{i}^{V}right]{F}_{k}$$
(5a)
$$frac{1}{{F}_{k}}frac{d{F}_{k}}{{dt}}=eleft[{y}_{i}^{V},{y}_{i}^{H}right]vleft[{x}_{i}^{V},{y}_{i}^{V}right]hleft[{x}_{i}^{H},{y}_{i}^{H}right]{H}_{i}-{delta }_{k}$$
(5b)
where b is the mutualistic benefits to the host, v is the visitation rate, Si is a saturation constant, h is the costs of antagonism to the host and its benefits to the partner, and e is the partner’s conversion efficiency. The mutualistic and antagonistic interactions are assumed to contribute additively to host population growth and multiplicatively to partner population growth, assumptions that may be valid for many types of interactions, but will not apply universally. To prevent unbounded population growth in the model, the effects of mutualism on host population growth are assumed to saturates with increasing partner density.The function b[xiB] gives the mutualistic benefits of the partner as a function of host trait xiB:$$bleft[{x}_{i}^{B}right]={b}_{{max },i}left(frac{2}{1+{e}^{-{B}_{i}^{{prime} }{x}_{i}^{B}}}-1right)$$
(6a)
where bmax,i gives the maximum mutualistic benefits and ({B}_{i}^{{prime} }) is a saturation constant. The interaction is purely antagonistic when xiB = 0. As xiB increases, the mutualistic benefits b[xiB] increase towards bmax,i.The function v[xiV, yiV] gives visitation rate as a sigmoid function of host trait xiV and partner trait yiV:$$vleft[{x}_{i}^{V},{y}_{i}^{V}right]=frac{{v}_{{max },i}}{1+{e}^{-{V}_{i}^{{prime} }left({x}_{i}^{V}+{y}_{i}^{V}right)}}$$
(6b)
where vmax,i is the maximum visitation rate and ({V}_{i}^{{prime} }) determines how rapidly visitation rate changes as host and partner traits change. As xiV or yiV increase, the visitation rate increases and approaches vmax,i when xiV + yiV → ∞. As xiV or yiV decrease, the visitation rate decreases and approaches zero when xiV + yiV → −∞. Negative values of xiV indicate that the host species is reducing its attraction of the partner species.The function h[xiH, yiH] gives the costs of antagonism to the host and its benefits to the partner, which is described via a sigmoid function of the difference between host trait xiH and partner trait yiH:$$hleft[{x}_{i}^{H},{y}_{i}^{H}right]=frac{{h}_{{max },i}}{1+{e}^{{H}_{i}^{{prime} }left({x}_{i}^{H}-{y}_{i}^{H}right)}}$$
(6c)
where hmax,i gives the maximum antagonism and ({H}_{i}^{{prime} }) determines how antagonism changes as the difference between host and partner traits increases. When xiH > yiH, antagonism declines and approaches zero when xiH – yiH → ∞, while when xiH < yiH, antagonism increases and approaches hmax,i when xiH – yiH → -∞ (unlike xiV, xiH cannot be negative).Partner traits yiV and yiH trade off with conversion efficiency via the function e[yiV, yiH] as defined by:$$eleft[{y}_{i}^{V},{y}_{i}^{H}right]={e}_{{max },i}{e}^{-left({c}_{I,i}^{V}{left({y}_{i}^{V}right)}^{2}+{c}_{I,i}^{H}{left({y}_{i}^{H}right)}^{2}right)}$$ (6d) where emax,i is the maximum conversion efficiency when interacting with host i (when yiV = yiH = 0), and cI,iV and cI,iH determine how rapidly conversion efficiency declines as yiV or yiH increase, thus quantifying the costliness of traits yiV and yiH, respectively. This trade-off shape was chosen because it is unimodal and constrains conversion efficiency to always be positive. Host trade-offs are defined below (Eq. 8c).Host-partner coevolutionary dynamicsWe model coevolution via the adaptive dynamics framework17,18. Coevolution of a mutant host trait ximut and partner trait yimut (for any general traits xi and yi) is given by:$$frac{d{{x}_{i}}^{{mut}}}{dtau }={mu }_{x}{left.frac{partial {W}_{H}left({{x}_{i}}^{{mut}},{x}_{i},{y}_{i}right)}{partial {{x}_{i}}^{{mut}}}right|}_{{{x}_{i}}^{{mut}}={x}_{i}}$$ (7a) $$frac{d{{y}_{i}}^{{mut}}}{dtau }={mu }_{y}{left.frac{partial {W}_{F}left({{y}_{i}}^{{mut}},{y}_{i},{x}_{i}right)}{partial {{y}_{i}}^{{mut}}}right|}_{{{y}_{i}}^{{mut}}={y}_{i}}$$ (7b) where τ is the evolutionary timescale, μx and μy give, respectively, the rates of host and partner evolution, and WH(ximut,xi,yi) and WF(yimut,yi,xi) are the invasion fitness (per capita growth rate when rare) of a mutant host and partner species with trait ximut and yimut in a resident community with trait xi and yi, respectively. The partial derivatives ({left.partial {W}_{H}/partial {{x}_{i}}^{{mut}}right|}_{{{x}_{i}}^{{mut}}={x}_{i}}) and ({left.partial {W}_{F}/partial {{y}_{i}}^{{mut}}right|}_{{{y}_{i}}^{{mut}}={y}_{i}}) are the selection gradients.We model coevolution of mutualistic benefits from the focal partner species (via b), attraction (via v), and defense (via h). The invasion fitness of the mutant host and a mutant partner are given by:$${W}_{H}={g}_{i}left(1-qleft[{{x}_{i}}^{{mut}},{x}_{i}right]{{H}_{i}}^{ast }right)+frac{bleft[{x}_{i}^{B,{mut}}right]vleft[{x}_{i}^{V,{mut}},{y}_{i}^{V}right]{{F}_{k}}^{ast }}{{S}_{i}+vleft[{x}_{i}^{V,{mut}},{y}_{i}^{V}right]{{F}_{k}}^{ast }}-hleft[{x}_{i}^{H,{mut}},{y}_{i}^{H}right]{vleft[{x}_{i}^{V,{mut}},{y}_{i}^{V}right]{F}_{k}}^{ast }$$ (8a) $${W}_{F}=eleft[{y}_{i}^{V,{mut}},{y}_{i}^{H,{mut}}right]vleft[{x}_{i}^{V},{y}_{i}^{V,{mut}}right]hleft[{x}_{i}^{H},{y}_{i}^{H,{mut}}right]{{H}_{i}}^{ast }-{delta }_{k}$$ (8b) where Hi* and Fk* are species’ densities at the ecological equilibrium (of Eq. 5). The functions b, v, h, and e are given by Eq. (6a–d), respectively, where xi and yi are replaced with ximut in Eq. (8a) and yimut in Eq. (8b). The function q[ximut,xi] describes trade-offs between mutant host traits and mutant host competitive ability as defined by:$$qleft[{{x}_{i}}^{{mut}},{x}_{i}right]=1+{c}_{H,i}^{B}left({left({{x}_{i}}^{B,{mut}}right)}^{{s}_{i}^{B}}-{left({{x}_{i}}^{B}right)}^{{s}_{i}^{B}}right)+{c}_{H,i}^{V}left({left({{x}_{i}}^{V,{mut}}right)}^{{s}_{i}^{V}}-{left({{x}_{i}}^{V}right)}^{{s}_{i}^{V}}right)+{c}_{H,i}^{H}left({left({{x}_{i}}^{H,{mut}}right)}^{{s}_{i}^{H}}-{left({{x}_{i}}^{H}right)}^{{s}_{i}^{H}}right)$$ (8c) If ximut > xi for any trait, the competitive effect experienced by the mutant host is increased by an amount taken to be proportional (for simplicity) to the difference between the trait values, ximut – xi, whereas if ximut < xi, the competitive effect experienced by the mutant host is decreased by that amount. The coefficients cH,iB, cH,iV, and cH,iH measure the costs associated with the trade-off for each trait, while the shape parameters siB, siV, and siH define whether the trade-offs are linear (si = 1), concave (si < 1), or convex (si > 1).Mutualism can evolve via the COA for all trade-off shapes (Supplementary Fig. 3). Parameter space plots show that the interaction transitions from antagonism to net mutualism when the costs associated with host traits underlying attraction (cH,iV) and defense (cH,iH) are within a range beyond which there is evolutionary purging of the partner (Supplementary Fig. 3a–c). Only with convex trade-offs can the net antagonism persist. The coevolution of mutualism also requires that the costs associated with partner traits underlying visitation (cFV) and antagonism (cFH) exceed a threshold (Supplementary Fig. 3d–f) below which there is evolutionary purging of the partner (linear or convex trade-offs) or the net antagonism persists (linear or concave trade-offs). Coevolution of mutualism occurs across greater parameter ranges when the trade-offs are linear or slightly concave because costs increase less rapidly than with convex trade-offs.Ecological model of plant-insect interactionsWe tailor the general model (Eq. 4) to model populations of D. wrightii (density Pw) and D. discolor (density Pd) interacting with M. sexta. We scale the model so that Pi = 1 in the absence of M. sexta: thus, Pi >1 indicates that pollination benefits exceed herbivory costs, and Pi < 1 indicates that herbivory costs exceed pollination benefits. The Datura species do not rely obligately on M. sexta and, consistent with ecology of the natural community (Box 1), the model incorporates the alternative host plant, Proboscidea parviflora (density Pp), and the alternative nectar source, Agave palmeri. The ecological dynamics of this community (without evolution) are given by:$$frac{1}{{P}_{i}}frac{d{P}_{i}}{{dt}}=left(1-{P}_{i}right)+frac{{b}_{i}{v}_{i}A}{H+{v}_{i}A}-{h}_{i}{L}_{i}$$ (9a) $$frac{d{L}_{i}}{{dt}}=varepsilon {e}_{i}{v}_{i}{P}_{i}A-{m}_{i}{h}_{i}{L}_{i}-{d}_{i}{L}_{i}$$ (9b) $$frac{{dA}}{{dt}}=mathop{sum}limits_{i}{rho }_{i}{m}_{i}{h}_{i}{L}_{i}-{d}_{A}A$$ (9c) Equation (9a) describes the population dynamics of plant species i (D. wrightii, D. discolor, or P. parviflora). Equation (9b,c) give the dynamics of M. sexta: Li gives the larvae density on plant species i, which recruit into the adult population, A. Pollination is described by the term biviA/(H+viA), where bi is the per capita growth of plant species i due to pollination by the antagonist, vi is the visitation rate to plant species i per antagonist adult, and H is the saturation constant for pollination. Oviposition is given by εeiviPiA, where ei is the oviposition efficiency (number of eggs laid per floral visit) and ε is the fractional increase in egg production due to nectar-feeding at A. palmeri. Floral visits lead to both pollination and oviposition because these behaviors have been shown to be tightly linked in M. sexta19. Pollination and oviposition are given by saturating and linear functions, respectively, based on our data (Supplementary Data 1). Herbivory damage is given by the term hiLi, where hi is the herbivory rate per larvae on plant species i. Larvae mature at rate mihiLi, where mi is the maturation efficiency (fraction of larvae maturing on plant species i). Larval mortality on plant species i is di, adult mortality is dA, and ρi is pupae survival (due to data constraints, we include pupae survival in our estimates of maturation mi, set ρi = 1, and drop ρi from equations hereafter). Equation (9a) gives the dynamics of the alternative larval host plant, P. parviflora (bp = 0 and cannot evolve), which can coevolve attraction and defense. The alternative nectar source, A. palmeri, is incorporated within the model via the parameter ε.Model scalingWithout the antagonist, plant population growth is given by gi (1 – qiPi), where gi is the per capita growth rate of plant species i due to autonomous self-pollination or pollination by other species and qi is plant self-limitation. As qi is very difficult to quantify in nature, we scale the model so that Pi = 1 without the antagonist. We scale plant density ((hat{{P}_{i}}={q}_{i}{P}_{i})), larvae density ((hat{{L}_{i}}={q}_{i}{L}_{i})), herbivory rate ((hat{{h}_{i}}={h}_{i}/{q}_{i})), maturation efficiency ((hat{{m}_{i}}={q}_{i}{m}_{i})), and survival of pupae ((hat{{rho }_{i}}={rho }_{i}/{q}_{i})); where the hats denote scaled quantities and are dropped elsewhere for clarity. Thus, the model is scaled for parameterization, but is not non-dimensionalized. We then scale gi to 1 such that pollination benefits, bi, are estimated by the ratio of the seed set of moth-pollinated flowers to autonomously self-pollinated flowers. Parameter estimates are for scaled quantities.Interaction breakdown boundary for ancestral interaction in a one-plant species communityFor the ancestral insect to persist, its per capita growth rate must be positive when it is rare (i.e., at Pi* = 1, Li* = 0, A* = 0). In stage-structured models, the per capita growth rate is given by the dominant eigenvalue (λD) of the matrix:$$left[begin{array}{cc}-{m}_{i}{h}_{i}-{d}_{i} & varepsilon {e}_{i}{v}_{i}{{P}_{i}}^*\ {m}_{i}{h}_{i} & -{d}_{A}end{array}right]$$which is given by:$${lambda }_{D}=frac{1}{2}left(-{d}_{A}-{d}_{i}-{m}_{i}{h}_{i}+sqrt{{({d}_{A}+{d}_{i}+{m}_{i}{h}_{i})}^{2}-4({d}_{A}left({d}_{i}+{m}_{i}{h}_{i}right)-varepsilon {e}_{i}{v}_{i}{m}_{i}{h}_{i})}right).$$For the insect to persist, λD must have a positive real part, which occurs only when the second term in the square root of λD is negative; i.e., ({d}_{A}left({d}_{i}+{m}_{i}{h}_{i}right)-varepsilon {e}_{i}{v}_{i}{m}_{i}{h}_{i} , , 1). Applying ({f}_{i}=frac{varepsilon {e}_{i}{v}_{i}}{{d}_{A}}) and ({s}_{i}=frac{{m}_{i}{h}_{i}}{{{m}_{i}{h}_{i}+d}_{i}}), where fi is insect lifetime fecundity and si is the larval success (probability of larvae maturing rather than dying), yields Eq. (1).Interaction transition boundary in a one-plant species communityFor the interaction to transition from antagonism to mutualism, equilibrium plant density, Pi* must exceed one (see “Model scaling”). Setting Eq. (9b) to zero and solving for Pi* yields:({{P}_{i}}^{ast }=frac{{{m}_{i}{h}_{i}+d}_{i}}{varepsilon {e}_{i}{v}_{i}}left(frac{{{L}_{i}}^{ast }}{{A}^{ast }}right)). Setting Eq. (9c) to zero and rearranging terms then yields: (frac{{{L}_{i}}^{ast }}{{A}^{ast }}=frac{{d}_{A}}{{m}_{i}{h}_{i}}). Thus, ({{P}_{i}}^{ast }=frac{{{m}_{i}{h}_{i}+d}_{i}}{varepsilon {e}_{i}{v}_{i}}left(frac{{d}_{A}}{{m}_{i}{h}_{i}}right)) and (rearranging slightly) the condition for mutualism to arise is: ({{P}_{i}}^{ast }=left(frac{{d}_{A}}{varepsilon {e}_{i}{v}_{i}}right)left(frac{{{m}_{i}{h}_{i}+d}_{i}}{{m}_{i}{h}_{i}}right) , > , 1). Rearranging and applying ({f}_{i}=frac{varepsilon {e}_{i}{v}_{i}}{{d}_{A}}) and ({s}_{i}=frac{{m}_{i}{h}_{i}}{{{m}_{i}{h}_{i}+d}_{i}}) yields Eq. (2).Interaction breakdown boundary in a one-plant species communityIn the ancestral interaction, insect persistence is evaluated by whether or not it can increase from low density, which yields Eq. (1). Within the net mutualistic region, however, the insect cannot increase from very low density because it cannot buoy plant density sufficiently to maintain a positive per capita growth rate (mathematically, Eq. 1 cannot hold when Eq. 2 is satisfied). The mutualistic region is thus characterized by bistability (see Supplementary Figure 1), and the interaction breakdown boundary is determined by the conditions for the coexistence equilibrium to exist. At the coexistence equilibrium, the larval and adult densities are: ({{L}_{i}}^{ast }=frac{-B+sqrt{{B}^{2}-4{A}_{L}{C}_{L}}}{2{A}_{L}}) and ({A}^{ast }=frac{-B+sqrt{{B}^{2}-4{A}_{A}{C}_{A}}}{{2A}_{A}}), where ({A}_{L}=varepsilon {e}_{i}{{v}_{i}}^{2}{h}_{i}{d}_{A}), ({A}_{A}=varepsilon {e}_{i}{{v}_{i}}^{2}{m}_{i}{{h}_{i}}^{2}), (B=varepsilon {e}_{i}{{v}_{i}}^{2}{m}_{i}{h}_{i}left(frac{1}{{f}_{i}{s}_{i}}+frac{H}{{v}_{i}{m}_{i}}-left(1+{b}_{i}right)right)), ({C}_{L}=varepsilon {e}_{i}{v}_{i}H{d}_{A}left(frac{1}{{f}_{i}{s}_{i}}-1right)), and ({C}_{A}=varepsilon {e}_{i}{v}_{i}{m}_{i}{h}_{i}Hleft(frac{1}{{f}_{i}{s}_{i}}-1right)). For the coexistence equilibrium to exist, either CL and CA must be negative or B must be negative and Li* and A* must be real. CL and CA are negative when fi si > 1, which is Eq. (1) and cannot hold within the mutualistic region because Eq. (2) must be satisfied. However, B is negative when ({f}_{i}{s}_{i}left(left(1+{b}_{i}right)-frac{H}{{v}_{i}{m}_{i}}right) , > , 1), which is approximated by Eq. (3) when the last term is assumed to be small. For Li* and A* to be real, B2 – 4ALCL > 0 and B2 – 4AACA > 0. Assuming that the pollination saturation constant is small (i.e., H ≈ 0) yields CL ≈ CA ≈ 0 such that ({{L}_{i}}^{ast }approx frac{-B}{{A}_{L}}approx frac{{m}_{i}}{{d}_{A}{f}_{i}{s}_{i}}left({f}_{i}{s}_{i}left(1+{b}_{i}right)-1right)) and ({A}^{ast }approx frac{-B}{{A}_{A}}approx frac{1}{{h}_{i}}left({f}_{i}{s}_{i}left(1+{b}_{i}right)-1right)), which are both positive when fi si (1 + bi) > 1 as approximated by Eq. (3).Interaction transition and breakdown boundaries in a two-plant species communityThese boundaries are analytically intractable and are estimated by simulation (see codes provided online).Coevolutionary dynamics of plants and insectThe effects of plant traits xi and insect traits yi on the ecological dynamics of the interactions are given by:$$frac{1}{{P}_{i}}frac{d{P}_{i}}{{dt}}=left(1-{P}_{i}right)+frac{bleft[{x}_{i}^{B}right]vleft[{x}_{i}^{V},{y}_{i}^{V}right]A}{H+vleft[{x}_{i}^{V},{y}_{i}^{V}right]A}-hleft[{x}_{i}^{H},{y}_{i}^{H}right]{L}_{i}$$
(10a)
$$frac{d{L}_{i}}{{dt}}=varepsilon eleft[{y}_{i}^{V},{y}_{i}^{H}right]vleft[{x}_{i}^{V},{y}_{i}^{V}right]{P}_{i}A-{m}_{i}hleft[{x}_{i}^{H},{y}_{i}^{H}right]{L}_{i}-{d}_{i}{L}_{i}$$
(10b)
$$frac{{dA}}{{dt}}=mathop{sum}limits_{i}{m}_{i}hleft[{x}_{i}^{H},{y}_{i}^{H}right]{L}_{i}-{d}_{A}A$$
(10c)
We model coevolution of plant-insect interactions using the adaptive dynamics framework17,18 to link population dynamics and trait coevolution. The coevolution of mutant plant trait xmut and insect trait ymut (for general traits x and y) is given by Equation (7). We model the coevolution of pollination benefits from the antagonist, bi (via mutant plant trait xiB,mut), attraction (via mutant plant trait xiV,mut and mutant insect trait yiV,mut), and defense (via mutant plant trait xiH,mut and mutant insect trait yiH,mut). The invasion fitness of a mutant plant is given by:$${W}_{P,i}left({x}_{i}^{{mut}},{x}_{i},{y}_{i}right)=left(1-qleft[{x}_{i}^{{mut}},{x}_{i}right]{{P}_{i}}^{ast }right)+frac{bleft[{x}_{i}^{B,{mut}}right]vleft[{x}_{i}^{V,{mut}},{y}_{i}^{V}right]{A}^{ast }}{H+vleft[{x}_{i}^{V,{mut}},{y}_{i}^{V}right]{A}^{ast }}-hleft[{x}_{i}^{H,{mut}},{y}_{i}^{H}right]{{L}_{i}}^{ast }$$
(11a)
where Pi*, Li*, and A* are the densities of the plant, insect larvae per plant, and insect adults, respectively, at the ecological equilibrium (of Eq. 10). The functions (b[x_{i}^{B,mut}], v[x_{i}^{V,mut} , y_{i}^{V}],) and (h[x_{i}^{H,mut}, y_{i}^{H}]), describe the effects of mutant plant traits (x_{i}^{B,mut}), (x_{i}^{V,mut}), (x_{i}^{H,mut}), and (x_{i}^{H,mut}) on pollination benefits, attraction, and defense, respectively, which are defined by Eq. (6a–c), where xi is replaced with ximut (where the plant is the host species and the insect is the partner species). The function q[x,mut, xi] defines the trade-offs between mutant plant traits and the competitive ability of mutant plants, which is given by Eq. (8c) (with si = 1). At a coESS, ximut = xi for all traits such that q[ximut, xi] = 1 and the original definition of Pi >1 indicating that pollination benefits exceed herbivory costs is retained when pollination benefits evolve.Invasion fitness of a mutant insect is given by the dominant eigenvalue of its system of equations evaluated at the resident equilibrium. In a one-plant species community, the insect invasion fitness is:$${W}_{I,i}=frac{1}{2}left(-{d}_{A}-{d}_{i}-{m}_{i}{h}_{i}^{{mut}}+sqrt{{left({d}_{A}+{d}_{i}+{m}_{i}{h}_{i}^{{mut}}right)}^{2}-4left({d}_{A}left({d}_{i}+{m}_{i}{h}_{i}^{{mut}}right)-frac{varepsilon {e}_{i}^{{mut}}{v}_{i}^{{mut}}{h}_{i}^{{mut}}{d}_{A}left({d}_{i}+{m}_{i}{h}_{i}right)}{{e}_{i}{v}_{i}{h}_{i}}right)}right)$$
(11b)
where vimut, himut, and eimut are functions describing the effects of mutant insect traits on attraction, defense, and mutant oviposition efficiency, respectively, which are given by Eq. (6b–d), where yi is replaced with yimut. Invasion fitness of a mutant insect in a two-plant species community is given by the dominant eigenvalue of its system of equations evaluated at the resident equilibrium, which is analytically tractable, but sufficiently complicated that we do not include it here (see codes provided online).The curves where the selection gradients (see Eqs. 7) become zero give the evolutionary isoclines for the coevolutionary system. The points where the isoclines intersect give the coevolutionary singularities, which are coevolutionary stable states (coESSs) when they are stable for both plants and the insect. For tractability, the local stability of the coevolutionary singularities was assessed by carefully inspecting the selection gradient of each trait in the neighborhood of its coESS with all other traits held at their coESS as well as by simulating coevolutionary dynamics. Importantly, all three plant traits (xiB, xiV, and xiH) and both insect traits (yiV and yiH) all coevolve simultaneously in the model.Coevolution of the ancestral antagonistic interactionIn the ancestral interaction, pollination by the antagonist is impossible (bi = 0) and thus visitation only contributes to oviposition. From the plant perspective, the selection gradients for attraction and defense in the ancestral interaction are given by:$${left.frac{partial {W}_{P,i}}{partial {x}_{i}^{V,{mut}}}right|}_{{x}_{i}^{{mut}}={x}_{i}}=-{c}_{P,i}^{V}{{P}_{i}}^{ast }$$
(12a)
$${left.frac{partial {W}_{P,i}}{partial {x}_{i}^{H,{mut}}}right|}_{{x}_{i}^{{mut}}={x}_{i}}=frac{{h}_{{max },i}{H}_{i}^{{prime} }{e}^{{H}_{i}^{{prime} }left({x}_{i}^{H}-{y}_{i}^{H}right)}}{{left(1+{e}^{{H}_{i}^{{prime} }left({x}_{i}^{H}-{y}_{i}^{H}right)}right)}^{2}}{{L}_{i}}^{ast }-{c}_{P,i}^{H}{{P}_{i}}^{ast }$$
(12b)
Equation (12a) predicts that selection favors plant traits that reduce attracting the antagonist (e.g., reduced production of volatiles) and lower costs associated with competitive ability. We constrain xiV to be non-negative in the ancestral interaction so that xiV = 0 at the coESS; otherwise, xiV → –∞ and the plant always purges the insect given this model parameterization. Selection balances reduced herbivory damage (first term of Eq. 12b) with costs of reduced competitive ability (second term of Eq. 12b). Selection gradients for insect traits are sufficiently complicated that we do not include them here (see codes provided online); however, selection balances traits that increase visitation and overcome plant defenses with the costs associated with reduced oviposition. The ancestral coESSs are given in Supplementary Table 3.Coevolution of pollination benefits, attraction, and defenseThe evolution of mutant plant traits that allow the antagonist to pollinate it (bimut > 0) initiates the evolution of pollination benefits from the antagonist. The selection gradient for pollination benefits from the antagonist is given by:$${left.frac{partial {W}_{P,i}}{partial {x}_{i}^{B,{mut}}}right|}_{{x}_{i}^{{mut}}={x}_{i}}=frac{2{b}_{{max },i}{B}_{i}^{{prime} }{e}^{-{B}_{i}^{{prime} }{x}_{i}^{B}}vleft[{x}_{i}^{V},{y}_{i}^{V}right]{A}^{ast }}{{left(1+{e}^{-{B}_{i}^{{prime} }{x}_{i}^{B}}right)}^{2}left(H+vleft[{x}_{i}^{V},{y}_{i}^{V}right]{A}^{ast }right)}-{c}_{P,i}^{B}{{P}_{i}}^{ast }$$
(13a)
Equation (13a) shows that plants evolve traits to benefit from floral visits by the antagonist when selection for increased pollination benefits (first term of Eq. 13a) exceeds the costs associated with reduced competitive ability (second term of Eq. 13a).In the model, pollination benefits from the antagonist evolve via Eq. (13a) simultaneously with plant and insect traits affecting attraction and defense. The plant selection gradient for attraction is now:$${left.frac{partial {W}_{P,i}}{partial {x}_{i}^{V,{mut}}}right|}_{{x}_{i}^{{mut}}={x}_{i}}=frac{bleft[{x}_{i}^{B}right]{v}_{{max },i}{V}_{i}^{{prime} }{e}^{-{V}_{i}^{{prime} }left({x}_{i}^{V}+{y}_{i}^{V}right)}H{A}^{ast }}{{left(Hleft(1+{e}^{-{V}_{i}^{{prime} }left({x}_{i}^{V}+{y}_{i}^{V}right)}right)+{v}_{{max },i}{A}^{ast }right)}^{2}}-{c}_{P,i}^{V}{{P}_{i}}^{ast }$$
(13b)
The co-option of the antagonist has fundamentally changed selection on attraction (Eq. 13b vs. Equation 12a), which now balances traits affecting attraction (first term of Eq. 13b) with the costs of reduced competitive ability (second term of Eq. 13b). Co-option of the antagonist also modifies selection on defense (which is still given by Eq. 12b) by changing both trait values and equilibrium densities.Model parameterizationAll ecological parameters are estimated from empirical data. Here we parameterize the saturation constant H, maturation efficiency mi, larval mortality di, and adult mortality dA as well as the parameters for the alternative larval host plant and the alternative nectar source (see “Model validation” for other parameters).We cannot fit the saturation constant H to data because seed set saturates with even a single floral visit. We therefore estimate H as follows: D. wrightii flowers have a 91% chance of setting fruit30; thus, ({v}_{w}A/(H+{v}_{w}A))= 0.91 for a single visit (({v}_{w})A = 1). Solving (1/(H+1))= 0.91 for H yields: H = 0.1. H is assumed to be the same for D. discolor as pollination benefits saturate with a single visit for D. discolor. For maturation efficiency mi, only 0.5% of M. sexta larvae on D. wrightii survive through the final larval instar in nature34; thus, mw = 0.005. As M. sexta suffers 40% lower larval survival on D. discolor (5/8 larvae surviving to pupation) than on D. wrightii (10/10 larvae surviving to pupation) in our experiment19, we estimate that maturation efficiency is ~40% lower on D. discolor than on D. wrightii; i.e., md = (1 – 0.4)mw = 0.003. To estimate larval mortality, we note that larval survival is given by: ({m}_{i}={e}^{-{d}_{i}{D}_{i}}), where Di is development time. M. sexta has a larval stage of ~20 days on D. wrightii35 and there is no difference in development on D. wrightii and D. discolor, at least to the 5th instar19. Solving for di yields: dw ≈ 0.25 and dd ≈ 0.3. Finally, adults live ~5 days in the wild36. Assuming adult mortality is roughly the inverse of the lifespan: dA ≈ 0.2.For the alternative larval host plant, females lay similar numbers of eggs on D. wrightii and P. parviflora34; thus, visitation rate and oviposition efficiency are assumed to be the same as with D. wrightii; i.e., vp = vw and ep = ew. Because P. parviflora plants are of similar size and architecture as D. wrightii34, we assume that herbivory rate on P. parviflora is the same as on D. wrightii; i.e., hp = hw. (see “Model validation” for estimates of vw, ew, and hw). Only 1% of M. sexta larvae on P. parviflora survive through the final larval stage34; thus, mp = 0.01. As larvae have roughly the same development time on P. parviflora as on D. wrightii (~20 days37), solving ({m}_{p}={e}^{-{d}_{p}{D}_{p}}) yields an estimate of larval mortality on P. parviflora of: dp ≈ 0.25.For the alternative nectar source, A. palmeri provides M. sexta with copious amounts of nectar that females likely utilize for egg production38. M. sexta females lay 100–300 eggs/night39. If females foraging exclusively on D. wrightii lay the minimum 100 eggs/night and females that also forage at A. palmeri lay the maximum 300 eggs/night, then A. palmerii is estimated to increase oviposition by a factor of: ε = 3.Model validationPollination benefits (bi), visitation rate (vi), herbivory rate (hi), and oviposition efficiency (ei) all evolve simultaneously in the model. We independently validate the coESSs predicted by the models whenever possible by estimating these parameters using data that were not used to parameterize the models. We estimate bi via the ratio of the seed set of moth-pollinated flowers to autonomously self-pollinated flowers (autonomously self-pollinated seeds germinate as readily as do outcrossed seeds;30). Pollinated D. wrightii and D. discolor flowers set bw = 4.6 ± 0.2 and bd = 3.6 ± 0.1 times more seeds, respectively, than do autonomously self-pollinated flowers (D. wrightii: n = 21 fruit; D. discolor: n = 85 fruit). Moths averaged vw = 4.3 ± 0.6 floral visits to D. wrightii (n = 89 plants) and vd = 2.4 ± 0.4 floral visits to D. discolor (n = 33 plants) in our experiment19. Estimating the herbivory rate is very difficult in nature; however, we can make cursory estimates based on our data. A single M. sexta larvae can consume 1400–1900 cm2 of leaves, which is more than many D. wrightii plants in nature30. Assuming that an average D. wrightii plant supplies larvae with 1400 cm2 of leaves, the variation in leaf consumption (500 cm2) represents ~0.4 plants (=500/1400). Thus, M. sexta larvae are estimated to consume: hw ≈ 1 ± 0.4 D. wrightii plants. M. sexta larvae consumed roughly two times more D. discolor leaf biomass than D. wrightii leaf biomass based on our cursory estimates from our experiments; thus, hd = 2hw ≈ 2 ± 0.8. We estimate oviposition efficiency by the slope of a linear regression of the number of eggs versus the number of floral visits that each plant received from each female moth in our experiments19, which yields: ew = 0.6 ± 0.1 (n = 34 plants) and ed = 0.6 ± 0.2 (n = 24 plants) (Supplementary Data 1).Estimating evolutionary model parametersDirectly estimating evolutionary parameters with data is not possible. We therefore use theory to predict how key parameters affect eco-coevolutionary outcomes and to select reasonable parameter estimates. Our approach is as follows. We set the rates of plant and insect evolution to one (μx = μy = 1); these rates affect the speed of evolution, but not the coESSs. For each trait, we need to estimate the maximum value (bmax,i, vmax,i, hmax,i, and emax,i), the coefficient (({R}_{i}^{{prime} }), ({V}_{i}^{{prime} }), and ({H}_{i}^{{prime} })), and the associated costs (cP,iB, cP,iV, and cP,iH for plant i and cI,iV and cI,iH for the insect). Maximum trait values were chosen to constrain coevolution to a realistic range. We set the coefficients ({R}_{i}^{{prime} }), ({V}_{i}^{{prime} }), and ({H}_{i}^{{prime} }) to one for simplicity because the exact value of any trait x and y are themselves somewhat arbitrary. The costs associated with the traits therefore largely determine the coevolutionary outcomes in the model.We estimate the costs of each trait by systematically varying the costs of plant traits in the one-plant species community given reasonable values for the insect costs and then systematically varying the costs of insect traits while holding plant costs constant at their chosen values (Fig. 5). Parameter space plots show that the interactions transition from antagonism to net mutualism provided that the costs associated with insect traits underlying visitation (cI,iV) exceed a threshold below which the plant and insect engage in an evolutionary arms-race that results in the evolutionary purging of the antagonist (Fig. 5a, b). Only very rarely does the net antagonism persist. We assigned all insect traits a cost of 0.5 (black points in Fig. 5a, b) and then systematically vary the costs of plant traits associated with attraction and defense.Parameter space plots show that interactions transition from antagonism to net mutualism when the costs associated with defense are high relative to the costs associated with attraction (cP,iH > cP,iV); otherwise, coevolution drives evolutionary purging of the antagonist (Fig. 5c, d). When the costs associated with attraction and defense are both fairly high, the net antagonism persists. We assigned values of cP,iH and cP,iV to D. wrightii and D. discolor such that the parameters for D. discolor are closer to the threshold at which evolutionary purging occurs than are those of D. wrightii (Fig. 5d vs. 5c), reflecting the smaller range of ecological parameters over which M. sexta can persist with D. discolor versus with D. wrightii (Fig. 2b vs. 2a). Finally, the costs associated with pollination benefits from the antagonist (cP,iB) must be very high for the net antagonism to persist and we never observed evolutionary purging of the insect within the range of values used (see codes provided online). We assigned values of cP,iB so that pollination benefits to D. wrightii and D. discolor are well below their maximum values. Our estimates of evolutionary parameters are reported in Supplementary Table 2. Evolutionary parameters for P. parviflora are set equal to D. discolor because, in the absence of more information, both species are annual plants that may face broadly similar evolutionary constraints, at least relative to the perennial D. wrightii.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.IPCC Special Report on the Ocean and Cryosphere in a Changing Climate (eds Pörtner, H.-O. et al.) (IPCC, 2019).2.Urban, M. C. Accelerating extinction risk from climate change. Science 348, 571–573 (2015).CAS 
Article 

Google Scholar 
3.McCauley, D. J. & Pinsky, M. L. Marine defaunation: animal loss in the global ocean. Science 347, 1255641 (2015).Article 
CAS 

Google Scholar 
4.Hoffmann, A. A. & Sgrò, C. M. Climate change and evolutionary adaptation. Nature 470, 479–485 (2011).CAS 

Google Scholar 
5.Somero, G. N. The physiology of climate change: how potentials for acclimatization and genetic adaptation will determine ‘winners’ and ‘losers’. J. Exp. Biol. 213, 912–920 (2010).CAS 
Article 

Google Scholar 
6.Kearney, M. & Porter, W. Mechanistic niche modelling: combining physiological and spatial data to predict species’ ranges. Ecol. Lett. 12, 334–350 (2009).Article 

Google Scholar 
7.Foden, W. B. et al. Identifying the world’s most climate change vulnerable species: a systematic trait-based assessment of all birds, amphibians and corals. PLoS ONE 8, e65427 (2013).CAS 
Article 

Google Scholar 
8.West, J. M. & Salm, R. V. Resistance and resilience to coral bleaching: implications for coral reef conservation and management. Conserv. Biol. 17, 956–967 (2003).Article 

Google Scholar 
9.Baskett, M. L., Nisbet, R. M., Kappel, C. V., Mumby, P. J. & Gaines, S. D. Conservation management approaches to protecting the capacity for corals to respond to climate change: a theoretical comparison. Glob. Change Biol. 16, 1229–1246 (2010).Article 

Google Scholar 
10.Beyer, H. L. et al. Risk-sensitive planning for conserving coral reefs under rapid climate change. Conserv. Lett. 11, e12587 (2018).Article 

Google Scholar 
11.Walsworth, T. E. et al. Management for network diversity speeds evolutionary adaptation to climate change. Nat. Clim. Change 9, 632–636 (2019).Article 

Google Scholar 
12.Hughes, T. P. et al. Spatial and temporal patterns of mass bleaching of corals in the Anthropocene. Science 359, 80–83 (2018).CAS 
Article 

Google Scholar 
13.Donner, S. D., Skirving, W. J., Little, C. M., Oppenheimer, M. & Hoegh-Guldberg, O. Global assessment of coral bleaching and required rates of adaptation under climate change. Glob. Change Biol. 11, 2251–2265 (2005).Article 

Google Scholar 
14.Frieler, K. et al. Limiting global warming to 2 °C is unlikely to save most coral reefs. Nat. Clim. Change 3, 165–170 (2012).Article 

Google Scholar 
15.Van Hooidonk, R., Maynard, J. A., Manzello, D. & Planes, S. Opposite latitudinal gradients in projected ocean acidification and bleaching impacts on coral reefs. Glob. Change Biol. 20, 103–112 (2014).Article 

Google Scholar 
16.Logan, C. A., Dunne, J. P., Eakin, C. M. & Donner, S. D. Incorporating adaptive responses into future projections of coral bleaching. Glob. Change Biol. 20, 125–139 (2014).Article 

Google Scholar 
17.Bay, R. A., Rose, N. H., Logan, C. A. & Palumbi, S. R. Genomic models predict successful coral adaptation if future ocean warming rates are reduced. Sci. Adv. 3, e1701413 (2017).Article 

Google Scholar 
18.Matz, M. V., Treml, E. A. & Haller, B. C. Estimating the potential for coral adaptation to global warming across the Indo-West Pacific. Glob. Change Biol. 26, 3473–3481 (2020).Article 

Google Scholar 
19.Baskett, M. L., Gaines, S. D. & Nisbet, R. M. Symbiont diversity may help coral reefs survive moderate climate change. Ecol. Appl. 19, 3–17 (2009).Article 

Google Scholar 
20.Matz, M. V., Treml, E. A., Aglyamova, G. V. & Bay, L. K. Potential and limits for rapid genetic adaptation to warming in a Great Barrier Reef coral. PLoS Genet. 14, e1007220 (2018).Article 
CAS 

Google Scholar 
21.Muscatine, L., Falkowski, P. G., Porter, J. W. & Dubinsky, Z. Fate of photosynthetic fixed carbon in light- and shade-adapted colonies of the symbiotic coral Stylophora pistillata. Proc. R. Soc. Lond. B. 222, 181–202 (1984).CAS 
Article 

Google Scholar 
22.Csaszar, N. B., Ralph, P. J., Frankham, R., Berkelmans, R. & van Oppen, M. J. Estimating the potential for adaptation of corals to climate warming. PLoS ONE 5, e9751 (2010).Article 
CAS 

Google Scholar 
23.Howells, E. J. et al. Coral thermal tolerance shaped by local adaptation of photosymbionts. Nat. Clim. Change 2, 116–120 (2012).Article 

Google Scholar 
24.Buerger, P. et al. Heat-evolved microalgal symbionts increase coral bleaching tolerance. Sci. Adv. 6, eaba2498 (2020).CAS 
Article 

Google Scholar 
25.Baker, A. C. Flexibility and specificity in coral–algal symbiosis: diversity, ecology, and biogeography of Symbiodinium. Annu. Rev. Ecol. Evol. Syst. 34, 661–689 (2003).26.Berkelmans, R. & van van Oppen, M. J. H. The role of zooxanthellae in the thermal tolerance of corals: a ‘nugget of hope’ for coral reefs in an era of climate change. Proc. R. Soc. Lond. B 273, 2305–2312 (2006).
Google Scholar 
27.National Academies of Sciences, Engineering, and Medicine A research Review of Interventions to Increase the Persistence and Resilience of Coral Reefs (National Academies Press, 2019).28.National Academies of Sciences, Engineering, and Medicine A Decision Framework for Interventions to Increase the Persistence and Resilience of Coral Reefs (National Academies Press, 2019).29.Darling, E. S., Alvarez-Filip, L., Oliver, T. A., McClanahan, T. R. & Côté, I. M. Evaluating life-history strategies of reef corals from species traits. Ecol. Lett. 15, 1378–1386 (2012).Article 

Google Scholar 
30.Chan, N. C. S. & Connolly, S. R. Sensitivity of coral calcification to ocean acidification: a meta-analysis. Glob. Change Biol. 19, 282–290 (2013).Article 

Google Scholar 
31.Hughes, T. P. et al. Global warming transforms coral reef assemblages. Nature 556, 492–496 (2018).CAS 
Article 

Google Scholar 
32.Darling, E. S. et al. Relationships between structural complexity, coral traits, and reef fish assemblages. Coral Reefs 36, 561–575 (2017).Article 

Google Scholar 
33.Howells, E. J. et al. Corals in the hottest reefs in the world exhibit symbiont fidelity not flexibility. Mol. Ecol. 29, 899–911 (2020).CAS 
Article 

Google Scholar 
34.Madin, J. S., Hughes, T. P. & Connolly, S. R. Calcification, storm damage and population resilience of tabular corals under climate change. PLoS ONE 7, e46637 (2012).CAS 
Article 

Google Scholar 
35.Hoegh-Guldberg, O. et al. in Global Warming of 1.5°C (eds Masson-Delmotte, V. et al.) Ch. 3 (IPCC, 2018).36.Darling, E. S. et al. Social–environmental drivers inform strategic management of coral reefs in the Anthropocene. Nat. Ecol. Evol. 3, 1341–1350 (2019).Article 

Google Scholar 
37.Wilkinson, C. R. Global and local threats to coral reef functioning and existence: review and predictions. Mar. Freshw. Res. 50, 867–878 (1999).
Google Scholar 
38.Palumbi, S. R., Barshis, D. J., Traylor-Knowles, N. & Bay, R. A. Mechanisms of reef coral resistance to future climate change. Science 344, 895–898 (2014).CAS 
Article 

Google Scholar 
39.Kleypas, J. A. et al. Larval connectivity across temperature gradients and its potential effect on heat tolerance in coral populations. Glob. Change Biol. 22, 3539–3549 (2016).Article 

Google Scholar 
40.Heron, S. F. et al. Validation of reef-scale thermal stress satellite products for coral bleaching monitoring. Remote Sens. 8, 59 (2016).Article 

Google Scholar 
41.Safaie, A. et al. High frequency temperature variability reduces the risk of coral bleaching. Nat. Commun. 9, 1671 (2018).Article 
CAS 

Google Scholar 
42.Forster, P. M. et al. Evaluating adjusted forcing and model spread for historical and future scenarios in the CMIP5 generation of climate models. J. Geophys. Res. Atmos. 118, 1139–1150 (2013).Article 

Google Scholar 
43.Ziegler, M., Eguíluz, V. & Duarte, C. et al. Rare symbionts may contribute to the resilience of coral–algal assemblages. ISME J 12, 161–172 (2018).Article 

Google Scholar 
44.Sampayo, E. M., Ridgway, T., Bongaerts, P. & Hoegh-Guldberg, O. Bleaching susceptibility and mortality of corals are determined by fine-scale differences in symbiont type. Proc. Natl Acad. Sci. USA 105, 10444–10449 (2008).CAS 
Article 

Google Scholar 
45.Thornhill, D. J., Xiang, Y. U., Fitt, W. K. & Santos, S. R. Reef endemism, host specificity and temporal stability in populations of symbiotic dinoflagellates from two ecologically dominant Caribbean corals. PLoS ONE 4, e6262 (2009).Article 
CAS 

Google Scholar 
46.Stat, M., Loh, W. K. W., LaJeunesse, T. C., Hoegh-Guldberg, O. & Carter, D. A. Stability of coral–endosymbiont associations during and after a thermal stress event in the southern Great Barrier Reef. Coral Reefs 28, 709–713 (2009).Article 

Google Scholar 
47.Chakravarti, L. J., Beltran, V. H. & van Oppen, M. J. Rapid thermal adaptation in photosymbionts of reef-building corals. Glob. Change Biol. 23, 4675–4688 (2017).Article 

Google Scholar 
48.Schulte, P. M., Healy, T. M. & Fangue, N. A. Thermal performance curves, phenotypic plasticity, and the time scales of temperature exposure. Integr. Comp. Biol. 51, 691–702 (2011).Article 

Google Scholar 
49.Loya, Y. et al. Coral bleaching: the winners and the losers. Ecol. Lett. 4, 122–131 (2001).Article 

Google Scholar 
50.Langmead, O. & Sheppard, C. Coral reef community dynamics and disturbance: a simulation model. Ecol. Modell. 175, 271–290 (2004).Article 

Google Scholar 
51.Chancerelle, Y. Methods to estimate actual surface areas of scleractinian coral at the colony- and community-scale. Oceanol. Acta 23, 211–219 (2000).Article 

Google Scholar 
52.Falkowski, P. G., Dubinsky, Z., Muscatine, L. & Porter, J. W. Light and the bioenergetics of a symbiotic coral. Bioscience 34, 705–709 (1984).CAS 
Article 

Google Scholar 
53.Huston, M. Variation in coral growth rates with depth at Discovery Bay, Jamaica. Coral Reefs 4, 19–25 (1985).Article 

Google Scholar 
54.Hoogenboom, M., Beraud, E. & Ferrier-Pagès, C. Relationship between symbiont density and photosynthetic carbon acquisition in the temperate coral Cladocora caespitosa. Coral Reefs 29, 21–29 (2010).Article 

Google Scholar 
55.Cunning, R. & Baker, A. C. Not just who, but how many: the importance of partner abundance in reef coral symbioses. Front. Microbiol. 5, 400 (2014).Article 

Google Scholar 
56.McClanahan, T., Muthiga, N. & Mangi, S. Coral and algal changes after the 1998 coral bleaching: interaction with reef management and herbivores on Kenyan reefs. Coral Reefs 19, 380–391 (2001).Article 

Google Scholar 
57.Fitt, W. K., McFarland, F. K., Warner, M. E. & Chilcoat, G. C. Seasonal patterns of tissue biomass and densities of symbiotic dinoflagellates in reef corals and relation to coral bleaching. Limnol. Oceanogr. 45, 677–685 (2000).CAS 
Article 

Google Scholar 
58.Eppley, R. W. Temperature and phytoplankton growth in the sea. Fish. Bull. 70, 1063–1085 (1972).
Google Scholar 
59.Jon, N. Biodiversity and ecosystem functioning: a complex adaptive systems approach. Limnol. Oceanogr. 49, 1269–1277 (2004).Article 

Google Scholar 
60.Mousseau, T. A. & Roff, D. A. Natural selection and the heritability of fitness components. Heredity 59, 181–197 (1987).Article 

Google Scholar 
61.Lynch, M. The rate of polygenic mutation. Genet. Res. 51, 137–148 (1988).CAS 
Article 

Google Scholar 
62.Donner, S. D., Rickbeil, G. J. & Heron, S. F. A new, high-resolution global mass coral bleaching database. PLoS ONE 12, e0175490 (2017).Article 
CAS 

Google Scholar 
63.Cunning, R., Gillette, P., Capo, T., Galvez, K. & Baker, A. C. Growth tradeoffs associated with thermotolerant symbionts in the coral Pocillopora damicornis are lost in warmer oceans. Coral Reefs 34, 155–160 (2015).Article 

Google Scholar 
64.Silverstein, R. N., Cunning, R. & Baker, A. C. Change in algal symbiont communities after bleaching, not prior heat exposure, increases heat tolerance of reef corals. Glob. Change Biol. 21, 236–249 (2015).Article 

Google Scholar 
65.Dunne, J. P. et al. GFDL’s ESM2 global coupled climate–carbon Earth system models part I: physical formulation and baseline simulation characteristics. J. Clim. 25, 6646–6665 (2012).Article 

Google Scholar 
66.Dunne, J. P. et al. GFDL’s ESM2 global coupled climate–carbon Earth system models Part II: carbon system formulation and baseline simulation characteristics. J. Clim. 26, 2247–2267 (2012).Article 

Google Scholar 
67.Lough, J. M. & Barnes, D. J. Environmental controls on growth of the massive coral Porites. J. Exp. Mar. Biol. Ecol. 245, 225–243 (2000).CAS 
Article 

Google Scholar 
68.UNEP-WCMC, WorldFish Centre, WRI & TNC. Global Distribution of Warm-water Coral Reefs, Compiled From Multiple Sources Including the Millennium Coral Reef Mapping Project. Version 4.1 (UN Environment World Conservation Monitoring Centre. Data, 2021); https://doi.org/10.34892/t2wk-5t3469.van Hooidonk, R., Maynard, J. A. & Planes, S. Temporary refugia for coral reefs in a warming world. Nat. Clim. Change 3, 508–511 (2013).Article 
CAS 

Google Scholar 
70.Fitt, W., Brown, B., Warner, M. & Dunne, R. Coral bleaching: interpretation of thermal tolerance limits and thermal thresholds in tropical corals. Coral Reefs 20, 51–65 (2001).Article 

Google Scholar 
71.González-Espinosa, P. C. & Donner, S. D. Predicting cold-water bleaching in corals: role of temperature, and potential integration of light exposure. Mar. Ecol. Prog. Ser. 642, 133–146 (2020).Article 

Google Scholar  More

1.Lamichhaney, S. et al. Evolution of Darwin’s finches and their beaks revealed by genome sequencing. Nature 518, 371–375 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Simões, M. et al. The evolving theory of evolutionary radiations. Trends Ecol. Evol. 31, 27–34 (2016).PubMed 
Article 

Google Scholar 
3.Arnegard, M. E. et al. Genetics of ecological divergence during speciation. Nature 511, 307–311 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
4.Hoffmann, A. A. & Sgrò, C. M. Climate change and evolutionary adaptation. Nature 470, 479–485 (2011).CAS 
PubMed 
Article 

Google Scholar 
5.Olsen, J. L. et al. The genome of the seagrass Zostera marina reveals angiosperm adaptation to the sea. Nature 530, 331–335 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
6.Becerra, J. X., Nogeb, K. & Venable, D. L. Macroevolutionary chemical escalation in an ancient plant–herbivore arms race. Proc. Natl Acad. Sci. USA 106, 18062–18066 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
7.Edger, P. P. et al. The butterfly plant arms-race escalated by gene and genome duplications. Proc. Natl Acad. Sci. USA 112, 8362–8366 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Adler, L. S. & Bronstein, J. L. Attracting antagonists: does floral nectar increase leaf herbivory? Ecology 85, 1519–1526 (2004).Article 

Google Scholar 
9.McCall, A. C. & Irwin, R. E. Florivory: the intersection of pollination and herbivory. Ecol. Lett. 9, 1351–1365 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
10.Cook, J. M. & Rasplus, J.-Y. Mutualists with attitude: coevolving fig wasps and figs. Trends Ecol. Evol. 18, 241–248 (2003).Article 

Google Scholar 
11.Zhang, X. et al. Genomes of the banyan tree and pollinator wasp provide insights into fig–wasp coevolution. Cell 183, 875–889 (2020).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
12.Herre, E. A., Jandér, K. C. & Machado, C. A. Evolutionary ecology of figs and their associates: recent progress and outstanding puzzles. Annu. Rev. Ecol. Evol. Syst. 39, 439–458 (2008).Article 

Google Scholar 
13.Souza, C. D. et al. Diversity of fig glands is associated with nursery mutualism in fig trees. Am. J. Bot. 102, 1564–1577 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
14.Souto-Vilarós, D. et al. Pollination along an elevational gradient mediated both by floral scent and pollinator compatibility in the fig and fig–wasp mutualism. J. Ecol. Evol. 106, 2256–2273 (2018).
Google Scholar 
15.Wang, R. et al. Loss of top-down biotic interactions changes the relative benefits for obligate mutualists. Proc. R. Soc. B 286, 20182501 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
16.Mori, K. et al. Identification of RAN1 orthologue associated with sex determination through whole genome sequencing analysis in fig (Ficus carica L.). Sci. Rep. 7, 41124 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
17.Proffit, M. et al. Chemical signal is in the blend: bases of plant–pollinator encounter in a highly specialized interaction. Sci. Rep. 10, 10071 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Chen, C. et al. Private channel: a single unusual compound assures specific pollinator attraction in Ficus semicordata. Funct. Ecol. 23, 941–950 (2009).Article 

Google Scholar 
19.Wang, G., Cannon, C. H. & Chen, J. Pollinator sharing and gene flow among closely related sympatric dioecious fig taxa. Proc. R. Soc. B 283, 20152963 (2016).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
20.Yu, H. et al. De novo transcriptome sequencing in Ficus hirta Vahl. (Moraceae) to investigate gene regulation involved in the biosynthesis of pollinator attracting volatiles. Tree Genet. Genomes 11, 91 (2015).Article 

Google Scholar 
21.Soler, C. C. L., Proffit, M., Bessière, J.-M., Hossaert -McKey, M. & Schatz, B. Evidence for intersexual chemical mimicry in a dioicous plant. Ecol. Lett. 15, 978–985 (2012).PubMed 
Article 

Google Scholar 
22.Volf, M. et al. Community structure of insect herbivores is driven by conservatism, escalation and divergence of defensive traits in Ficus. Ecol. Lett. 21, 83–92 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
23.Martinson, E. O., Hackett, J. D., Machado, C. A. & Arnold, A. E. Metatranscriptome analysis of fig flowers provides insights into potential mechanisms for mutualism stability and gall induction. PLoS ONE 10, e0130745 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
24.Zhang, H. et al. Leaf-mining by Phyllonorycter blancardella reprograms the host–leaf transcriptome to modulate phytohormones associated with nutrient mobilization and plant defense. J. Insect Physiol. 84, 114–127 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
25.Schultz, J. C., Edger, P. P., Body, M. & Appel, H. M. A galling insect activates plant reproductive programs during gall development. Sci. Rep. 9, 1833 (2019).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
26.The Nasonia Genome Working Group Functional and evolutionary insights from the genomes of three parasitoid Nasonia species. Science 327, 343–348 (2010).PubMed Central 
Article 
CAS 

Google Scholar 
27.Xiao, J.-H. et al. Obligate mutualism within a host drives the extreme specialization of a fig wasp genome. Genome Biol. 14, R141 (2013).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
28.Ohri, D. & Khoshoo, T. N. Nuclear DNA contents in the genus Ficus (Moraceae). Plant Syst. Evol. 156, 1–4 (1987).Article 

Google Scholar 
29.Chen, Y., Compton, S. G., Liu, M. & Chen, X.-Y. Fig trees at the northern limit of their range: the distributions of cryptic pollinators indicate multiple glacial refugia. Mol. Ecol. 21, 1687–1701 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
30.Chen, L.-G. et al. Binding affinity characterization of an antennae-enriched chemosensory protein from the white-backed planthopper, Sogatella furcifera (Horváth), with host plant volatiles. Pestic. Biochem. Phys. 152, 1–7 (2018).Article 
CAS 

Google Scholar 
31.Gu, S.-H. et al. Functional characterization and immunolocalization of odorant binding protein 1 in the lucerne plant bug, Adelphocoris lineolatus (GOEZE). Arch. Insect Biochem. Physiol. 77, 81–99 (2011).CAS 
PubMed 
Article 

Google Scholar 
32.Leal, W. S. et al. Reverse and conventional chemical ecology approaches for the development of oviposition attractants for Culex mosquitoes. PLoS ONE 3, e3045 (2008).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
33.Rizzo, W. B. et al. Fatty aldehyde and fatty alcohol metabolism: review and importance for epidermal structure and function. Biochim. Biophys. Acta 1841, 377–389 (2014).CAS 
PubMed 
Article 

Google Scholar 
34.Schwab, W., Davidovich‐Rikanati, R. & Lewinsohn, E. Biosynthesis of plant‐derived flavor compounds. Plant J. 54, 712–732 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
35.Capella, M., Ribone, P. A., Arce, A. L. & Chan, R. L. Arabidopsis thaliana HomeoBox 1 (AtHB1), a Homedomain-Leucine Zipper I (HD-Zip I) transcription factor, is regulated by PHYTOCHROME-INTERACTING FACTOR 1 to promote hypocotyl elongation. New Phytol. 207, 669–682 (2015).CAS 
PubMed 
Article 

Google Scholar 
36.Jiang, W. et al. Two transcription factors TaPpm1 and TaPpb1 co-regulate anthocyanin biosynthesis in purple pericarps of wheat. J. Exp. Bot. 69, 2555–2567 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Guan, R. et al. Draft genome of the living fossil Ginkgo biloba. GigaScience 5, 49 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
38.Mithöfer, A. & Boland, W. Plant defense against herbivores: chemical aspects. Annu. Rev. Plant Biol. 63, 431–450 (2012).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
39.Salazar, D. et al. Origin and maintenance of chemical diversity in a species-rich tropical tree lineage. Nat. Ecol. Evol. 2, 983–990 (2018).PubMed 
Article 

Google Scholar 
40.Després, L., David, J. P. & Gallet, C. The evolutionary ecology of insect resistance to plant chemicals. Trends Ecol. Evol. 22, 298–307 (2007).PubMed 
Article 

Google Scholar 
41.Segar, S. T., Volf, M., Sisol, M., Pardikes, N. & Souto-Vilarós, A. D. Chemical cues and genetic divergence in insects on plants: conceptual cross pollination between mutualistic and antagonistic systems. Curr. Opin. Insect Sci. 32, 83–90 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
42.Cook, J. M. & Segar, S. T. Speciation in fig wasps. Ecol. Entomol. 35, 54–66 (2010).Article 

Google Scholar 
43.Cruaud, A. et al. An extreme case of plant–insect codiversification: figs and fig-pollinating wasps. Syst. Biol. 61, 1029–1047 (2012).PubMed 
PubMed Central 
Article 

Google Scholar 
44.Yu, H. et al. Multiple parapatric pollinators have radiated across a continental fig tree displaying clinal genetic variation. Mol. Ecol. 28, 2391–2405 (2019).PubMed 
Article 

Google Scholar 
45.Satler, J. D. et al. Inferring processes of coevolutionary diversification in a community of Panamanian strangler figs and associated pollinating wasps. Evolution 73, 2295–2311 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
46.Wang, G. et al. Genomic evidence of prevalent hybridization throughout the evolutionary history of the fig–wasp pollination mutualism. Nat. Commun. 12, 718 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
47.Hoballah, M. E. et al. Single gene-mediated shift in pollinator attraction in Petunia. Plant Cell 19, 779–790 (2007).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
48.Potts, S. G. et al. Global pollinator declines: trends, impacts and drivers. Trends Ecol. Evol. 25, 345–353 (2010).PubMed 
Article 

Google Scholar 
49.Kiers, E. T., Palmer, T. M., Ives, A. R., Bruno, J. F. & Bronstein, J. L. Mutualisms in a changing world: an evolutionary perspective. Ecol. Lett. 13, 1459–1474 (2010).Article 

Google Scholar 
50.Segar, S. T. et al. The role of evolution in shaping ecological networks. Trends Ecol. Evol. 35, 454–466 (2020).PubMed 
Article 

Google Scholar 
51.Stoy, K. S., Gibson, A. K., Gerardo, N. M. & Morran, L. T. A need to consider the evolutionary genetics of host–symbiont mutualisms. J. Evol. Biol. 33, 1656–1668 (2020).PubMed 
Article 

Google Scholar 
52.Marçais, G. & Kingsford, C. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics 27, 764–770 (2011).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
53.Xiao, C.-L. et al. MECAT: fast mapping, error correction, and de novo assembly for single-molecule sequencing reads. Nat. Methods 14, 1072–1074 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
54.Walker, B. J. et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS ONE 9, e112963 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
55.Pryszcz, L. P. & Gabaldón, T. Redundans: an assembly pipeline for highly heterozygous genomes. Nucleic Acids Res. 44, e113 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
56.Zhang, Z., Schwartz, S., Wagner, L. & Miller, W. A greedy algorithm for aligning DNA sequences. J. Comput. Biol. 7, 203–214 (2000).CAS 
PubMed 
Article 

Google Scholar 
57.English, A. C. et al. Mind the gap: upgrading genomes with Pacific Biosciences RS long-read sequencing technology. PLoS ONE 7, e47768 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
58.Sahlin, K., Chikhi & Arvestad, R. L. Assembly scaffolding with PE-contaminated mate-pair libraries. Bioinformatics 32, 1925–1932 (2016).PubMed 
Article 

Google Scholar 
59.Belton, J. M. et al. Hi-C: a comprehensive technique to capture the conformation of genomes. Methods 58, 268–276 (2012).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
60.Servant, N. et al. HiC-Pro: an optimized and flexible pipeline for Hi-C data processing. Genome Biol. 16, 259 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
61.Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows–Wheeler Transform. Bioinformatics 25, 1754–1760 (2009).CAS 
PubMed 
PubMed Central 

Google Scholar 
62.Durand, N. C. et al. Juicer provides a one-click system for analyzing loop-resolution Hi-C experiments. Cell Syst. 3, 95–98 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
63.Dudchenko, O. et al. De novo assembly of the Aedes aegypti genome using Hi-C yields chromosome-length scaffolds. Science 356, 92–95 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
64.Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, E. M. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210–3212 (2015).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
65.Wu, T. D. & Watanabe, C. K. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics 21, 1859–1875 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
66.Benson, G. Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acids Res. 27, 573–580 (1999).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
67.Bao, W., Kojima, K. K. & Kohany, O. Repbase Update, a database of repetitive elements in eukaryotic genomes. Mob. DNA 6, 11 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
68.Xu, Z. & Wang, H. LTR_FINDER: an efficient tool for the prediction of full-length LTR retrotransposons. Nucleic Acids Res. 35, W265–W268 (2007).PubMed 
PubMed Central 
Article 

Google Scholar 
69.Edgar, R. C. & Myers, E. W. PILER: identification and classification of genomic repeats. Bioinformatics 21, i152–i158 (2005).CAS 
PubMed 
Article 

Google Scholar 
70.Price, A. L., Jones, N. C. & Pevzner, P. A. De novo identification of repeat families in large genomes. Bioinformatics 21, i351–i358 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
71.Lowe, T. M. & Eddy, S. R. tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 25, 955–964 (1997).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
72.Nawrocki, E. P. & Eddy, S. R. Infernal 1.1: 100-fold faster RNA homology searches. Bioinformatics 29, 2933–2935 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
73.Nawrocki, E. P. et al. Rfam 12.0: updates to the RNA families database. Nucleic Acids Res. 43, D130–D137 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
74.Holt, C. & Yandell, M. MAKER2: an annotation pipeline and genome-database management tool for second-generation genome projects. BMC Bioinf. 12, 491 (2011).Article 

Google Scholar 
75.Kim, D., Langmead, B. & Salzberg, S. L. HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
76.Pertea, M. et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
77.Johnson, A. D. et al. SNAP: a web-based tool for identification and annotation of proxy SNPs using HapMap. Bioinformatics 24, 2938–2939 (2008).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
78.Stanke, M. & Morgenstern, B. AUGUSTUS: a web server for gene prediction in eukaryotes that allows user-defined constraints. Nucleic Acids Res. 33, W465–W467 (2005).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
79.Buels, R. et al. JBrowse: a dynamic web platform for genome visualization and analysis. Genome Biol. 17, 66 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
80.Buchfink, B., Xie, C. & Huson, D. H. Fast and sensitive protein alignment using DIAMOND. Nat. Methods 12, 59–60 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
81.Bairoch, A. & Apweiler, R. The SWISS-PROT protein sequence database and its supplement TrEMBL in 2000. Nucleic Acids Res. 28, 45–48 (2000).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
82.Kanehisa, M., Sato, Y., Kawashima, M., Furumichi, M. & Tanabe, M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 44, D457–D462 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
83.Jones, P. et al. InterProScan 5: genome-scale protein function classification. Bioinformatics 30, 1236–1240 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
84.Conesa, A. et al. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
85.Li, L., Stoeckert, C. J. Jr & Roos, D. S. OrthoMCL: identification of ortholog groups for eukaryotic genomes. Genome Res. 13, 2178–2189 (2003).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
86.Li, H. et al. TreeFam: a curated database of phylogenetic trees of animal gene families. Nucleic Acids Res. 34, D572–D580 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
87.Edgar, R. C. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
88.Capella-Gutierrez, S., Silla-Martinez, J. M. & Gabaldon, T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25, 1972–1973 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
89.Guindon, S., Delsuc, F., Dufayard, J. F. & Gascuel, O. Estimating maximum likelihood phylogenies with PhyML. Methods Mol. Biol. 537, 113–137 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
90.Yang, Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, 1586–1591 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
91.De Bie, T., Cristianini, N., Demuth, J. P. & Hahn, M. W. CAFE: a computational tool for the study of gene family evolution. Bioinformatics 22, 1269–1271 (2006).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
92.Tang, H. et al. Unraveling ancient hexaploidy through multiply-aligned angiosperm gene maps. Genome Res. 18, 1944–1954 (2008).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
93.Schwartz, S. et al. Human–mouse alignments with BLASTZ. Genome Res. 13, 103–107 (2003).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
94.Bailey, J. A., Yavor, A. M., Massa, H. F., Trask, B. J. & Eichler, E. E. Segmental duplications: organization and impact within the current human genome project assembly. Genome Res. 11, 1005–1017 (2001).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
95.Birney, E., Clamp, M. & Durbin, R. GeneWise and Genomewise. Genome Res. 14, 988–995 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
96.Ruan, J., Li, H., Chen, Z. & Coghlan, A. TreeFam: 2008 update. Nucleic Acids Res. 36, D735–D740 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
97.Tholl, D. et al. Practical approaches to plant volatile analysis. Plant J. 45, 540–560 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
98.Wen, B., Mei, Z., Zeng, C. & Liu, S. metaX: a flexible and comprehensive software for processing metabolomics data. BMC Bioinform. 18, 183 (2017).Article 
CAS 

Google Scholar 
99.Wen, P. et al. The sex pheromone of a globally invasive honey bee predator, the Asian eusocial hornet, Vespa velutina. Sci. Rep. 7, 12956 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
100.Chen, Y. et al. SOAPnuke: a MapReduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data. Gigascience 7, 1–6 (2018).PubMed 
PubMed Central 

Google Scholar 
101.Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
102.Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from RNA-seq data with or without a reference genome. BMC Bioinform. 12, 323 (2011).CAS 
Article 

Google Scholar 
103.Tian, X., Chen, L., Wang, J., Qiao, J. & Zhang, W. Quantitative proteomics reveals dynamic responses of Synechocystis sp. PCC 6803 to next-generation biofuel butanol. J. Proteom. 78, 326–345 (2013).CAS 
Article 

Google Scholar 
104.Wen, B. et al. IQuant: an automated pipeline for quantitative proteomics based upon isobaric tags. Proteomics 14, 2280–2285 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
105.Brosch, M., Yu, L., Hubbard, T. & Choudhary, J. Accurate and sensitive peptide identification with Mascot Percolator. J. Proteome Res. 8, 3176–3181 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
106.Savitski, M. M., Wilhelm, M., Hahne, H., Kuster, B. & Bantscheff, M. A scalable approach for protein false discovery rate estimation in large proteomic data sets. Mol. Cell Proteomics 14, 2394–2404 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
107.Cox, J. & Mann, M. MaxQuant enables high peptide identification rates, individualized p.p.b.-range mass accuracies and proteome-wide protein quantification. Nat. Biotechnol. 26, 1367–1372 (2008).CAS 
Article 

Google Scholar 
108.Bruderer, R. et al. Extending the limits of quantitative proteome profiling with data-independent acquisition and application to acetaminophen-treated three-dimensional liver microtissues. Mol. Cell Proteomics 14, 1400–1410 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
109.Dunn, W. B. et al. Procedures for large-scale metabolic profiling of serum and plasma using gas chromatography and liquid chromatography coupled to mass spectrometry. Nat. Protoc. 6, 1060–1083 (2011).CAS 
PubMed 
Article 

Google Scholar 
110.Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
111.Choi, M. et al. MSstats: an R package for statistical analysis of quantitative mass spectrometry-based proteomic experiments. Bioinformatics 30, 2524–2526 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
112.Wen, X.-L., Wen, P., Dahlsjö, C. A. L., Sillam-Dussès, D. & Šobotník, J. Breaking the cipher: ant eavesdropping on the variational trail pheromone of its termite prey. Proc. R. Soc. B 284, 20170121 (2017).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
113.Bailey, T. L. & Elkan, C. Fitting a mixture model by expectation maximization to discover motifs in biopolymers. Proc. Int. Conf. Intell. Syst. Mol. Biol. 2, 28–36 (1994).CAS 
PubMed 
PubMed Central 

Google Scholar 
114.Lescot, M. et al. PlantCARE, a database of plant cis-acting regulatory elements and a portal to tools for in silico analysis of promoter sequences. Nucleic Acids Res. 30, 325–327 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
115.Langfelder, P. & Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 9, 559 (2008).Article 
CAS 

Google Scholar 
116.Yip, A. M. & Horvath, S. Gene network interconnectedness and the generalized topological overlap measure. BMC Bioinform. 8, 22 (2007).Article 
CAS 

Google Scholar 
117.Langfelder, P., Zhang, B. & Horvath, S. Defining clusters from a hierarchical cluster tree: the Dynamic Tree Cut package for R. Bioinformatics 24, 719–720 (2008).CAS 
PubMed 
Article 

Google Scholar 
118.Gendrel, A. V., Lippman, Z., Martienssen, R. & Colot, V. Profiling histone modification patterns in plants using genomic tiling microarrays. Nat. Methods 2, 213–218 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar  More

1.Partecke, J., Gwinner, E. & Bensch, S. Is urbanisation of European blackbirds (Turdus merula) associated with genetic differentiation?. J. Ornithol. 147, 549–552 (2006).Article 

Google Scholar 
2.Perrier, C. et al. Great tits and the city: Distribution of genomic diversity and gene–environment associations along an urbanization gradient. Evol. Appl. 11, 593–613 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
3.Chace, J. F. & Walsh, J. J. Urban effects on native avifauna: A review. Landsc. Urban Plan. 74, 46–69 (2006).Article 

Google Scholar 
4.Evans, K. L. et al. Independent colonization of multiple urban centres by a formerly forest specialist bird species. Proc. R. Soc. B 276(1666), 2403–2410 (2009).PubMed 
Article 
PubMed Central 

Google Scholar 
5.Björklund, M., Ruiz, I. & Senar, J. C. Genetic differentiation in the urban habitat: The great tits (Parus major) of the parks of Barcelona city. Biol. J. Linn. Soc. 99, 9–19 (2010).Article 

Google Scholar 
6.Crooks, K. R. & Sanjayan, M. A. Connectivity conservation: Maintaining connections for nature. In Connectivity Conservation (eds Crooks, K. R. & Sanjayan, M. A.) 1–20 (Cambridge University Press, Cambridge, 2006).
Google Scholar 
7.Evans, K. L., Chamberlain, D. E., Hatchwell, B. J., Gregory, R. D. & Gaston, K. J. What makes an urban bird?. Glob. Change Biol. 17, 32–44 (2011).ADS 
Article 

Google Scholar 
8.Seress, G. & Liker, A. Habitat urbanization and its effects on birds. Acta Zool. Acad. Sci. Hung. 61(4), 373–408 (2015).Article 

Google Scholar 
9.Miles, L. S., Rivkin, L. R., Johnson, M. T. J., Munshi-South, J. & Verrelli, B. C. Gene flow and genetic drift in urban environments. Mol. Ecol. 28, 4138–4151 (2019).PubMed 
Article 

Google Scholar 
10.Shochat, E., Warren, P. S., Faeth, S. H., McIntyre, N. E. & Hope, D. From patterns to emerging processes in mechanistic urban ecology. Trends Ecol. Evol. 21, 186–191 (2006).PubMed 
Article 

Google Scholar 
11.Chamberlain, D. E. et al. Avian productivity in urban landscapes: A review and meta-analysis. Ibis 151, 1–18 (2009).Article 

Google Scholar 
12.Delaney, K. S., Riley, S. P. D. & Fisher, R. N. A rapid, strong, and convergent genetic response to urban habitat fragmentation in four divergent and widespread vertebrates. PLoS ONE 5(9), e12767 (2010).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
13.Unfried, T. M., Hauser, L. & Marzluff, J. M. Effects of urbanization on Song Sparrow (Melospiza melodia) population connectivity. Conserv. Genet. 14(1), 41–53 (2013).Article 

Google Scholar 
14.Cureton, J. C. et al. Effects of urbanization on genetic diversity, gene flow, and population structure in the ornate box turtle (Terrapene ornato). Amphib-Reptil. 35, 87–97 (2014).Article 

Google Scholar 
15.Indykiewicz, P., Podlaszczuk, P., Janiszewska, A. & Minias, P. Extensive gene flow along the urban-rural gradient in a migratory colonial bird. J. Avian Biol. 49(6), e01723 (2018).Article 

Google Scholar 
16.Hurtado, G. & Mabry, K. E. Genetic structure of an abundant small mammal is influenced by low intensity urbanization. Conserv. Genet. 20, 705–715 (2019).CAS 
Article 

Google Scholar 
17.Khimoun, A. et al. Urbanization without isolation: The absence of genetic structure among cities and forests in the tiny acorn ant Temnothorax nylanderi. Biol. Lett. 16, 20190741 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Munshi-South, J., Zolnik, C. P. & Harris, S. E. Population genomics of the Anthropocene: Urbanization is negatively associated with genome-wide variation in white-footed mouse populations. Evol. Appl. 9, 546–564 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
19.Brewer, V. N., Lane, S. J., Sewall, K. B. & Mabry, K. E. Effects of low-density urbanization on genetic structure in the Song Sparrow. PLoS ONE 15(6), e0234008 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
20.Slatkin, M. Gene flow in natural populations. Annu. Rev. Ecol. Syst. 16, 393–430 (1985).Article 

Google Scholar 
21.Balloux, F. & Lugon-Moulin, N. The estimation of population differentiation with microsatellite markers. Mol. Ecol. 11, 155–165 (2002).PubMed 
Article 

Google Scholar 
22.Vangestel, C., Mergeay, J., Dawson, D. A., Vandomme, V. & Lens, L. Spatial heterogeneity in genetic relatedness among house sparrows along an urban—rural gradient as revealed by individual-based analysis. Mol. Ecol. 20, 4643–4653 (2011).PubMed 
Article 

Google Scholar 
23.Barnett, J. R., Ruiz-Gutierrez, V., Coulon, A. & Lovette, I. J. Weak genetic structuring indicates ongoing gene flow across White-ruffed Manakin (Corapipo altera) populations in a highly fragmented Costa Rica landscape. Conserv. Genet. 9, 1403–1412 (2008).Article 

Google Scholar 
24.Riegert, J., Fainová, D. & Bystrická, D. Genetic variability, body characteristics and reproductive parameters of neighbouring rural and urban common kestrel (Falco tinnuculus) populations. Popul. Ecol. 52, 73–79 (2009).Article 

Google Scholar 
25.MacDougall-Shackleton, E. A., Clinchy, M., Zanette, L. & Neff, B. D. Songbird genetic diversity is lower in anthropogenically versus naturally fragmented landscapes. Conserv. Genet. 12, 1195–1203 (2011).Article 

Google Scholar 
26.Caizergues, A. E. et al. Testing for parallel genomic and epigenomic footprints of adaptation to urban life in a passerine bird. bioRxiv. https://doi.org/10.1101/2021.02.10.43045227.Schmidt, C., Domaratzki, M., Kinnunen, R. P., Bowman, J. & Garroway, C. J. Continent-wide effects of urbanization on bird and mammal genetic diversity. Proc. R. Soc. B. 287, 20192497 (2020).CAS 
PubMed 
Article 

Google Scholar 
28.Cramp, S. & Perrins, C. M. The Birds of the Western Palearctic Vol. 7 (Oxford University Press, 1993).
Google Scholar 
29.Dauwe, T. et al. Great and Blue tit feathers as biomonitors for heavy metal pollution. Ecol. Indic. 1, 227–234 (2002).CAS 
Article 

Google Scholar 
30.Bańbura, J. & Bańbura, M. Blue tits Cyanistes caeruleus and great tits Parus major as urban habitat breeders. Inter Stud. Sparrows 36, 66–72 (2012).Article 

Google Scholar 
31.Charmantier, A., Doutrelant, C., Dubuc-Messier, G., Fargevieille, A. & Szulkin, M. Mediterranean blue tits as a case study of local adaptation. Evol. Appl. 9, 135–152 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
32.Lemoine, M. et al. Low but contrasting neutral genetic differentiation shaped by winter temperature in European Great Tits. Biol. J. Linn. Soc. 118, 668–685 (2016).Article 

Google Scholar 
33.Porlier, M. Garant, D. Perret, P. and Charmantier, A. Habitat-linked population genetic differentiation in the Blue tit Cyanistes caeruleus. J. Hered. 103, 781–791 (2012).34.Szulkin, M., Gagnaire, P. A., Bierne, N. & Charmantier, A. Population genomic footprints of fine-scale differentiation between habitats in Mediterranean blue tits. Mol. Ecol. 25, 542–558 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
35.Dubuc-Messier, G. et al. Gene flow does not prevent personality and morphological differentiation between two blue tit populations. J. Evol. Biol. 31, 1127–1137 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
36.Postma, E. D., Tex, R.-J., Noordwijk, A. J. & Mateman, A. C. Neutral markers mirror small-scale quantitative genetic differentiation in an avian island population. Biol. J. Linn. Soc. 97, 867–875 (2009).Article 

Google Scholar 
37.Salmón, P. et al. Repeated genomic signature of adaptation to urbanisation in a songbird across Europe. bioRxiv. https://doi.org/10.1101/2020.05.05.078568 (2020).38.Dhondt, A. A. Effects of competition on great and blue tit reproduction: Intensity and importance in relation to habitat quality. J. Anim. Ecol. 79, 257–265 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
39.Nilsson, A. L. K., Lindström, Å., Jonzén, N., Nilsson, S. G. & Karlsson, L. The effect of climate change on partial migration: The blue tit paradox. Glob. Change Biol. 12, 2014–2022. https://doi.org/10.1111/j.1365-2486.2006.01237.x (2006).ADS 
Article 

Google Scholar 
40.Nilsson, A. L. K., Alerstam, T. & Nilsson, J. Å. Diffuse, short and slow migration among Blue Tits. J. Ornithol. 149, 365–373. https://doi.org/10.1007/s10336-008-0280-3 (2008).Article 

Google Scholar 
41.Bańbura, J. et al. Spatial and temporal variation in heterophil-to-lymphocyte ratios of nestling passerine birds: Comparison of blue tits and great tits. PLoS ONE 8(9), e74226 (2013).42.Adamou, A.-E., Bańbura, M. & Bańbura, J. Subtle differences in breeding performance between Great Tits Parus major and Afrocanarian Blue Tits Cyanistes teneriffae in the peripheral zone of the species geographic ranges in NE Algeria. Eur. Zool. J. 87, 263–271 (2020).Article 

Google Scholar 
43.Dhondt, A. A. & Eyckerman, R. Competition between the great tit and the blue tit outside the breeding season in field experiments. Ecology 61, 1291–1296 (1980).Article 

Google Scholar 
44.Ortego, J., Garcia-Navas, V., Ferrer, E. S. & Sanz, J. J. Genetic structure reflects natal dispersal movements at different spatial scales in the blue tit Cyanistes caeruleus. Anim. Behav. 82, 131–137 (2011).Article 

Google Scholar 
45.Langin, K. M. et al. Characterizing range-wide divergence in an alpine-endemic bird: A comparison of genetic and genomic approaches. Conserv. Genet. 19(6), 1471–1485 (2018).CAS 
Article 

Google Scholar 
46.Roques, S., Chancerel, E., Boury, C., Pierre, M. & Acolas, M. L. From microsatellites to single nucleotide polymorphisms for the genetic monitoring of a critically endangered sturgeon. Ecol. Evol. 9(12), 7017–7029 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
47.Zimmerman, S. J., Aldridge, C. L. & Oyler-McCance, S. J. An empirical comparison of population genetic analyses using microsatellite and SNP data for a species of conservation concern. BMC Genom. 21, 1–16 (2020).Article 
CAS 

Google Scholar 
48.Markowski, M. et al. Effects of experimental lead exposure on physiological indices of nestling great tits Parus major: Haematocrit and heterophile-to-lymphocyte ratio. Conserv. Physiol. 7, coz067 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Bańbura, J. et al. Habitat and year-to-year variation in haemoglobin concentration in nestling blue tits Cyanistes caeruleus. Comp. Biochem. Phys. A 148, 572–577 (2007).Article 
CAS 

Google Scholar 
50.Kiedrzyński, M. The impact of forest management on the flora and vegetation of old oak-stands (an example from The Spała Forests, central Poland). Nat. Conserv. 65, 51–62 (2008).
Google Scholar 
51.Glądalski, M. et al. Effects of human-related disturbance on breeding success of urban and non-urban blue tits (Cyanistes caeruleus). Urban Ecosyst. 19, 1325–1334 (2016).Article 

Google Scholar 
52.Markowski, M. et al. Spatial and temporal variation of lead, cadmium, and zinc in feathers of great tit and blue tit nestlings in Central Poland. Arch. Environ. Contam. Toxicol. 67, 507–518 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Richard, M. & Thorpe, R. S. Highly polymorphic microsatellites in the lacertid Gallotia Gallowi from the western Canary Islands. Mol. Ecol. 9, 1919–1952 (2000).CAS 
PubMed 
Article 

Google Scholar 
54.Saladin, V., Bonfils, D., Binz, T. & Richner, H. Isolation and characterization of 16 microsatellite loci in the Great Tit Parus major. Mol. Ecol. Notes 3, 520–522 (2003).CAS 
Article 

Google Scholar 
55.Dawson, D. A., Hanotte, O., Greig, C., Stewart, I. R. K. & Burke, T. Polymorphic microsatellites in the blue tit Parus caeruleus and their cross-species utility in 20 songbird families. Mol. Ecol. 9, 1941–1944 (2000).CAS 
PubMed 
Article 

Google Scholar 
56.van Oosterhout, C., Hutchinson, W. F., Wills, D. P. & Shipley, P. MICRO-CHECKER: Software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Res. 4, 535–538 (2004).
Google Scholar 
57.Guo, S. W. & Thompson, E. A. Performing the exact test of Hardy-Weinberg proportion for multiple alleles. Biometrics 48, 361–372 (1992).CAS 
MATH 
Article 

Google Scholar 
58.Excoffier, L. & Lischer, H. E. L. Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Res. 10, 564–567 (2010).Article 

Google Scholar 
59.Goudet, J. FSTAT (version 12): A computer program to calculate F-statistics. J. Hered. 86, 485–486 (1995).Article 

Google Scholar 
60.Kalinowski, S. T., Taper, M. L. & Marshall, T. C. Revising how the computer program CERVUS accommodates genotyping error increases success in paternity assignment. Mol. Ecol. 16, 1099–1106 (2007).PubMed 
Article 

Google Scholar 
61.Peakall, R. & Smouse, P. E. GENALEX 6: Genetic analysis in Excel: Population genetic software for teaching and research. Mol. Ecol. Notes 6, 288–295 (2006).Article 

Google Scholar 
62.Peakall, R. & Smouse, P. E. GENALEX 6.5: Genetic analysis in Excel. Population genetic software for teaching and research: An update. Bioinformatics 28, 2537–2539 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
63.Slatkin, M. A measure of subdivision based on microsatellite allele frequencies. Genetics 139, 457–462 (1995).CAS 
PubMed 
PubMed Central 

Google Scholar 
64.Hardy, O. J., Charbonnel, N., Fréville, H. & Heuertz, M. Microsatellite allele sizes: A simple test to assess their significance on genetic differentiation. Genetics 163, 1467–1482 (2003).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
65.Hardy, O. J. & Vekemans, X. SPAGeDi: A versatile computer program to analyse spatial genetic structure at the individual or population levels. Mol. Ecol. Notes 2, 618–620 (2002).Article 
CAS 

Google Scholar 
66.Hedrick, P. W. A standardized genetic differentiation measure. Evolution 59, 1633–1638 (2005).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
67.Meirmans, P. G. Using the AMOVA framework to estimate a standardized genetic differentiation measure. Evolution 60, 2399–2402 (2006).Article 

Google Scholar 
68.Nei, M. Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics 89, 583–590 (1978).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
69.Dray, S. & Dufour, A. The ade4 Package: Implementing the duality diagram for ecologists. J. Stat. Softw. 22(4), 1–20. https://doi.org/10.18637/jss.v022.i04 (2007).Article 

Google Scholar 
70.Jombart, T. adegenet: A R package for the multivariate analysis of genetic markers. Bioinformatics 24(11), 1403–1405 (2008).CAS 
Article 

Google Scholar 
71.TIBCO Software Inc. Statistica (Data Analysis Software System), Version 13. http://statistica.io. (2017). More

1.Bontemps, J. D. & Bouriaud, O. Predictive approaches to forest site productivity: Recent trends, challenges and future perspectives. Forestry 87, 109–128 (2014).Article 

Google Scholar 
2.Skovsgaard, J. P. & Vanclay, J. K. Forest site productivity: a review of the evolution of dendrometric concepts for even-aged stands. Forestry 81, 13–31 (2008).Article 

Google Scholar 
3.Albert, M. & Schmidt, M. Climate-sensitive modelling of site-productivity relationships for Norway spruce (Picea abies (L.) Karst) and common beech (Fagus sylvatica L.). For. Ecol. Manag. 259, 739–749 (2010).Article 

Google Scholar 
4.Véga, C. & St-Onge, B. Mapping site index and age by linking a time series of canopy height models with growth curves. For. Ecol. Manag. 257, 951–959 (2009).Article 

Google Scholar 
5.Hägglund, B. & Lundmark, J. E. Site index estimation by means of site properties of Scots pine and Norway spruce in Sweden. Stud. For. Suec. 138, 5–38 (1977).
Google Scholar 
6.Johansson, T. Site index curves for common alder and grey alder growing on different types of forest soil in Sweden. Scand. J. For. Res. 14, 441–453 (1999).Article 

Google Scholar 
7.Raulier, F., Lambert, M.-C., Pothier, D. & Ung, C.-H. Impact of dominant tree dynamics on site index curves. For. Ecol. Manag. 184, 65–78 (2003).Article 

Google Scholar 
8.Corral Rivas, J. J., Álvarez González, J. G., Ruíz González, A. D. & von Gadow, K. Compatible height and site index models for five pine species in El Salto, Durango (Mexico). For. Ecol. Manag. 201, 145–160 (2004).Article 

Google Scholar 
9.Tewari, V. P., Rivas, J. J. C., VilČko, F. & Von Gadow, K. Height-Growth and site index equations for social forestry plantations of acacia nilotica and eucalyptus hybrid in gujarat state of India. For. Trees Livelihoods 17, 125–140 (2007).Article 

Google Scholar 
10.Monserud, R. A. & Rehfeldt, G. E. Genetic and environmental components of variation of site index in inland douglas-fir. For. Sci. 36, 1–9 (1990).
Google Scholar 
11.Alvarez-González, J. G., Ruiz-González, A. D., Rodríguez-Soalleiro, R. & Barrio-Anta, M. Ecorregional site index models for Pinus pinaster in Galicia (northwestern Spain). Ann. For. Sci. 62, 115–127 (2005).Article 

Google Scholar 
12.Bravo-Oviedo, A., Tomé, M., Bravo, F., Montero, G. & del Río, M. Dominant height growth equations including site attributes in the generalised algebraic difference approach. Can. J. For. Res. 38, 2348–2358 (2008).Article 

Google Scholar 
13.Johansson, T. Site index curves for Norway spruce plantations on farmland with different soil types. Stud. For. Suec. 198, 1–19 (1995).
Google Scholar 
14.Adams, J. P., Matney, T. G., Land, S. B. Jr., Belli, K. L. & Duzan, H. W. Jr. Incorporating genetic parameters into a loblolly pine growth-and-yield model. Can. J. For. Res. 36, 1959–1967 (2006).Article 

Google Scholar 
15.Buford, M. A. & Burkhart, H. E. Genetic improvement effects on growth and yield of loblolly pine plantations. For. Sci. 33, 707–724 (1987).
Google Scholar 
16.Monserud, R. A. Height growth and site index curves for inland Douglas-fir based on stem analysis data and forest habitat type. For. Sci. 30, 943–965 (1984).
Google Scholar 
17.García, O. Dynamical implications of the variability representation in site-index modelling. Eur. J. For. Res. 130, 671–675 (2010).Article 

Google Scholar 
18.Calama, R., Cañadas, N. & Montero, G. Inter-regional variability in site index models for even-aged stands of stone pine (Pinus pinea L.) in Spain. Ann. For. Sci. 60, 259–269 (2003).Article 

Google Scholar 
19.Adame, P., Hynynen, J., Cañellas, I. & del Río, M. Individual-tree diameter growth model for rebollo oak (Quercus pyrenaica Willd.) coppices. For. Ecol. Manag. 255, 1011–1022 (2008).Article 

Google Scholar 
20.Bravo-Oviedo, A., del Río, M. & Montero, G. Geographic variation and parameter assessment in generalised algebraic difference site index modelling. For. Ecol. Manag. 247, 107–119 (2007).Article 

Google Scholar 
21.Bontemps, J. & Bouriaud, O. Predictive approaches to forest site productivity: recent trends, challenges and future perspectives. Forestry 87, 1–20 (2013).
Google Scholar 
22.Claessens, H., Pauwels, D., Thibaut, A. & Rondeux, J. Site index curves and autecology of ash, sycamore and cherry in Wallonia (Southern Belgium). Forestry 72, 171–182 (1999).Article 

Google Scholar 
23.Karlsson, K. Height growth patterns of Scots pine and Norway spruce in the coastal areas of western Finland. For. Ecol. Manag. 135, 205–216 (2000).Article 

Google Scholar 
24.Martín-Benito, D., Gea-Izquierdo, G., del Río, M. & Cañellas, I. Long-term trends in dominant-height growth of black pine using dynamic models. For. Ecol. Manag. 256, 1230–1238 (2008).Article 

Google Scholar 
25.Dyrekcja Generalna Lasów Państwowych, Raport o stanie lasów w Polsce (2018).26.Socha, J. Long-term effect of wetland drainage on the productivity of Scots pine stands in Poland. For. Ecol. Manag. 274, 172–180 (2012).Article 

Google Scholar 
27.Bravo, F. & Montero, G. Site index estimation in Scots pine (Pinus sylvestris L.) stands in the High Ebro Basin (northern Spain) using soil attributes. Forestry 74, 395–406 (2001).Article 

Google Scholar 
28.Gadow, K.; Hui, G. Modelling Forest Development, Vol. 57, Forestry Sciences (Springer Netherlands, Dordrecht, 1999). ISBN 978–1–4020–0276–2.29.Schwappach, A. Ertragstafeln der Wichtigeren Holzarten (Druckerei Merkur, Prag, 1943).
Google Scholar 
30.Szymkiewicz, B. Niektóre zagadnienia dotyczące tablic zasobności drzewostanów sosnowych. Pr. IBL Seria A 67 (1948)31.Szymkiewicz, B. Tablice Zasobności i Przyrostu Drzewostanów (Państwowe Wydawnictwo Rolnicze i Leśne, Warszawa, 2001). ISBN 83–09–01745–6.32.Pretzsch, H., Biber, P., Schütze, G., Uhl, E. & Rötzer, T. Forest stand growth dynamics in Central Europe have accelerated since 1870. Nat. Commun. 5, 4967 (2014).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
33.Socha, J. & Orzeł, S. Dynamic site index curves for Scots pine (Pinus sylvestris L.) in southern Poland. Sylwan 157, 26–38 (2013).
Google Scholar 
34.Socha, J., Tymińska-Czabańska, L., Grabska, E. & Orzeł, S. Site index models for main forest-forming tree species in Poland. Forests 11, 8–10 (2020).
Google Scholar 
35.Esri Inc. ArcGIS Pro (Version 2.2.0). Esri Inc. https://www.esri.com/en-us/arcgis/products/arcgis-pro/overview (2020).36.Caudullo, G., Welk, E. & San-Miguel-Ayanz, J. Chorological maps for the main European woody species. Data Br. 12, 662–666 (2017).Article 

Google Scholar 
37.Zielony, R., Kliczkowska, A. Regionalizacja Przyrodniczo-Leśna Polski 2010 (2012). ISBN 9788361633624.38.Carmean, W. H. Site index curves for upland oaks in the central states. For. Sci. 18, 109–120 (1972).
Google Scholar 
39.Mehtätalo, L. & Lappi, J. Biometry for Forestry and Environmental Data: with Examples in R 1st edn. (Chapman and Hall/CRC, Boca Raton, 2020).Book 

Google Scholar 
40.Nigh, G. Engelmann spruce site index models: a comparison of model functions and parameterisations. PLoS ONE 10, e0124079 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
41.Pinheiro J., Bates D., & DebRoy S. S. D. nlme: Linear and Nonlinear Mixed Effects Models (2020).42.Wang, M., Borders, B. E. & Zhao, D. An empirical comparison of two subject-specific approaches to dominant heights modeling: the dummy variable method and the mixed model method. For. Ecol. Manag. 255, 2659–2669 (2008).Article 

Google Scholar 
43.Bronisz, K. & Mehtätalo, L. Mixed-effects generalised height–diameter model for young silver birch stands on post-agricultural lands. For. Ecol. Manage. 460, 117901 (2020).Article 

Google Scholar 
44.Wang, M., Bhatti, J., Wang, Y. & Varem-Sanders, T. Examining the gain in model prediction accuracy using serial autocorrelation for dominant height prediction. For. Sci. 57, 241–251 (2011).
Google Scholar 
45.R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/ (2020).46.Bailey, R. L. & Clutter, J. L. Base-age invariant polymorphic site curves. For. Sci. 20, 155–159 (1974).
Google Scholar 
47.Cieszewski, J. & Bailey, L. Generalized algebraic difference approach: theory based derivation of dynamic site equations with polymorphism and variable asymptotes. For. Sci. 46, 116–126 (2000).
Google Scholar 
48.Cieszewski, C. J. Three methods of deriving advanced dynamic site equations demonstrated on inland Douglas-fir site curves. Can. J. For. Res. For. 31, 165–173 (2001).Article 

Google Scholar 
49.Krumland, B. & Eng, H. Site index systems for major young-growth forest and Woodland species in North California. Calif. For. 4, 1–220 (2005).
Google Scholar 
50.Cieszewski, C. J. Developing a well-behaved dynamic site equation using a modified hossfeld IV function Y 3 = (axm)/(c + x m–1), a simplified mixed-model and scant subalpine fir data. For. Sci. 49, 539–554 (2003).
Google Scholar 
51.Sharma, R. P., Brunner, A., Eid, T. & Øyen, B.-H. Modelling dominant height growth from national forest inventory individual tree data with short time series and large age errors. For. Ecol. Manag. 262, 2162–2175 (2011).Article 

Google Scholar 
52.Anta, M. B. et al. Development of a basal area growth system for maritime pine in northwestern Spain using the generalised algebraic difference approach. Can. J. For. Res. 36, 1461–1474 (2006).Article 

Google Scholar 
53.Cieszewski, C. J., Harrison, M. & Martin, S. W. Examples of practical methods for unbiased parameter estimation in self-referencing functions. In Proceedings of the First International Conference on Measurements and Quantitative Methods and Management and The 1999 Southern Mensur (ed. Cieszewski, C. J.) (2000).54.Nunes, L., Patrício, M., Tomé, J. & Tomé, M. Modeling dominant height growth of maritime pine in Portugal using GADA methodology with parameters depending on soil and climate variables. Ann. For. Sci. 68, 311–323 (2011).Article 

Google Scholar 
55.Cieszewski, C. J. Comparing properties of self-referencing models based on nonlinear-fixed-effects versus nonlinear-mixed-effects modeling approaches. Math. Comput. For. Nat. Sci. 10, 46–57 (2018).
Google Scholar 
56.Seynave, I. et al. Picea abies site index prediction by environmental factors and understorey vegetation: a two-scale approach based on survey databases. Can. J. For. Res. 35, 10 (2005).Article 

Google Scholar 
57.Brandl, S. et al. Static site indices from different national forest inventories: harmonisation and prediction from site conditions. Ann. For. Sci. 75, 1–17 (2018).Article 

Google Scholar 
58.Breda, N., Huc, R. & André Granier, E. D. Temperate forest trees and stands under severe drought: a review of ecophysiological responses, adaptation processes and long-term consequences. Ann. For. Sci. 63, 625–644 (2006).Article 

Google Scholar 
59.Loik, M. E., Breshears, D. D., Lauenroth, W. K. & Belnap, J. A multi-scale perspective of water pulses in dryland ecosystems: climatology and ecohydrology of the western USA. Oecologia 141, 269–281 (2004).ADS 
PubMed 
Article 

Google Scholar 
60.Kundzewicz, Z. W., Hov, Ø., Okruszko, T. Zmiany klimatu i ich wpływ na wybrane sektory w Polsce. 257 (2017).61.Barrio Anta, M. & Dieguez-Aranda, U. Site quality of pedunculate oak (Quercus robur L.) stands in Galicia (northwest Spain). Eur. J. For. Res. 124, 19–28 (2005).Article 

Google Scholar 
62.Bowman, D. M. J. S., Brienen, R. J. W., Gloor, E., Phillips, O. L. & Prior, L. D. Detecting trends in tree growth: not so simple. Trends Plant Sci. 18, 11–17 (2013).CAS 
PubMed 
Article 

Google Scholar 
63.Isaev, A., Korovin, G., Zamolodkchikov, D., Utkin, A. & Pryaznikov, A. Carbon stock and depostion in phytomass of the Russian forests 247–256 (Kluwer, Amsterdam, 1995).
Google Scholar 
64.Fang, J. Y. & Wang, Z. M. Forest biomass estimation at regional and global levels, with special reference to China’s forest biomass. Ecol. Res. 16, 587–592 (2001).CAS 
Article 

Google Scholar 
65.Schroeder, P., Brown, S., Mo, J., Birdsey, R. & Cieszewski, C. Biomass estimation for temperate broadleaf forests of the United States using inventory data. For. Sci. 43, 424–434 (1997).
Google Scholar 
66.Bravo-Oviedo, A., Gallardo-Andrés, C., del Río, M. & Montero, G. Regional changes of Pinus pinaster site index in Spain using a climate-based dominant height model. Can. J. For. Res. 40, 2036–2048 (2010).Article 

Google Scholar 
67.Woodbury, P. B., Smith, J. E., Weinstein, D. & Laurence, J. Assessing potential climate change effects on loblolly pine growth: a probabilistic regional modeling approach. For. Ecol. Manag. 107, 99–116 (1998).Article 

Google Scholar 
68.Latta, G., Temesgen, H. & Barrett, T. M. B. M. Mapping and imputing potential productivity of Pacific Northwest forests using climate variables. Can. J. For. Res. 39, 1197–1207 (2009).69.Coops, N. C., Hember, R. & Waring, R. H. Assessing the impact of current and projected climates on Douglas-Fir productivity in British Columbia, Canada, using a process-based model (3-PG). Can. J. For. Res. 40, 511–524 (2010).Article 

Google Scholar 
70.Coops, N. C., Coggins, S. B. & Kurz, W. a Mapping the environmental limitations to growth of coastal Douglas-fir stands on Vancouver Island, British Columbia. Tree Physiol. 27, 805–815 (2007).PubMed 
Article 
PubMed Central 

Google Scholar  More

1.Latrubesse, E. M. et al. Damming the rivers of the Amazon basin. Nature 546, 363–369 (2017).ADS 
CAS 
PubMed 
Article 

Google Scholar 
2.Encalada, A. C. et al. A global perspective on tropical montane rivers. Science 365, 1124–1129 (2019).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
3.Winemiller, K. O. et al. Development and environment. Balancing hydropower and biodiversity in the Amazon, Congo, and Mekong. Science 351, 128–129 (2016).ADS 
CAS 
PubMed 
Article 

Google Scholar 
4.McIntyre, P. B., Reidy Liermann, C. A. & Revenga, C. Linking freshwater fishery management to global food security and biodiversity conservation. Proc. Natl. Acad. Sci. U. S. A. 113, 12880–12885 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.IPBES. Global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (eds Brondizio, E. S., Settele, J., Díaz, S. & Ngo, H. T.) (Bonn, Germany, 2019).6.Allen, D. J., Smith, K. G. & Darwall, W. R. T. The Status and Distribution of Freshwater Biodiversity in Indo-Burma (IUCN, 2012).
Google Scholar 
7.Barlow, J. et al. The future of hyperdiverse tropical ecosystems. Nature 559, 517–526 (2018).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Dudgeon, D. Multiple threats imperil freshwater biodiversity in the Anthropocene. Curr. Biol. 29, 960–967 (2019).Article 
CAS 

Google Scholar 
9.Ziv, G., Baran, E., Nam, S., Rodríguez-Iturbe, I. & Levin, S. A. Trading-off fish biodiversity, food security, and hydropower in the Mekong River Basin. Proc. Natl. Acad. Sci. U. S. A. 109, 5609–5614 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
10.Grill, G. et al. Mapping the world’s free-flowing rivers. Nature 569, 215–221 (2019).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Dudgeon, D. et al. Freshwater biodiversity: importance, threats, status and conservation challenges. Biol. Rev. Camb. Philos. Soc. 81, 163–182 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
12.Itakura, H. et al. Environmental DNA analysis reveals the spatial distribution, abundance, and biomass of Japanese eels at the river-basin scale. Aquat. Conserv. 29, 361–373 (2019).Article 

Google Scholar 
13.Wallace, A. R. The Geographical Distribution of Animals: With a Study of the Relations of Living and Extinct Faunas as Elucidating the Past Changes of the Earth’s Surface (Macmillan and Co, 1876).
Google Scholar 
14.Kreft, H. & Jetz, W. A framework for delineating biogeographical regions based on species distributions: Global quantitative biogeographical regionalizations. J. Biogeogr. 37, 2029–2053 (2010).Article 

Google Scholar 
15.Holt, B. G. et al. An update of Wallace’s zoogeographic regions of the world. Science 339, 74–78 (2013).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
16.Costello, M. J. et al. Marine biogeographic realms and species endemicity. Nat. Commun. 8, 1057 (2017).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
17.Leray, M. & Knowlton, N. DNA barcoding and metabarcoding of standardized samples reveal patterns of marine benthic diversity. Proc. Natl. Acad. Sci. U. S. A. 112, 2076–2081 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Deiner, K., Fronhofer, E. A., Mächler, E., Walser, J.-C. & Altermatt, F. Environmental DNA reveals that rivers are conveyer belts of biodiversity information. Nat. Commun. 7, 12544 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
19.Bush, A. et al. DNA metabarcoding reveals metacommunity dynamics in a threatened boreal wetland wilderness. Proc. Natl. Acad. Sci. U. S. A. 117, 8539–8545 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
20.Pawlowski, J., Apothéloz-Perret-Gentil, L. & Altermatt, F. Environmental DNA: What’s behind the term? Clarifying the terminology and recommendations for its future use in biomonitoring. Mol. Ecol. 29, 4258–4264 (2020).PubMed 
Article 
PubMed Central 

Google Scholar 
21.Hänfling, B. et al. Environmental DNA metabarcoding of lake fish communities reflects long-term data from established survey methods. Mol. Ecol. 25, 3101–3119 (2016).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
22.Li, J. et al. Ground-truthing of a fish-based environmental DNA metabarcoding method for assessing the quality of lakes. J. Appl. Ecol. 56, 1232–1244 (2019).CAS 
Article 

Google Scholar 
23.Li, J. et al. Limited dispersion and quick degradation of environmental DNA in fish ponds inferred by metabarcoding. Environ. DNA 1, 238–250 (2019).Article 

Google Scholar 
24.Pont, D. et al. Environmental DNA reveals quantitative patterns of fish biodiversity in large rivers despite its downstream transportation. Sci. Rep. 8, 10361 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
25.Olds, B. P. et al. Estimating species richness using environmental DNA. Ecol. Evol. 6, 4214–4226 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
26.Jerde, C. L. Can we manage fisheries with the inherent uncertainty from eDNA? J. Fish Biol. 98(2), 341–353 (2019).27.Bellemain, E. et al. Trails of river monsters: detecting critically endangered Mekong giant catfish Pangasianodon gigas using environmental DNA. Glob. Ecol. Conserv. 7, 148–156 (2016).Article 

Google Scholar 
28.Sakata, M. K., Maki, N., Sugiyama, H. & Minamoto, T. Identifying a breeding habitat of a critically endangered fish, Acheilognathus typus, in a natural river in Japan. Naturwissenschaften 104, 100 (2017).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
29.Mizumoto, H., Urabe, H., Kanbe, T., Fukushima, M. & Araki, H. Establishing an environmental DNA method to detect and estimate the biomass of Sakhalin taimen, a critically endangered Asian salmonid. Limnology 19, 219–227 (2018).CAS 
Article 

Google Scholar 
30.Cilleros, K. et al. Unlocking biodiversity and conservation studies in high-diversity environments using environmental DNA (eDNA): a test with Guianese freshwater fishes. Mol. Ecol. Resour. 19, 27–46 (2019).CAS 
PubMed 
Article 

Google Scholar 
31.Cantera, I. et al. Optimizing environmental DNA sampling effort for fish inventories in tropical streams and rivers. Sci. Rep. 9, 3085 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
32.Doble, C. J. et al. Testing the performance of environmental DNA metabarcoding for surveying highly diverse tropical fish communities: a case study from Lake Tanganyika. Environ. DNA 2, 24–41 (2020).Article 

Google Scholar 
33.Altermatt, F. et al. Uncovering the complete biodiversity structure in spatial networks: the example of riverine systems. Oikos 129, 607–618 (2020).Article 

Google Scholar 
34.Carraro, L., Mächler, E., Wüthrich, R. & Altermatt, F. Environmental DNA allows upscaling spatial patterns of biodiversity in freshwater ecosystems. Nat. Commun. 11, 3585 (2020).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
35.Zarfl, C. et al. Future large hydropower dams impact global freshwater megafauna. Sci. Rep. 9, 18531 (2019).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Kelly, R. P., Port, J. A., Yamahara, K. M. & Crowder, L. B. Using environmental DNA to census marine fishes in a large mesocosm. PLoS ONE 9, e86175 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
37.Riaz, T. et al. ecoPrimers: inference of new DNA barcode markers from whole genome sequence analysis. Nucleic Acids Res 39, e145 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
38.Miya, M. et al. MiFish, a set of universal PCR primers for metabarcoding environmental DNA from fishes: detection of more than 230 subtropical marine species. R. Soc. Open Sci. 2, 150088 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
39.Baselga, A. & Orme, C. D. L. betapart : an R package for the study of beta diversity : betapart package. Methods Ecol. Evol. 3, 808–812 (2012).Article 

Google Scholar 
40.Altermatt, F. Diversity in riverine metacommunities: a network perspective. Aquat. Ecol. 47, 365–377 (2013).Article 

Google Scholar 
41.Tonkin, J. D. et al. The role of dispersal in river network metacommunities: patterns, processes, and pathways. Freshw. Biol. 63, 141–163 (2018).Article 

Google Scholar 
42.Muneepeerakul, R. et al. Neutral metacommunity models predict fish diversity patterns in Mississippi-Missouri basin. Nature 453, 220–222 (2008).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
43.Azaele, S., Muneepeerakul, R., Maritan, A., Rinaldo, A. & Rodriguez-Iturbe, I. Predicting spatial similarity of freshwater fish biodiversity. Proc. Natl. Acad. Sci. U. S. A. 106, 7058–7062 (2009).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
44.Carrara, F., Altermatt, F., Rodriguez-Iturbe, I. & Rinaldo, A. Dendritic connectivity controls biodiversity patterns in experimental metacommunities. Proc. Natl. Acad. Sci. U. S. A. 109, 5761–5766 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
45.Muneepeerakul, R., Bertuzzo, E., Rinaldo, A. & Rodriguez-Iturbe, I. Evolving biodiversity patterns in changing river networks. J. Theor. Biol. 2019(462), 418–424 (2019).MATH 
Article 

Google Scholar 
46.Kang, B., Huang, X., Yan, Y., Yan, Y. & Lin, H. Continental-scale analysis of taxonomic and functional fish diversity in the Yangtze river. Glob. Ecol. Conserv. 15, e00442 (2018).Article 

Google Scholar 
47.Mächler, E. et al. Assessing different components of diversity across a river network using eDNA. Environ. DNA 1, 290–301 (2019).Article 

Google Scholar 
48.Carraro, L., Hartikainen, H., Jokela, J., Bertuzzo, E. & Rinaldo, A. Estimating species distribution and abundance in river networks using environmental DNA. Proc. Natl. Acad. Sci. U. S. A. 115, 11724–11729 (2018).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Roux, D. J. et al. Designing protected areas to conserve riverine biodiversity: lessons from a hypothetical redesign of the Kruger National Park. Biol. Conserv. 141, 100–117 (2008).Article 

Google Scholar 
50.Troia, M. J. & McManamay, R. A. Biogeographic classification of streams using fish community– and trait–environment relationships. Divers Distrib 26, 108–125 (2020).Article 

Google Scholar 
51.Vannote, R. L., Minshall, G. W., Cummins, K. W., Sedell, J. R. & Cushing, C. E. The river continuum concept. Can. J. Fish Aquat. Sci. 37, 130–137 (1980).Article 

Google Scholar 
52.He, Y., Wang, J., Lek, S., Cao, W. & Lek-Ang, S. Structure of endemic fish assemblages in the upper Yangtze River Basin. River Res Appl 27, 59–75 (2011).Article 

Google Scholar 
53.Lawson Handley, L. et al. Temporal and spatial variation in distribution of fish environmental DNA in England’s largest lake. Environ. DNA 1, 26–39 (2019).Article 

Google Scholar 
54.Taberlet, P., Bonin, A., Zinger, L. & Coissac, E. Environmental DNA: For Biodiversity Research and Monitoring (Oxford University Press, 2018).Book 

Google Scholar 
55.Monkolprasit, S., Sontirat, S., Vimollohakarn, S. & Songsirikul, T. Checklist of Fishes in Thailand: OEPP Biodiversity Series Vol. 4 (Office of Environmental Policy and Planning, 1997).
Google Scholar 
56.Weigand, H. et al. DNA barcode reference libraries for the monitoring of aquatic biota in Europe: gap-analysis and recommendations for future work. Sci. Total Environ. 678, 499–524 (2019).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
57.Blackman, R. et al. Detection of a new non-native freshwater species by DNA metabarcoding of environmental samples—first record of Gammarus fossarum in the UK. Aquat. Invasions 12, 177–189 (2017).Article 

Google Scholar 
58.Csárdi, G. & Nepusz, T. The igraph software package for complex network research. Inter Journal, Complex Systems 1695: 1–9 (2006). R package version 1.2.5. Available from https://cran.r-project.org/web/packages/igraph/index.html. Accessed 27 June 2020.59.Oksanen, J. et al. vegan: Community Ecology Package 2.5-6. Available from https://CRAN.R-project.org/package=vegan. Accessed 27 June 2020.60.Kindt, R. & Coe, R. Tree diversity analysis. A manual and software for common statistical methods for ecological and biodiversity studies. World Agroforestry Centre (ICRAF), Nairobi. ISBN 92-9059-179-X. Accessed 23 March 2021.61.Brock, G., Pihur, V., Datta, S. & Datta, S. clValid: An R Package for Cluster Validation. J Stat Softw 25: 1–22 (2008). R package version 0.6-9. Available from https://cran.r-project.org/web/packages/clValid/index.html. Accessed 27 June 2020.62.R Studio Team. RStudio: Integrated Development for R. RStudio, Inc., Boston, MA. (2019).63.Kassambara, A. & Mundt, F. factoextra: Extract and Visualize the Results of Multivariate Data Analyses. (2020) R package version 1.0.7. Available from: https://cran.r-project.org/web/packages/factoextra/index.html. Accessed 27 June 2020. More

Sample collectionThe 267 samples were collected between 2009 and 2018 from 41 locations expanding from the subarctic region to the tropics (Fig. 1, Auxillary Table S1)10 and processed using the same analytical pipeline (Fig. 2). The majority of the samples were collected using a depth-discrete Limnos tube-sampler (Limnos, Poland), with the exception of the samples from La Plata reservoir (Puerto Rico), which were collected using horizontal Van Dorn sampler (5 L capacity) and samples from Lake Loclat, which were collected using a deployed PVC-inlet connected to a peristaltic pump via tubing. Of all the lakes, 29 were sampled during the open water season and the majority of the lakes were sampled once. For 12 of the lakes only surface samples taken during the ice-covered period in winter were available, and one of the Swedish lakes (Lake Lomtjärnan) was sampled twice during the ice-covered period. Moreover, a total of 5 samples (one depth profile) from the time series of the Swiss lake (Loclat) were taken from under the ice. Time series samples were taken for Lake Loclat (seven time points, Auxillary Table S1)10 and for Lake Mekkojärvi (22 time points, see Saarenheimo et al.11 for details). For most lakes and ponds, samples were collected from multiple depths, including samples from the oxic surface layer (epilimnion), the layer with steepest change in oxygen concentration and temperature (metalimnion) and from the layer where oxygen levels were below the detection limit (hypolimnion). The exception to this were the 12 Swedish lakes sampled during ice-covered period, and five shallow ponds in Canada, for which only one sample from the oxic surface layer was taken (see Auxillary Table S1)10.Fig. 2Overview of the workflow from sample collection to mOTUs.Full size imageFrom two of the lakes, Lake Lomtjärnan in Sweden and Lake Alinen Mustajärvi in Finland, samples were collected also for single cell sorting. From both locations samples were preserved in glycerol-TE (gly-TE) and from Lomtjärnan samples were preserved also using phosphate buffered saline (PBS). For both preservants, the samples were flash frozen in liquid nitrogen after first incubating for 1 minute at ambient temperature.Simultaneous to collection of the DNA samples, also samples for environmental variables were taken. Variables included temperature, pH, conductivity, oxygen, total and dissolved nutrients (P and N species), gases (CO2 or dissolved inorganic carbon and methane (CH4)), total or dissolved organic carbon, iron, sulfate and chlorophyll a (Auxillary Table S1 and Auxillary Table S210 for the methods). As the samples were collected during multiple years and by different research groups, there was some variation for the procedures between the different sampling occasions, leading to variation in the final set of environmental data across the samples.DNA extraction and metagenome sequencingMost of the DNA samples were collected on 0.2 µm Sterivex filters (Millipore), except for the time-series samples collected from Loclat, which were collected by vacuum filtration onto 47 mm polycarbonate membrane filters with 0.2 μm pore size, and time series samples from Finnish Lake Mekkojärvi, for which the water for DNA extraction was collected from epilimnion (0–0.5 m), metalimnion (0.5–1 m) and hypolimnion (1–3 m) and pooled samples from each stratum were stored in 100 ml plastic containers and frozen at −20 °C and eventually freeze-dried (Alpha 1–4 LD plus, Christ). For all filter samples, water was filtered until the filter clogged. All filters were stored frozen (−20 to −80 °C) until the extraction of DNA. For all samples, DNA was extracted using PowerSoil DNA extraction kit (MoBio, Carlsbad, CA, USA) following the manufacturer’s instructions and the DNA concentrations were measured using Qubit dsDNA HS kit (Thermo Fisher Scientific Inc.).Sequencing libraries were prepared from 10 or 20 ng of DNA using the ThruPLEX DNA-seq Prep Kit according to the manufacturer’s preparation guide. Briefly, the DNA was fragmented using a Covaris E220 system, aiming at 400 bp fragments. The ends of the fragments were end-repaired and stem-loop adapters were ligated to the 5′ ends of the fragments. The 3′ end of the stem loop were subsequently extended to close the nick. Finally, the fragments were amplified and unique index sequences were introduced using 7 cycles of PCR followed by purification using AMPure XP beads (Beckman Coulter).The quality of the libraries was evaluated using the Agilent Fragment Analyzer system and the DNF-910-kit. The adapter-ligated fragments were quantified by qPCR using the Library quantification kit for Illumina (KAPA Biosystems/Roche) on a CFX384Touch instrument (BioRad) prior to cluster generation and sequencing.The sequencing libraries were pooled and subjected to cluster generation and paired-end sequencing with 150 bp read length S2/S4 flow-cells and the NovaSeq 6000 system (Illumina Inc.) using the v1 chemistry according to the manufacturer’s protocols. Negative controls were included to the sequencing as well as 1% of PhiX control library as a positive control.Base calling was done on the instrument by RTA (v3.3.3, 3.3.5, 3.4.4) and the resulting.bcl files were demultiplexed and converted to fastq format with tools provided by Illumina Inc., allowing for one mismatch in the index sequence. Additional statistics on sequence quality were compiled with an in-house script from the fastq-files, RTA and CASAVA output files. Sequencing was performed by the SNP&SEQ Technology Platform in Uppsala, Sweden.Single-cell sorting and DNA amplificationAll Gly-TE cryopreserved samples were thawed and diluted in 1 xPBS if needed while all plates with PBS were UV-treated with a dose of 2 J prior to sorting. Samples collected from both lakes were sorted, and then screened for organisms belonging to candidate phyla radiation. Samples collected from Lake Lomtjärnan were additionally subjected to sorting based on autofluorescence to identify and sequence cells belonging to lineage Chlorobia.For obtaining SAGs from representatives of the candidate phyla radiation (CPR), samples were first stained with 1 x SYBR Green I for approximate 30 minutes. Subsequent single cell sorting was performed with a MoFlo Astrios EQ (Beckman Coulter, USA) cell sorter using a 488 nm laser for excitation, 70 µm nozzle, sheath pressure of 60 psi and 0.1 µm sterile filtered 1x PBS as sheath fluid. Individual cells were deposited into empty 384-well plates (Biorad, CA USA) UVed at 2 Joules using a CyCloneTM robotic arm and the most stringent single cell sort settings (single mode, 0.5 drop envelope). Green fluorescence (488–530/40) was used as trigger and sort decisions were made based on combined gates of 488–530/40 Height log vs 488–530/40 Area log and 488–530/40 Height log vs SSC with increasing side scatter divided up in three different regions. Flow sorting data was interpreted and displayed using the associated software Summit v 6.3.1. Next, individual cells were subject to lysis, neutralization and whole genome amplification using MDA based on the protocol and workflow described by Rinke et al.12 but with several modifications. Reagent mastermixes were added using the MANTIS liquid dispenser (Formulatrix) and the LV or HV silicone chips. The lysozyme, D2 buffer, stop solution and MDA-mastermix were each dispensed with its own chip. Most MDA-reactions were run using the phi29 from ThermoFisher but a few were run with a more heat-stable phi29, EquiPhi also provided by ThermoFisher. The MDA reaction was carried out in a total volume of 5.2 µl. Thawed, sorted cells were first pre-treated with 400 nl/well of 12 U/µl of Ready-Lyse™ Lysozyme Solution (R1804M, Lucigen) at room temperature for 15 minutes before adding 400 nl Qiagen lysis buffer D2 followed by incubation at 95 °C for 10 seconds and 10 minutes on ice. Reactions were neutralized by adding 400 nl Qiagen Stop solution. Four µl of MDA mix containing 1x reaction buffer, 0.4 mM dNTP, 0.05 mM exonuclease-resistant Hexamers, 10 mM DTT, 1.7 U phi29 DNA polymerase (ThermoFisher Scientific) and 0.5 µM Syto13 was added to a final reaction volume of 5.2 µl. All reagents except SYTO13 were UV decontaminated at 2 Joules in a UV crosslinker. The whole genome amplification was run at 30 °C for 7 or 10 h followed by an inactivation step at 65 °C for 5 min. The reaction was monitored in real time by detection of SYTO13 fluorescence every 15 minutes using a FLUOstar® Omega plate reader (BMG Labtech, Germany) or a qPCR instrument. The EquiPhi protocol was run as previously described for ThermoFisher phi29 with the following exceptions; the EquiPhi polymerase was added in 1U/reaction, reaction buffer included with the polymerase was used and the reaction was carried out at 45 °C. The single amplified genome (SAG) DNA was stored at −20 °C until further PCR screening, library preparation and Illumina sequencing.The CPR SAGs were screened using the bacterial PCR primers targeting the 16 S rRNA gene, Bact_341 F and Bact_805 R13. The reactions were run in a LightCycler 480 PCR machine (ROCHE, MA USA) in 10 µl and a final concentration of 1 x LightCycler480 SYBR Green I Master mix, 0.25 µM of each primer and 2 µl of 60 to 80 times diluted SAGs. Following a 3 min denaturation at 95 °C, targets were amplified for 40 cycles of 95 °C for 10 s, 55 °C for 20 s, 72 °C for 30 s and a final 10 min extension at 72 °C followed by melting curve analysis. The products were purified using the NucleoSpin Gel and PCR clean-up purification kit (Macherey-Nagel, Germany), quantified using the Quant-iT TM PicoGreen® dsDNA assay kit (Invitrogen, MA USA) in a FLUOstar® Omega microplate reader (BMG Labtech, Germany) and submitted for identification by Sanger sequencing at Eurofin Genomics. All SAGs were further screened using the newly designed primers targeting the phylum Parcubacteria 684F-OD1 (3′ GTAGKRRTRAAATSCGTT 5′) and 784 R (5′ TAMNVGGGTATCTAATCC -3′). These primers target with good specificity 67% of Parcubacteria in the SILVA database14. Parcu-PCR was run at 3 min at 95 °C, 40 cycles of 95 °C for 10 s, 55 °C for 20 s, 72 °C for 30 s and a final 10 min extension at 72 °C followed by melting curve analysis. The products were purified using the NucleoSpin Gel and PCR clean-up purification kit (Macherey-Nagel, Germany), quantified using the Quant-iT TM PicoGreen® dsDNA assay kit (Invitrogen, MA USA) in a FLUOstar® Omega microplate reader (BMG Labtech, Germany) and submitted for identification by Sanger sequencing at Eurofin Genomics.To recover Chlorobia single amplified genomes, sorting was done in 2016 on a MoFlo™ Astrios EQ sorter (Beckman Coulter, USA) using a 488 and 532 nm laser for excitation, a 70 μm nozzle, a sheath pressure of 60 psi, and 0.1 μm filtered 1x PBS as sheath fluid. An ND filter ND = 1 and the masks M1 and M2 were used. The trigger channel was set to the forward scatter (FSC) at a threshold of 0.025% and sort regions were defined on autofluorescence using laser 532 nm and band pass filters 710/45 and 664/22. Three populations were sorted based on differences in autofluorescence signals. The sort mode was set to single cell with a drop envelope of 0.5. The target populations were sorted at approximately 400 events per second into 96-well plates containing 1 µl 1x PBS per well with either 1 or 10 cells (positive control) deposited. A few wells remained empty (no cell sorted) were kept as negative controls. Sorted plates were stored frozen at −80 °C.The subsequent whole genome amplification was performed in 2018 using the REPLI-g Single Cell kit (QIAGEN) following the instructions provided by the manufacturer but with total reaction volume reduced to 12.5 µl. The denaturation reagent D2, stop solution, water, and reagent tubes and strips were UV-treated at 2.5 J. The lysis was changed slightly to 10 min at 65 °C, followed by 5 min on ice before adding the stop solution. To the master mix containing water, reaction buffer, and the DNA the polymerase we added SYTO 13 (Invitrogen) at a final concentration of 0.5 µM. The amplification was performed at 30 °C for 8 hours in a plate reader with fluorescence readings every 15 min. The reaction was stopped by incubating it for 5 min at 65 °C. The plate was stored for less than a week at −20 °C. Amplified DNA was mixed thoroughly by pipetting up and down 20 times before diluting it 50x and 100x in nuclease-free water. The DNA was screened for bacterial 16 S rRNA applying the primers Bact_341 F (5′- CCTACGGGNGGCWGCAG- 3′) and Bact_805 R (5′- GACTACHVGGGTATCTAATCC-3′)13 using the LightCycler® 480 SYBR Green I Master (Roche) kit. The PCR mix contained 1.5 µl diluted amplified DNA, 1x the LightCycler® 480 SYBR Green I Master mix, 0.25 µM of each primer, and nuclease-free water in a total reaction volume of 10 µl. The PCR cycling (5 min at 95 °C, followed by 40 cycles of 10 sec at 95 °C, 20 sec at 60 °C, 30 sec at 72 °C) was followed by meltcurve analysis on the LightCycler® 480 Instrument (Roche). DNA of confirmed Chlorobia was sent to sequencing as outlined below.Library preparation and Illumina sequencing of the single cellsFor the CPR-targeted analysis, Illumina libraries were prepared from sixty SAGs mainly selected from the screening procedure in a PCR-free workflow using the sparQ DNA Frag & Library Prep Kit (Quantabio) and IDT for Illumina TruSeq UD Indexes (Illumina). Libraries were prepared from 50–250 ng of MDA-products in 25% of the recommended reaction volumes according to manufacturer’s instructions. The MDA-products were fragmented for 7 minutes (5 minutes for 4 samples) without using the DNA Frag Enhancer Solution. Library insert sizes were determined using Bioanalyzer High Sensitivity DNA Kit (Agilent). Each library was quantified using the KAPA Library Quantification kit (Roche) in 5 µl reaction volumes in a 384-well plate run on LightCycler 480 (Roche) to allow equimolar pooling before sequencing on Illumina HiSeqX v2.5 PE 2 × 150 bp including negative and positive (PhiX) controls.For the Chlorobia-targeted sequencing, amplified DNA from 23 SAGs were quantified individually with Qubit dsDNA HS assay kit (ThermoFisher Scientific) and diluted to 0.2 ng/ul in nuclease free water. Sequencing libraries were prepared with Nextera XT DNA Library Preparation Kit and combinatorial combinations of molecular identifiers in the Nextera XT Index Kit (Illumina, CA USA) according to manufacturer’s instructions. Libraries with an average length of 1200 bp were quantified with Qubit dsDNA HS assay kit to allow pooling of equal amounts of the libraries based on mass. The libraries were sequenced on an Illumina MiSeq v3 PE 2 × 300 bp including negative and positive (PhiX) controls.Data processing of the metagenome and single cell sequencesThe metagenome sequencing resulted in a total of ~107 paired-end reads of length 2 × 150 bp, amounting to a total of total 3 Tbp. The raw data was trimmed using Trimmomatic (version 0.36; parameters: ILLUMINACLIP:TruSeq 3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36)15 (Auxillary Table S3)10. The trimmed data was assembled using Megahit (version 1.1.13)16 with default settings. Two types of assemblies were done, single sample assemblies for all the samples individually and a total of 53, mainly lake-wise, co-assemblies (see Auxillary Table S4)10, some samples of the Canadian ponds have also been coassembled with previously sequenced libraries of the same sample (see Auxillary Table S5)10. The relevant quality controlled reads were mapped to all the assemblies using BBmap17 with default settings and the mapping results were used to bin the contigs using Metabat (version 2.12.1, parameters –maxP 93 –minS 50 -m 1500 -s 10000)18. Genes of obtained bins were predicted and annotated using Prokka (version 1.13.3)19 using standard parameters except for the bin containing all the unbinned contigs where the –metagenome flag was used. Single-cell libraries were processed similarly to the metagenomes, but without the binning step, and using the single-cell variant of the SPAdes20 assembler instead of Megahit.Prokaryotic completeness and redundancy of all bins from Metabat and for all assembled single cells were computed using CheckM (version 1.0.13)21 (Auxillary Tables S6 and S7 for MAGs and SAGs, respectively)10. Average Nucleotide Identity (ANI) for all bin-pairs was computed with fastANI (version 1.3)22. The bins were clustered into metagenomic Operational Taxonomic Units (mOTUs) starting with 40% complete genomes with less than 5% contamination. Genome pairs with ANI above 95% were clustered into connected components. Additionally, less complete genomes were recruited to the mOTU if its ANI similarity was above 95%. Bins were taxonomically annotated in a two-step process. GTDB-Tk (version 102 with database release 89)23 was used first with default settings. Using this classification an lca database for SourMASH (version 1.0)24 was made. This database as well as one based on the GTDB release 89 was then used with SourMASH’s lca classifier for a second round of classification of bins that were not annotated with GTDB-tk (Auxillary Table S8)10.The taxonomic diversity of the bacterial (Fig. 3) and archaeal (Fig. 4) mOTUs, respectively, were visualized in a tree format. The trees were computed using GTDB-tk with one representative MAG per mOTU of the stratfreshDB, and one random representative genome per family of the GTDB. Trees were visualized using anvi’o25.Fig. 3Bacterial diversity of the stratfreshDB27. The insert illustrates the quality of the MAGs and SAGs included in the tree. Interactive version of the tree with more information available at https://anvi-server.org/moritzbuck/bacterial_diversity_of_the_stratfreshdb.Full size imageFig. 4Archaeal diversity of the stratfreshDB27. The insert illustrates the quality of the MAGs included in the tree. Interactive version of the tree with more information available at https://anvi-server.org/moritzbuck/archaeal_diversity_of_the_stratfreshdb.Full size image More