Philippa Kaur

Godfray, H. C. J. et al. Science 327, 812–818 (2010).ADS 
CAS 
Article 

Google Scholar 
Pimm, S. L. et al. Science 344, 1246752 (2014).CAS 
Article 

Google Scholar 
Chung, M. G. & Liu, J. Nat. Food https://doi.org/10.1038/s43016-022-00499-7 (2022).Myers, N., Mittermeier, R. A., Mittermeier, C. G., Da Fonseca, G. A. & Kent, J. Nature 403, 853–858 (2000).ADS 
CAS 
Article 

Google Scholar 
A complex prairie ecosystem. National Park Service https://www.nps.gov/tapr/learn/nature/a-complex-prairie-ecosystem.htm (2022)Davalos, L. M. et al. Environ. Sci. Technol. 45, 1219–1227 (2011).ADS 
CAS 
Article 

Google Scholar 
Vijay, V., Pimm, S. L., Jenkins, C. N. & Smith, S. J. PLoS ONE 11, e0159668 (2016).Article 

Google Scholar 
Liu, J. et al. Ecol. Soc. 18, 26 (2013).CAS 
Article 

Google Scholar 
Liu, J. Consumption patterns and biodiversity. The Royal Society https://go.nature.com/3M19vup (2020).Xu, Z. et al. Nat. Sustain. 3, 964–971 (2020).Article 

Google Scholar 
Dou, Y., da Silva, R. F. B., Yang, H. & Liu, J. J. Geogr. Sci. 28, 1715–1732 (2018).Article 

Google Scholar  More

Siński, E. & Welc-Falęciak, R. Risk of Infections Transmitted by Ticks in Forest Ecosystems of Poland. (Zarządzanie Ochroną Przyrody w Lasach 6, 2012).Kantar Public. Zwierzęta w polskich domach, 2017. http://www.tnsglobal.pl/archiwumraportow/files/2017/05/K.021_Zwierzeta_domowe_O04a-17.pdf.Maia, C. et al. Bacterial and protozoal agents of feline vector-borne diseases in domestic and stray cats from southern Portugal. Parasit. Vectors. 7, 115. https://doi.org/10.1186/1756-3305-7-115 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
Maia, C. et al. Bacterial and protozoal agents of canine vector-borne diseases in the blood of domestic and stray dogs from southern Portugal. Parasit. Vectors. 23(8), 138. https://doi.org/10.1186/s13071-015-0759-8 (2015).Article 

Google Scholar 
Baturo, I.M. Parki narodowe i krajobrazowe, rezerwaty przyrody. (Departament Turystyki, Sportu, Promocji Urzędu Marszłkowskiego Województwa Małopolskiego, 2010).Dulias, R. & Hibszer, A. Województwo śląskie—przyroda, gospodarka, dziedzictwo kulturowe (Wydawnictwo Kubajak, 2004).
Google Scholar 
Siuda, K. Kleszcze Polski Acari Ixodida). Część II. Systematyka i Rozmieszczenie (Polskie Towarzystwo Parazytologiczne, 1993).
Google Scholar 
Nowak-Chmura, M. Fauna kleszczy (Ixodida) Europy Środkowej (Wydawnictwo Naukowe Uniwersytetu Pedagogicznego, 2013).
Google Scholar 
Rijpkema, S., Golubić, D., Molkenboer, M., Verbeek-De Kruif, N. & Schellekens, J. Identification of four genomic groups of Borrelia burgdorferi sensu lato in Ixodes ricinus ticks collected in a Lyme borreliosis endemic region of northern Croatia. Exp. Appl. Acarol. 20, 23–30 (1996).CAS 
Article 

Google Scholar 
Wójcik-Fatla, A., Szymańska, J., Wdowiak, L., Buczek, A. & Dutkiewicz, J. Coincidence of three pathogens (Borrelia burgdorferi sensu lato, Anaplasma phagocytophilum and Babesia microti) in Ixodes ricinus ticks in the Lublin makroregion. Ann. Agric. Environ. Med. 16(1), 151–158 (2009).PubMed 

Google Scholar 
Massung, R. F. et al. Nested PCR assay for detection of granulocytic ehrlichiae. J. Clin. Microbiol. 36(4), 1090–1095 (1998).CAS 
Article 

Google Scholar 
Persing, D. H. et al. Detection of Babesia microti by polymerase chain reaction. J. Clin. Microbiol. 30, 2097–2103 (1992).CAS 
Article 

Google Scholar 
Sroka, J., Szymańska, J. & Wójcik-Fatla, A. The occurence of Toxoplasma gondii and Borrelia burgdorferi sensu lato in Ixodes ricinus ticks from east Poland with the use of PCR. Ann. Agric. Environ. Med. 16(2), 313–319 (2009).PubMed 

Google Scholar 
Siuda, K., Nowak, M., Gierczak, M., Wierzbowska, I. & Faber, M. Kleszcze (Acari: Ixodida) pasożytujące na psach i kotach domowych w Polsce. Wiad. Parazytol. 53, 155 (2007).
Google Scholar 
Zajkowska, P. Ticks (Acari:Ixodida) attacking domestic dogs in the Malopolska voivodeship, Poland. In Arthropods: In the contemporary world (eds Buczek, A. & Błaszak, C. Z.) 87–99 (Koliber, 2015).Chapter 

Google Scholar 
Szymański, S. Przypadek masowego rozwoju kleszcza Rhipicephalus sanguineus (Latreile, 1806) w warszawskim mieszkaniu. Wiad. Parazytol. 25, 453–458 (1979).PubMed 

Google Scholar 
Król, N., Obiegala, A., Pfeffer, M., Lonc, E. & Kiewra, D. Detection of selected pathogens in ticks collected from cats and dogs in the Wrocław Agglomeration South-West Poland. Parasit. Vectors. 9(1), 351. https://doi.org/10.1186/s13071-016-1632-0 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Kocoń, A., Nowak-Chmura, M., Kłyś, M. & Siuda, K. Ticks (Acari: Ixodida) attacking domestic cats (Felis catus L.) in southern Poland. In Arthropods in Urban and Suburban Environments (eds Buczek, A. & Błaszak, C.) 51–61 (Koliber, 2017).
Google Scholar 
Roczeń-Karczmarz, M. et al. Comparison of the occurrence of tick-borne diseases in ticks collected from vegetation and animals in the same area. Med. Weter. 74(8), 484–488. https://doi.org/10.21521/mw.6107 (2018).Article 

Google Scholar 
Cuber, P., Asman, M., Solarz, K., Szilman, E. & Szilman, P. Pierwsze stwierdzenia obecności wybranych patogenów chorób transmisyjnych w kleszczach Ixodes ricinus (Acari: Ixididae) zebranych w okolicach zbiorników wodnych w Rogoźniku (województwo śląskie) in Stawonogi. Ekologiczne i patologiczne aspekty układu pasożyt – żywiciel (eds. Buczek, C. & Błaszak, Cz.). 155-164 (Akapit, Lublin, 2010).Asman, M. et al. The risk of exposure to Anaplasma phagocytophilum, Borrelia burgdorferi sensu lato, Babesia sp. and co-infections in Ixodes ricinus ticks on the territory of Niepołomice Forest (southern Poland). Ann. Parasitol. 59(1), 13–19 (2013).PubMed 

Google Scholar 
Pawełczyk, O. et al. The PCR detection of Anaplasma phagocytophilum, Babesia microti and Borrelia burgdorferi sensu lato in ticks and fleas collected from pets in the Będzin district area (Upper Silesia, Poland) – the preliminary studies in Stawonogi: zagrożenie zdrowia człowieka i zwierząt (eds. Buczek, C. & Błaszak, Cz.). 111–119 (Koliber, Lublin, 2014).Strzelczyk, J. K. et al. Prevalence of Borrelia burgdorferi sensu lato in Ixodes ricinus ticks collected from southern Poland. Acta Parasitol. 60(4), 666–674. https://doi.org/10.1515/ap-2015-0095 (2015).CAS 
Article 
PubMed 

Google Scholar 
Zygner, W. & Wędrychowicz, H. Occurrence of hard ticks in dogs from Warsaw area. Ann. Agric. Environ. Med. 13(2), 355–359 (2006).PubMed 

Google Scholar 
Kilar, P. Ticks attacking domestic dogs in the area of the Rymanów district, Subcarpathian province Poland. Wiad. Parazytol. 57(3), 189–1991 (2011).PubMed 

Google Scholar 
Claerebout, E. et al. Ticks and associated pathogens collected from dogs and cats in Belgium. Parasit. Vectors. 6, 183. https://doi.org/10.1186/1756-3305-6-183 (2013).Article 
PubMed 
PubMed Central 

Google Scholar 
Schreiber, C. et al. Pathogens in ticks collected from dogs in Berlin/Brandenburg, Germany. Parasit. Vectors. 7, 535. https://doi.org/10.1186/s13071-014-0535-1 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
Eichenberger, R. M., Deplazes, P. & Mathis, A. Ticks on dogs and cats: A pet owner-based survey in a rural town in northeastern Switzerland. Ticks Tick-borne Dis. 6, 267–271. https://doi.org/10.1016/j.ttbdis.2015.01.007 (2015).Article 
PubMed 

Google Scholar 
Michalski, M. M. Skład gatunkowy kleszczy psów (Acari: Ixodida) z terenu aglomeracji miejskiej w cyklu wieloletnim. Med. Weter. 73(11), 698–701 (2017).
Google Scholar 
Geurden, T. et al. Detection of tick-borne pathogens in ticks from dogs and cats in different European countries. Ticks. Tick. Borne. Dis. 9(6), 1431–1436. https://doi.org/10.1016/j.ttbdis.2018.06.013 (2018).Article 
PubMed 

Google Scholar 
Namina, A. et al. Tick-borne pathogens in ticks collected from dogs, Latvia, 2011–2016. BMC Vet. Res. 15, 398. https://doi.org/10.1186/s12917-019-2149-5 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Bhide, M., Travnicek, M., Curlik, J. & Stefancikova, A. The importance of dogs in eco-epidemiology of Lyme borreliosis: a review. Vet. Med. Czech 49(4), 135–142 (2004).Article 

Google Scholar 
Burgess, E. C. Experimentally induced infection of cats with Borrelia burgdorferi. Am. J. Vet. Res. 53, 1507–1511 (1992).CAS 
PubMed 

Google Scholar 
Skotarczak, B. & Wodecka, B. Identification of Borrelia burgdorferi genospecies inducing Lyme disease in dogs from western Poland. Acta Vet. Hung. 53(1), 13–21 (2005).Article 

Google Scholar 
Skotarczak, B. et al. Prevalence of DNA and antibodies to Borrelia burgdorferi sensu lato in dogs suspected of borreliosis. Ann. Agric. Environm. Med. 12(2), 199–205 (2005).CAS 

Google Scholar 
Adaszek, Ł, Winiarczyk, S., Kutrzeba, J., Puchalski, A. & Dębiak, P. Przypadki boreliozy u psow na Lubelszczyźnie. Życie Wet. 83, 311–313 (2008).
Google Scholar 
Hovius, K. E. Borreliosis. In Arthropod-borne Infectious Diseases of the Dog and Cat (eds Shaw, S. E. & Day, M. J.) 100–109 (Manson Publishing, 2005).Chapter 

Google Scholar 
Zygner, W., Jaros, S. & Wędrychowicz, H. Prevalence of Babesia canis, Borrelia afzelii, and Anaplasma phagocytophilum infection in hard ticks removed from dogs in Warsaw (central Poland). Vet. Parasitol. 153, 139–142. https://doi.org/10.1016/j.vetpar.2008.01.036 (2008).Article 
PubMed 

Google Scholar 
Welc-Falęciak, R., Rodo, A., Siński, E. & Bajer, A. Babesia canis and other tick-borne infections in dogs in Central Poland. Vet. Parasitol. 166(3–4), 191–198. https://doi.org/10.1016/j.vetpar.2009.09.038 (2009).Article 
PubMed 

Google Scholar 
Michalski, M. M., Kubiak, K., Szczotko, M., Chajęcka, M. & Dmitryjuk, M. Molecular Detection of Borrelia burgdorferi Sensu Lato and Anaplasma phagocytophilum in Ticks Collected from Dogs in Urban Areas of North-Eastern Poland. Pathogens. 9(6), 455. https://doi.org/10.3390/pathogens9060455 (2020).CAS 
Article 
PubMed Central 

Google Scholar 
Nijhof, A. M. et al. Ticks and associated pathogens collected from domestic animals in the Netherlands. Vector. Borne. Zoonot. Dis. 7, 585–595. https://doi.org/10.1089/vbz.2007.0130 (2007).Article 

Google Scholar 
Adaszek, Ł, Martinez, A. C. & Winiarczyk, S. The factors affecting the distribution of babesiosis in dogs in Poland. Vet. Parasitol. 181, 160–165. https://doi.org/10.1016/j.vetpar.2011.03.059 (2011).Article 
PubMed 

Google Scholar 
Adaszek, Ł, Łukaszewska, J., Winiarczyk, S. & Kunkel, M. Pierwszy przypadek babeszjozy u kota w Polsce. Życie Wet. 83(8), 668–670 (2008).
Google Scholar 
Kocoń, A. et al. Molecular detection of tick-borne pathogens in ticks collected from pets in selected mountainous areas of Tatra County (Tatra Mountains, Poland). Sci. Rep. 10, 15865. https://doi.org/10.1038/s41598-020-72981-w (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Asman, M. et al. Detection of protozoans Babesia microti and Toxoplasma gondii and their co-existence in ticks (Acari: Ixodida) collected in Tarnogórski district (Upper Silesia, Poland). Ann. Agric. Environ. Med. 22(1), 80–83. https://doi.org/10.5604/12321966.1141373 (2015).Article 
PubMed 

