Philippa Kaur

Spillover events, in which a pathogen that originates in animals jumps into people, have probably triggered every viral pandemic that’s occurred since the start of the twentieth century1. What’s more, an August 2021 analysis of disease outbreaks over the past four centuries indicates that the yearly probability of pandemics could increase several-fold in the coming decades, largely because of human-induced environmental changes2.Fortunately, for around US$20 billion per year, the likelihood of spillover could be greatly reduced3. This is the amount needed to halve global deforestation in hotspots for emerging infectious diseases; drastically curtail and regulate trade in wildlife; and greatly improve the ability to detect and control infectious diseases in farmed animals.That is a small investment compared with the millions of lives lost and trillions of dollars spent in the COVID-19 pandemic. The cost is also one-twentieth of the statistical value of the lives lost each year to viral diseases that have spilled over from animals since 1918 (see ‘Spillovers: a growing threat’), and less than one-tenth of the economic productivity erased per year1.

Source: Ref. 1

Yet many of the international efforts to better defend the world from future outbreaks, prompted by the COVID-19 pandemic, still fail to prioritize the prevention of spillover. Take, for example, the Independent Panel for Pandemic Preparedness and Response, established by the World Health Organization (WHO). The panel was convened in September 2020, in part to ensure that any future infectious-disease outbreak does not become another pandemic. In its 86-page report released last May, wildlife is mentioned twice; deforestation once.We urge the decision-makers currently developing three landmark international endeavours to make the prevention of spillover central to each.First, the G20 group of the world’s 20 largest economies provisionally agreed last month to create a global fund for pandemics. If realized, this could provide funding at levels that infectious-disease experts have been recommending for decades — around $5 per person per year globally (see go.nature.com/3yjitwx). Second, an agreement to improve global approaches to pandemics is under discussion by the World Health Assembly (WHA), the decision-making body of the WHO. Third, a draft framework for biodiversity conservation — the post-2020 global biodiversity framework — is being negotiated by parties to the Convention on Biological Diversity.Designed in the right way, these three international endeavours could foster a more proactive global approach to infectious diseases. This opportunity — to finally address the factors that drive major disease outbreaks, many of which also contribute to climate change and biodiversity loss — might not present itself again until the world faces another pandemic.Four actions The risk of spillover is greater when there are more opportunities for animals and humans to make contact, for instance in the trade of wildlife, in animal farming or when forests are cleared for mining, farming or roads. It is also more likely to happen under conditions that increase the likelihood of infected animals shedding viruses – when they are housed in cramped conditions, say, or not fed properly.Decades of research from epidemiology, ecology and genetics suggest that an effective global strategy to reduce the risk of spillover should focus on four actions1,3.First, tropical and subtropical forests must be protected. Various studies show that changes in the way land is used, particularly tropical and subtropical forests, might be the largest driver of emerging infectious diseases of zoonotic origin globally4. Wildlife that survives forest clearance or degradation tends to include species that can live alongside people, and that often host pathogens capable of infecting humans5. For example, in Bangladesh, bats that carry Nipah virus — which can kill 40–75% of people infected — now roost in areas of high human population density because their forest habitat has been almost entirely cleared6.Furthermore, the loss of forests is driving climate change. This could in itself aid spillover by pushing animals, such as bats, out of regions that have become inhospitable and into areas where many people live7.Yet forests can be protected even while agricultural productivity is increased — as long as there is enough political will and resources8. This was demonstrated by the 70% reduction in deforestation in the Amazon during 2004–12, largely through better monitoring, law enforcement and the provision of financial incentives to farmers. (Deforestation rates began increasing in 2013 due to changes in environmental legislation, and have risen sharply since 2019 during Jair Bolsonaro’s presidency.)Second, commercial markets and trade of live wild animals that pose a public-health risk must be banned or strictly regulated, both domestically and internationally.Doing this would be consistent with the call made by the WHO and other organizations in 2021 for countries to temporarily suspend the trade in live caught wild mammals, and to close sections of markets selling such animals. Several countries have already acted along these lines. In China, the trade and consumption of most terrestrial wildlife has been banned in response to COVID-19. Similarly, Gabon has prohibited the sale of certain mammal species as food in markets.

A worker in a crowded chicken farm in Anhui province, China.Credit: Jianan Yu/Reuters