Google Scholar 
Stensvold, C. R. et al. Babesia spp. and other pathogens in ticks recovered from domestic dogs in Denmark. Parasit. Vectors. 8(8), 262. https://doi.org/10.1186/s13071-015-0843-0 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Abdullah, S., Helps, C., Tasker, S., Newbury, H. & Wall, R. Prevalence and distribution of Borrelia and Babesia species in ticks feeding on dogs in the UK. Med. Vet. Entomol. 32(1), 14–22. https://doi.org/10.1111/mve.12257 (2018).CAS 
Article 
PubMed 

Google Scholar 
Bjoersdorff, A., Svendenius, L., Owens, J. H. & Massung, R. F. Feline granulocytic ehrlichiosis– a report of a new clinical entity and characterisation of the infectious agent. J. Small. Anim. Pract. 40(1), 20–24. https://doi.org/10.1111/j.1748-5827.1999.tb03249 (1999).CAS 
Article 
PubMed 

Google Scholar 
Lappin, M. R. et al. Molecular and serologic evidence of Anaplasma phagocytophilum infection in cats in North America. J. Am. Vet. Med. Assoc. 225(6), 893–896. https://doi.org/10.2460/javma.2004.225.893 (2004).Article 
PubMed 

Google Scholar 
Shaw, S. E. et al. Molecular evidence of tick-transmitted infections in dogs and cats in the United Kingdom. Vet. Rec. 157(21), 645–648. https://doi.org/10.1136/vr.157.21.645 (2005).CAS 
Article 
PubMed 

Google Scholar 
Tarello, W. Microscopic and clinical evidence for Anaplasma (Ehrlichia) phagocytophilum infection in Italian cats. Vet. Rec. 156(24), 772–774. https://doi.org/10.1136/vr.156.24.772 (2005).CAS 
Article 
PubMed 

Google Scholar 
Schaarschmidt-Kiener, D., Graf, F., von Loewenich, F. D. & Muller, W. Anaplasma phagocytophilum infection in a cat in Switzerland. Schweiz. Arch. Tierheilkd. 151(7), 336–341. https://doi.org/10.1024/0036-7281.151.7.336 (2009).CAS 
Article 
PubMed 

Google Scholar 
Heikkila, H. M., Bondarenko, A., Mihalkov, A., Pfister, K. & Spillmann, T. Anaplasma phagocytophilum infection in a domestic cat in Finland. Acta. Vet. Scand. 52(1), 62. https://doi.org/10.1186/1751-0147-52-62 (2010).Article 
PubMed 
PubMed Central 

Google Scholar 
Hamel, D., Bondarenko, A., Silaghi, C., Nolte, I. & Pfister, K. Seroprevalence and bacteremia of Anaplasma phagocytophilum in cats from Bavaria and Lower Saxony (Germany). Berl. Munch. Tierarztl. Wochenschr. 125(3–4), 163–167 (2012).PubMed 

Google Scholar 
Morgenthal, D. et al. Prevalence of haemotropic Mycoplasma spp., Bartonella spp. and Anaplasma phagocytophilum in cats in Berlin/Brandenburg (Northeast Germany). Berl. Munch Tierarztl. Wochenschr. 125(9–10), 418–427 (2012).PubMed 

Google Scholar 
Adaszek, Ł, Winiarczyk, S. & Łukaszewska, J. A first case of ehrlichiosis in a horse in Poland. Dtsch. Tierarztl. Wchschr. 116(9), 330–334 (2009).
Google Scholar 
Adaszek, Ł, Policht, K., Gorna, M., Kutrzuba, J. & Winiarczyk, S. Pierwszy w Polsce przypadek anaplazmozy (erlichiozy) granulocytarnej u kota. Życie Wet. 86, 132–135 (2011).
Google Scholar 
Adaszek, Ł, Kotowicz, W., Klimiuk, P., Gorna, M. & Winiarczyk, S. Ostry przebieg anaplazmozy granulocytarnej u psa—przypadek własny. Wet. w Praktyce 9, 59–62 (2011).
Google Scholar 
Adaszek, Ł et al. Three clinical cases of Anaplasma phagocytophilum infection in cats in Poland. J. Feline Med. Surg. 15, 333–337. https://doi.org/10.1177/1098612X12466552 (2013).Article 
PubMed 

Google Scholar 
Pusterla, N. et al. Seroprevalence of Ehrlichia canis and of canine granulocytic ehrlichia infection in dogs in Switzerland. J. Clin. Microbiol. 36, 3460–3462. https://doi.org/10.1128/JCM.36.12.3460-3462.1998 (1998).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Egenvall, A. et al. Detection of granulocytic Ehrlichia species DNA by PCR in persistently infected dogs. Vet. Rec. 146(7), 186–190. https://doi.org/10.1136/vr.146.7.186 (2000).CAS 
Article 
PubMed 

Google Scholar 
Shaw, S. E. et al. Review of exotic infectious dise-ases in small animals entering the United Kingdom from abroad diagnosed by PCR. Vet. Rec. 152(6), 176–177. https://doi.org/10.1136/vr.152.6.176 (2003).CAS 
Article 
PubMed 

Google Scholar 
Skotarczak, B., Adamska, M. & Supron, M. Blood DNA analysis for Ehrlichia (Anaplasma) phagocytophila and Babesia spp in dogs from Northern Poland. Acta Vet. Brno. 73, 347–351. https://doi.org/10.1136/vr.152.6.176 (2004).Article 

Google Scholar 
Adaszek, Ł. Wybrane Aspekty Epidemiologii Babeszjozy, Boreliozy i Erlichiozy Psów (Praca doktorska, 2007).
Google Scholar 
Kybicová, K. et al. Detection of Anaplasma phagocytophilum and Borrelia burgdorferi sensu lato in dogs in the Czech Republic. Vec. Born Zoon Dis. 9(6), 655–661. https://doi.org/10.1089/vbz.2008.0127 (2009).Article 

Google Scholar  More

Gordon, D. M. From division of labor to the collective behavior of social insects. Behav. Ecol. Sociobiol. 70, 1101–1108 (2016).PubMed 
Article 

Google Scholar 
Gordon, D. M. The organization of work in social insect colonies. Nature 380, 121–124 (1996).CAS 
Article 

Google Scholar 
Bonabeau, E., Theraulaz, G. & Deneubourg, J.-L. Quantitative study of the fixed threshold model for the regulation of division of labour in insect societies. Proc. R. Soc. Lond. Ser. B Biol. Sci. 263, 1565–1569 (1996).Article 

Google Scholar 
Pankiw, T. & Page, R. E. Jr. The effect of genotype, age, sex, and caste on response thresholds to sucrose and foraging behavior of honey bees (Apis mellifera L.). J. Comp. Physiol. A 185, 207–213 (1999).CAS 
PubMed 
Article 

Google Scholar 
Bonabeau, E., Sobkowski, A., Theraulaz, G. & Deneubourg, J.-L. Adaptive task allocation inspired by a model of division of labor in social insects. In BCEC 36–45 (1997).Robinson, G. E. & Page, R. E. J. Genetic basis for division of labor in an insect society. In The Genetics of Social Evolution (ed. Breed, R. P.) 61–80 (Westview Press, 1989).
Google Scholar 
Hogeweg, P. & Hesper, B. The ontogeny of the interaction structure in bumble bee colonies: A MIRROR model. Behav. Ecol. Sociobiol. 12, 271–283 (1983).Article 

Google Scholar 
Theraulaz, G., Bonabeau, E. & Denuebourg, J. N. Response threshold reinforcements and division of labour in insect societies. Proc. R. Soc. Lond. Ser. B Biol. Sci. 265, 327–332 (1998).Article 

Google Scholar 
Robinson, G. E. Labor in insect societies. Annu. Rev. Entomol. 37, 637–665 (1992).CAS 
PubMed 
Article 

Google Scholar 
Hölldobler, B. & Wilson, E. O. The Superorganism: The Beauty, Elegance, and Strangeness of Insect Societies (WW Norton & Company, 2009).
Google Scholar 
Gordon, D. M. The organization of work in social insect colonies. Complexity 8, 43–46 (2002).Article 

Google Scholar 
Bourke, A. F. G. & Franks, N. R. Social Evolution in Ants (Princeton University Press, 1995).
Google Scholar 
Robinson, E. J. H., Feinerman, O. & Franks, N. R. Flexible task allocation and the organization of work in ants. Proc. R. Soc. B Biol. Sci. 276, 4373–4380 (2009).Article 

Google Scholar 
Pinter-Wollman, N., Hubler, J., Holley, J.-A., Franks, N. R. & Dornhaus, A. How is activity distributed among and within tasks in Temnothorax ants?. Behav. Ecol. Sociobiol. 66, 1407–1420 (2012).Article 

Google Scholar 
Ishii, Y. & Hasgeawa, E. The mechanism underlying the regulation of work-related behaviors in the monomorphic ant, Myrmica kotokui. J. Ethol. 31, 61–69 (2013).Article 

Google Scholar 
Baudier, K. M. et al. Changing of the guard: Mixed specialization and flexibility in nest defense (Tetragonisca angustula). Behav. Ecol. 30, 1041–1049 (2019).Article 

Google Scholar 
Schmid-Hempel, P. & Schmid-Hempel, R. Life duration and turnover of foragers in the antcataglyphis bicolor (hymenoptera, formicidae). Insectes Soc. 31, 345–360 (1984).Article 

Google Scholar 
O’Donnell, S. & Jeanne, R. L. Lifelong patterns of forager behaviour in a tropical swarm-founding wasp: Effects of specialization and activity level on longevity. Anim. Behav. 44, 1021–1027 (1992).Article 

Google Scholar 
Calabi, P. Behavioral flexibility in Hymenoptera: a re-examination of the concept of caste. In Advances in Myrmecology (ed. J. C. Trager) 237–258 (Leiden,1988).Gordon, D. M. Dynamics of task switching in harvester ants. Anim. Behav. 38, 194–204 (1989).Article 

Google Scholar 
Giray, T. & Robinson, G. E. Effects of intracolony variability in behavioral development on plasticity of division of labor in honey bee colonies. Behav. Ecol. Sociobiol. 35, 13–20 (1994).Article 

Google Scholar 
Cartar, R. V. Adjustment of foraging effort and task switching in energy-manipulated wild bumblebee colonies. Anim. Behav. 44, 75–87 (1992).Article 

Google Scholar 
Huang, Z. Y. & Robinson, G. E. Regulation of honey bee division of labor by colony age demography. Behav. Ecol. Sociobiol. 39, 147–158 (1996).Article 

Google Scholar 
Gordon, D. M. The dynamics of the daily round of the harvester ant colony (Pogonomyrmex barbatus). Anim. Behav. 34, 1402–1419 (1986).Article 

Google Scholar 
Wilson, E. O. Caste and division of labor in leaf-cutter ants (Hymenoptera: Formicidae: Atta): III. Ergonomic resiliency in foraging by A. cephalotes. Behav. Ecol. Sociobiol. 14, 47–54 (1983).Article 

Google Scholar 
Middleton, E. J. T. & Latty, T. Resilience in social insect infrastructure systems. J. R. Soc. Interface 13, 20151022 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Nalepa, C. A. Origin of termite eusociality: Trophallaxis integrates the social, nutritional, and microbial environments. Ecol. Entomol. 40, 323–335 (2015).Article 

Google Scholar 
McMahan, E. A. Mound repair and foraging polyethism in workers of Nasutitermes exitiosus (Hill):(Isoptera: Termitidae). Insectes Soc. 24, 225–232 (1977).Article 

Google Scholar 
Watson, J. A. L. & McMahan, E. A. Polyethism in the Australian harvester Termite Drepanotermes (Isoptera, Termitinae). Insectes Soc. 25, 53–62 (1978).Article 

Google Scholar 
Du, H., Chouvenc, T. & Su, N.-Y. Development of age polyethism with colony maturity in Coptotermes formosanus (Isoptera: Rhinotermitidae). Environ. Entomol. 46, 311–318 (2017).PubMed 

Google Scholar 
Gerber, C., Badertscher, S. & Leuthold, R. H. Polyethism in Macrotermes bellicosus (Isoptera). Insectes Soc. 35, 226–240 (1988).Article 

Google Scholar 
Rosengaus, R. B. & Traniello, J. F. A. Temporal polyethism in incipient colonies of the primitive termite Zootermopsis angusticollis: A single multiage caste. J. Insect Behav. 6, 237–252 (1993).Article 

Google Scholar 
Crosland, M. W. J., Lok, C. M., Wong, T. C., Shakarad, M. & Traniello, J. F. A. Division of labour in a lower termite: The majority of tasks are performed by older workers. Anim. Behav. 54, 999–1012 (1997).CAS 
PubMed 
Article 

Google Scholar 
Miura, T. & Matsumoto, T. Foraging organization of the open-air processional lichen-feeding termite Hospitalitermes (Isoptera, Termitidae) in Borneo. Insectes Soc. 45, 17–32 (1998).Article 

Google Scholar 
Hinze, B. & Leuthold, R. H. Age related polyethism and activity rhythms in the nest of the termite Macrotermes bellicosus (Isoptera, Termitidae). Insectes Soc. 46, 392–397 (1999).Article 

Google Scholar 
Konate, S., Leuthold, R., Hari, M. & Veivers, P. Colour variation and polyethism of the soldier caste in the termite Macrotermes bellicosus. Entomol. Exp. Appl. 94, 51–55 (2000).Article 