Restrictions on urban and peri-urban commercial markets and trade must not infringe on the rights and needs of Indigenous peoples and local communities, who often rely on wildlife for food security, livelihoods and cultural practices. There are already different rules for hunting depending on the community in many countries, including Brazil, Canada and the United States.Third, biosecurity must be improved when dealing with farmed animals. Among other measures, this could be achieved through better veterinary care, enhanced surveillance for animal disease, improvements to feeding and housing animals, and quarantines to limit pathogen spread.Poor health among farmed animals increases their risk of becoming infected with pathogens — and of spreading them. And nearly 80% of livestock pathogens can infect multiple host species, including wildlife and humans9.Fourth, particularly in hotspots for the emergence of infectious diseases, people’s health and economic security should be improved.People in poor health — such as those who have malnutrition or uncontrolled HIV infection — can be more susceptible to zoonotic pathogens. And, particularly in immunosuppressed individuals such as these, pathogens can mutate before being passed on to others10.What’s more, some communities — especially those in rural areas — use natural resources to produce commodities or generate income in a way that brings them into contact with wildlife or wildlife by-products. In Bangladesh, for example, date palm sap, which is consumed as a drink in various forms, is often collected in pots attached to palm trees. These can become contaminated with bodily substances from bats. A 2016 investigation linked this practice to 14 Nipah virus infections in humans that caused 8 deaths11.Providing communities with both education and tools to reduce the risk of harm is crucial. Tools can be something as simple as pot covers to prevent contamination of date palm sap, in the case of the Bangladesh example.In fact, providing educational opportunities alongside health-care services and training in alternative livelihood skills, such as organic agriculture, can help both people and the environment. For instance, the non-governmental organization Health in Harmony in Portland, Oregon, has invested in community-designed interventions in Indonesian Borneo. During 2007–17, these contributed to a 90% reduction in the number of households that were reliant on illegal logging for their main livelihood. This, in turn, reduced local rainforest loss by 70%. Infant mortality also fell by 67% in the programme’s catchment area12.Systems-oriented interventions of this type need to be better understood, and the most effective ones scaled up.Wise investmentSuch strategies to prevent spillover would reduce our dependence on containment measures, such as human disease surveillance, contact tracing, lockdowns, vaccines and therapeutics. These interventions are crucial, but are often expensive and implemented too late — in short, they are insufficient when used alone to deal with emerging infectious diseases.The COVID-19 pandemic has exposed the real-world limitations of these reactive measures — particularly in an age of disinformation and rising populism. For example, despite the US federal government spending more than $3.7 trillion on its pandemic response as of the end of March, nearly one million people in the United States — or around one in 330 — have died from COVID-19 (see go.nature.com/39jtdfh and go.nature.com/38urqvc). Globally, between 15 million and 21 million lives are estimated to have been lost during the COVID-19 pandemic beyond what would be expected under non-pandemic conditions (known as excess deaths; see Nature https://doi.org/htd6; 2022). And a 2021 model indicates that, by 2025, $157 billion will have been spent on COVID-19 vaccines alone (see go.nature.com/3jqds76).

A farmer in Myanmar gathers sap from a palm tree to make wine. Contamination of the collection pots with excretions from bats can spread diseases to humans.Credit: Wolfgang Kaehler/LightRocket via Getty

Preventing spillover also protects people, domesticated animals and wildlife in the places that can least afford harm — making it more equitable than containment. For example, almost 18 months since COVID-19 vaccines first became publicly available, only 21% of the total population of Africa has received at least one dose. In the United States and Canada, the figure is nearly 80% (see go.nature.com/3vrdpfo). Meanwhile, Pfizer’s total drug sales rose from $43 billion in 2020 to $72 billion in 2021, largely because of the company’s COVID-19 vaccine, the best-selling drug of 202113.Lastly, unlike containment measures, actions to prevent spillover also help to stop spillback, in which zoonotic pathogens move back from humans to animals and then jump again into people. Selection pressures can differ across species, making such jumps a potential source of new variants that can evade existing immunity. Some researchers have suggested that spillback was possibly responsible for the emergence of the Omicron variant of SARS-CoV-2 (see Nature 602, 26–28; 2022).Seize the dayOver the past year, the administration of US President Joe Biden and two international panels (one established in 2020 by the WHO and the other in 2021 by the G20) have released guidance on how to improve approaches to pandemics. All recommendations released so far acknowledge spillover as the predominant cause of emerging infectious diseases. None adequately discusses how that risk might be mitigated. Likewise, a PubMed search for the spike protein of SARS-CoV-2 yields thousands of papers, yet only a handful of studies investigate coronavirus dynamics in bats, from which SARS-CoV-2 is likely to have originated14.Spillover prevention is probably being overlooked for several reasons. Upstream animal and environmental sources of pathogens might be being neglected by biomedical researchers and their funders because they are part of complex systems — research into which does not tend to lead to tangible, profitable outputs. Also, most people working in public health and biomedical sciences have limited training in ecology, wildlife biology, conservation and anthropology.There is growing recognition of the importance of cross-sectoral collaboration, including soaring advocacy for the ‘One Health’ approach — an integrated view of health that recognizes links between the environment, animals and humans. But, in general, this has yet to translate into action to prevent pandemics.Another challenge is that it can take decades to realize the benefits of preventing spillover, instead of weeks or months for containment measures. Benefits can be harder to quantify for spillover prevention, no matter how much time passes, because, if measures are successful, no outbreak occurs. Prevention also runs counter to individual, societal and political tendencies to wait for a catastrophe before taking action.The global pandemic fund, the WHA pandemic agreement and the post-2020 global biodiversity framework all present fresh chances to shift this mindset and put in place a coordinated global effort to reduce the risk of spillover alongside crucial pandemic preparedness efforts.Global fund for pandemicsFirst and foremost, a global fund for pandemics will be key to ensuring that the wealth of evidence on spillover prevention is translated into action. Funding for spillover prevention should not be folded into existing conservation funds, nor draw on any other existing funding streams.Investments must be targeted to those regions and practices where the risk of spillover is greatest, from southeast Asia and Central Africa to the Amazon Basin and beyond. Actions to prevent spillover in these areas, particularly by reducing deforestation, would also help to mitigate climate change and reduce loss of biodiversity. But conservation is itself drastically underfunded. As an example, natural solutions (such as conservation, restoration and improved management of forests, wetlands and grasslands) represent more than one-third of the climate mitigation needed by 2030 to stabilize warming to well below 2 °C15. Yet these approaches receive less than 2% of global funds for climate mitigation16. (Energy systems receive more than half.)In short, the decision-makers backing the global fund for pandemics must not assume that existing funds are dealing with the threat of spillover — they are not. The loss of primary tropical forest was 12% higher in 2020 than in 2019, despite the economic downturn triggered by COVID-19. This underscores the continuing threat to forests.Funding must be sustained for decades to ensure that efforts to reduce the risk of spillover are in place long enough to yield results.WHA pandemic agreementIn 2020, the president of the European Council, Charles Michel, called for a treaty to enable a more coordinated global response to major epidemics and pandemics. Last year, more than 20 world leaders began echoing this call, and the WHA launched the negotiation of an agreement (potentially, a treaty or other international instrument) to “strengthen pandemic prevention, preparedness, and response” at the end of 2021.Such a multilateral agreement could help to ensure more-equitable international action around the transfer of scientific knowledge, medical supplies, vaccines and therapeutics. It could also address some of the constraints currently imposed on the WHO, and define more clearly the conditions under which governments must notify others of a potential disease threat. The COVID-19 pandemic exposed the shortcomings of the International Health Regulations on many of these fronts17. (This legal framework defines countries’ rights and obligations in the handling of public-health events and emergencies that could cross borders.)We urge negotiators to ensure that the four actions to prevent spillover outlined here are prioritized in the WHA pandemic agreement. For instance, it could require countries to create national action plans for pandemics that include reducing deforestation and closing or strictly regulating live wildlife markets. A reporting mechanism should also be developed to evaluate progress in implementing the agreement. This could build on experience from existing schemes, such as the WHO Joint External Evaluation process (used to assess countries’ capacities to handle public-health risks) and the verification regime of the Chemical Weapons Convention.Commitments to expand pathogen surveillance at interfaces between humans, domesticated animals and wildlife — from US mink farms and Asian wet markets to areas of high deforestation in South America — should also be wrapped into the WHA agreement. Surveillance will not prevent spillover, but it could enable earlier detection and better control of zoonotic outbreaks, and provide a better understanding of the conditions that cause them. Disease surveillance would improve simply through investing in clinical care for both people and animals in emerging infectious-disease hotspots.Convention on Biological DiversityWe are in the midst of the sixth mass extinction, and activities that drive the loss of biodiversity, such as deforestation, also contribute to the emergence of infectious disease. Meanwhile, epidemics and pandemics resulting from the exploitation of nature can lead to further conservation setbacks — because of economic damage from lost tourism and staff shortages affecting management of protected areas, among other factors18. Also, pathogens that infect people can be transmitted to other animals and decimate those populations. For instance, an Ebola outbreak in the Republic of Congo in 2002–03 is thought to have killed 5,000 gorillas19.Yet the global biodiversity framework currently being negotiated by the Convention on Biological Diversity fails to explicitly address the negative feedback cycle between environmental degradation, wildlife exploitation and the emergence of pathogens. The first draft made no mention of pandemics. Text about spillover prevention was proposed in March, but it has yet to be agreed on.Again, this omission stems largely from the siloing of disciplines and expertise. Just as the specialists relied on for the WHA pandemic agreement tend to be those in the health sector, those informing the Convention on Biological Diversity tend to be specialists in environmental science and conservation.The global biodiversity framework, scheduled to be agreed at the Conference of the Parties later this year, must strongly reflect the environment–health connection. This means explicitly including spillover prevention in any text relating to the exploitation of wildlife and nature’s contributions to people. Failing to connect these dots weakens the ability of the convention to achieve its own objectives around conservation and the sustainable use of resources.Preventive health careA reactive response to catastrophe need not be the norm. In many countries, preventive health care for chronic diseases is widely embraced because of its obvious health and economic benefits. For instance, dozens of colorectal cancer deaths are averted for every 1,000 people screened using colonoscopies or other methods20. A preventive approach does not detract from the importance of treating diseases when they occur.With all the stressors now being placed on the biosphere — and the negative implications this has for human health — leaders urgently need to apply this way of thinking to pandemics. 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

Cockroach strainsAll cockroaches were maintained on rodent diet (Purina 5001, PMI Nutrition International, St. Louis, MO) and distilled water at 27 °C, ~40% RH, and a 12:12 h L:D cycle. The WT colony (Orlando Normal) was collected in Florida in 1947 and has served as a standard insecticide-susceptible strain. The GA colony (T-164) was collected in 1989, also in Florida, and shown to be aversive to glucose; continued artificial selection with glucose-containing toxic bait fixed the homozygous GA trait in this population (approximately 150 generations as of 2020).Generating recombinant lines and life history dataTo homogenize the genetic backgrounds of the WT and GA strains, two recombinant colonies were initiated in 2013 by crossing 10 pairs of WT♂ × GA♀ and 10 pairs of GA♂ × WT♀ (Fig. 3a). At the F8 generation (free bulk mating without selection), 400 cockroaches were tested in two-choice feeding assays (see below) that assessed their initial response to tastants, as described in previous studies11,26. The cockroaches were separated into glucose-accepting and glucose-rejecting groups by the rapid Acceptance-Rejection assay (described in Feeding Bioassays). These colonies were bred for three more generations, and 200 cockroaches from each group were assayed in the F11 generation and backcrossed to obtain homozygous glucose-accepting (aa) and glucose-averse (AA) lines. Similar results were obtained in both directions of the cross, confirming previous findings of no sex linkage of the GA trait27. These two lines were defined as WT_aa (homozygotes, glucose-accepting) and GA_AA (homozygotes, glucose-averse). To obtain heterozygous GA cockroaches, GA_Aa, a single intercross group was generated from crosses of 10 pairs of WT_aa♂ × GA_AA♀ and 10 pairs of GA_AA♂ × WT_aa♀.The GA trait follows Mendelian inheritance. Therefore, we used backcrosses, guided by two-choice feeding assays and feeding responses in Acceptance-rejection assays, to determine the homozygosity of WT and GA cockroaches. The cross of WT♂ × WT♀ produced homozygous F1 cockroaches showing maximal glucose-acceptance. The cross of GA♂ × GA♀ produced homozygous F1 cockroaches showing maximal glucose-aversion. The cross of WT × GA produced F1 heterozygotes with intermediate glucose-aversion. When the F1 heterozygotes were backcrossed with WT cockroaches, they produced F2 cockroaches with a 1:1 ratio of WT and GA phenotypes.The two-choice feeding assay assessed whether cockroaches accepted or rejected glucose (binary: yes-no). Insects were held for 24 h without water, or starved without food and water. Either 10 adults or 2 day-old first instar siblings (30–40) were placed in a Petri dish (either 90 mm or 60 mm diameter × 15 mm height). Each Petri dish contained two agar discs: one disc contained 1% agar and 1 mmol l−1 red food dye (Allura Red AC), and the second disc contained 1% agar, 0.5 mmol l−1 blue food dye (Erioglaucine disodium salt) and either 1000 mmol l−1 or 3000 mmol l−1 glucose. The assay duration was 2 h during the dark phase of the insects’ L:D cycle. After each assay, the color of the abdomen of each cockroach was visually inspected under a microscope to infer the genotype.We assessed whether the recombinant colonies had different traits from the parental WT and GA lines. We paired single newly eclosed females (day 0) with single 10–12 days-old males of the same line in a Petri dish (90 mm diameter, 15 mm height) with fresh distilled water in a 1.5 ml microcentrifuge tube and a pellet of rodent food, and monitored when they mated. When females formed egg cases, each gravid female was placed individually in a container (95 × 95 × 80 mm) with food and water until the eggs hatched. After removing the female, her offspring were monitored until adult emergence. We recorded the time to egg hatch, first appearance of each nymphal stage, first appearance of adults and the end of adult emergence. The first instar nymphs and adults in each cohort were counted to obtain measures of survivorship. Although there were significant differences in some of these parameters across all four strains, we found no significant differences between the two recombinant lines, except mating success, which was significantly lower in GA_AA♀ than WT_aa♀ (Supplementary Table 11).Mating bioassaysAll mating sequences were recorded using an infra-red-sensitive camera (Polestar II EQ610, Everfocus Electronics, New Taipei City, Taiwan) coupled to a data acquisition board and analyzed by searchable and frame-by-frame capable software (NV3000, AverMedia Information) at 27 °C, ~40% RH and a 12:12 h L:D cycle. For behavioral analysis, tested pairs were classified into two groups: mated (successful courtship) and not-mated (failed courtship). Four distinct behavioral events (Fig. 1c, Contact, Wing raising, Nuptial feeding, and Copulation) were analyzed using seven behavioral parameters as shown in Supplementary Table 2.We extracted behavioral data from successful courtship sequences, defined as courtship that led to Copulation. For failed courtship sequences, we extracted the behavioral data from the first courtship of both mated and not-mated groups, because most pairs in both groups failed to copulate in their first encounter, and there were no significant differences in behavioral parameters between the two groups.To assay female choice, we conducted two-choice mating assays (Fig. 1a). A single focal WT♀ or GA♀ and two males, one WT and one GA, were placed in a Petri dish (90 mm diameter, 15 mm height) with fresh distilled water in a 1.5 ml microcentrifuge tube and a pellet of rodent food (n = 25 WT♀ and 27 GA♀). To assay male choice, a single focal WT♂ or GA♂ was given a choice of two females, one WT♀ and one GA♀ (n = 27 WT♂ and 18 GA♂). Experiments were started using 0 day-old sexually unreceptive females and 10–12 days-old sexually mature males. Newly emerged (0 day-old) females were used to avoid the disruption of introducing a sexually mature female into the bioassay. B. germanica females become sexually receptive at 5–7 days of age, so the mating behavior of the focal insect was video-recorded for several days until they mated. Fertility of mated females was evaluated by the number of offspring produced. We assessed the gustatory phenotype of nymphs (either WT-type or GA-type) to determine which of the two adult cockroaches mated with the focal insect. Each gravid female was maintained individually in a container (95 × 95 × 80 mm) with food and water until the eggs hatched. Two day-old first instar nymphs were starved for one day without water and food, and then they were tested in Two-choice feeding assays using 1000 mmol l−1 glucose-containing agar with 0.5 mmol l−1 blue food dye vs. plain sugar-free agar with 1 mmol l−1 red food dye. If all the nymphs chose the glucose-containing agar, their parents were considered WT♂ and WT♀. When all the nymphs showed glucose-aversion, they were raised to the adult stage. Newly emerged adults were backcrossed with WT cockroaches, and their offspring were tested in the Two-choice assay. When the parents were both GA, 100% of the offspring exhibited glucose-aversion. When the parents were WT and GA, the offspring showed a 1:1 ratio of glucose-accepting and glucose-aversive behavior. Mate choice, mating success ratio and the number of offspring were analyzed statistically.We conducted no-choice mating assay using the WT and GA strains (Fig. 1b, d). A female and a male were placed in a Petri dish with fresh water and a piece of rodent food and video-recorded for 24 h. The females were 5–7 days-old and males were 10–12 days-old. Four treatment pairs were tested: WT♂ × WT♀ (n = 20, 18 and 14 pairs for 5, 6 and 7 day-old females, respectively); GA♂ × GA♀ (n = 23, 22 and 35 pairs); GA♂ × WT♀ (n = 21, 14 and 17 pairs); and WT♂ × GA♀ (n = 33, 19 and 15 pairs).To confirm that gustatory stimuli guide nuptial feeding, we artificially augmented the male nuptial secretion and assessed whether the duration of nuptial feeding and mating success of GA♀ were affected (Fig. 2c). Before starting the mating assay with 5 day-old GA♀, 10–12 days-old WT♂ were separated into three groups: A control group did not receive any augmentation; A water control group received distilled water with 1 mmol l−1 blue dye (+Blue); A fructose group received 3000 mmol l−1 fructose solution with blue dye (+Blue+Fru). Approximately 50 nl of the test solution was placed into the tergal gland reservoirs using a glass microcapillary. No-choice mating assays were carried out for 24 h. n = 20–25 pairs for each treatment.We evaluated the association of short nuptial feeding (Fig. 1c) and the GA trait we conducted no-choice mating assays using females from the recombinant lines (Fig. 3c). Before starting each mating assay with 4 day-old females from the WT, GA and recombinant lines (WT_aa, GA_AA and GA_Aa), the EC50 for glucose was obtained by the instantaneous Acceptance-Rejection assay using 0, 10, 30, 100, 300, 1000 and 3000 mmol l−1 glucose (WT♀ and WT_aa♀, non-starved; GA♀, GA_AA♀ and GA_Aa♀, 1-day starved). After the Acceptance-Rejection assay, GA_Aa♀ were separated into two groups according to their sensitivity for rejecting glucose; the GA_Aa_high sensitivity group rejected glucose at 100 and 300 mmol l−1, whereas the GA_Aa_low sensitivity group rejected glucose at 1000 and 3000 mmol l−1. We paired these females with 10–12 days-old WT♂ (n = 15 WT_aa♀, n = 20 GA_AA♀, n = 20 GA_Aa_high♀ and n = 17 GA_Aa_low♀).Feeding bioassayWe conducted two feeding assays: Acceptance-Rejection assay and Consumption assay. The Acceptance-Rejection assay assessed the instantaneous initial responses (binary: yes-no) of cockroaches to tastants, as previously described7,22,27. Briefly, acceptance means that the cockroach started drinking. Rejection means that the cockroach never initiated drinking. The percentage of positive responders was defined as the Number of insects accepting tastants/Total number of insects tested. The effective concentration (EC50) for each tastant was obtained from dose-response curves using this assay. The Consumption assay was previously described27. Briefly, we quantified the amount of test solution females ingested after they started drinking. Females were observed until they stopped drinking, and we considered this a single feeding bout.We used the Acceptance-Rejection assay and Consumption assay, respectively, to assess the sensitivity of 5 day-old WT♀ and GA♀ for accepting and consuming the WT♂ nuptial secretion (Fig. 2a, b). The secretion was diluted with HPLC-grade water to 0.001, 0.01, 0.03, 0.1, 0.3 and 1 male-equivalents/µl (n = 20 non-starved females each). The amount of nuptial secretion consumed was tested at 0.1 male-equivalents/µl in the Consumption assay (n = 10 each).The Acceptance-Rejection assay was used to calculate the effective concentration (EC50) of glucose for females in the WT, GA and recombinant lines (Fig. 3a, b). A glucose concentration series of 0.1, 1, 10, 100 and 1000 mmol l−1 was tested with one-day starved 4-day old females (n = 65 GA_Aa♀, n = 50 GA_AA♀ and n = 50 GA♀) and non-starved females (n = 50 WT_aa♀ and n = 16 WT♀).The effects of female saliva on feeding responses of 5 day-old WT♀ and GA♀ were tested using the Acceptance-Rejection assay (Fig. 4a). Freshly collected saliva of WT♀ and GA♀ was immediately used in experiments. Assays were prepared as follows: 3 µl of 200 mmol l−1 maltose or maltotriose were mixed with 3 µl of either HPLC-grade water or saliva of WT♀ or GA♀. The final concentration of each sugar was 100 mmol l−1 in a total volume of 6 µl. This concentration represented approximately the acceptance EC70 for WT♀ and GA♀27. Nuptial secretion (1 µl representing 10 male-equivalents) was mixed with 1 µl of either HPLC-grade water or saliva from WT♀ or GA♀, and 8 µl of HPLC-grade water was added to the mix. The final concentration of the nuptial secretion was 1 male-equivalent/µl in a total volume of 10 µl. This concentration also represented approximately the acceptance EC70 for WT♀ and GA♀ (Fig. 2a). The mix of saliva and either sugar or nuptial secretion was incubated for 300 s at 25 °C. Additionally, we tested the effect of only saliva in the Acceptance-Rejection assay. Either 1-day starved or non-starved females were tested with water only and then a 1:1 mixture of saliva and water. Saliva alone did not affect acceptance or rejection of stimuli. n = 20–33 females from each strain.To evaluate whether salivary enzymes are involved in the hydrolysis of oligosaccharides, the contribution of salivary glucosidases was tested using the glucosidase inhibitor acarbose in the Acceptance-Rejection assay (Fig. 4b), as previously described27. We first confirmed that the range of 0–125 mmol l−1 acarbose in HPLC-grade water did not disrupt the acceptance and rejection of tastants. Test solutions were prepared as follows: 2 µl of either HPLC-grade water or saliva of GA♀ was mixed with 1 µl of either 250 µmol l−1 of acarbose or HPLC-grade water, then the mixture was added to 1 µl of 400 mmol l−1 of either maltose or maltotriose solution. The total volume was 4 µl, with the final concentration of sugar being 100 mmol l−1. For assays with nuptial secretion, 1 µl of either HPLC-grade water or saliva from 5 day-old GA♀ was mixed with 0.5 µl of either 250 µmol l−1 of acarbose or HPLC-grade water. This mixture was added to 0.5 µl of 10 male-equivalents of nuptial secretion (i.e., 20 male-equivalents/µl). HPLC-grade water was added for a total volume of 10 µl and a final concentration of 1 male-equivalent/µl. The mix of saliva and either sugars or nuptial secretion was incubated for 5 min at 25 °C. All test solutions contained blue food dye. Test subjects were 5 day-old GA♀ and 20–25 females were tested in each assay.Nuptial secretion and saliva collectionsThe nuptial secretion of WT♂ was collected by the following method: Five 10–12 days-old males were placed in a container (95 × 95 × 80 mm) with 5 day-old GA♀. After the males displayed wing-raising courtship behavior toward the females, individual males were immediately decapitated and the nuptial secretion in their tergal gland reservoirs was drawn into a calibrated borosilicate glass capillary (76 × 1.5 mm) under the microscope. The nuptial secretions from 30 males were pooled in a capillary and stored at −20 °C until use. Saliva from 5 day-old WT♀ and GA♀ was collected by the following method: individual females were briefly anesthetized with carbon dioxide under the microscope and the side of the thorax was gently squeezed. A droplet of saliva that accumulated on the mouthparts was then collected into a microcapillary (10 µl, Kimble Glass). Fresh saliva was immediately used in experiments.GC-MS procedures for analysis of sugarsStandards of D-( + )-glucose (Sigma-Aldrich), D-( + )-maltose (Fisher Scientific) and maltotriose (Sigma-Aldrich) were diluted in HPLC-grade water (Fisher Scientific) at 10, 50, 100, 500 and 1000 ng/µl to generate calibration curves. Samples were vortexed for 20 s and a 10 μl aliquot of each sample was transferred to a Pyrex reaction vial containing a 10 μl solution of 5 ng/μl sorbitol (≥98%) in HPLC-grade water as internal standard and dried under a gentle flow of N2 for 20 min.Samples containing degradation products from nuptial secretions were prepared by adding 15 μl of HPLC-water to each sample in a 1.5 ml Eppendorf tube, vortexed for 30 s and centrifuged at 8000 rpm (5223 RCF) for 5 min to separate lipids from the water layer. The water phase was transferred to a reaction vial using a glass capillary. This procedure was repeated with the remaining lipid layer and the water layers were combined in the same reaction vial containing 10 μl of a solution of 5 ng/μl sorbitol and dried under N2 for 20 min.For derivatization of sugars and samples, each reaction vial received 12 μl of anhydrous pyridine under a constant N2 flow, then vortexed and incubated at 90 °C for 5 min. Three μl of N-methyl-N-(trimethylsilyl)trifluoroacetamide (MSTFA; Sigma-Aldrich) was added to each reaction vial and centrifuged at 1000 rpm (118 RCF) for 2 min. Vials were incubated in a heat block at 90 °C for 1.5 hr and vortexed every 10 min for the first 30 min of incubation.The total volume of sample was ~10 μl, and 1 μl was injected into the GC-MS (6890 GC coupled to a 5975 MS, Agilent Technologies, Palo Alto, CA). The inlet was operated in splitless mode (17.5 psi) at 290 °C. The GC was equipped with a DB-5 column (30 m, 0.25 mm, 0.25 μm, Agilent), and helium was used as the carrier gas at an average velocity of 50 cm/s. The oven temperature program started at 80 °C for 1 min, increased at 10 °C/min to 180 °C, then increased at 5 °C/min to 300 °C, and held for 10 min. The transfer line was set at 250 °C for 24 min, ramped at 5 °C/min to 300 °C and held until the end of program. The ion source operated at 70 eV and 230 °C, while the MS quadrupole was maintained at 200 °C. The MSD was operated in scan mode, starting after 9 min (solvent delay time) with a mass range of 33–650 AMU.For GC-MS data analysis, the sorbitol peak area was obtained from the extracted ion chromatograms with m/z = 205, the sorbitol base peak. The area of peaks of glucose, maltose and maltotriose were obtained from the extracted ion chromatograms using m/z = 204, the base peak of the three sugars. The most abundant peaks of each sugar were selected for quantification36, and these peaks did not coelute with other peaks. Then, the peak areas of the three sugars were divided by the area of the respective sorbitol peak in each sample to normalize the data and to correct technical variability during sample processing. This procedure was performed to obtain the calibration curves and quantification of sugars in our experiments.The results of sugar analysis using GC-MS are reported in Supplementary Figs. 1–4.Analysis of nuptial secretionsWe focused the GC-MS analysis on glucose, maltose and maltotriose in WT♂ nuptial secretion (Fig. 4c). To quantify the time-course of saliva-catalyzed hydrolysis of WT♂ nuptial secretion to glucose, 1 µl of GA♀ saliva was mixed with 1 µl of 10 male-equivalents/µl. We incubated the mixtures for 0, 5, 10 and 300 s at 25 °C, and added 4 µl of methanol to stop the enzyme activity (n = 5 each treatment). Each sample contained the nuptial secretions of 5 males to obtain enough detectable amount of sugars. For the statistical analysis, the amounts of sugars were divided by 5 to obtain the amount of sugars in 1 male (1 male-equivalent). These amounts were also used for generating Fig. 4c and Supplementary Table 9. In calculations of the concentration of the three sugars (mmol l−1), the mass and volume of the nuptial secretion were measured using 70–130 male-equivalents of undiluted secretion of each strain (n = 3). The mass and volume of the nuptial secretion/male, including both lipid and aqueous layers, were approximately 30–50 µg and 40–50 nl. Because it was difficult to separate the lipid layer from the water layer at this small scale, we roughly estimated that the tergal reservoirs of the four cockroach lines had 30 nl of aqueous layer that contained sugars.To quantify the time-course of saliva-catalyzed hydrolysis of maltose and maltotriose to glucose, 1 µl of GA♀ saliva was mixed with 1 µl of 200 mmol l−1 of either maltose or maltotriose (Fig. 4d, e). Incubation time points were 0, 5, 10 and 300 s at 25 °C and methanol was used to stop the enzyme activity. Controls without saliva were also prepared using HPLC-grade water instead of saliva and 300 s incubations. n = 5 for each treatment.PhotomicroscopyThe photographs of the tergal glands and mouthparts (Fig. 5) were obtained using an Olympus Digital camera attached to an Olympus CX41 microscope (Olympus America, Center Valley, PA).Statistics and reproducibilityThe sample size and number of replicates for each experiment are noted in the respective section describing the experimental details. In summary, the samples sizes were: Mating bioassays, n = 18–80; Feeding assays, n = 16–65; Sugar analysis, n = 5; Life history parameters, n  > 14. All statistical analyses were conducted in R Statistical Software (v4.1.0; R Core Team 2021) and JMP Pro 15.2 software (SAS Institute Inc., Carey, NC). For bioassay data and sugar analysis data, we calculated the means and standard errors, and we used the Chi-square test with Holm’s method for post hoc comparisons, t-test, and ANOVA followed by Tukey’s HSD test (all α = 0.05), as noted in each section describing the experimental details, results, and in Supplementary Tables 1–11.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Study areaHarbin is located in the centre of Northeast Asia, between 44°04’–46° 40′ N and 125° 42′–130° 10′ E24,26. The site has a mid-temperate continental monsoon climate, with an average annual temperature of 3.6° C and an average annual precipitation is 569.1 mm. The main precipitation months being from June to September, accounting for about 60% of the annual precipitation, the main snow months are from November to January24,25. The overall topography is high in the east and low in the west, with mountains and hills predominating in the east and plains predominating in the west27. In this study, we identified the central district of Harbin, where urban construction activities are frequent and the population is dense, as the study area. According to the “Harbin City Urban Master Plan (2011–2020)” (revised draft in 2017), the specific scope includes Daoli District, Daowai District, Nangang District, Xiangfang District, Pingfang District, Songbei District’s administrative district, Hulan District, and Acheng District part of the area, with a total area of 4187 km2 (Fig. 2). The blue–green space in this study included woodland, grassland, cultivated land, wetland and water that permeate inside and outside the construction sites. They all have integrated functions such as ecology, supply, beautification, culture, and disaster prevention and avoidance, and have a decisive influence on the urban ecological environment.Figure 2Schematic of study area. The Figure is created using ArcGIS ver.10.2 (https://www.esri.com/).Full size imageData sourcesThe data used in this research included the following: land-cover date (30 m × 30 m) of two periods (2011, 2020) spported by the China Geographic National Conditions Data Cloud Platform (http://www.dsac.cn/), Meteorological datasets (1 km × 1 km) were obtained from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http:∥www.resdc.cn/), including air temperature, precipitation, and surface runoff. ASTER GDFM elevation data (30 m × 30 m) came from the Geospatial Data Cloud (http:∥www.gscloud.cn), from which the slope was extracted. Soil data (1 km × 1 km) were from the World Soil Database (HWSD) China Soil Data Set (v1.1). The normalized difference vegetation index (NDVI) and modified normalized difference water index (MNDWI) data (30 m × 30 m) came from the National Comprehensive Earth Observation Data Sharing Platform (http://www.chinageoss.org/), ET datasets (30 m × 30 m) were drawn from the NASA-USGS (https://lpdaac.usgs.gov/). Social and economic data were mainly obtained through the Harbin statistical yearbook and the Harbin social and economic bulletin.Framework of urban blue–green space LEH assessmentUrban blue–green space is a politically defined man-land coupling region composed of ecological, economic, and social systems, which is greatly disturbed by human activities11. The essence of urban blue–green space LEH is that the landscape ecological function sustainably meets human needs28,29. The landscape ecological function reflects the value orientation of human beings to blue–green space, and to a large extent affects the blue–green landscape ecological pattern and process. The interaction between the blue–green landscape ecological pattern and process drives the overall dynamics of blue–green space. Meanwhile, presenting certain landscape ecological function characteristics, which provide ecological support for various human activities30,31,32. While the pattern and process of blue–green space both profoundly influence and are influenced by human activities33,34. This influence is long-term, the standard of LEH should not be fixed in real-time health, but should fully consider the sustainability of the health state.In summary, the landscape ecological pattern, process, function, and sustainability are not separate, but a complex of mutual integration, and organic unity. In this study, we constructed an integrated assessment framework of blue–green space LEH that included four units: pattern, process, service, and sustainability (Fig. 3). In the assessment framework, the LEH of urban blue–green space involves two dimensions: the first is the health status of the urban blue–green space itself, emphasizing the maintenance of the ecological conditions, thereby potentially satisfying a series of diversity goals. The other is that urban blue–green space, as a part of social and economic development, could sustainably provide the ability to meet (subject) needs and goals.Figure 3Key units, interactions of urban blue–green space LEH.Full size imageLandscape ecological patternThe landscape ecological pattern of urban blue–green space is a spatial mosaic combination of landscape elements at different levels or the same level. Affected by human activities interference31, the landscape ecological pattern shows the changing trend of landscape structure complexity, landscape type diversification, and landscape fragmentation. The assessment of urban landscape ecological pattern should be a comprehensive reflection of this changing trend1. Landscape pattern indexes are the most frequently applied which could reflect the structural composition and spatial configuration characteristics of the landscape4,35. This study took landscape ecology as the entry point and selected the landscape pattern indexes that can quantitatively reflect the change characteristics of landscape structural composition and spatial configuration under the disturbance. In this way, the landscape disturbance index (U), landscape connectivity index (CON), and landscape adaptability index (LAI) were used as the indexes for the assessment of landscape ecological pattern health.

(1)

Landscape disturbance index (U)

There are two kinds of relationships between the landscape ecological pattern and the external disturbance: compatibility and conflict. As the landscape ecological pattern has accommodating characteristics, the disturbance beyond the accommodating capacity will degrade the landscape ecological pattern36,37. The landscape disturbance index (U) could characterize the degree of fragmentation, dispersion, and morphological changes in landscape pattern38. The index is a comprehensive index that can reflect the health of the landscape pattern by quantifying the ability of ecosystems to accommodate external disturbances. It consists of the landscape fragmentation index, the inverse of the fractional dimension, and the dominance index. They measure the response of the landscape pattern to external disturbance from the perspective of different landscape types, the same landscape type, and landscape diversity, respectively36,38, and their weights were determined by the entropy weight method. The formula is as follows:$$ U = alpha N_{{{Fi}}} + bD_{{{Fi}}} + cD_{{{Oi}}} $$
(1)
where NFi is the landscape fragmentation index, DFi is the inverse of the fractional dimension, DOi is the dominance index, and a, b, and c are the corresponding weights, which were 0.20, 0.5, and 0.3 in this study, respectively.

(2)

Landscape connectivity index (CON)

The most direct result of landscape ecological pattern degradation caused by external disturbance is that the flow of energy, material, and information among ecological patches is reduced or even blocked, ultimately the stability of the landscape pattern is decreased. The connectivity could characterize the ability of landscape ecological pattern to mitigate risk transmission, which is significant for the dynamic stability of landscape ecological pattern39,40. The landscape connectivity index (CON) could measure the connectivity between ecosystem components through the aggregation or dispersion trend of patches41. The better the connectivity, the stronger the stability of landscape ecological pattern. The formula is as follows:$$ CON = frac{{100sumlimits_{s = 1}^{q} {sumlimits_{h ne l}^{p} {C_{{{shl}}} } } }}{{sumlimits_{s = 1}^{s} {left[ {q_{{s}} (q_{{s}} – 1)/2} right]} }} $$
(2)
where qs is the number of plaques of patch type s, Cshl is the link between patch h and patch l in s within the delimited distance.

(3)

Landscape Restorability Index (LRI)

The ability to recover to its original structure when subjected to disturbances is an important criterion for the landscape ecological pattern42. Research confirmed that the restorability of the landscape ecological pattern is closely related to the structure, function, diversity, and uniformity of distribution. The landscape restorability index (LRI) combines the above landscape information and could indicate the restorability of the landscape ecological pattern in response to disturbance43. The index consists of the patch density, Shannon diversity index, and the landscape evenness, the patch density is the number of patches per square kilometer. The Shannon diversity index reflects the change in the proportion of landscape types. The landscape evenness index shows the distribution evenness of patches in terms of area. The larger the LRI index, the more complex and evenly distributed the structure is, and the more recovery ability of the landscape pattern against disturbance is. The formula is as follows:$$ LRI = PD times SHDI times SHEI $$
(3)
where PD is the patch density, SHDI is the Shannon diversity index, and SHEI is the landscape evenness index.Landscape ecological processThe landscape ecological process of urban blue–green space is extremely complex for it involves multiple factors such as natural ecology, economy, and culture. Landscape ecological process assessment is the measure of the self-organized capacity and the efficiency of ecological processes within and among patches44. A blue–green space with a healthy landscape ecological process should have the ability to adapt to conventional land use under human management and maintain physiological integrity while maintaining the balance of ecological components. Specifically, the landscape ecological process could quickly restore its balance after severe disturbances, with strong organization, suitability, recoverability, and low sensitivity45,46. A single model hardly to gets good research on landscape ecological process under the urban scale. The comprehensive application of multidisciplinary methods is effective means to solve the problem. Regarding this, we selected ecological indexes and models from four aspects: organization, suitability, restoration, and sensitivity to assess the landscape ecological process of urban blue–green space.

(1)

Organization index (O)

The organization of the landscape ecological process is the maintenance ability of stable and orderly material cycling and energy flow within and between landscapes47. The normalized vegetation index (NDVI) and the modified normalized difference water index (MNDWI) could reflect the efficiency and order of ecological processes. Such as accumulation of organic matter, fixation of solar energy, nutrient cycling, regeneration, and metabolism13. The indexes are the external performance of the internal dynamics and organizational capabilities of the ecological process. In recent years, it has been widely used in the assessment of related to landscape ecological process. The formulas are as follows:$$ NDVI = frac{NIR – R}{{NIR + R}} $$$$ MNDWI = frac{p(green) – p(MIR)}{{p(green) + p(MIR)}} $$
(4)
where (NDVI) is the normalized vegetation index, (MNDWI) is the modified water body index, (NIR) is the reflectance value in the near-infrared band, (R) is the reflectance value in the visible channel, (p(green)) and (p(MIR)) are the normalized values in the green and mid-infrared bands.

(2)

Suitability index (Q)

The suitability of the landscape ecological process is a measurement of the self-regulating ability of the landscape ecosystem. That is, to effectively maintain the ecological process in a state of being protected from disturbance during the occasional changes caused by the external environment2. The water conservation amount index (Q) can measure the operating capacity of ecosystems to maintain ecological balance, water conservation, climate regulation, and other ecological processes by integrating the water balance of rainfall, surface runoff, and evaporation41. It could reflect the suitability of landscape ecological process to regional environment and developmental conditions. The formula is as follows:$$ Q = R – J – ET $$
(5)
where Q is the water conservation amount, R is the annual rainfall, J is the surface runoff, ET is the evapotranspiration.

(3)

Recoverability index (ECO)

The recoverability of the landscape ecological process refers to the ability of an ecosystem to return to its original operating state after being subjected to external impacts. Land-use types play an essential role in landscape ecological recoverability48. The ecological recoverability index (ECO) uses the resilience coefficients of land-use types to reflect the level of ecosystem resilience38. Based on previous studies, the resilience coefficient of land-use types was assigned (Table 1).

(4)

Sensitivity index(A)

Table 1 Resilience coefficients of different land use types.Full size tableThe sensitivity index (A) could be used to indicate landscape ecological process formation, change, and vulnerability to disturbance31. We started from the physical effects of blue–green space on sand production, water confluence, and sediment transport, introduced the Soil Erosion Modulus to characterize the sensitivity of landscape ecological processes to disturbance. The index effectively combines landscape ecology, erosion mechanics, soil science, and sediment dynamics49. The formula is as follows:$$ begin{gathered} A = R_{{i}} cdot K cdot LS cdot C cdot P hfill \ L = (l/22.1)^{m} hfill \ S = left{ begin{gathered} 10.8sin theta + 0.03,theta < 5^{ circ } hfill \ 16.8sin theta - 0.50,5^{ circ } le theta < 10^{ circ } hfill \ 21.9sin theta - 0.96,theta ge 10^{ circ } hfill \ end{gathered} right. hfill \ C = left{ begin{gathered} 1,c = 0 hfill \ 0.6508 - 0.3436lg c,0 < c le 78.3% hfill \ 0,c > 78.3% hfill \ end{gathered} right. hfill \ end{gathered} $$
(6)
where A is the soil erosion modulus. Ri is the rainfall erosion factor, K is the soil erosion factor, L and S are slope the length factor and the slope factor respectively, C is the vegetation coverage and management factor, P is the soil and water conservation factor, l is the slope length value, m is the slope length index, and θ the is slope value.Landscape ecological functionThe landscape ecological function determines the ability of ecological service50,51,52, the ecological service of urban blue–green space depends on the human value orientation48. It includes four categories: supply, support, regulation, and culture. Based on Maslow’s Hierarchy of Needs and Alderfer’s ERG theory, scholars have summarized the three major needs of human beings for urban blue–green space. Namely, securing the living environment to meet the survival needs, improving social relationships to meet the interaction needs, and cultivating cultural cultivation to meet the development needs53. Specifically corresponding to the landscape ecological function of urban blue–green space, supply is not the main function, only plays a subsidiary role, support is the basic guarantee, regulation is the basic need for urban environmental construction, and culture is an important element of high-quality social life. Ecosystem service value (ESV) can realize the measurement of ecological service function by calculating the specific value of life support products and services produced by the ecosystem54,55,56. Considering the human value orientation of the urban blue–green space landscape ecological function, the weights were given by consulting 16 experts, with supply, regulation, support, and culture weights of 0.2, 0.3, 0.3, 0.2, respectively. The formula is as follows:$$ ESV = sumlimits_{k = 1}^{n} {S_{k} times V_{k}^{{}} } $$
(7)
where Sk is the area of landscape type k, Vk is the value coefficient of the ecosystem service function of landscape type k .Landscape ecological sustainabilityWu (2013) proposed a research framework for landscape sustainability based on a summary of related studies, stating that landscape ecological sustainability is the ability to provide ecosystem services in a long-term and stable manner34. The framework emphasized that landscape sustainability should focus on the analysis of ecosystem service trade-offs effect34,57. In the process of dynamic change of urban blue–green space ecosystem, there are complex trade-offs among various ecosystem services. This is important for promoting the optimal overall benefits of various ecosystem services and achieving sustainable development of urban ecology58. In addition, as a special type of human-centered ecosystem developed by humans based on nature, human well-being is also very important for the landscape ecological sustainability of urban blue–green space. For this reason, we introduced ecosystem service trade-offs (EST) and ecological construction input (IEC) as assessment indexes of landscape ecological sustainability.

(1)

Ecosystem service trade-offs (EST)

This study applied the root mean square deviation of ecological services to quantify ecosystem service trade-offs (EST). The index could effectively measure the average difference in standard deviation between individual ecosystem services and the average ecosystem services. It is a simple and effective way to evaluate the trade-offs among ecosystem services. The formula is as follows:$$ EST = sqrt {frac{1}{n – 1}sumnolimits_{i = 1}^{n} {(ES_{std} – overline{ES}_{std} } } )^{2} $$
(8)
where ESstd is the normalized ecosystem services, n is the number of ecosystem services , and (overline{ES}_{std}) is the mean value of normalized ecosystem services.

(2)

Ecological construction input (ECI)

Human well-being is a premise for the landscape ecological sustainability of urban blue–green spaces, it is closely related to government investment in ecological construction planning34. From the perspective of economics, this study assessed the human well-being obtained by urban blue–green space with the ratio of urban ecological construction investment to GDP, that is, the ecological construction input (ECI). The formula is as follows:$$ ECI = EI/G $$
(9)
where EI is the amount of ecological construction investment, and G is the gross regional product.Evaluation methodThe index weight determines its relative importance in the index system, and the selection of the weight calculation method in the decision-making of multi-attribute problems has an important impact on the assessment results21. Traditional weighting methods can be divided into two categories, subjective weighting method and objective weighting method21,38. The subjective weighting method is represented by the analytic hierarchy process (AHP), Delphi method, and so on. It has the advantage of simplicity, but the disadvantage is too subjective and randomness because it was completely dependent on the knowledge and experience of decision makers. The objective weighting method is represented by the entropy weighting method (EWM), principal component analysis, variation coefficient method, and so on. And it has been widely recognized for reflecting the variability of assessment results18, but the values of indexes have significant influence and the calculation results are not stable. Considering the limitations of the single weighting method, the weights of each assessment index in this study were determined by the combination of subjective weight and objective weight. Among them, the subjective weighting selected the AHP, and the objective weighting selected the EWM (Table 2). The formula is as follows:$$ w_{{j}} = alpha w_{{j}}^{{{AHP}}} + (1 – alpha )w_{{j}}^{{{EWM}}} $$
(10)
$$ w_{{j}}^{{{EWM}}} = d_{{j}} /sumlimits_{i = 1}^{m} {d_{{j}} } $$
(11)
$$ d_{{j}} = 1 – e_{{j}} $$
(12)
$$ e_{{j}} = – ksumlimits_{i = 1}^{n} {f_{{{ij}}} ln (f_{{{ij}}} )} ,;k = 1/ln (n) $$
(13)
$$ f_{{{ij}}} = X^{prime}_{{{ij}}} /sumlimits_{i = 1}^{n} {X^{prime}_{{{ij}}} } $$
(14)
where (W_{{j}}^{{}}) is the combined weight. (W_{{j}}^{{_{AHP} }}) is the weight of the j-th index of the AHP, (W_{{j}}^{{{EWM}}}) is the weight of the j-th index of the EWM, dj is the information entropy of the j-th index, ej is the entropy value of the j-th index, (f_{{{ij}}}) is the proportion of the index value of the j-th sample under the i-th indexm, (X^{prime}_{{{ij}}}) is the standardized value of the i-th sample of the j-th index, m is the number of index, n is the number of samples, and (alpha) was taken as 0.5.Table 2 Weight of assessment index.Full size tableSince the dimensions of indexes are different, it is necessary to unify the dimensions of the index to avoid the errors caused by direct calculation to make the evaluation results inaccurate. The range standardization was used to normalize the index data and bound its value in the interval [0, 1], the range standardization can be expressed as follows15,23:$$ {text{Positive indicator}}left( + right):A_{{{ij}}} = (X_{{{ij}}} – X_{{{jmin}}} )/(X_{{{jmax}}} – X_{{{jmin}}} ) $$
(15)
$$ {text{Negative indicator}}left( – right):A_{{{ij}}} = (X_{{{jmax}}} – X_{ij} )/(X_{{{jmax}}} – X_{{{jmin}}} ) $$
(16)
Additionally, we divided the LEH index into five levels from high to low using an equal-interval approach as follows40: [1–0.8) healthy, [0.8–0.6) sub-healthy, [0.6–0.4) moderately healthy, [0.4–0.2) unhealthy, [0.2–0] pathological, corresponding level I–V. And the level transfer of LEH in different periods was divided into three types: improvement type, degradation type, and stabilization type. For example, III-II means that the transfer from level III to level II is the improvement type.Spatial autocorrelation analysisSpatial autocorrelation analysis is one of the basic methods in theoretical geography. It could deeply investigate the spatial correlation characteristics of data, including global spatial autocorrelation and local spatial autocorrelation23. The global spatial autocorrelation uses global Moran’s I to evaluate the degree of their spatial agglomeration or differentiation of an attribute value in the study area. The local spatial autocorrelation is a decomposed form of the global spatial autocorrelation18,21, including four types: HH(High-High), LL(Low-Low), HL(High-Low), LH(Low–High). In this study, spatial autocorrelation analysis was applied to study the spatial correlation characteristics of blue–green space LEH. The calculation formulas are as follows:$$ I = frac{{Nsumlimits_{i} {sumlimits_{v} {W_{iv} (Y_{i} – overline{Y} )(Y_{v} – overline{Y} )} } }}{{(sumlimits_{i} {sumlimits_{v} {W_{iv} } } )sumlimits_{i} {(Y_{i} – overline{Y} )} }} $$
(17)
$$ I_{i} = frac{{Y_{i} – overline{Y} }}{{S_{x}^{2} }}sumlimits_{v} {left[ {W_{iv} (Y_{i} – overline{Y} )} right]} $$
(18)
where N is the number of space units, (W_{iv}) is the spatial weight, (Y_{i} ,Y_{v}) are the variable attribute values of the area (i,v), (overline{Y}) is the variable mean, (S_{x}^{2}) is the variance, (I) is the global Moran’s I index, and (I_{i}) is the local Moran’s I index. 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

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

1500-year stable carbon and oxygen isotopes in larch tree-ring celluloseThe δ13Ccell (Fig. 1a, Fig. S2) and δ18Ocell (Fig. 1b, Fig. S3) records span 516–2016 CE, at annual resolution. The δ13Ccell timeseries shows mostly increasing trends during the first millennium of the Common Era (516–1120 CE), and similarly at the end of the last millennium (1720–2016 CE). The maximum δ13Ccell value occurs in 2016 CE (−19.6‰; + 3.2σ), while the minimum occurs in 686 CE (−24.7‰, −3.6σ) relative to the average for the period 516–2016 CE (−22.04‰) (Table S2, Fig. S2). The standard error (SE) for the whole analysed period is 0.02.Figure 1Annually resolved δ13Ccell (a) and δ18O cell (b) in Siberian larch tree-ring cellulose chronologies for the period from 516 to 2016 CE. Chronologies are smoothed by a 101-year Hamming window to highlight a centennial scale. The dotted and dashed lines indicate the number of trees analysed.Full size imageThe δ18Ocell timeseries (Fig. 1b, Fig. S3) showed two positive and one negative extreme over the past 1500 years, with the minimum value (19.9‰; −6.3σ), occurring in 536 CE, and maximum values (31.9‰; + 3.8σ and 32.2‰; + 4.4σ), occurring in 1266 and 2008 CE, respectively (Table S2, Fig. S3). The SE for the whole analysed period is 0.03. The δ18Ocell data has higher standard deviation (SD) (1.15) than δ13Ccell (0.75).Less than 1% of values in the δ18Ocell record are classified as extreme, with the standard deviation ≥  ± 3σ. The δ13Ccell and δ18Ocell records are significantly correlated (r = 0.1, p = 0.0001, n = 1500).Local climate signals preserved in δ13Ccell and δ18Ocell recordsWe used weather observations from the local Mugur-Aksy weather station (50°N, 90°E, 1850 m asl) (Table S1) to derive quantitative paleoclimatic reconstructions from our δ13Ccell and δ18Ocell timeseries. A multiple linear regression analysis revealed significant correlations between δ13Ccell and July precipitation (r = −0.58; p  More

Virus propagation and enumerationEchovirus-11 (E11, Gregory strain, ATCC VR737) and Coxsackievirus-A9 (CVA9, environmental strain from sewage, kindly provided by the Finnish National Institute for Health and Welfare) stocks were produced by infecting sub-confluent monolayers of BGMK cells as described previously [7]. Viruses were released from infected cells by freezing and thawing the culture flasks three times. To eliminate cell debris, the suspensions were centrifuged at 3000 × g for 5 min. Each stock solution was stored at −20 °C until use. Infectious virus concentrations were enumerated by a most probable number (MPN) infectivity assay as described in the Supplementary Information. The assay limit of detection (LoD), defined as the concentration corresponding to one positive cytopathic effect in the lowest dilution of the MPN assay under the experimental conditions used, corresponding to 2 MPN/mL.Inactivation of enteroviruses by bacterial consortia from lake waterTo study the inactivation of CVA9 and E11 by a bacterial consortium from lake water, four surface water samples were collected from Lake Geneva (Ecublens, Switzerland) during the summer 2021. Each sampling event was conducted on warm and sunny days, to minimize biological variation. Immediately after sampling, large particles of the sample were removed by filtering 500 mL of water on a 8 μm nitrocellulose filter membrane (Merck Millipore, Cork, Ireland). The sample was then filtered through a 0.8 μm nitrocellulose filter membrane (Merck Millipore) to remove large microorganisms such as protists. The resulting water sample corresponds to the bacterial fraction used to study virus inactivation.For inactivation experiments, each virus was spiked into individual 1 mL aliquots of fractionated lake water to a final concentration of 106 MPN/mL, and samples were incubated for 48 h at 30 °C without shaking. Duplicate experiments were conducted for each virus and each lake water sample. Experiments to control for thermal inactivation were conducted using the same procedure but by replacing the fractionated lake water with sterile milliQ water. Viral infectivity at times 0 h and 48 h was determined by MPN as described above. Virus decay was calculated as log10 (C/C0), where C is the residual titer after 48 h of incubation, and C0 is the initial titer. The experimental LoD was approximately 5-log10.These same experiments were conducted for three new water samples in the presence of four protease inhibitors with the following final concentrations: E64—10 μM (E3132, Sigma–Aldrich, Saint-Louis, MO, USA), GM6001—4 μM (CC1010, Sigma–Aldrich), Chymostatin—100 μM (C7268, Sigma–Aldrich), and PMSF—100 μM (P7626, Sigma–Aldrich). Each inhibitor was added to 1 mL of fractionated lake water, vortexed for 30 seconds, and incubated at room temperature for 15 min, before adding the two viral strains under the same conditions as described above.Bacterial isolation, cultivation, and storageBacteria were isolated from two water samples from Lake Geneva’s Ecublens beach, taken in November 2019 (Fall, 89 isolates) and May 2020 (Spring, 47 isolates). Bacteria recovery was performed on R2A agar plate (BD Difco, Franklin Lakes, NJ, USA) as described previously [15]. Briefly, successive dilutions from 10−1 to 10−5 were carried out in sterile water for each sample. For each dilution, a volume of 1 mL was deposited on three separate R2A plates, before being incubated at 22, 30, and 37 °C. After 5 days of incubation, each colony was picked and enriched on a new R2A plate. To ensure purity, each isolate was successively plated five times on R2A plate and incubated at the same temperature as the initial isolation. Each purified isolate was cryopreserved in R2A / 20% glycerol at −80 °C. The isolates were named based on the water body (Lake (L)), isolation temperature, and the isolation order (L-T°C-number).Bacterial identificationThe identification of each isolate was performed by 16 S rRNA gene sequencing using the pair of primers 27 F (5’- AGA GTT TGA TCM TGG CTC AG- 3’, Microsynth AG, Balgach, Switzerland) / 786 R (5’- CTA CCA GGG TAT CTA ATC – 3’, Microsynth AG), following a methodology previously described [15]. The thermocycling conditions and the purification of PCR products are described in the Supplementary Information. The complete list of isolated bacteria and associated accession numbers is given in Supplementary Table 1.Phylogenetic inference and metadata visualizationThe consensus from 16 S rRNA gene sequences of the 136 isolates was aligned using the MUSCLE algorithm [16]. The phylogenetic analysis of 566 bp aligned sequences from the V2-V4 16 S rRNA gene regions (Positions: 152–717) was performed using Molecular Evolutionary Genetics Analysis X software [17]. Phylogeny was inferred by maximum likelihood, with 1000 bootstrap iterations to test the robustness of the nodes. The resulting tree was uploaded and formatted using iTOL [18].Virus incubation with bacterial isolatesFor the preparation of the bacteria before co-incubation, each one was first cultured on R2A agar for 48 h at their initial isolation temperature. Overnight suspensions of each bacterial isolate were grown in R2A broth at room temperature under constant agitation (180 rpm). For co-incubation experiments, 200 μL of each bacterial suspension were mixed with 100 μL of a 105 MPN/mL stock of E11 or CVA9. Then, each condition was supplemented with 600 μL of R2A broth. Incubation was carried out for 96 h at room temperature, without shaking. At the end of the co-incubation, each tube was centrifuged for 15 min at 9000 × g (4 °C) to eliminate bacteria, and the residual infectious viral titer was enumerated by MPN assay as described above [7]. Each co-incubation experiment was carried out in triplicate. Control experiments were performed under the same conditions but using sterile R2A. Virus decay was quantified as log10 (Cexp/Cctrl), where Cexp is the residual titer after a co-incubation for 96 h, and Cctrl is the titer after incubation of the virus in sterile R2A for 96 h. The experimental LoD was 3-log10.Protease activity measurement using casein and gelatin agar platesCasein agar was prepared as follows: 20 g of skim milk (BD Difco), supplemented with 1 g glucose were reconstituted with 200 mL of distilled water. Likewise, a 10% bacteriological agar solution was prepared in a final volume of 200 mL. Finally, a solution consisting of 0.8% NaCl, 0.02% KCl, 0.144% Na2HPO4, and 0.024% KH2PO4 was reconstituted in 600 mL of water. All solutions were autoclaved for 15 min at 110 °C. The solutions were mixed, and 25 mL were poured into each Petri dish. Gelatin agar was composed of 0.4% peptone, 0.1% yeast extract, 1.5% gelatin and 1.5% bacteriological agar. The mixture was autoclaved 15 min at 120 °C, and 25 mL of medium was poured into each Petri dish.For each isolate, an overnight suspension was performed in R2A broth at room temperature, before spotting 15 μL of each suspension at the center of both gelatin and casein agar plates. Each plate was incubated at 22, 30, or 37 °C for 72 h, depending on the initial isolation temperature of the bacteria. Casein-degrading activity (cas), which is exerted by many different protease classes, and gelatin-degrading activity (gel), which is mostly caused by MMPs, were revealed by a hydrolysis halo around the producing bacteria. Hydrolysis diameters were measured in millimeters (mm) to report the extent of the proteolytic effect of each strain on both substrates.Protease activity quantification in cell-free supernatantUsing the same bacterial suspensions as for bacterial/virus co-incubation, 200 μL of each suspension was inoculated into 600 μL of R2A broth and incubated without shaking for 96 h at room temperature. Each culture was centrifuged for 15 min at 9000 × g at 4 °C. The resulting cell-free supernatants (CFS) were stored at −20 °C until use. For each CFS, protease activity was measured using the Protease Activity Assay Kit (ab112152, Abcam, Cambridge, UK), which measures general protease activity (pgen) except MMPs, and the MMP Activity Assay Kit (ab112146, Abcam), which selectively measures MMP activity (mmp). Briefly, for the Protease Activity Assay kit, 50 μL of the substrate was added into each well of a dark-bottom plate containing 50 μL of each CFS. Standard trypsin provided by the kit was used as a positive control. For the MMP Activity Assay kit, 50 μL of each CFS was incubated with 50 μL of 2 mM APMA for 3 h at 37 °C, prior to the activity test. Collagenase I (C0130, Sigma–Aldrich) was used as a positive control. R2A broth was used as a negative control for each assay. Protease activity was measured at time 0 and after 60 min, using a Synergy MX fluorescence reader (BioTek). The excitation and emission wavelengths were set to 485 and 530 nm, respectively. The emitted fluorescence, generated by proteolytic cleavage of the substrate of each kit, was calculated as follows: ∆RFU = RFU (60 min) − RFU (0 min). Proteolytic activity was calculated in mmol/min/μL based on the emitted fluorescence measured for trypsin and collagenase I at known proteolytic activities.Data analysisStatistical analyses to compare inactivation data were performed by one-way t-test or one-way ANOVA with Dunnett’s post-hoc test in GraphPad Prism v.9. An alpha value of 0.05 was used as a threshold for statistical significance. For each dataset we confirmed that data were normally distributed.To analyze a potential correlation between protease activity and viral decay, the decay values for each virus strain was related to the four protease activity tests of this study using a scatterplot combined with a Kernel density estimation. The analyses were performed with R v.3.6.1 using the SmoothScatter function of the R Base package.A Left-Censored Tobit model (CTM) with mixed effects was chosen to investigate interactions between protease activity and the decay measured for each virus strain. Briefly, the CTM with mixed effect was chosen for three reasons: (1) The protocol used to measure viral decay had a limit of quantification of −3-log10, and 152 measurement points reached the detection limit, requiring the use of this value as the left-censored value of the model; (2) The two virus strains used in the study showed distinct responses after exposure to environmental bacteria, preventing the use of a multiple linear regression model; (3) Among biological replicates of co-incubation experiments, inactivation variability was observed, suggesting the concomitant action of random biological effects (e.g., production of other compounds than proteases by bacteria, or differences in protease production rate between replicates for each bacterial isolate). The resulting statistical model was then formulated as follows:$$log left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right) = ; beta _0 + beta _1;{rm I}_{{{{{{{{mathrm{virus}}}}}}}}_i = 2} + beta _2sqrt {left[ {pgen} right]_i} + beta _3sqrt {left[ {mmp} right]_i} + beta _4sqrt {left[ {cas} right]_i} \ + beta _5sqrt {left[ {gel} right]_i} + beta _6I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {pgen} right]_i} + beta _7I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {mmp} right]_i} \ + beta _8I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {cas} right]_i} + beta _9I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {gel} right]_i} + alpha _{{{{{{{{mathrm{id}}}}}}}}_i} + varepsilon _i$$$${{{mbox{where}}}}; log left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right) = left{ {begin{array}{*{20}{c}} { – 3} & {{{{{{{{mathrm{if}}}}}}}};{{{{{{{mathrm{log}}}}}}}}left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right) le – 3} \ {{{{{{{{mathrm{log}}}}}}}}left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right)} & {{{{{{{{mathrm{otherwise}}}}}}}}} end{array}} right.$$$$alpha _{{{{{{{{mathrm{id}}}}}}}}_i}sim {{{{{{{mathrm{i}}}}}}}}.{{{{{{{mathrm{i}}}}}}}}.;{{{{{{{mathrm{d}}}}}}}}.;{rm N}left( {0,;sigma _{{{{{{{{mathrm{id}}}}}}}}}^2} right)$$$${{{{{{{mathrm{for}}}}}}}};i in left{ {1,2, ldots } right}$$for which β0 defines the model intercept, (beta _1{rm I}_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}) corresponds to the main effect of the virus factor on the viral decay, (beta _2,;beta _3,;beta _4,;{{{{{{{mathrm{and}}}}}}}};beta _5) corresponds to the main effects of the different protease activity measurements on viral decay, (beta _6I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2},;beta _7I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2},;beta _8I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2},{{{{{{{mathrm{and}}}}}}}};beta _9I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}) corresponds to the interaction effects between each of these variables and the viral decay, (alpha _{{{{{{{{mathrm{id}}}}}}}}_i}) corresponds to the mixed effect of the model and (varepsilon _i) corresponds to the error term of the model. The selection of the model is further detailed in the Supplementary Information (Supplementary Material and Figs. S1 and S2).The full dataset included in the correlation analysis and the CTM is provided in Supplementary Table 2. A description of the variables used is given in the Supplementary Information. The dataset was analyzed using the censReg package in R [19]. The R code is given in the Supplementary Information. More