Google Scholar 
Yang, R.-L., Su, N.-Y. & Bardunias, P. Individual task load in tunnel excavation by the Formosan subterranean termite (Isoptera: Rhinotermitidae). Ann. Entomol. Soc. Am. 102, 906–910 (2009).Article 

Google Scholar 
Yanagihara, S., Suehiro, W., Mitaka, Y. & Matsuura, K. Age-based soldier polyethism: Old termite soldiers take more risks than young soldiers. Biol. Lett. 14, 20180025 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Su, N. & Scheffrahn, R. H. Foraging population and territory of the Formosan subterranean termite (Isoptera, Rhinotermitidae) in an urban-environment. Sociobiology 14, 353–360 (1988).
Google Scholar 
King, E. G. & Spink, W. T. Foraging galleries of the Formosan subterranean termite, Coptotermes formosanus, in Louisiana. Ann. Entomol. Soc. Am. 62, 536–542 (1969).Article 

Google Scholar 
Abe, T. Evolution of life types in termites. In Evolution and Coadaptation in Biotic Communities (eds. J.H. Connell and J. Hidaka) 125-148 (University of Tokyo Press, 1987)Shellman-Reeve, J. S. The Spectrum of Eusociality in Termites. The Evolution of Social Behavior in Insects and Arachnids (Cambridge University Press, 1997).
Google Scholar 
Legendre, F. et al. The phylogeny of termites (Dictyoptera: Isoptera) based on mitochondrial and nuclear markers: Implications for the evolution of the worker and pseudergate castes, and foraging behaviors. Mol. Phylogenet. Evol. 48, 615–627 (2008).CAS 
PubMed 
Article 

Google Scholar 
Du, H., Chouvenc, T., Osbrink, W. L. A. & Su, N. Y. Heterogeneous distribution of castes/instars and behaviors in the nest of Coptotermes formosanus Shiraki. Insectes Soc. 64, 103–112 (2017).Article 

Google Scholar 
Su, N. Y., Osbrink, W., Kakkar, G., Mullins, A. & Chouvenc, T. Foraging distance and population size of juvenile colonies of the Formosan subterranean termite (Isoptera: Rhinotermitidae) in laboratory extended arenas. J. Econ. Entomol. 110, 1728–1735 (2017).PubMed 
Article 

Google Scholar 
Osbrink, W. L. A., Cornelius, M. L. & Lax, A. R. Effects of flooding on field populations of Formosan subterranean termites (Isoptera: Rhinotermitidae) in New Orleans, Louisiana. J. Econ. Entomol. 101, 1367–1372 (2008).PubMed 
Article 

Google Scholar 
Tuma, J., Eggleton, P. & Fayle, T. M. Ant-termite interactions: An important but under-explored ecological linkage. Biol. Rev. 95, 555–572 (2020).PubMed 
Article 

Google Scholar 
Rust, M. K. & Su, N.-Y. Managing social insects of urban importance. Annu. Rev. Entomol. 57, 355–375 (2012).CAS 
PubMed 
Article 

Google Scholar 
Evans, T. A., Forschler, B. T. & Grace, J. K. Biology of invasive termites: A worldwide review. Annu. Rev. Entomol. 58, 455–474 (2013).CAS 
PubMed 
Article 

Google Scholar 
Beverly, B. D., McLendon, H., Nacu, S., Holmes, S. & Gordon, D. M. How site fidelity leads to individual differences in the foraging activity of harvester ants. Behav. Ecol. 20, 633–638 (2009).Article 

Google Scholar 
Tenczar, P., Lutz, C. C., Rao, V. D., Goldenfeld, N. & Robinson, G. E. Automated monitoring reveals extreme interindividual variation and plasticity in honeybee foraging activity levels. Anim. Behav. 95, 41–48 (2014).Article 

Google Scholar 
O’Donnell, S. Effects of experimental forager removals on division of labour in the primitively eusocial wasp Polistes instabilis (Hymenoptera: Vespidae). Behaviour 135, 173–193 (1998).Article 

Google Scholar 
Crall, J. D. et al. Spatial fidelity of workers predicts collective response to disturbance in a social insect. Nat. Commun. 9, 1–13 (2018).Article 

Google Scholar 
Charbonneau, D. & Dornhaus, A. When doing nothing is something. How task allocation strategies compromise between flexibility, efficiency, and inactive agents. J. Bioeconomics 17, 217–242 (2015).Article 

Google Scholar 
Gordon, D. M. The regulation of foraging activity in red harvester ant colonies. Am. Nat. 159, 509–518 (2002).PubMed 
Article 

Google Scholar 
O’Donnell, S. Polybia wasp biting interactions recruit foragers following experimental worker removals. Anim. Behav. 71, 709–715 (2006).Article 

Google Scholar 
Gentry, J. B. Response to predation by colonies of the Florida harvester ant, Pogonomyrmex badius. Ecology 55, 1328–1338 (1974).Article 

Google Scholar 
Schafer, R. J., Holmes, S. & Gordon, D. M. Forager activation and food availability in harvester ants. Anim. Behav. 71, 815–822 (2006).Article 

Google Scholar 
Tschinkel, W. R. Biomantling and bioturbation by colonies of the Florida harvester ant, Pogonomyrmex badius. PLoS ONE 10, e0120407 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
Kwapich, C. L. & Tschinkel, W. R. Demography, demand, death, and the seasonal allocation of labor in the Florida harvester ant (Pogonomyrmex badius). Behav. Ecol. Sociobiol. 67, 2011–2027 (2013).Article 

Google Scholar 
Perry, C. J., Søvik, E., Myerscough, M. R. & Barron, A. B. Rapid behavioral maturation accelerates failure of stressed honey bee colonies. Proc. Natl. Acad. Sci. 112, 3427–3432 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Vance, J. T., Williams, J. B., Elekonich, M. M. & Roberts, S. P. The effects of age and behavioral development on honey bee (Apis mellifera) flight performance. J. Exp. Biol. 212, 2604–2611 (2009).PubMed 
PubMed Central 
Article 

Google Scholar 
Nalepa, C. A. Body size and termite evolution. Evol. Biol. 38, 243–257 (2011).Article 

Google Scholar 
Chouvenc, T. & Su, N. Y. Colony age-dependent pathway in caste development of Coptotermes formosanus Shiraki. Insectes Soc. 61, 171–182 (2014).Article 

Google Scholar 
Robinson, G. E., Page, R. E. Jr. & Huang, Z. Y. Temporal polyethism in social insects is a developmental process. Anim. Behav. 48, 467–469 (1994).Article 

Google Scholar 
Kakkar, G., Chouvenc, T., Osbrink, W. & Su, N. Y. Temporal assessment of molting in workers of Formosan subterranean termites (Isoptera: Rhinotermitidae). J. Econ. Entomol. 109, 2175–2181 (2016).CAS 
PubMed 
Article 

Google Scholar 
Kakkar, G., Osbrink, W., Mullins, A. & Su, N. Y. Molting site fidelity in workers of Formosan subterranean termites (Isoptera: Rhinotermitidae). J. Econ. Entomol. https://doi.org/10.1093/jee/tox246 (2017).Article 
PubMed 

Google Scholar 
Raina, A., Park, Y. I. & Gelman, D. Molting in workers of the Formosan subterranean termite Coptotermes formosanus. J. Insect Physiol. 54, 155–161 (2008).CAS 
PubMed 
Article 

Google Scholar 
Lee, S.-B., Chouvenc, T. & Su, N.-Y. Differential time allocation of foraging workers in the subterranean termite. Front. Zool. 18, 1–8 (2021).Article 

Google Scholar 
Lee, S.-B., Chouvenc, T. & Su, N.-Y. A reproductives excluder for subterranean termites in laboratory experiments. J. Econ. Entomol. 112, 2882–2887 (2019).PubMed 
Article 

Google Scholar 
Team, R. C. R: A language and environment for statistical computing. (2022). More

Growth in harvests of crops meant for exports, processing and industrial use, together with their higher yields and faster yield gains, stands out globally; at a more granular level, this was driven by specific global regions that are getting increasingly specialized in harvesting crops for these usages.Changes in global-level harvested areasAt the global scale, we find that crops harvested for direct food utilization have the highest area and have been relatively stable over the study period (Fig. 1a). However, as the total harvested hectares have increased globally (Supplementary Table 1), this has translated into decreasing fractions of crops harvested for direct food utilization, from ~51% in the 1960s (average over 1964 to 1968) to ~37% in the 2010s (average over 2009 to 2013), with a similar reduction in feed crop harvests (Table 1). Conversely, there has been a substantial increase in crops for processing, exports and industrial use (Fig. 1a, Table 1 and Supplementary Table 1). The increase in industrial crop harvests occurred after year 2000. Around the same time, harvested hectares for exported crops ramped up and by the 2010s had surpassed those of crops harvested for feed use (Fig. 1a). Crops harvested for seed usage and losses are relatively minor, and we will not discuss them further. If the global trends observed in the past 20 years continue (Fig. 1a), by 2030, crops harvested for exports, processing and industrial use will account for ~ 23%, 17% and 8% of overall harvested hectares, whereas those for food will decrease to ~29% (Table 1).Fig. 1: Sector-based global crop-utilization trends.a–d, Observed total harvested ha (a), average yield in kcal ha−1 per year (b), average yield in protein ha−1 per year (c) and average yield in fat ha−1 per year (d) in the seven sectors of food, feed, processing, export, other uses (non-food/industrial), seed and losses from 1964 to 2013, annually, and projections to 2030 based on the past 20 years. The shading shows the 90% confidence interval for the significant linear model projections.Full size imageTable 1 Sector-based global crop-utilization changesFull size tableChanges in global-level crop yieldsWe find that crops harvested for direct food usage generally have had lower yields than all other sectors at the global scale over the time period of the study (Fig. 1b–d). This is not a new phenomenon, as crops harvested for direct food utilization have always had lower yields relative to other sectors (Supplementary Table 1). What has changed, however, is the ramping up (steeper positive slopes) of industrial, export and processing crop yields (Fig. 1b–d and Table 1). At these rates, caloric yields of industrial-use crops could increase by 28% from the 2010s to 2030 compared with 24% and 21% yield increases of crops harvested for directly consumed food and for feed use (Fig. 1b). Given that caloric yields of industrial-use crops are already substantially higher than food and feed crops (2× and 1.4×, respectively, in the 2010s), the faster caloric yield increases for industrial-use crops will widen this gap (2.1× and 1.5×, respectively). Yield measurements in other units of protein and fat show similar results (Table 1, Fig. 1c,d and Supplementary Table 1).Changes in the spatial patterns of harvested areas and productionWithin country-level information on harvested areas and productivity based on utilization categories is required for developing more locally effective agricultural policies. Over the course of the study time period 1964 to 2013 (Fig. 2a,b and Supplementary Video 1), we find changes in all continents when spatially analysed at the grid-cell level, except for most parts of Africa. Even in Africa, there are locations with fractional reductions in food crop harvests over the study period, such as parts of Angola, Ghana, Nigeria and South Africa. Within these and other countries, the exact location, magnitude and direction of the change varies from one region to the next (that is, compare Fig. 2a with Fig. 2b).Fig. 2: Sector-based spatial changes in crop harvests.a,b, The fraction of a grid cell in one of seven categories—food, feed, processing, export, other (non-food/industrial use), seed and losses—in each period, 1964–1968 (a) and 2009–2013 (b).Full size imageCrops harvested for direct food utilization have been prevalent in Asia, though much has changed since the 1960s (Fig. 2a,b and Supplementary Video 1). In China, there appears to be an imaginary belt, north and west of which harvests of crops used as directly consumed food decreased between the 1960s (Fig. 2a) and 2010s (Fig. 2b), while those for other uses increased. This belt appears to roughly extend from the northern half of Jiangsu (a province on the Yellow Sea in the east), curving westwards and southwards through northern Anhui, southern Henan, central Hubei and the northern tip of Hunan, and then turning sharply south and splitting Guangdong (a province on the South China Sea) through the middle. The sector gaining from the 10–20% fractional food harvest reduction varies. The increase in crops for feed, processing and industrial usage increases as one moves northward, especially north of Jiangsu and Anhui (Fig. 2a,b and Supplementary Video 1).Similarly, in India, there is a north–south zone encompassing eastern Haryana in the north, moving southwards through eastern Rajasthan, western Madhya Pradesh to eastern Maharashtra in the south, where there was a drastic reduction in crops harvested for direct food utilization over the study period (Fig. 2 and Supplementary Video 1); crops harvested for processing primarily increased. Changes in South and Southeast Asia over the study period are primarily away from once-dominant harvests of directly consumed food crops to feed crops, followed by processing crops, export crops and industrial-use crops, as in Myanmar and Thailand. In Malaysia, the growth was in export and industrial-usage crops, whereas in Indonesia, it was export crops and smaller increases in industrial-utilization crops. Central Asian states, especially Kazakhstan and some parts of Russia, witnessed a large reduction in crops harvested for direct food use over the study period, replaced by the crops destined for exports between the two periods (Fig. 2 and Supplementary Video 1).In Australia in the 1960s, food crops were harvested everywhere, accounting for ~10% of the total, which declined to ~5% by the 2010s. This was accompanied by small reductions in crops harvested for feed and export and balanced mainly by increases in crops for processing and industrial utilization (Fig. 2 and Supplementary Video 1).In Europe in the 1960s, crops were dominantly harvested for food and feed, but by the 2010s, this changed to include crops harvested for processing (Fig. 2 and Supplementary Video 1). In France, major reductions in feed crops have been balanced by growth in processing, export and industrial-use crops. In Spain, the primary change is from crops harvested for direct food to those of feed. In Germany, crops harvested for export have replaced those for direct food utilization.Latin America used to dominantly harvest food crops (as in Mexico) or food and feed crops (as in Brazil and Argentina) (Fig. 2 and Supplementary Video 1). Midwestern Brazil used to harvest only food crops, and feed and processing crop harvests were restricted to the Atlantic states (the 1960s; Fig. 2a), but by the 2010s (Fig. 2b), harvests of food crops had become a negligible fraction in Midwestern Brazil (as in Mato Grosso), and crops harvested for processing and exports are dominant now. In the Atlantic states of Brazil, one of the major changes is the increased proportion of harvests for industrial crops. In Argentina, over the study period, the proportion of crops harvested for food and feed has decreased, and this utilization has been mainly replaced by crops harvested for processing; crops harvested for exports changed, but the direction of change was spatially heterogeneous across Argentina (Fig. 2 and Supplementary Video 1). In Mexico, the primary change is the reduction in the fraction of crops harvested for direct food consumption and the increased harvests of crops for feed.Crops harvested for food and feed are also on the decline proportionally in North America. The United States has experienced a change from the dominance of food and feed crops in the 1960s to processing and industrial-usage crops in the 2010s. Detailed changes in the United States and Canada vary from one location to the next (Fig. 2), though the major change is the lower fraction of crops harvested for direct food consumption.Results are similar when viewed through the lens of calories, protein and fat with local-level differences as yields vary based on the measurement units (Supplementary Fig. 1). Further dramatic changes can be expected if observed linear trends from 1994 to 2013 at each grid cell continued until 2030 (Supplementary Fig. 2).Calories harvested in 2030 and achieving UN SDG 2We compare the extra food calories that will potentially be harvested in 2030 (Fig. 3a and Supplementary Data 2) to those required for both the projected extra population and feeding the projected undernourished population in each country (Fig. 3b and Supplementary Data 2). As an extreme case, we also compared whether total calories (all seven utilization sectors) would be sufficient (Fig. 3c and Supplementary Data 2). Altogether, we evaluated 156 countries, of which 86 had reported undernourished populations (Supplementary Data 2). On the basis of the minimum dietary energy requirement (MDER), we find that countries with reported undernourished populations will have a shortfall of ~675.4 trillion kcal per year to nourish the increased population and the expected undernourished from their extra harvested food calories. However, compared with the more realistic average dietary energy requirement (ADER), this shortfall will be ~993.9 trillion kcal per year (or ~70% from requirements) in 2030 (15 additional scenarios of undernourished populations in 2030 (provided in Supplementary Data 3) show global calorie shortfalls may similarly range from ~587.2 trillion kcal per year to ~1,269.3 trillion kcal per year based on the MDER level of nutrition requirement, and ~880.7 trillion kcal per year to ~1,755.6 trillion kcal per year based on the more realistic ADER level of nutrition requirement in 2030).Fig. 3: Meeting UN SDG goal 2 in 2030.a, Same as Fig. 2 but for the projected kcal ha−1 per year in 2030 per utilization sector and then mapping the fraction of total kcal ha−1 per year projected as harvested. b, Shortfall or gap from kcal per year harvested in 2030 as crops for direct food use and those to plug the gap from population growth and/or undernourished population. Computed based on the 2018 to 2020 ADER number for the country. c, Same as b but the kcal per year harvested used for computation is the total across all the seven sectors and shortfall is from whether the total calories harvested were used for direct food consumption (little to no processing).Full size imageCountries reporting undernourishment can, however, meet their requirement of extra calories in 2030 for both population change and those for the undernourished if calories from other utilization sectors are diverted and consumed directly as food calories (Fig. 3c and Supplementary Data 2 and 3). Though at the global scale, it appears that countries with high levels of undernourishment in 2030 can divert just a portion of their total harvested calories and meet some of the requirements of UN’s SDG 2 (ref. 4). In reality, many of the individual countries concentrated in sub-Saharan Africa have limited scope of diversion of calories from other sectors such as feed, processing or exports as crops for direct food use, as they already harvest most crops for direct food consumption (Fig. 2 and Supplementary Figs. 1 and 2). As such, many countries in this region may see deepening reliance on food imports. Note that the UN’s second SDG goal is broader in scope, including efforts to end malnutrition and increase agricultural productivity, among other goals4. Reconfiguration planning19 can use our spatially detailed information (Figs. 2 and 3, Supplementary Figs. 1 and 2 and Supplementary Data 2 and 3) in conjunction with policies that incentivize increased food crop harvests globally and ensure their equitable distribution to undernourished regions when local production is not sufficient20,21. This will require supply chain management22,23 and detailed analysis of optimization scenarios24 with our maps and tables as an important step linking specific production regions with the initial use of that production. More

LookingGlass design and optimizationDataset generationThe taxonomic organization of representative Bacterial and Archaeal genomes was determined from the Genome Taxonomy Database, GTDB51 (release 89.0). The complete genome sequences were downloaded via the NCBI Genbank ftp52. This resulted in 24,706 genomes, comprising 23,458 Bacterial and 1248 Archaeal genomes.Each genome was split into read-length chunks. To determine the distribution of realistic read lengths produced by next-generation short-read sequencing machines, we obtained the BioSample IDs52 for each genome, where they existed, and downloaded their sequencing metadata from the MetaSeek53 database using the MetaSeek API. We excluded samples with average read lengths less than 60 or greater than 300 base pairs. This procedure resulted in 7909 BioSample IDs. The average read lengths for these sequencing samples produced the read-length distribution (Supplementary Fig. 1) with a mean read length of 136 bp. Each genome was split into read-length chunks (with zero overlap in order to maximize information density and reduce data redundancy in the dataset): a sequence length was randomly selected with replacement from the read-length distribution and a sequence fragment of that length was subset from the genome, with a 50% chance that the reverse complement was used. The next sequence fragment was chosen from the genome starting at the end point of the previous read-length chunk, using a new randomly selected read length, and so on. These data were partitioned into a training set used for optimization of the model; a validation set used to evaluate model performance during parameter tuning and as a benchmark to avoid overfitting during training; and a test set used for final evaluation of model performance. To ensure that genomes in the training, validation, and test sets had low sequence similarity, the sets were split along taxonomic branches such that genomes from the Actinomycetales, Rhodobacterales, Thermoplasmata, and Bathyarchaeia were partitioned into the validation set; genomes from the Bacteroidales, Rhizobiales, Methanosarcinales, and Nitrososphaerales were partitioned into the test set; and the remaining genomes remained in the training set. This resulted in 529,578,444 sequences in the training set, 57,977,217 sequences in the validation set, and 66,185,518 sequences in the test set. We term this set of reads the GTDB representative set (Table 1).Table 1 Summary table of datasets used.Full size tableThe amount of data needed for training was also evaluated (Supplementary Fig. 2). Progressively larger amounts of data were tested by selecting at random 1, 10, 100, or 500 read-length chunks from each of the GTDB representative genomes in the GTDB representative training set. Additionally, the performance of smaller but more carefully selected datasets, representing the diversity of the microbial tree of life, were tested by selecting for training one genome at random from each taxonomic class or order in the GTDB taxonomy tree. In general, better accuracy was achieved in fewer epochs with a greater amount of sequencing data (Supplementary Fig. 2); however, a much smaller amount of data performed better if a representative genome was selected from each GTDB taxonomy class.The final LookingGlass model was trained on this class-level partition of the microbial tree of life. We term this dataset the GTDB class set (Table 1). The training, validation, and test sets were split such that no classes overlapped across sets: the validation set included 8 genomes from each of the classes Actinobacteria, Alphaproteobacteria, Thermoplasmata, and Bathyarchaeia (32 total genomes); the test set included 8 genomes from each of the classes Bacteroidia, Clostridia, Methanosarcinia, and Nitrososphaeria (32 total genomes); and the training set included 1 genome from each of the remaining classes (32 archaeal genomes and 298 bacterial genomes for a total of 330 genomes). This resulted in a total of 6,641,723 read-length sequences in the training set, 949,511 in the validation set, and 632,388 in the test set (Supplementary Data 1).Architecture design and trainingRecurrent neural networks (RNNs) are a type of neural network designed to take advantage of the context dependence of sequential data (such as text, video, audio, or biological sequences), by passing information from previous items in a sequence to the current item in a sequence54. Long short-term memory networks (LSTMs)55 are an extension of RNNs, which better learn long-term dependencies by handling the RNN tendency to “forget” information farther away in a sequence56. LSTMs maintain a cell state which contains the “memory” of the information in the previous items in the sequence. LSTMs learn additional parameters which decide at each step in the sequence which information in the cell state to “forget” or “update”.LookingGlass uses a three-layer LSTM encoder model with 1152 units in each hidden layer and an embedding size of 104 based on the results of hyperparameter tuning (see below). It divides the sequence into characters using a kmer size of 1 and a stride of 1, i.e., is a character-level language model. LookingGlass is trained in a self-supervised manner to predict a masked nucleotide, given the context of the preceding nucleotides in the sequence. For each read in the training sequence, multiple training inputs are considered, shifting the nucleotide that is masked along the length of the sequence from the second position to the final position in the sequence. Because it is a character-level model, a linear decoder predicts the next nucleotide in the sequence from the possible vocabulary items “A”, “C”, “G”, and “T”, with special tokens for “beginning of read”, “unknown nucleotide” (for the case of ambiguous sequences), “end of read” (only “beginning of read” was tokenized during LookingGlass training), and a “padding” token (used for classification only).Regularization and optimization of LSTMs require special approaches to dropout and gradient descent for best performance57. The fastai library58 offers default implementations of these approaches for natural language text, and so we adopt the fastai library for all training presented in this paper. We provide the open source fastBio python package59 which extends the fastai library for use with biological sequences.LookingGlass was trained on a Pascal P100 GPU with 16GB memory on Microsoft Azure, using a batch size of 512, a back propagation through time (bptt) window of 100 base pairs, the Adam optimizer60, and utilizing a Cross Entropy loss function (Supplementary Table 1). Dropout was applied at variable rates across the model (Supplementary Table 1). LookingGlass was trained for a total of 12 days for 75 epochs, with progressively decreasing learning rates based on the results of hyperparameter optimization (see below): for 15 epochs at a learning rate of 1e−2, for 15 epochs at a learning rate of 2e−3, and for 45 epochs at a learning rate of 1e−3.Hyperparameter optimizationHyperparameters used for the final training of LookingGlass were tuned using a randomized search of hyperparameter settings. The tuned hyperparameters included kmer size, stride, number of LSTM layers, number of hidden nodes per layer, dropout rate, weight decay, momentum, embedding size, bptt size, learning rate, and batch size. An abbreviated dataset consisting of ten randomly selected read-length chunks from the GTDB representative set was created for testing many parameter settings rapidly. A language model was trained for two epochs for each randomly selected hyperparameter combination, and those conditions with the maximum performance were accepted. The hyperparameter combinations tested and the selected settings are described in the associated Github repository61.LookingGlass validation and analysis of embeddingsFunctional relevanceDataset generation
In order to assess the ability of the LookingGlass embeddings to inform the molecular function of sequences, metagenomic sequences from a diverse set of environments were downloaded from the Sequence Read Archive (SRA)62. We used MetaSeek53 to choose ten metagenomes at random from each of the environmental packages defined by the MIxS metadata standards63: built environment, host-associated, human gut, microbial mat/biofilm, miscellaneous, plant-associated, sediment, soil, wastewater/sludge, and water, for a total of 100 metagenomes. The SRA IDs used are available in (Supplementary Table 2). The raw DNA reads for these 100 metagenomes were downloaded from the SRA with the NCBI e-utilities. These 100 metagenomes were annotated with the mi-faser tool27 with the read-map option to generate predicted functional annotation labels (to the fourth digit of the Enzyme Commission (EC) number), out of 1247 possible EC labels, for each annotatable read in each metagenome. These reads were then split 80%/20% into training/validation candidate sets of reads. To ensure that there was minimal overlap in sequence similarity between the training and validation set, we compared the validation candidate sets of each EC annotation to the training set for that EC number with CD-HIT64, and filtered out any reads with >80% DNA sequence similarity to the reads of that EC number in the training set (the minimum CD-HIT DNA sequence similarity cutoff). In order to balance EC classes in the training set, overrepresented ECs in the training set were downsampled to the mean count of read annotations (52,353 reads) before filtering with CD-HIT. After CD-HIT processing, any underrepresented EC numbers in the training set were oversampled to the mean count of read annotations (52,353 reads). The validation set was left unbalanced to retain a distribution more realistic to environmental settings. The final training set contained 61,378,672 reads, while the validation set contained 2,706,869 reads. We term this set of reads and their annotations the mi-faser functional set (Table 1).
As an external test set, we used a smaller number of DNA sequences from genes with experimentally validated molecular functions. We linked the manually curated entries of Bacterial or Archaeal proteins from the Swiss-Prot database65 corresponding to the 1247 EC labels in the mi-faser functional set with their corresponding genes in the EMBL database66. We downloaded the DNA sequences, and selected ten read-length chunks at random per CDS. This resulted in 1,414,342 read-length sequences in the test set. We term this set of reads and their annotations the Swiss-Prot functional set (Table 1).

Fine-tuning procedure
We fine-tuned the LookingGlass language model to predict the functional annotation of DNA reads, to demonstrate the speed with which an accurate model can be trained using our pretrained LookingGlass language model. The architecture of the model retained the 3-layer LSTM encoder and the weights of the LookingGlass language model encoder, but replaced the language model decoder with a new multiclass classification layer with pooling (with randomly initialized weights). This pooling classification layer is a sequential model consisting of the following layers: a layer concatenating the output of the LookingGlass encoder with min, max, and average pooling of the outputs (for a total dimension of 104*3 = 312), a batch normalization67 layer with dropout, a linear layer taking the 312-dimensional output of the batch norm layer and producing a 50-dimensional output, another batch normalization layer with dropout, and finally a linear classification layer that is passed through the log(Softmax(x)) function to output the predicted functional annotation of a read as a probability distribution of the 1247 possible mi-faser EC annotation labels. We then trained the functional classifier on the mi-faser functional set described above. Because the >61 million reads in the training set were too many to fit into memory, training was done in 13 chunks of ~5-million reads each until one total epoch was completed. Hyperparameter settings for the functional classifier training are seen in Supplementary Table 1.

Encoder embeddings and MANOVA test
To test whether the LookingGlass language model embeddings (before fine-tuning, above) are distinct across functional annotations, a random subset of ten reads per functional annotation was selected from each of the 100 SRA metagenomes (or the maximum number of reads present in that metagenome for that annotation, whichever was greater). This also ensured that reads were evenly distributed across environments. The corresponding fixed-length embedding vectors for each read was produced by saving the output from the LookingGlass encoder (before the embedding vector is passed to the language model decoder) for the final nucleotide in the sequence. This vector represents a contextually relevant embedding for the overall sequence. The statistical significance of the difference between embedding vectors across all functional annotation groups was tested with a MANOVA test using the R stats package68.
Evolutionary relevance
Dataset generation
The OrthoDB database69 provides orthologous groups (OGs) of proteins at various levels of taxonomic distance. For instance, the OrthoDB group “77at2284” corresponds to proteins belonging to “Glucan 1,3-alpha-glucosidase at the Sulfolobus level”, where “2284” is the NCBI taxonomy ID for the genus Sulfolobus.
We tested whether embedding similarity of homologous sequences (sequences within the same OG) is higher than that of nonhomologous sequences (sequences from different OGs). We tested this in OGs at multiple levels of taxonomic distance—genus, family, order, class, and phylum. At each taxonomic level, ten individual taxa at that level were chosen from across the prokaryotic tree of life (e.g., for the genus level, Acinetobacter, Enterococcus, Methanosarcina, Pseudomonas, Sulfolobus, Bacillus, Lactobacillus, Mycobacterium, Streptomyces, and Thermococcus were chosen). For each taxon, 1000 randomly selected OGs corresponding to that taxon were chosen; for each of these OGs, five randomly chosen genes within this OG were chosen.
OrthoDB cross-references OGs to UniProt65 IDs of the corresponding proteins. We mapped these to the corresponding EMBL CDS IDs66 via the UniProt database API65; DNA sequences of these EMBL CDSs were downloaded via the EMBL database API. For each of these sequences, we generated LookingGlass embedding vectors.

Homologous and nonhomologous sequence pairs
To create a balanced dataset of homologous and nonhomologous sequence pairs, we compared all homologous pairs of the five sequences in an OG (total of ten homologous pairs) to an equal number of randomly selected out-of-OG comparisons for the same sequences; i.e., each of the five OG sequences was compared to 2 other randomly selected sequences from any other randomly selected OG (total of ten nonhomologous pairs). We term this set of sequences, and their corresponding LookingGlass embeddings, the OG homolog set (Table 1).

Embedding and sequence similarity
For each sequence pair, the sequence and embedding similarity were determined. The embedding similarity was calculated as the cosine similarity between embedding vectors. The sequence similarity was calculated as the Smith-Waterman alignment score using the BioPython70 pairwise2 package, with a gap open penalty of −10 and a gap extension penalty of −1. The IDs of chosen OGs, the cosine similarities of the embedding vectors, and sequence similarities of the DNA sequences are available in the associated Github repository61.

Comparison to HMM-based domain searches for distant homology detection
Distantly related homologous sequences that share, e.g., Pfam71, domains can be identified using HMM-based search methods. We used hmmscan25 (e-val threshold = 1e−10) to compare homologous (at the phylum level) sequences in the OG homolog set, for which the alignment score was less than 50 bits and the embedding similarity was greater than 0.62 (total: 21,376 gene pairs). Specifically, we identified Pfam domains in each sequence and compared whether the most significant (lowest e-value) domain for each sequence was identified in common for each homologous pair.
Environmental relevance
Encoder embeddings and MANOVA test
The LookingGlass embeddings and the environment of origin for each read in the mi-faser functional set were used to test the significance of the difference between the embedding vectors across environmental contexts. The statistical significance of this difference was evaluated with a MANOVA test using the R stats package68.
Oxidoreductase classifier
Dataset generation
The manually curated, reviewed entries of the Swiss-Prot database65 were downloaded (June 2, 2020). Of these, 23,653 entries were oxidoreductases (EC number 1.-.-.-) of Archaeal or Bacterial origin (988 unique ECs). We mapped their UniProt IDs to both their EMBL CDS IDs and their UniRef50 IDs via the UniProt database mapper API. Uniref50 IDs identify clusters of sequences with >50% amino acid identity. This cross-reference identified 28,149 EMBL CDS IDs corresponding to prokaryotic oxidoreductases, belonging to 5451 unique UniRef50 clusters. We split this data into training, validation, and test sets such that each UniRef50 cluster was contained in only one of the sets, i.e., there was no overlap in EMBL CDS IDs corresponding to the same UniRef50 cluster across sets. This ensures that the oxidoreductase sequences in the validation and test sets are dissimilar to those seen during training. The DNA sequences for each EMBL CDS ID were downloaded via the EMBL database API. These data generation process were repeated for a random selection of non-oxidoreductase UniRef50 clusters, which resulted in 28,149 non-oxidoreductase EMBL CDS IDs from 13,248 unique UniRef50 clusters.
Approximately 50 nucleotide read-length chunks (selected from the representative read-length distribution, as above) were selected from each EMBL CDS DNA sequence, with randomly selected start positions on the gene and a 50% chance of selecting the reverse complement, such that an even number of read-length sequences with “oxidoreductase” and “not oxidoreductase” labels were generated for the final dataset. This procedure produced a balanced dataset with 2,372,200 read-length sequences in the training set, 279,200 sequences in the validation set, and 141,801 sequences in the test set. We term this set of reads and their annotations the oxidoreductase model set (Table 1). In order to compare the oxidoreductase classifier performance to an HMM-based method, reads with “oxidoreductase” labels in the oxidoreductase model test set (71,451 reads) were 6-frame translated and searched against the Swiss-Prot protein database using phmmer25 (reporting e-val threshold = 0.05, using all other defaults).

Fine-tuning procedure
Since our functional annotation classifier addresses a closer classification task to the oxidoreductase classifier than LookingGlass itself, the architecture of the oxidoreductase classifier was fine-tuned starting from the functional annotation classifier, replacing the decoder with a new pooling classification layer (as described above for the functional annotation classifier) and with a final output size of 2 to predict “oxidoreductase” or “not oxidoreductase”. Fine tuning of the oxidoreductase classifier layers was done successively, training later layers in isolation and then progressively including earlier layers into training, using discriminative learning rates ranging from 1e−2 to 5e−4, as previously described72. The fine-tuned model was trained for 30 epochs, over 18 h, on a single P100 GPU node with 16GB memory.

Model performance in metagenomes
Sixteen marine metagenomes from the surface (SRF, ~5 meters) and mesopelagic (MES, 175–800 meters) from eight stations sampled as part of the TARA expedition37 were downloaded from the SRA62 (Supplementary Table 3, SRA accession numbers ERR598981, ERR599063, ERR599115, ERR599052, ERR599020, ERR599039, ERR599076, ERR598989, ERR599048, ERR599105, ERR598964, ERR598963, ERR599125, ERR599176, ERR3589593, and ERR3589586). Metagenomes were chosen from a latitudinal gradient spanning polar, temperate, and tropical regions and ranging from −62 to 76 degrees latitude. Mesopelagic depths from four out of the eight stations were sampled from oxygen minimum zones (OMZs, where oxygen More

Smale, D. A. et al. Marine heatwaves threaten global biodiversity and the provision of ecosystem services. Nat. Clim. Change https://doi.org/10.1038/s41558-019-0412-1 (2019).Article 

Google Scholar 
Smith, K. E. et al. Socioeconomic impacts of marine heatwaves: Global issues and opportunities. Science 374, eabj3593 (2021).CAS 
Article 

Google Scholar 
Frolicher, T. L., Fischer, E. M. & Gruber, N. Marine heatwaves under global warming. Nature 560, 360–364. https://doi.org/10.1038/s41586-018-0383-9 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Oliver, et al. Longer and more frequent marine heatwaves over the past century. Nat. Commun. https://doi.org/10.1038/s41467-018-03732-9 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Oliver, et al. Projected marine heatwaves in the 21st century and the potential for ecological impact. Front. Mar. Sci. https://doi.org/10.3389/fmars.2019.00734 (2019).Article 

Google Scholar 
Hobday, A. J. et al. A hierarchical approach to defining marine heatwaves. Prog. Oceanogr. 141, 227–238 (2016).ADS 
Article 

Google Scholar 
Banzon, V., Smith, T. M., Chin, T. M., Liu, C. & Hankins, W. A long-term record of blended satellite and in situ sea-surface temperature for climate monitoring, modeling and environmental studies. Earth Syst. Sci. Data 8, 165–176 (2016).ADS 
Article 

Google Scholar 
Wernberg, T. et al. An extreme climatic event alters marine ecosystem structure in a global biodiversity hotspot. Nat. Clim. Change 3, 78–82. https://doi.org/10.1038/nclimate1627 (2013).ADS 
Article 

Google Scholar 
Arias-Ortiz, A. et al. A marine heatwave drives massive losses from the world’s largest seagrass carbon stocks. Nat. Clim. Change https://doi.org/10.1038/s41558-018-0096-y (2018).Article 

Google Scholar 
Smale, D. A. Impacts of ocean warming on kelp forest ecosystems. New Phytol. 225, 1447–1454 (2020).Article 

Google Scholar 
Couch, C. S. et al. Mass coral bleaching due to unprecedented marine heatwave in Papahānaumokuākea Marine National Monument (Northwestern Hawaiian Islands). PLoS ONE 12, e0185121. https://doi.org/10.1371/journal.pone.0185121 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Oliver, E. C. J. et al. The unprecedented 2015/16 Tasman Sea marine heatwave. Nat. Commun. https://doi.org/10.1038/ncomms16101 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Montie, S., Thomsen, M. S., Rack, W. & Broady, P. A. Extreme summer marine heatwaves increase chlorophyll a in the Southern Ocean. Antarct. Sci. 32, 508–509 (2020).ADS 
Article 

Google Scholar 
Gupta, A. S. et al. Drivers and impacts of the most extreme marine heatwaves events. Sci. Rep. 10, 1–15 (2020).ADS 
Article 

Google Scholar 
Holbrook, N. J. et al. A global assessment of marine heatwaves and their drivers. Nat. Commun. 10, 1–13 (2019).CAS 
Article 

Google Scholar 
La Sorte, F. A., Johnston, A. & Ault, T. R. Global trends in the frequency and duration of temperature extremes. Clim. Change 166, 1. https://doi.org/10.1007/s10584-021-03094-0 (2021).ADS 
Article 

Google Scholar 
Thomsen, et al. Local extinction of bull kelp (Durvillaea spp.) due to a marine heatwave. Front. Mar. Sci. https://doi.org/10.3389/fmars.2019.00084 (2019).Article 

Google Scholar 
Strydom, S. et al. Too hot to handle: Unprecedented seagrass death driven by marine heatwave in a World Heritage Area. Glob. Change Biol. 26, 3525–3538. https://doi.org/10.1111/gcb.15065 (2020).ADS 
Article 

Google Scholar 
Leggat, W. P. et al. Rapid coral decay is associated with marine heatwave mortality events on reefs. Curr. Biol. 29, 2723-2730.e2724. https://doi.org/10.1016/j.cub.2019.06.077 (2019).CAS 
Article 
PubMed 

Google Scholar 
Wernberg, T. et al. Climate-driven regime shift of a temperate marine ecosystem. Science 353, 169–172. https://doi.org/10.1126/science.aad8745 (2016).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Thomsen, M. S. & McGlathery, K. Facilitation of macroalgae by the sedimentary tube forming polychaete Diopatra cuprea. Estuar. Coast. Shelf Sci. 62, 63–73. https://doi.org/10.1016/j.ecss.2004.08.007 (2005).ADS 
Article 

Google Scholar 
Spalding, M. D. et al. Marine ecoregions of the world: A bioregionalization of coastal and shelf areas. Bioscience 57, 573–583 (2007).Article 

Google Scholar 
Costello, M. J. & Chaudhary, C. Marine biodiversity, biogeography, deep-sea gradients, and conservation. Curr. Biol. 27, R511–R527. https://doi.org/10.1016/j.cub.2017.04.060 (2017).CAS 
Article 
PubMed 

Google Scholar 
Tait, L. W., Thoral, F., Pinkerton, M. H., Thomsen, M. S. & Schiel, D. R. Loss of giant kelp, Macrocystis pyrifera, driven by marine heatwaves and exacerbated by poor water clarity in New Zealand. Front. Mar. Sci. https://doi.org/10.3389/fmars.2021.721087 (2021).Article 

Google Scholar 
Marin, M., Feng, M., Phillips, H. E. & Bindoff, N. L. A global, multiproduct analysis of coastal marine heatwaves: Distribution, characteristics, and long-term trends. J. Geophys. Res. Oceans 126, e2020JC016708. https://doi.org/10.1029/2020JC016708 (2021).ADS 
Article 

Google Scholar 
Kain, J. M. The seasons in the subtidal. Brit. Phycol. J. 24, 203–215 (1989).Article 

Google Scholar 
Atkinson, J., King, N. G., Wilmes, S. B. & Moore, P. J. Summer and winter marine heatwaves favor an invasive over native seaweeds. J. Phycol. 56, 1591–1600. https://doi.org/10.1111/jpy.13051 (2020).CAS 
Article 
PubMed 

Google Scholar 
Salinger, M. J. et al. The unprecedented coupled ocean-atmosphere summer heatwave in the New Zealand region 2017/18: Drivers, mechanisms and impacts. Environ. Res. Lett. 14, 044023 (2019).ADS 
Article 

Google Scholar 
Amaya, D. J., Miller, A. J., Xie, S.-P. & Kosaka, Y. Physical drivers of the summer 2019 North Pacific marine heatwave. Nat. Commun. 11, 1–9 (2020).Article 

Google Scholar 
Di Lorenzo, E. & Mantua, N. Multi-year persistence of the 2014/15 North Pacific marine heatwave. Nat. Clim. Change 6, 1042–1047. https://doi.org/10.1038/nclimate3082 (2016).ADS 
Article 

Google Scholar 
Cayan, D. R. Large-scale relationships between sea surface temperature and surface air temperature. Mon. Weather Rev. 108, 1293–1301 (1980).ADS 
Article 

Google Scholar 
Hipel, K. W. & McLeod, A. I. Time Series Modelling of Water Resources and Environmental Systems (Elsevier, 1994).
Google Scholar 
trend: non-parametric trend tests and changepoint detection.–R package ver. 1.1. 2 (2020).Costanza, R. et al. The value of the world’s ecosystem services and natural capital. Nature 387, 253–260 (1997).ADS 
CAS 
Article 

Google Scholar 
Halpern, B. S. et al. A global map of human impact on marine ecosystems. Science 319, 948–952 (2008).ADS 
CAS 
Article 

Google Scholar 
Harley, C. D. et al. The impacts of climate change in coastal marine systems. Ecol. Lett. 9, 228–241. https://doi.org/10.1111/j.1461-0248.2005.00871.x (2006).ADS 
Article 
PubMed 

Google Scholar 
Thomsen, M. S. & South, P. M. Communities and attachment networks associated with primary, secondary and alternative foundation species; a case study of stressed and disturbed stands of southern bull kelp. Diversity 11, 56. https://doi.org/10.3390/d11040056 (2019).Article 

Google Scholar 
Smale, D. A. & Wernberg, T. Extreme climatic event drives range contraction of a habitat-forming species. Proc. R. Soc. B Biol. Sci. 280, 20122829 (2013).Article 

Google Scholar 
Thomsen, M. S. et al. Cascading impacts of earthquakes and extreme heatwaves have destroyed populations of an iconic marine foundation species. Divers. Distrib. (2021).Rogers-Bennett, L. & Catton, C. Marine heat wave and multiple stressors tip bull kelp forest to sea urchin barrens. Sci. Rep. 9, 1–9 (2019).CAS 
Article 

Google Scholar 
Filbee-Dexter, K. et al. Marine heatwaves and the collapse of marginal North Atlantic kelp forests. Sci. Rep. 10, 1–11 (2020).Article 

Google Scholar 
Thomson, J. A. et al. Extreme temperatures, foundation species, and abrupt ecosystem change: An example from an iconic seagrass ecosystem. Glob. Change Biol. 21, 1463–1474. https://doi.org/10.1111/gcb.12694 (2015).ADS 
Article 

Google Scholar 
Hughes, T. P. et al. Global warming and recurrent mass bleaching of corals. Nature 543, 373–377. https://doi.org/10.1038/nature21707 (2017).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Le Nohaïc, M. et al. Marine heatwave causes unprecedented regional mass bleaching of thermally resistant corals in northwestern Australia. Sci. Rep. 7, 14999. https://doi.org/10.1038/s41598-017-14794-y (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Smale, D. A., Wernberg, T. & Vanderklift, M. A. Regional-scale variability in the response of benthic macroinvertebrate assemblages to a marine heatwave. Mar. Ecol. Prog. Ser. 568, 17–30. https://doi.org/10.3354/meps12080 (2017).ADS 
Article 

Google Scholar 
Cavole, L. et al. Biological impacts of the 2013–2015 warm-water anomaly in the Northeast Pacific: Winners, losers, and the future. Oceanography (Washington D.C.) https://doi.org/10.5670/oceanog.2016.32 (2016).Article 

Google Scholar 
Coleman, M. A., Minne, A. J. P., Vranken, S. & Wernberg, T. Genetic tropicalisation following a marine heatwave. Sci. Rep. 10, 12726. https://doi.org/10.1038/s41598-020-69665-w (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Collette, B. B. in Reproduction and sexuality in marine fishes 21–64 (University of California Press, 2010).Yatsu, A. & Shimada, H. Distributions of Epipelagic Fishes, Squids, Marine Mammals. Bulletin 53 North Pacific Commission, 111–146.Hirst, A., Roff, J. & Lampitt, R. A synthesis of growth rates in marine epipelagic invertebrate zooplankton. Adv. Mar. Biol. 44, 1–142 (2003).CAS 
Article 

Google Scholar 
Smale, D. A. & Wernberg, T. Satellite-derived SST data as a proxy for water temperature in nearshore benthic ecology. Mar. Ecol. Prog. Ser. 387, 27–37 (2009).ADS 
Article 

Google Scholar 
Bernardello, R., Serrano, E., Coma, R., Ribes, M. & Bahamon, N. A comparison of remote-sensing SST and in situ seawater temperature in near-shore habitats in the western Mediterranean Sea. Mar. Ecol. Prog. Ser. 559, 21–34 (2016).ADS 
Article 

Google Scholar 
Brewin, R. J. et al. Evaluating operational AVHRR sea surface temperature data at the coastline using benthic temperature loggers. Remote Sens. 10, 925 (2018).ADS 
Article 

Google Scholar 
Smit, A. J. et al. A coastal seawater temperature dataset for biogeographical studies: Large biases between in situ and remotely-sensed data sets around the Coast of South Africa. PLoS ONE 8, e81944. https://doi.org/10.1371/journal.pone.0081944 (2013).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Marin, M., Bindoff, N. L., Feng, M. & Phillips, H. E. Slower long-term coastal warming drives dampened trends in coastal marine heatwave exposure. J. Geophys. Res. Oceans https://doi.org/10.1029/2021jc017930 (2021).Article 

Google Scholar 
Lourenço, C. R. et al. Upwelling areas as climate change refugia for the distribution and genetic diversity of a marine macroalga. J. Biogeogr. 43, 1595–1607. https://doi.org/10.1111/jbi.12744 (2016).Article 

Google Scholar 
Riegl, B. & Piller, W. E. Possible refugia for reefs in times of environmental stress. Int. J. Earth Sci. 92, 520–531. https://doi.org/10.1007/s00531-003-0328-9 (2003).Article 

Google Scholar 
El Glynn, P. W. Nino-southern oscillation 1982–1983: Nearshore population, community, and ecosystem responses. Annu. Rev. Ecol. Syst. 19, 309–346. https://doi.org/10.1146/annurev.es.19.110188.001521 (1988).Article 

Google Scholar 
Glynn, P. W. & D’Croz, L. Experimental evidence for high temperature stress as the cause of El Niño-coincident coral mortality. Coral Reefs 8, 181–191. https://doi.org/10.1007/bf00265009 (1990).ADS 
Article 

Google Scholar 
Glynn, P. W., Maté, J. L., Baker, A. C. & Calderón, M. O. Coral bleaching and mortality in panama and ecuador during the 1997–1998 El Niño-Southern Oscillation Event: Spatial/temporal patterns and comparisons with the 1982–1983 event. Bull. Mar. Sci. 69, 79–109 (2001).
Google Scholar 
Podestá, G. P. & Glynn, P. W. The 1997–98 El Niño event in Panama and Galápagos: An update of thermal stress indices relative to coral bleaching. Bull. Mar. Sci. 69, 43–59 (2001).
Google Scholar 
Ladah, L. B. & Zertuche-Gonzalez, J. A. Giant kelp (Macrocystis pyrifera) survival in deep water (25–40 m) during El Nino of 1997–1998 in Baja California, Mexico. Bot. Marina 47, 367–372. https://doi.org/10.1515/bot.2004.054 (2004).Article 

Google Scholar 
Kayanne, H. Validation of degree heating weeks as a coral bleaching index in the northwestern Pacific. Coral Reefs 36, 63–70 (2017).ADS 
Article 

Google Scholar 
Le Nohaïc, M. et al. Marine heatwave causes unprecedented regional mass bleaching of thermally resistant corals in northwestern Australia. Sci. Rep. 7, 1–11 (2017).ADS 
Article 

Google Scholar 
Marba, N. & Duarte, C. M. Mediterranean warming triggers seagrass (Posidonia oceanica) shoot mortality. Glob. Change Biol. 16, 2366–2375. https://doi.org/10.1111/j.1365-2486.2009.02130.x (2010).ADS 
Article 

Google Scholar 
Bennett, S., Wernberg, T., Arackal Joy, B., de Bettignies, T. & Campbell, A. H. Central and rear-edge populations can be equally vulnerable to warming. Nat. Commun. 6, 10280. https://doi.org/10.1038/ncomms10280 (2015).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Filbee-Dexter, K. et al. Marine heatwaves and the collapse of marginal North Atlantic kelp forests. Sci. Rep. 10, 13388. https://doi.org/10.1038/s41598-020-70273-x (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Snover, M. L. Ontogenetic habitat shifts in marine organisms: Influencing factors and the impact of climate variability. Bull. Mar. Sci. 83, 53–67 (2008).
Google Scholar 
Harley, C. D. Climate change, keystone predation, and biodiversity loss. Science 334, 1124–1127 (2011).ADS 
CAS 
Article 

Google Scholar 
Kelaher, B. P., Coleman, M. A. & Bishop, M. J. Ocean warming, but not acidification, accelerates seagrass decomposition under near-future climate scenarios. Mar. Ecol. Prog. Ser. 605, 103–110 (2018).ADS 
CAS 
Article 

Google Scholar 
De Senerpont Domis, L. N. et al. Plankton dynamics under different climatic conditions in space and time. Freshw. Biol. 58, 463–482 (2013).Article 

Google Scholar 
Fossheim, M. et al. Recent warming leads to a rapid borealization of fish communities in the Arctic. Nat. Clim. Change 5, 673–677. https://doi.org/10.1038/nclimate2647 (2015).ADS 
Article 

Google Scholar 
Morales-Nin, B. & Panfili, J. Seasonality in the deep sea and tropics revisited: What can otoliths tell us?. Mar. Freshw. Res. 56, 585–598 (2005).Article 

Google Scholar 
Alongi, D. M. Ecology of tropical soft-bottom benthos: A review with emphasis on emerging concepts. Rev. Biol. Trop. 37, 85–100 (1989).
Google Scholar 
Hobday, A. J., Spillman, C. M., Paige Eveson, J. & Hartog, J. R. Seasonal forecasting for decision support in marine fisheries and aquaculture. Fish. Oceanogr. 25, 45–56. https://doi.org/10.1111/fog.12083 (2016).Article 

Google Scholar 
Spillman, C. M., Smith, G. A., Hobday, A. J. & Hartog, J. R. Onset and decline rates of marine heatwaves: Global trends, seasonal forecasts and marine management. Front. Clim. https://doi.org/10.3389/fclim.2021.801217 (2021).Article 

Google Scholar 
Schlegel, R. W., Oliver, E. C. J., Wernberg, T. & Smit, A. J. Nearshore and offshore co-occurrence of marine heatwaves and cold-spells. Prog. Oceanogr. 151, 189–205. https://doi.org/10.1016/j.pocean.2017.01.004 (2017).ADS 
Article 

Google Scholar 
Huang, B. et al. Improvements of the daily optimum interpolation sea surface temperature (DOISST) Version 2.1. J. Clim. 34, 2923–2939. https://doi.org/10.1175/jcli-d-20-0166.1 (2021).ADS 
Article 

Google Scholar 
OBPG, N. & Stumpf, R. P. Distance to Nearest Coastline: 0.01-Degree Grid. Distributed by the Pacific Islands Ocean Observing System (PacIOOS). http://pacioos.org/metadata/dist2coast_1deg.html and https://data.noaa.gov/dataset/dataset/distance-to-nearest-coastline-0-01-degree-grid2http://www.pacioos.hawaii.edu/metadata/dist2coast_1deg.html (2012).Schlegel, R. W. & Smit, A. J. heatwaveR: A central algorithm for the detection of heatwaves and cold-spells. J. Open Source Softw. 3(27), 821 (2018).ADS 
Article 

Google Scholar 
Sakurai, T., Yukio, K. & Kuragano, T. in Proceedings. 2005 IEEE International Geoscience and Remote Sensing Symposium, 2005. IGARSS’05. 2606–2608 (IEEE). More

Foster, G. & Rahmstorf, S. Global temperature evolution 1979–2010. Environ. Res. Lett. 6, 044022 (2011).ADS 
Article 

Google Scholar 
Cahill, N., Rahmstorf, S. & Parnell, A. C. Change points of global temperature. Environ. Res. Lett. 10, 084002 (2015).ADS 
Article 

Google Scholar 
Yan, X. H. et al. The global warming hiatus: Slowdown or redistribution?. Earth’s Future 4, 472–482 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Karl, T. R. et al. Possible artifacts of data biases in the recent global surface warming hiatus. Science 348, 5632 (2015).Article 

Google Scholar 
Cohen, J. L., Furtado, J. C., Barlow, M., Alexeev, V. A. & Cherry, J. E. Asymmetric seasonal temperature trends. Geophys. Res. Lett. 39, 04705. https://doi.org/10.1029/2011GL050582 (2012).ADS 
Article 

Google Scholar 
Pepin, N. C. & Lundquist, J. D. Temperature trends at high elevations: patterns across the globe. Geophys. Res. Lett. 35, 14 (2008).Article 

Google Scholar 
Rangwala, I. & Miller, J. R. Climate change in mountains: a review of elevation-dependent warming and its possible causes. Clim. Change 114, 527–547 (2012).ADS 
Article 

Google Scholar 
Wang, Q., Fan, X. & Wang, M. Recent warming amplification over high elevation regions across the globe. Clim. Dyn. 43, 87–101 (2014).CAS 
Article 

Google Scholar 
Fan, X., Wang, Q., Wang, M. & Jiménez, C. V. Warming amplification of minimum and maximum temperatures over high-elevation regions across the globe. PLoS ONE 10, e0140213 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
Pepin, N. et al. Elevation-dependent warming in mountain regions of the world. Nat. Clim. Chang. 5, 424 (2015).ADS 
Article 

Google Scholar 
Piccarreta, M., Lazzari, M. & Pasini, A. Trends in daily temperature extremes over the Basilicata region (southern Italy) from 1951 to 2010 in a Mediterranean climatic context. Int. J. Climatol. 35, 1964–1975 (2015).Article 

Google Scholar 
Gonzalez-Hidalgo, J. C., Peña-Angulo, D., Brunetti, M. & Cortesi, N. Recent trend in temperature evolution in Spanish mainland (1951–2010): from warming to hiatus. Int. J. Climatol. 36, 2405–2416 (2016).Article 

Google Scholar 
McCullough, I. M. et al. High and dry: high elevations disproportionately exposed to regional climate change in Mediterranean-climate landscapes. Landsc. Ecol. 31, 1063–1075 (2016).Article 

Google Scholar 
Sanz-Elorza, M., Dana, E. D., González, A. & Sobrino, E. Changes in the high-mountain vegetation of the central Iberian Peninsula as a probable sign of global warming. Ann. Bot. 92, 273–280 (2003).PubMed 
PubMed Central 
Article 

Google Scholar 
Peñuelas, J. & Boada, M. A global change induced biome shift in the Montseny mountains (NE Spain). Glob. Change Biol. 9, 131–140 (2003).ADS 
Article 

Google Scholar 
Linares, J. C. & Tíscar, P. A. Climate change impacts and vulnerability of the southern populations of Pinus nigra subsp. salzmannii. Tree Physiol. 30, 795–806 (2010).PubMed 
Article 

Google Scholar 
Giorgi, F., Hurrell, J. W., Marinucci, M. R. & Beniston, M. Elevation dependency of the surface climate change signal: a model study. J. Clim. 10, 288–296 (1997).ADS 
Article 

Google Scholar 
Palazzi, E., Mortarini, L., Terzago, S. & Von Hardenberg, J. Elevation-dependent warming in global climate model simulations at high spatial resolution. Clim. Dyn. 52, 2685–2702 (2019).Article 

Google Scholar 
Poyatos, R., Latron, J. & Llorens, P. Land use and land cover change after agricultural abandonment. Mt. Res. Dev. 23, 362–368 (2003).Article 

Google Scholar 
Mouillot, F., Ratte, J. P., Joffre, R., Mouillot, D. & Rambal, S. Long-term forest dynamic after land abandonment in a fire prone Mediterranean landscape (central Corsica, France). Landsc. Ecol. 20, 101–112 (2005).Article 

Google Scholar 
Zellweger, F. et al. Forest microclimate dynamics drive plant responses to warming. Science 368, 772–775 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Ríos-Cornejo, D., Penas, Á., Álvarez-Esteban, R. & Del Río, S. Links between teleconnection patterns and mean temperature in Spain. Theor. Appl. Climatol. 122, 1–18 (2015).ADS 
Article 

Google Scholar 
Nogués-Bravo, D., Araújo, M. B., Errea, M. P. & Martinez-Rica, J. P. Exposure of global mountain systems to climate warming during the 21st Century. Glob. Environ. Chang. 17, 420–428 (2007).Article 

Google Scholar 
Vicente-Serrano, S. M., Beguería, S., López-Moreno, J. I., El Kenawy, A. M. & Angulo-Martínez, M. Daily atmospheric circulation events and extreme precipitation risk in northeast Spain: Role of the North Atlantic Oscillation, the Western Mediterranean Oscillation, and the Mediterranean Oscillation. J. Geophys. Res. Atmos. 114, D8 (2009).Article 

Google Scholar 
Guzman-Morales, J. & Gershunov, A. Climate change suppresses Santa Ana winds of southern California and sharpens their seasonality. Geophys. Res. Lett. 46, 2772–2780. https://doi.org/10.1029/2018GL080261 (2019).ADS 
Article 

Google Scholar 
Yu, M. & Ruggieri, E. Change point analysis of global temperature records. Int. J. Climatol. 39, 3679–3688 (2019).Article 

Google Scholar 
Giorgi, F. Climate change hot-spots. Geophys. Res. Lett. 33, 08707. https://doi.org/10.1029/2006GL025734 (2006).ADS 
Article 

Google Scholar 
García, M. J. L. Recent warming in the Balearic Sea and Spanish Mediterranean coast: Towards an earlier and longer summer. Atmósfera 28, 149–160 (2015).Article 

Google Scholar 
Toreti, A., Desiato, F., Fioravanti, G. & Perconti, W. Seasonal temperatures over Italy and their relationship with low-frequency atmospheric circulation patterns. Clim. Change 99, 211–227 (2010).ADS 
Article 

Google Scholar 
Scorzini, A. R. & Leopardi, M. Precipitation and temperature trends over central Italy (Abruzzo Region): 1951–2012. Theor. Appl. Climatol. 135, 959–977 (2019).ADS 
Article 

Google Scholar 
Lee, X. et al. Observed increase in local cooling effect of deforestation at higher latitudes. Nature 479, 384–387 (2011).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Juang, J.-Y., Katul, G., Siqueira, M., Stoy, P. & Novick, K. Separating the effects of albedo from eco-physiological changes on surface temperature along a successional chronosequence in the southeastern United States. Geophys. Res. Lett. 34, 21408. https://doi.org/10.1029/2007.GL03129 (2007).ADS 
Article 

Google Scholar 
Boulant, N., Kunstler, G., Rambal, S. & Lepart, J. Seed supply, drought, and grazing determine spatio-temporal patterns of recruitment for native and introduced invasive pines in grasslands. Divers. Distrib. 14, 862–874 (2008).Article 

Google Scholar 
Améztegui, A. Land-use changes as major drivers of mountain pine (Pinus uncinata Ram.) expansion in the Pyrenees. Glob. Ecol. Biogeogr. 19, 632–641 (2010).
Google Scholar 
Rambal, S. Relations entre couverts végétaux des parcours et cycle de l’eau. In L’eau des troupeaux en alpages et sur parcours: une ressource à gérer, aménager, partager (ed. Lepart, J.) 25–37 (Association Française de Pastoralisme et Cardère éditeur, 2015).
Google Scholar 
Fonderflick, J., Lepart, J., Caplat, P., Debussche, M. & Marty, P. Managing agricultural change for biodiversity conservation in a Mediterranean upland. Biol. Conserv. 143, 737–746 (2010).Article 

Google Scholar 
Abadie, J. et al. Forest recovery since 1860 in a Mediterranean region: Drivers and implications for land use and land cover spatial distribution. Landsc. Ecol. 33, 289–305 (2018).Article 

Google Scholar 
Cervera, T., Pino, J., Marull, J., Padró, R. & Tello, E. Understanding the long-term dynamics of forest transition: From deforestation to afforestation in a Mediterranean landscape (Catalonia, 1868–2005). Land Use Policy 80, 318–331 (2019).Article 

Google Scholar 
Wolpert, F., Quintas-Soriano, C. & Plieninger, T. Exploring land-use histories of tree-crop landscapes: a cross-site comparison in the Mediterranean Basin. Sustain. Sci. 15, 1267–1283 (2020).Article 

Google Scholar 
Lasanta-Martínez, T., Vicente-Serrano, S. M. & Cuadrat-Prats, J. M. Mountain Mediterranean landscape evolution caused by the abandonment of traditional primary activities: A study of the Spanish Central Pyrenees. Appl. Geogr. 25, 47–65 (2005).Article 

Google Scholar 
Malandra, F., Vitali, A., Urbinati, C., Weisberg, P. J. & Garbarino, M. Patterns and drivers of forest landscape change in the Apennines range, Italy. Reg. Environ. Change 19, 1973–1985 (2019).Article 

Google Scholar 
Zhang, Q. et al. Reforestation and surface cooling in temperate zones: Mechanisms and implications. Glob. Change Biol. 26, 3384–3401 (2020).ADS 
Article 

Google Scholar 
Rambal, S., Lacaze, B. & Winkel, T. Testing an area-weighted model for albedo or surface temperature of mixed pixels in Mediterranean woodlands. Int. J. Remote Sens. 11, 1495–1499 (1990).Article 

Google Scholar 
Luyssaert, S. et al. Land management and land-cover change have impacts of similar magnitude on surface temperature. Nat. Clim. Change 4, 389–393. https://doi.org/10.1038/nclimate2196 (2014).ADS 
Article 

Google Scholar 
Novick, K. A. & Katul, G. G. The duality of reforestation impacts on surface and air temperature. J. Geophys. Res. Biogeosci. 125, e05543 (2020).Article 

Google Scholar 
Davy, R. & Esau, I. Differences in the efficacy of climate forcings explained by variations in atmospheric boundary layer depth. Nat. Commun. 7, 1–8 (2016).Article 

Google Scholar 
Serafin, S. et al. Exchange processes in the atmospheric boundary layer over mountainous terrain. Atmosphere 9, 102. https://doi.org/10.3390/atmos9030102 (2018).ADS 
CAS 
Article 

Google Scholar 
Perugini, L. et al. Biophysical effects on temperature and precipitation due to land cover change. Environ. Res. Lett. 12, 053002 (2017).ADS 
Article 

Google Scholar 
Visbeck, M. H., Hurrell, J. W., Polvani, L. & Cullen, H. M. The North Atlantic oscillation: Past, present, and future. Proc. Natl. Acad. Sci. 98, 12876–12877 (2001).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hurrell, J. W. Decadal trends in the North Atlantic oscillation: Regional temperatures and precipitation. Science 269, 676–679 (1995).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Martín, P., Sabatés, A., Lloret, J. & Martin-Vide, J. Climate modulation of fish populations: the role of the Western Mediterranean Oscillation (WeMO) in sardine (Sardina pilchardus) and anchovy (Engraulis encrasicolus) production in the north-western Mediterranean. Clim. Change 110, 925–939 (2012).ADS 
Article 

Google Scholar 
Schwingshackl, C., Hirschi, M. & Seneviratne, S. I. Global contributions of incoming radiation and land surface conditions to maximum near surface air temperature variability and trend. Geophys. Res. Lett. 45, 5034–5044 (2018).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Philipona, R., Behrens, K. & Ruckstuhl, C. How declining aerosols and rising greenhouse gases forced rapid warming in Europe since the 1980s. Geophys. Res. Lett. 36, L02806. https://doi.org/10.1029/2008GL036350 (2009).ADS 
Article 

Google Scholar 
Schwarz, M., Folini, D., Yang, S., Allan, R. P. & Wild, M. Changes in atmospheric shortwave absorption as important driver of dimming and brightening. Nat. Geosci. 13, 110–115 (2020).ADS 
CAS 
Article 

Google Scholar 
Norris, J. R. & Wild, M. Trends in aerosol radiative effects over Europe inferred from observed cloud cover, solar “dimming”, and solar “brightening”. J. Geophys. Res. Atmos. 112, D08214. https://doi.org/10.1029/2006JD007794 (2007).ADS 
Article 

Google Scholar 
Mateos, D. et al. Quantifying the respective roles of aerosols and clouds in the strong brightening since the early 2000s over the Iberian Peninsula. J. Geophys. Res. Atmos. 119, 10–382 (2014).Article 

Google Scholar 
Sanchez-Lorenzo, A. et al. Reassessment and update of long-term trends in downward surface shortwave radiation over Europe (1939–2012). J. Geophys. Res. Atmos. 120, 9555–9569 (2015).ADS 
Article 

Google Scholar 
Kambezidis, H. D., Kaskaoutis, D. G., Kalliampakos, G. K., Rashki, A. & Wild, M. The solar dimming/brightening effect over the Mediterranean Basin in the period 1979–2012. J. Atmos. Solar Terr. Phys. 150, 31–46 (2016).ADS 
Article 

Google Scholar 
Chiacchio, M. & Wild, M. Influence of NAO and clouds on long-term seasonal variations of surface solar radiation in Europe. J. Geophys. Res. Atmos. 115, 0022. https://doi.org/10.1029/2009JD012182 (2010).Article 

Google Scholar 
Wild, M. Decadal changes in radiative fluxes at land and ocean surfaces and their relevance for global warming. Wiley Interdiscipl. Rev. Clim. Change 7, 91–107 (2016).Article 

Google Scholar 
Held, I. M. & Soden, B. J. Water vapor feedback and global warming. Annu. Rev. Energy Environ. 25, 441–475 (2000).Article 

Google Scholar 
Dessler, A. E. & Sherwood, S. C. A matter of humidity. Science 323, 1020–1021 (2009).CAS 
PubMed 
Article 

Google Scholar 
Ruckstuhl, C., Philipona, R., Morland, J. & Ohmura, A. Observed relationship between surface specific humidity, integrated water vapor, and longwave downward radiation at different altitudes. J. Geophys. Res. Atmos. 112(D03302), 2007. https://doi.org/10.1029/2006JD007850 (2007).Article 

Google Scholar 
Parras-Berrocal, I. M. et al. The climate change signal in the Mediterranean Sea in a regionally coupled atmosphere–ocean model. Ocean Sci. 16, 743–765. https://doi.org/10.5194/os-16-743-2020 (2020).ADS 
Article 

Google Scholar 
Reale, M. et al. The regional earth system model RegCM-ES: Evaluation of the Mediterranean climate and marine biogeochemistry. J. Adv. Model. Earth Syst. 12, e001812 (2020).Article 

Google Scholar 
Sen, P. K. Estimates of the regression coefficient based on Kendall’s tau. J. Am. Stat. Assoc. 63, 1379–1389 (1968).MathSciNet 
MATH 
Article 

Google Scholar 
Kelliher, F. M., Leuning, R. & Schulze, E. D. Evaporation and canopy characteristics of coniferous forests and grasslands. Oecologia 95, 153–163 (1993).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Linacre, E. T. Simpler empirical expression for actual evapotranspiration rates-a discussion. Agric. Meteorol. 11, 451–452 (1973).Article 

Google Scholar 
Jones, P. D., Jónsson, T. & Wheeler, D. Extension to the North Atlantic Oscillation using early instrumental pressure observations from Gibraltar and south-west Iceland. Int. J. Climatol. 17, 1433–1450 (1997).Article 

Google Scholar 
Palutikof, J. P. Analysis of Mediterranean climate data: measured and modelled. In Mediterranean Climate: Variability and Trends (ed. Bolle, H. J.) (Springer, 2003).
Google Scholar 
Martin-Vide, J. & Lopez-Bustins, J. A. The western Mediterranean oscillation and rainfall in the Iberian Peninsula. Int. J. Climatol. 26, 1455–1475 (2006).Article 

Google Scholar  More

We selected eleven major wheat production provinces of China for the study area, which comprise the largest winter wheat-sowing fraction: Henan, Shandong, Anhui, Jiangsu, Hebei, Hubei, Shanxi, Shaanxi, Sichuan, Xinjiang, and Gansu (Fig. 1). The wheat planting area is about 22 million ha in these provinces, accounting for more than 93% of the total wheat planting area. The total wheat production in these regions contributes more than 96% of the total wheat production in China, with more than 128 million tons in 201933.We developed a methodological framework for high-resolution AGB mapping. It mainly includes three parts: (1) Data acquisition and processing. (2) The WOFOST model parameterization and calibration. (3) Data assimilation (Fig. 1). Each part is explained in more detail below.Data acquisition and processingMeteorological dataChina Meteorological Forcing Dataset34,35 is used as weather driving data for the WOFOST model. The dataset is based on the internationally existing Princeton reanalysis data, Global Land Data Assimilation System data, Global Energy and Water Cycle Experiment-Surface Radiation Budget radiation data, and Tropical Rainfall Measuring Mission precipitation data. It is made by fusing the conventional meteorological observation data of the China Meteorological Administration. It includes seven elements: near-surface air temperature, air pressure, near-surface total humidity, wind speed, ground downward shortwave radiation, ground downward longwave radiation, and ground precipitation rate. The meteorological drive elements required for WOFOST are daily radiation, minimum temperature, maximum temperature, water vapor pressure, average wind speed, and precipitation. Details of these variables that participated in the WOFOST model can be referred to in Table S1.Soil characteristics measurements and phenology observationsSoil and phenology data were collected at 177 agricultural meteorological stations (AMS) from 2007 to 2015 (Fig. 1). Soil characteristics include soil moisture content at wilting points, field capacity, and saturation. To be consistent with the corresponding units in the crop model, the original data in weight water content was converted into volume water content through the corresponding soil bulk density measurements. Winter wheat phenology observations include the date of emergence (more than 50% of the wheat seedlings in the field show the first green leaves and reached about 2 cm), anthesis (the inner and outer glumes of the middle and upper florets of more than 50% of the wheat ears in the whole field are open, and the anthers loose powder), and maturity (more than 80% of the wheat grains turn yellow, the glumes and stems turn yellow, and only the upper first and second nodes are still slightly green). In most cases, the phenological stage “anthesis” is missing. The anthesis date was calculated by adding seven days to the observed heading date (when more than 50% of the wheat in the whole field exposes the tip of the ear from the sheath of the flag leaf).County-level yield statistics dataThe county-level yield data was collected from city statistical yearbooks of the study area from 2007 to 2015. Since most statistical yearbooks do not directly record per-unit yield data, the county-level yield was obtained by dividing the total yield and planting area. It is worth noting that all yields were calculated in units of metric kilograms per cultivated hectares (kg·ha−1).The winter wheat land cover dataWe used a winter wheat land cover product from a 1 km resolution dataset named ChinaCropArea1km36. This data was derived from GLASS leaf area index products and crop phenology from 2000 to 2015. This dataset is the base map of our data production.The MODIS LAI dataWe used the improved 8-days MODIS LAI products (i.e., 1 km) generated by Yuan et al.32 to assimilate the WOFOST model. The products applied the modified temporal-spatial filter and Savitzky-Golay filter to overcome the spatial-temporal discontinuity and inconsistence of raw MODIS LAI products, which makes them more applicable for the realm of land surface and climate modeling. The products can be accessed via the Land-Atmosphere Interaction Research Group website at Sun Yat-sen University (http://globalchange.bnu.edu.cn/research/lai).The WOFOST model parameterization and calibrationThe WOFOST model introductionThe WOFOST model was initially developed as a crop growth simulation model to evaluate the yield potential of various crops in tropical countries37. In this study, we chose the WOFOST model because the model reaches a trade-off of the complexity of the crop model and is suitable for large-scale simulations3. The WOFOST model is a typical crop growth model that explains crop growth based on underlying processes such as photosynthesis and respiration and their response to changing environmental conditions38. The WOFOST model estimates phenology, LAI, aboveground biomass, and storage organ biomass (i.e., grain yield) at a daily time step39 (Fig. 2).Fig. 2Schematic overview of the major processes implemented in WOFOST. The Astronomical module calculates day length, some variables relating to solar elevation, and the fraction of diffuse radiation.Full size imageZonal parameterizationWe first divided the study area covered by AMS into seamless Thiessen polygon zones. Each Thiessen polygon contains only a single AMS. These zones represent the whole areas where any location is closer to its associated AMS point than any other AMS point. Then, we assigned parameters to the entire zone based on the AMS data, including crop calendar (date of emergence) and soil water retention parameters (soil moisture content at wilting point, field capacity, and saturation). Besides, we also optimized two main crop parameters for controlling phenological development stages, namely TSUM1 (accumulated temperature required from emergence to anthesis) and TSUM2 (accumulated temperature required from anthesis to maturity), by minimizing the cost function of the observational and simulated date corresponding to anthesis and maturity.Parameter calibration within a single zoneWe implemented the calibration of parameters within every single zone, as illustrated in Fig. 3. We calculated the average statistical yield of each county within every single zone from 2007 to 2015, then ranked the counties in descending order and divided them into three groups, namely high, medium, and low-level yield counties, by the 33% quantile and 67% quantile of the average statistical yield. The three counties corresponding to 17% quantile, 50% quantile, and 83% quantile would be used for subsequent calibration and represent the corresponding three yield level groups. We used the statistical yields (converted to dry matter mass based on the standard moisture content of 12.5%) of the three counties for multiple years and a harvest index for each province to convert the county-level yield to AGB for calibration. The harvest index of each province was mainly estimated from the AMS’s dynamic growth records on the biomass composition of the dominant winter wheat varieties of the province and a published literature40. Besides, we collected the maximum LAI observations on all agrometeorological stations in all years in the study area, according to its histogram. We found that the histogram follows a normal distribution with a mean of 6.5 and a standard deviation of 1.5. Finally, we calibrated three sets of parameters corresponding to three yield level groups in each single zone according to the three selected counties.Fig. 3Flow chart of parameter calibration within a single zone.Full size imageWe designed a three-step calibration strategy for a specific yield level group. Firstly, as winter wheat varieties did not change significantly according to information recorded by agrometeorological stations from 2007 to 2015, we assumed the crop parameters of winter wheat remain unchanged every three years to combine three years of observational data to calibrate the parameters of the WOFOST model better. We maximized a log-likelihood function based on the maximum LAI statistics and every three-year county-level yield and AGB data mentioned to optimize selected crop parameters (see Table S2 in the Supplement Materials).The log-likelihood function was constructed as follows:$$log;{{rm{L}}}_{{rm{LAI}}}=-frac{1}{2}left[dlogleft(2pi right)+logleft(left|{Sigma }_{{rm{LAI}}}right|right)+{rm{MD}}{left({{bf{x}}}_{{rm{LAI}}};{mu }_{{rm{LAI}}},{Sigma }_{{rm{LAI}}}right)}^{2}right]$$
(1)
$$log;{{rm{L}}}_{{rm{TWSO}}}=-frac{1}{2}left[dlog(2pi )+logleft(left|{{boldsymbol{Sigma }}}_{{rm{TWSO}}}right|right)+{rm{MD}}{left({{bf{x}}}_{{rm{TWSO}}};{{boldsymbol{mu }}}_{{rm{TWSO}}},{{boldsymbol{Sigma }}}_{{rm{TWSO}}}right)}^{2}right]$$
(2)
$$log;{{rm{L}}}_{{rm{AGB}}}=-frac{1}{2}left[dlog(2pi )+logleft(left|{{boldsymbol{Sigma }}}_{{rm{AGB}}}right|right)+{rm{MD}}{left({{bf{x}}}_{{rm{AGB}}};{{boldsymbol{mu }}}_{{rm{AGB}}},{{boldsymbol{Sigma }}}_{{rm{AGB}}}right)}^{2}right]$$
(3)
$$log;{rm{L}}=log;{L}_{{rm{LAI}}}+log;{L}_{{rm{TWSO}}}+log;{L}_{{rm{AGB}}}$$
(4)
Where log L is the natural logarithm of the likelihood function, d is the dimension, that is, the number of years of joint calibration, which is set to 3 in this study xLAI is the vector composed of the maximum value of the 3-year LAI simulated by WOFOST, μLAI and ΣLAI are the mean vector and error covariance matrix of maximum LAI based on observation statistics. The annual maximum LAI was assumed to be independent, and the mean and standard deviation for each year was set the same as the result of Fig. 3. Similarly, xTWSO and xAGB are the yield vector and AGB vector at maturity of 3 years simulated by WOFOST, and μTWSO, μAGB are their corresponding county-level statistic vector, ΣTWSO and ΣAGB are their corresponding error covariance matrix. In this study, we assumed that the annual yield or AGB was independent, and their corresponding standard deviation was 10% of their statistical value. |Σ| is the determinant of Σ. The expression ({rm{MD}}{({bf{x}};{boldsymbol{mu }},{boldsymbol{Sigma }})}^{2}={({bf{x}}-{boldsymbol{mu }})}^{{rm{T}}}{{boldsymbol{Sigma }}}^{-1}({bf{x}}-{boldsymbol{mu }})), where MD is the Mahalanobis distance between the point x and the mean vector μ.Secondly, we optimized the inter-annual irrigation. We optimized two parameters every year: the critical value of soil moisture (denoted as SMc) and the amount of irrigation (denoted as V). When the soil moisture simulated by WOFOST is lower than SMc, an irrigation event will be triggered, and the irrigation amount is V cm. In this study, we combined three-year data for calibration with six parameters for optimization. The optimization strategy is the same as the previous step by maximizing the log-likelihood function. Finally, we fixed the optimized irrigation parameters and repeated the first step to calibrate the selected crop parameters and obtain the final optimal parameters.Data assimilationConsidering that MODIS LAI is relatively low compared to the actual LAI of winter wheat41, we select a weak-constraint cost function based on the least square of normalized observational and simulated LAI as shown in Eq. (5), which is assimilating the trend information of MODIS LAI into the crop growth model.$$J={sum }_{{rm{t}}=1}^{{rm{n}}}{left(frac{{{rm{LAI}}}_{{rm{MODIS}}}^{{rm{t}}}-{{rm{LAI}}}_{{rm{MODIS}}}^{min}}{{{rm{LAI}}}_{{rm{MODIS}}}^{max}-{{rm{LAI}}}_{{rm{MODIS}}}^{min}}-frac{{{rm{LAI}}}_{{rm{WOFOS}}}^{{rm{t}}}-{{rm{LAI}}}_{{rm{WOFOS}}}^{min}}{{{rm{LAI}}}_{{rm{WOFOS}}}^{max}-{{rm{LAI}}}_{{rm{WOFOS}}}^{min}}right)}^{2}$$
(5)
Where ({{rm{LAI}}}_{{rm{MODIS}}}^{{rm{t}}}) and .. are MODIS LAI and WOFOST simulated LAI of time t. ({{rm{LAI}}}_{{rm{MODIS}}}^{max}) and ({{rm{LAI}}}_{{rm{WOFOS}}}^{max}) are maximum of MODIS LAI and WOFOST simulated LAI. ({{rm{LAI}}}_{{rm{MODIS}}}^{min}) and ({{rm{LAI}}}_{{rm{WOFOS}}}^{min}) are minimum of MODIS LAI and WOFOST simulated LAI. J is the value of the cost function.We reinitialize the day of emergence (IDEM), the life span of leaves growing at 35 °C (SPAN), and thermal time from emergence to anthesis (TSUM1) in the WOFOST model on each 1 km winter wheat pixel according to cost function between WOFOST LAI and MODIS LAI. Besides, we applied the Subplex algorithm from the NLOPT library (https://github.com/stevengj/nlopt) for parameter optimization. More