Ecology

1.Perry, A. L., Low, P. J., Ellis, J. R. & Reynolds, J. D. Climate change and distribution shifts in marine fishes. Science 308, 1912–1915 (2005).ADS 
CAS 
PubMed 
Article 

Google Scholar 
2.Nye, J. A., Link, J. S., Hare, J. A. & Overholtz, W. J. Changing spatial distribution of fish stocks in relation to climate and population size on the Northeast United States continental shelf. Mar. Ecol. Prog. Ser. 393, 111–129 (2009).ADS 
Article 

Google Scholar 
3.Fodrie, F. J., Heck, K. L. Jr., Powers, S. P., Graham, W. M. & Robinson, K. L. Climate-related, decadal-scale assemblage changes of seagrass-associated fishes in the northern Gulf of Mexico. Glob. Change Biol. 16, 48–59 (2010).ADS 
Article 

Google Scholar 
4.Pinsky, M. L., Worm, B., Fogarty, M. J., Sarmiento, J. L. & Levin, S. A. Marine taxa track local climate velocities. Science 341, 1239–1242 (2013).ADS 
CAS 
PubMed 
Article 

Google Scholar 
5.Free, C. M. et al. Impacts of historical warming on marine fisheries production. Science 363, 979–983 (2019).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
6.Simpson, S. D. et al. Continental shelf-wide response of a fish assemblage to rapid warming of the sea. Curr. Biol. 21, 1565–1570 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
7.Pecl, G. T. et al. Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being. Science 355 (2017).8.Poloczanska, E. S. et al. Responses of marine organisms to climate change across oceans. Front. Mar. Sci. 3, 62 (2016).Article 

Google Scholar 
9.Oswald, S. A. & Arnold, J. M. Direct impacts of climatic warming on heat stress in endothermic species: seabirds as bioindicators of changing thermoregulatory constraints. Integr. Zool. 7, 121–136 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
10.Boyles, J. G., Seebacher, F., Smit, B. & McKechnie, A. E. Adaptive thermoregulation in endotherms may alter responses to climate change. Integr. Comp. Biol. 51, 676–690 (2011).PubMed 
Article 
PubMed Central 

Google Scholar 
11.Khaliq, I., Hof, C., Prinzinger, R., Böhning-Gaese, K. & Pfenninger, M. Global variation in thermal tolerances and vulnerability of endotherms to climate change. Proc. R. Soc. B Biol. Sci. 281, 20141097 (2014).Article 

Google Scholar 
12.Gibson-Reinemer, D. K., Sheldon, K. S. & Rahel, F. J. Climate change creates rapid species turnover in montane communities. Ecol. Evol. 5, 2340–2347 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
13.Pörtner, H.-O. et al. Climate induced temperature effects on growth performance, fecundity and recruitment in marine fish: developing a hypothesis for cause and effect relationships in Atlantic cod (Gadus morhua) and common eelpout (Zoarces viviparus). Cont. Shelf Res. 21, 1975–1997 (2001).ADS 
Article 

Google Scholar 
14.Neuheimer, A., Thresher, R., Lyle, J. & Semmens, J. Tolerance limit for fish growth exceeded by warming waters. Nat. Clim. Chang. 1, 110–113 (2011).ADS 
Article 

Google Scholar 
15.Pörtner, H.-O. Climate variations and the physiological basis of temperature dependent biogeography: systemic to molecular hierarchy of thermal tolerance in animals. Comp. Biochem. Physiol. A Mol. Integr. Physiol. 132, 739–761 (2002).PubMed 
Article 
PubMed Central 

Google Scholar 
16.Buckley, L. B., Hurlbert, A. H. & Jetz, W. Broad-scale ecological implications of ectothermy and endothermy in changing environments. Glob. Ecol. Biogeogr. 21, 873–885 (2012).Article 

Google Scholar 
17.Loarie, S. R. et al. The velocity of climate change. Nature 462, 1052–1055 (2009).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
18.Burrows, M. T. et al. The pace of shifting climate in marine and terrestrial ecosystems. Science 334, 652–655 (2011).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
19.Przeslawski, R., Falkner, I., Ashcroft, M. B. & Hutchings, P. Using rigorous selection criteria to investigate marine range shifts. Estuar. Coast. Shelf Sci. 113, 205–212 (2012).ADS 
Article 

Google Scholar 
20.Sunday, J. M. et al. Species traits and climate velocity explain geographic range shifts in an ocean-warming hotspot. Ecol. Lett. 18, 944–953 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
21.MacLeod, C. D. Global climate change, range changes and potential implications for the conservation of marine cetaceans: a review and synthesis. Endanger. Species Res. 7, 125–136 (2009).Article 

Google Scholar 
22.Sydeman, W., Poloczanska, E., Reed, T. & Thompson, S. Climate change and marine vertebrates. Science 350, 772–777 (2015).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
23.Bowen, W. Role of marine mammals in aquatic ecosystems. Mar. Ecol. Prog. Ser. 158, 74 (1997).Article 

Google Scholar 
24.Williams, T. M., Estes, J. A., Doak, D. F. & Springer, A. M. Killer appetites: assessing the role of predators in ecological communities. Ecology 85, 3373–3384 (2004).Article 

Google Scholar 
25.Roman, J. et al. Whales as marine ecosystem engineers. Front. Ecol. Environ. 12, 377–385 (2014).Article 

Google Scholar 
26.Neutel, A.-M., Heesterbeek, J. A. & de Ruiter, P. C. Stability in real food webs: weak links in long loops. Science 296, 1120–1123 (2002).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
27.Bascompte, J., Melián, C. J. & Sala, E. Interaction strength combinations and the overfishing of a marine food web. Proc. Natl. Acad. Sci. U.S.A. 102, 5443–5447 (2005).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
28.Macnab, B. K. The Physiological Ecology of Vertebrates: A View from Energetics (Cornell University Press, 2002).
Google Scholar 
29.Robinson, R. A. et al. Climate change and migratory species (2005).30.Worthy, G. A. & Edwards, E. F. Morphometric and biochemical factors affecting heat loss in a small temperate cetacean (Phocoena phocoena) and a small tropical cetacean (Stenella attenuata). Physiol. Zool., 432–442 (1990).31.Koopman, H. N. Phylogenetic, ecological, and ontogenetic factors influencing the biochemical structure of the blubber of odontocetes. Mar. Biol. 151, 277–291 (2007).Article 

Google Scholar 
32.Adamczak, S. K., Pabst, D. A., McLellan, W. A. & Thorne, L. H. Do bigger bodies require bigger radiators? Insights into thermal ecology from closely related marine mammal species and implications for ecogeographic rules. J. Biogeogr. 47, 1193–1206 (2020).Article 

Google Scholar 
33.Silber, G. K. et al. Projecting marine mammal distribution in a changing climate. Front. Mar. Sci. 4, 413 (2017).Article 

Google Scholar 
34.Kaschner, K., Tittensor, D. P., Ready, J., Gerrodette, T. & Worm, B. Current and future patterns of global marine mammal biodiversity. PLoS ONE 6, e19653 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Salvadeo, C. J., Lluch-Belda, D., Gómez-Gallardo, A., Urbán-Ramírez, J. & MacLeod, C. D. Climate change and a poleward shift in the distribution of the Pacific white-sided dolphin in the northeastern Pacific. Endanger. Species Res. 11, 13–19 (2010).Article 

Google Scholar 
36.Kovacs, K. M., Lydersen, C., Overland, J. E. & Moore, S. E. Impacts of changing sea-ice conditions on Arctic marine mammals. Mar. Biodivers. 41, 181–194 (2011).Article 

Google Scholar 
37.MacLeod, C. D. et al. Climate change and the cetacean community of north-west Scotland. Biol. Cons. 124, 477–483 (2005).Article 

Google Scholar 
38.Higdon, J. W. & Ferguson, S. H. Loss of Arctic sea ice causing punctuated change in sightings of killer whales (Orcinus orca) over the past century. Ecol. Appl. 19, 1365–1375 (2009).PubMed 
Article 

Google Scholar 
39.Evans, P. G. & Hammond, P. S. Monitoring cetaceans in European waters. Mammal Rev. 34, 131–156 (2004).Article 

Google Scholar 
40.Kaschner, K., Quick, N. J., Jewell, R., Williams, R. & Harris, C. M. Global coverage of cetacean line-transect surveys: status quo, data gaps and future challenges. PLoS ONE 7, e44075 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
41.Taylor, B. L., Martinez, M., Gerrodette, T., Barlow, J. & Hrovat, Y. N. Lessons from monitoring trends in abundance of marine mammals. Mar. Mamm. Sci. 23, 157–175 (2007).Article 

Google Scholar 
42.Pyenson, N. D. The high fidelity of the cetacean stranding record: insights into measuring diversity by integrating taphonomy and macroecology. Proc. R. Soc. B Biol. Sci. 278, 3608–3616 (2011).Article 

Google Scholar 
43.Leeney, R. H. et al. Spatio-temporal analysis of cetacean strandings and bycatch in a UK Wsheries hotspot. Biodivers. Conserv. 17, 2323–2338 (2008).Article 

Google Scholar 
44.Lambert, E. et al. Quantifying likely cetacean range shifts in response to global climatic change: implications for conservation strategies in a changing world. Endanger. Species Res. 15, 205–222 (2011).Article 

Google Scholar 
45.Poloczanska, E. S. et al. Global imprint of climate change on marine life. Nat. Clim. Chang. 3, 919–925 (2013).ADS 
Article 

Google Scholar 
46.Pershing, A. J. et al. Slow adaptation in the face of rapid warming leads to collapse of the Gulf of Maine cod fishery. Science 350, 809–812 (2015).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
47.Nawojchik, R., St. Aubin, D. J. & Johnson, A. Movements and dive behavior of two stranded, rehabilitated long-finned pilot whales (Globicephala melas) in the northwest Atlantic. Mar. Mammal Sci. 19, 232–239 (2003).Article 

Google Scholar 
48.Bloch, D. et al. Short-term movements of long-finned pilot whales Globicephala melas around the Faroe Islands. Wildl. Biol. 9, 47–8 (2003).Article 

Google Scholar 
49.Hayes, S., Josephson, E., Maze‐Foley, K. & Rosel, P. US Atlantic and Gulf of Mexico marine mammal stock assessments–2019. NOAA Tech Memo NMFS‐NE 264 (2020).50.Gannon, D., Read, A., Craddock, J., Fristrup, K. & Nicolas, J. Feeding ecology of long-finned pilot whales Globicephala melas in the western North Atlantic. Mar. Ecol. Prog. Ser. Oldendorf 148, 1–10 (1997).ADS 
Article 

Google Scholar 
51.Harden Jones, F. R. In Animal migration. Soc. Exp. Biol. Sem. Ser. 13 (ed. Aidley, D. J.) 139–165 (Cambridge Univ. Press, 1981).
Google Scholar 
52.Alerstam, T., Hedenström, A. & Åkesson, S. Long-distance migration: evolution and determinants. Oikos 103, 247–260 (2003).Article 

Google Scholar 
53.Heide-Jørgensen, M. P. et al. Diving behaviour of long-finned pilot whales Globicephala melas around the Faroe Islands. Wildl. Biol. 8, 307–313 (2002).Article 

Google Scholar 
54.Baird, R. W., Borsani, J. F., Hanson, M. B. & Tyack, P. L. Diving and night-time behavior of long-finned pilot whales in the Ligurian Sea. Mar. Ecol. Prog. Ser. 237, 301–305 (2002).ADS 
Article 

Google Scholar 
55.Adamczak, S. K., McLellan, W. A., Read, A. J., Wolfe, C. L. & Thorne, L. H. The impact of temperature at depth on estimates of thermal habitat for short‐finned pilot whales. Mar. Mammal Sci. (2020).56.Jorda, G. et al. Ocean warming compresses the three-dimensional habitat of marine life. Nat. Ecol. Evolut. 4, 109–114 (2020).Article 

Google Scholar 
57.Burrows, M. T. et al. Ocean community warming responses explained by thermal affinities and temperature gradients. Nat. Clim. Chang. 9, 959–963 (2019).ADS 
Article 

Google Scholar 
58.Kleisner, K. M. et al. The effects of sub-regional climate velocity on the distribution and spatial extent of marine species assemblages. PLoS ONE 11, e0149220 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
59.Kleisner, K. M. et al. Marine species distribution shifts on the US Northeast Continental Shelf under continued ocean warming. Prog. Oceanogr. 153, 24–36 (2017).ADS 
Article 

Google Scholar 
60.Kavanaugh, M. T., Rheuban, J. E., Luis, K. M. & Doney, S. C. Thirty-three years of ocean benthic warming along the US northeast continental shelf and slope: Patterns, drivers, and ecological consequences. J. Geophys. Res. Oceans 122, 9399–9414 (2017).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
61.Grady, J. M. et al. Metabolic asymmetry and the global diversity of marine predators. Science 363 (2019).62.Williams, T. M. et al. The diving physiology of bottlenose dolphins (Tursiops truncatus). III. Thermoregulation at depth. J. Exp. Biol. 202, 2763–2769 (1999).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
63.Pabst, D. A., Rommel, S. A. & McLELLAN, W. A. The emergence of whales 379–397 (Springer, 1998).Book 

Google Scholar 
64.McNab, B. K. Short-term energy conservation in endotherms in relation to body mass, habits, and environment. J. Therm. Biol 27, 459–466 (2002).Article 

Google Scholar 
65.Yeates, L. C. & Houser, D. S. Thermal tolerance in bottlenose dolphins (Tursiops truncatus). J. Exp. Biol. 211, 3249–3257 (2008).PubMed 
Article 
PubMed Central 

Google Scholar 
66.Saba, V. S. et al. Enhanced warming of the Northwest Atlantic Ocean under climate change. J. Geophys. Res. Oceans (2016).67.Kenney, R. D., Scott, G. P., Thompson, T. J. & Winn, H. E. Estimates of prey consumption and trophic impacts of cetaceans in the USA northeast continental shelf ecosystem. J. Northwest Atl. Fish. Sci. 22, 155–171 (1997).Article 

Google Scholar 
68.Read, A. J. & Brownstein, C. R. Considering other consumers: fisheries, predators, and Atlantic herring in the Gulf of Maine. Conserv. Ecol. 7, 2 (2003).
Google Scholar 
69.Overholtz, W. & Link, J. Consumption impacts by marine mammals, fish, and seabirds on the Gulf of Maine-Georges Bank Atlantic herring (Clupea harengus) complex during the years 1977–2002. ICES J. Mar. Sci. J. Conseil 64, 83–96 (2007).Article 

Google Scholar 
70.Smith, L. A., Link, J. S., Cadrin, S. X. & Palka, D. L. Consumption by marine mammals on the Northeast US continental shelf. Ecol. Appl. 25, 373–389 (2015).PubMed 
Article 

Google Scholar 
71.Estes, J. A., Tinker, M. T., Williams, T. M. & Doak, D. F. Killer whale predation on sea otters linking oceanic and nearshore ecosystems. Science 282, 473–476 (1998).ADS 
CAS 
PubMed 
Article 

Google Scholar 
72.Pace, M. L., Cole, J. J., Carpenter, S. R. & Kitchell, J. F. Trophic cascades revealed in diverse ecosystems. Trends Ecol. Evol. 14, 483–488 (1999).CAS 
PubMed 
Article 

Google Scholar 
73.Myers, R. A., Baum, J. K., Shepherd, T. D., Powers, S. P. & Peterson, C. H. Cascading effects of the loss of apex predatory sharks from a coastal ocean. Science 315, 1846–1850 (2007).ADS 
CAS 
PubMed 
Article 

Google Scholar 
74.Durant, J. M., Hjermann, D. Ø., Ottersen, G. & Stenseth, N. C. Climate and the match or mismatch between predator requirements and resource availability. Clim. Res. (CR) 33, 271–283 (2007).ADS 
Article 

Google Scholar 
75.Schweiger, O., Settele, J., Kudrna, O., Klotz, S. & Kühn, I. Climate change can cause spatial mismatch of trophically interacting species. Ecology 89, 3472–3479 (2008).PubMed 
Article 

Google Scholar 
76.Both, C., Van Asch, M., Bijlsma, R. G., Van Den Burg, A. B. & Visser, M. E. Climate change and unequal phenological changes across four trophic levels: constraints or adaptations?. J. Anim. Ecol. 78, 73–83 (2009).PubMed 
Article 

Google Scholar 
77.Evans, K. et al. Periodic variability in cetacean strandings: links to large-scale climate events. Biol. Let. 1, 147–150 (2005).CAS 
Article 

Google Scholar 
78.Overholtz, W., Hare, J. & Keith, C. Impacts of interannual environmental forcing and climate change on the distribution of Atlantic mackerel on the US Northeast continental shelf. Mar. Coastal Fish. 3, 219–232 (2011).Article 

Google Scholar 
79.Roper, C., Lu, C. & Vecchione, M. A revision of the systematics and distribution of Illex species (Cephalopoda: Ommastrephidae). Smithsonian Contrib. Zool., 405–424 (1998).80.Brodziak, J. & Hendrickson, L. An analysis of environmental effects on survey catches of squids Loligo pealei and Illex illecebrosus in the northwest Atlantic. Fish. Bull. 97, 9–24 (1999).
Google Scholar 
81.Henderson, M. E., Mills, K. E., Thomas, A. C., Pershing, A. J. & Nye, J. A. Effects of spring onset and summer duration on fish species distribution and biomass along the Northeast United States continental shelf. Rev. Fish Biol. Fisheries 27, 411–424 (2017).Article 

Google Scholar 
82.Sosebee, K. A. & Cadrin, S. X. A historical perspective on the abundance and biomass of northeast demersal complex stocks from NMFS and Massachusetts inshore bottom trawl surveys, 1963–2002. (2006).83.Brito-Morales, I. et al. Climate velocity can inform conservation in a warming world. Trends Ecol. Evol. 33, 441–457 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
84.Burrows, M. T. et al. Geographical limits to species-range shifts are suggested by climate velocity. Nature 507, 492–495 (2014).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
85.García Molinos, J., Schoeman, D. S., Brown, C. J. & Burrows, M. T. VoCC: an r package for calculating the velocity of climate change and related climatic metrics. Methods Ecol. Evolut. 10, 2195–2202 (2019).Article 

Google Scholar  More

1.Huai, J. J. Dynamics of resilience of wheat to drought in Australia from 1991–2010. Sci. Rep. 7, 9532 (2017).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
2.Li, M., Peterson, C. A., Tautges, N. E., Scow, K. M. & Gaudin, A. C. M. Yields and resilience outcomes of organic cover crop, and conventional practices in a Mediterranean climate. Sci. Rep. 9, 12283 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
3.Keersmaecker, W. D. et al. A model quantifying global vegetation resistance and resilience to short-term climate anomalies and their relationship with vegetation cover. Glob. Ecol. Biogeogr. 24, 539–548 (2015).Article 

Google Scholar 
4.Griffith, G. P. et al. Ecological resilience of Arctic marine food webs to climate change. Nat. Clim. Change 9, 868–872 (2019).ADS 
Article 

Google Scholar 
5.You, N. S., Meng, J. J. & Zhu, L. K. Sensitivity and resilience of ecosystems to climate variability in the semi-arid to hyper-arid areas of Northern China: a case study in the Heihe River Basin. Ecol. Res. 33, 161–174 (2018).Article 

Google Scholar 
6.Reijers, V. C. et al. Resilience of beach grasses along a biogeomorphic successive gradient: resource availability vs. clonal integration. Oceologia https://doi.org/10.1007/s00442-019-04568-w (2019).Article 

Google Scholar 
7.Chambers, J. C. et al. Resilience to stress and disturbance, and resistance to Bromus tectorum L. invasion in clod desert shrublands of western North America. Ecosystems 17, 360–375 (2014).CAS 
Article 

Google Scholar 
8.Driessen, M. M. Fire resilience of a rare, freshwater crustacean in a fire-prone ecosystem and the implications for fire management. Austral Ecol. 44, 1030–1039 (2019).Article 

Google Scholar 
9.Ren, H., Lu, H. F., Li, Y. D. & Wen, Y. G. Vegetation restoration and its research advancement in Southern China. J. Trop. Subtrop. Bot. 27(5), 469–480 (2019).
Google Scholar 
10.Yan, H. M., Zhan, J. Y. & Zhang, T. Review of ecosystem resilience research progress. Prog. Geogr. 31(3), 303–314 (2012).
Google Scholar 
11.Zhan, J. Y., Yan, H. M., Deng, X. Z. & Zhang, T. Assessment of forest ecosystem resilience in Lianhua County of Jiangxi Province. J. Nat. Resour. 27(8), 1304–1315 (2012).
Google Scholar 
12.Pérez-Girón, J. C., Álvarez-Álvarez, P., Díaz-Valera, E. R. & Lopes, D. M. M. Influence of climate variations on primary production indicators and on the resilience of forest ecosystems in a future scenario of climate change: application to sweet chestnut agroforestry systems in the Iberian Peninsula. Ecol. Indic. 113, 106199 (2020).Article 

Google Scholar 
13.Meng, Y. Y. et al. Analysis of ecological resilience to evaluate the inherent maintenance capacity of a forest ecosystem using a dense Landsat time series. Ecol. Inform. 57, 101064 (2020).Article 

Google Scholar 
14.Han, L. et al. Species composition, community structure, and floristic characteristics of desert riparian forest community along the mainstream of Tarim River. Plant Sci. J. 37(3), 324–336 (2019).
Google Scholar 
15.Zhou, H. H. et al. Climate change may accelerate the decline of desert riparian forest in the lower Tarim River, Northwestern China: evidence from tree-rings of Populus euphratica. Ecol. Indic. 111, 105997 (2020).Article 

Google Scholar 
16.Aini, A. et al. Analysis of stakeholders’ cognition on desert riparian forest ecosystem services in the lower reaches of Tarim River, China. Res. Soil Water Conserv. 23(1), 205–209 (2016).
Google Scholar 
17.Li, Y. Q., Chen, Y. N., Zhang, Y. Q. & Xia, Y. Rehabilitating China’s largest inland river. Conserv. Biol. 23(3), 531–536 (2009).PubMed 
Article 

Google Scholar 
18.Dai, J. S. Evaluation of eco-environment and socio-economic benefits on comprehensive reclamation projects on the Tarim River Basin. Doctoral Dissertation of Xinjiang Agricultural University (2015).19.Han, L., Wang, H. Z., Niu, J. L., Wang, J. Q. & Liu, W. Y. Response of Populus euphratica communities in a desert riparian forest to the groundwater level gradient in the Tarim River Basin. Acta Ecol. Sin. 37, 6836–6846 (2017).
Google Scholar 
20.Yang, G. & Guo, Y. P. The change and prospect of vegetation in the end of the lower reaches of Tarim River after ecological water delivery. J. Desert Res. 24(2), 167–172 (2004).
Google Scholar 
21.Yan, H. M., Zhan, J. Y. & Zhang, T. Resilience of forest ecosystems and its influencing factors. Procedia Environ. Sci. 10, 2201–2206 (2011).Article 

Google Scholar 
22.Abenayake, C. C., Mikami, Y., Matsuda, Y. & Jayasinghe, A. Ecosystem service-based composite indicator for assessing community resilience to floods. Environ. Dev. 27, 34–46 (2018).Article 

Google Scholar 
23.Maestas, J. D., Campbell, S. B., Chambers, J. C., Pellant, M. & Miller, R. F. Tapping soil survey information for rapid assessment of sagebrush ecosystem resilience and resistance. Rangelands 38(3), 120–128 (2016).Article 

Google Scholar 
24.Ponce-Campos, G. E. et al. Ecosystem resilience despite large-scale altered hydroclimatic conditions. Nature 494, 349–352 (2013).ADS 
CAS 
PubMed 
Article 

Google Scholar 
25.Frazier, A. E., Renschler, C. S. & Miles, S. B. Evaluating post-disaster ecosystem resilience using MODIS GPP data. Int. J. Appl. Earth Obs. Geoinform. 21, 43–52 (2013).ADS 
Article 

Google Scholar 
26.Kahiluoto, H. et al. Decline in climate resilience of European wheat. PNAS 116(1), 123–128 (2019).CAS 
PubMed 
Article 

Google Scholar 
27.Li, X. Y. et al. Temporal trade-off between gymnosperm resistance and resilience increases forest sensitivity to extreme drought. Nat. Ecol. Evolut. 4, 1075–1083 (2020).Article 

Google Scholar 
28.Li, C. H., Zhou, M., Wang, Y. T., Zhu, T. B., Sun, H., Yin, H. H., Cao, H. J., Han, H. Y. Inter-annual variations of vegetation net primary productivity and their spatial-temporal contribution and climate driving in arid Northwest China: a case study of Hexi Corridor. Chin. J. Ecol. (2020).29.Song, J. et al. A global database of plant production and carbon exchange from global change manipulative experiments. Sci. Data 7, 1–7 (2020).Article 
CAS 

Google Scholar 
30.Yang, G. et al. Research progress of ecosystem resilience assessment. Zhejiang Agric. Sci. 60(3), 508–513 (2019).
Google Scholar 
31.Liu, J. Z. & Chen, Y. N. Analysis on converse succession of plant communities at the lower reaches of Tarim River. Arid Land Geogr. 25(3), 231–236 (2002).
Google Scholar 
32.Chen, X., Bao, A. M., Wang, X. P., Guli, J. P. E. & Huang, Y. Recent ecological effectiveness assessment of integrated management projects in the Tarim River. Bull. Chin. Acad. Sci. 32(1), 20–28 (2017).
Google Scholar 
33.Zhao, H., Yan, L. & Ji, F. The dynamics of land utilization in the upper reaches of Tarim River. J. Arid Land Resour. Environ. 15(4), 40–43 (2001).
Google Scholar 
34.Sun, F., Wang, Y. & Chen, Y. N. Dynamics of desert-oasis ecotone and its influencing factors in the Tarim Basin. Chin. J. Ecol. 39(10), 1–11 (2020).
Google Scholar 
35.Xu, G. H. A genetic explanation of the recent changes of ecological environment in the Tarim River Basin, southern Xinjiang. Xinjiang Meteorol. 28–31 (2005).36.Kamkin, A. & Lozinsky, I. Mechanically Gated Channels and Their Regulation (Springer, 2012).Book 

Google Scholar 
37.Feyisa, K. et al. Effects of enclosure management on carbon sequestration, soil properties and vegetation attributes in East African rangelands. CATENA 159, 9–19 (2017).Article 

Google Scholar 
38.Wang, G. H., Ren, Y. J. & Gou, Q. Q. The changes of soil physical and chemical property during the enclosure process in a typical desert oasis ecotone of the Hexi Corridor in northwestern China. J. Desert Res. 40(2), 222–231 (2020).
Google Scholar 
39.Xu, H. L., Ye, M. & Li, J. M. Changes in groundwater levels and the response of natural vegetation to the transfer of water to the lower reaches of the Tarim River. J. Environ. Sci. 19(10), 1199–1207 (2007).Article 

Google Scholar 
40.Zhang, P. F., Guli, J., Bao, A. M., Meng, F. H. & Guo, H. Ecological effects evaluation for short term planning of the Tarim River. Arid Land Geogr. 40(1), 156–164 (2017).
Google Scholar 
41.Gulimire, H., Wang, G. Y., Zhang, Y., Liu, Q. Q. & Su, L. T. Influence mechanisms of intermittent ecological water conveyance on groundwater level and vegetation in arid land. Arid Land Geogr. 41(4), 726–733 (2018).
Google Scholar 
42.Guo, H. W., Xu, H. L. & Ling, H. B. Study of ecological water transfer mode and ecological compensation scheme of the Tarim River Basin in dry years. J. Nat. Resour. 32(10), 1705–1717 (2017).
Google Scholar 
43.Wu, T. Z., Ding, J., Guan, W. K., Ruan, C. J. & Guan, Y. Populus euphratica forest replacement and photosynthetic characteristics in Tarim Populus euphratica national nature reserve. Prot. For. Sci. Technol. 8, 1–4 (2020).
Google Scholar 
44.Zhu, C. G., Aikeremu, A., Li, W. H. & Zhou, H. H. Ecosystem restoration of Populus euphratica forest under the ecological water conveyance in the lower reaches of Tarim River. Arid Land Geography, 44(3), 629–636 (2021).
Google Scholar 
45.Chen, Y. N. Study on Eco-hydrological Problems of the Tarim River Basin in Xinjiang (Science Press, 2010).
Google Scholar 
46.Halik, U., Aishan, T., Betz, F., Kurban, A. & Rouzi, A. Effectiveness and challenges of ecological engineering for desert riparian forest restoration along China’s largest inland river. Ecol. Eng. 127, 11–22 (2019).Article 

Google Scholar 
47.Xinjiang Morning News. In the past three years, the area of the Populus euphratica forest reserve in the Tarim River Basin has increased by 569.95 km2. https://www.sohu.com/a/308626663_100034331?sec=wd (2019).48.China News Service. Ecological water transfer for desert vegetation in lower reaches of Konqi River in Xinjiang. https://news.sina.com.cn/o/2020-02-22/doc-iimxyqvz4945915.shtml (2020). More

Research cruisesThis dataset consists of sequence data from 4 separate cruises: ARK-XXVII/1 (PS80)—17th June to 9th July 2012; Stratiphyt-II— April to May 2011; ANT-XXIX/1 (PS81)—1st to 24th November 2012 and ANT-XXXII/2 (PS103)—16th December 2016 to 3rd February 2017 and covers a transect of the Atlantic Ocean from Greenland to the Weddell Sea (71.36°S to 79.09°N) (Supplementary Table 1). In order to study the composition, distribution and activity of microbial communities in the upper ocean across the broadest latitudinal ranges possible, samples have been collected during four field campaigns as shown in Fig. 1A. The first collection of samples was collected in the North Atlantic Ocean from April to May 2011 by Dr. Willem van de Poll of the University of Groningen, Netherlands and Dr. Klaas Timmermans of the Royal Netherlands Institute for Sea Research. The second set of samples was collected in the Arctic Ocean from June to July 2012, and the third set of samples was collected in the South Atlantic Ocean from October to November 2012. Both of which were collected by Dr. Katrin Schmidt of the University of East Anglia. The final set of samples was collected in the Antarctic Ocean from December 2016 to January 2017 by Dr. Allison Fong of the Alfred-Wegener Institute for Polar and Marine Research, Bremerhaven, Germany.SamplingWater samples from the Arctic Ocean and South Atlantic Ocean expeditions were collected using 12 L Niskin bottles (Rosette sampler with an attached Sonde (CTD, conductivity, temperature, depth) either at the chlorophyll maximum (10–110 m) and/or upper of the ocean (0–10 m). As soon as the rosette sampler was back on board, water samples were immediately transferred into plastic containers and transported to the laboratory. All samples were accompanied by measurements on salinity, temperature, sampling depth and silicate, nitrate, phosphate concentration (Supplementary Table 1). Water samples were pre-filtered with a 100 μm mesh to remove larger organisms and subsequently filtered onto 1.2 μm polycarbonate filters (Isopore membrane, Millipore, MA, USA). All filters were snap frozen in liquid nitrogen and stored at −80 °C until further analysis.Water samples from the North Atlantic Ocean cruise were also taken with 12 L Niskin bottles attached to a Rosette sampler with a Sonde. However, these samples were filtered onto 0.2 μm polycarbonate filters (Isopore membrane, Millipore, MA, USA) without pre-filtration but snap frozen in liquid nitrogen and stored at −80 °C as the other samples.Water samples from the Southern Ocean cruise were taken with 12 L Niskin bottles attached to an SBE911plus CTD system equipped with 24 Niskin samplers. These samples were filtered onto 1.2 μm polycarbonate membrane filters (Merck Millipore, Germany) in a container cooled to 4 °C and snap frozen in liquid nitrogen and stored at −80 °C as the other samples. Environmental data recorded at the time of sampling can be found in Supplementary Table 1.DNA extractions: Arctic Ocean and South Atlantic Ocean samplesDNA was extracted with the EasyDNA Kit (Invitrogen, Carlsbad, CA, USA) with modification to optimise DNA quantity and quality. Briefly, cells were washed off the filter with pre-heated (65 °C) Solution A and the supernatant was transferred into a new tube with one small spoon of glass beads (425–600 μm, acid washed) (Sigma-Aldrich, St. Louis, MO, USA). Samples were vortexed three times in intervals of 3 s to break the cells. RNase A was added to the samples and incubated for 30 min at 65 °C. The supernatant was transferred into a new tube and Solution B was added followed by a chloroform phase separation and an ethanol precipitation step. DNA was pelleted by centrifugation and washed several times with isopropanol, air dried and suspended in 100 μL TE buffer (10 mM Tris-HCl, pH 7.5, 1 mM EDTA, pH 8.0). Samples were snap frozen in liquid nitrogen and stored at −80 °C until sequencing.DNA extractions: North Atlantic Ocean samplesNorth Atlantic Ocean samples were extracted with the ZR-Duet™DNA/RNA MiniPrep kit (Zymo Research, Irvine, USA) allowing simultaneous extraction of DNA and RNA from one sample filter. Briefly, cells were washed from the filters with DNA/RNA Lysis Buffer and one spoon of glass beads (425–600 μm, Sigma-Aldrich, MO, USA) was added. Samples were vortexed quickly and loaded onto Zymno-Spin™IIIC columns. The columns were washed several times and DNA was eluted in 60 μmL, DNase-free water. Samples were snap frozen in liquid nitrogen and stored at −80 °C until sequencing.DNA extractions: Southern Ocean samplesDNA from the Southern Ocean samples was extracted with the NucleoSpin Soil DNA extraction kit (Macherey‐Nagel) following the manufacturer’s instructions. Briefly, cells were washed from the filters with DNA Lysis Buffer and into a lysis tube containing glass beads was added. Samples were disrupted by bead beating for 2 × 30 s interrupted by 1 min cooling on ice and loaded onto the NucleoSpin columns. The columns were washed three times and DNA was eluted in 50 μL, DNase-free water. Samples were stored at −20 °C until further processing.Amplicon sequencing of 16S and 18S rDNAAll extracted DNA samples were sequenced and pre-processed by the Joint Genome Institute (JGI) (Department of Energy, Berkeley, CA, USA). iTAG amplicon sequencing was performed at JGI with primers for the V4 region of the 16S (FW(515F): GTGCCAGCMGCCGCGGTAA; RV(806R): GGACTACNVGGGTWTCTAAT)49 and 18S (FW(565F): CCAGCASCYGCGGTAATTCC; RV(948R): ACTTTCGTTCTTGATYRA)50. (Supplementary Table 6) rRNA gene (on an Illumina MiSeq instrument with a 2 × 300 base pairs (bp) read configuration51. 18S sequences were pre-processed, this consisted of scanning for contamination with the tool Duk (US Department of Energy Joint Genome Institute (JGI), 2017,a) and quality trimming of reads with cutadapt52. Paired end reads were merged using FLASH53 with a max mismatch set to 0.3 and min overlap set to 20. A total of 54 18S samples passed quality control after sequencing. After read trimming, there was an average of 142,693 read pairs per 18S sample with an average length of 367 bp and 2.8 Gb of data over all samples.16S sequences were pre-processed, this consisted of merging the overlapping read pairs using USEARCH’s merge pairs54 with the parameter minimum number of differences (merge max diff pct) set to 15.0 into unpaired consensus sequences. Any reads that could not be merged are discarded. JGI then applied the tool USEARCH’s search oligodb tool with the parameters mean length (len mean) set to 292, length standard deviation (len stdev) set to 20, primer trimmed max difference (primer trim max diffs) set to 3, a list of primers and length filter max difference (len filter max diffs) set to 2.5 to ensure the Polymerase Chain Reaction (PCR) primers were located with the correct direction and inside the expected spacing. Reads that did not pass this quality control step were discarded. With a max expected error rate (max exp err rate) set to 0.02, JGI evaluated the quality score of the reads and those with too many expected errors were discarded. Any identical sequence was de-duplicated. These are then counted and sorted alphabetically for merging with other such files later. A total of 57 × 16S samples passed quality control after sequencing. There was an average 393,247 read pairs per sample and an average base length of 253 bp for each sequence with a total of 5.6 Gb.RNA extractions: Arctic Ocean and Atlantic samplesRNA from the Arctic and Atlantic Ocean samples was extracted using the Direct-zol RNA Miniprep Kit (Zymo Research, USA). Briefly, cells were washed off the filters with Trizol into a tube with one spoon of glass beads (425–600 μm, Sigma-Aldrich, MO, USA). Filters were removed and tubes bead beaten for 3 min. An equal volume of 95% ethanol was added, and the solution was transferred onto Zymo-Spin™ IICR Column and the manufacturer instructions were followed. Samples were treated with DNAse to remove DNA impurities, snap frozen in liquid nitrogen and stored at −80 °C until sequencing.RNA extractions: Southern OceanRNA from the Southern Ocean samples was extracted using the QIAGEN RNeasy Plant Mini Kit (QIAGEN, Germany) following the manufacturer’s instructions with on-column DNA digestion. Cells were broken by bead beating like for the DNA extractions before loading samples onto the columns. Elution was performed with 30 µm RNase-free water. Extracted samples were snap frozen in liquid nitrogen and stored at −80 °C until sequencing.Metatranscriptome sequencingAll samples were sequenced and pre-processed by the U.S. Department of Energy Joint Genome Institute (JGI). Metatranscriptome sequencing was performed on an Illumina HiSeq-2000 instrument27. A total of 79 samples passed quality control after sequencing with 19.87 Gb of sequence read data over all samples for analysis. This comprised a total of 34,241,890 contigs, with an average length of 503 and an average GC% of 51%. This resulted in 36354419 of non-redundant genes detected.JGI employed their suite of tools called BBTools55 for preprocessing the sequences. First, the sequences were cleaned using Duk a tool in the BBTools suite that performs various data quality procedures such as quality trimming and filtering by kmer matching. In our dataset, Duk identified and removed adaptor sequences, and also quality trimmed the raw reads to a phred score of Q10. In Duk the parameters were; kmer-trim (ktrim) was set to r, kmer (k) was set to 25, shorter kmers (mink) set to 12, quality trimming (qtrim) was set to r, trimming phred (trimq) set to 10, average quality below (maq) set to 10, maximum Ns (maxns) set to 3, minimum read length (minlen) set to 50, the flag “tpe” was set to t, so both reads are trimmed to the same length and the “tbo” flag was set to t, so to trim adaptors based on pair overlap detection. The reads were further filtered to remove process artefacts also using Duk with the kmer (k) parameter set to 16.BBMap55 is another a tool in the BBTools suite, that performs mapping of DNA and RNA reads to a database. BBMap aligns the reads by using a multi-kmer-seed-and-extend approach. To remove ribosomal RNA reads, the reads were aligned against a trimmed version of the SILVA database using BBMap with parameters set to; minratio (minid) set to 0.90, local alignment converter flag (local) set to t and fast flag (fast) set to t. Also, any human reads identified were removed using BBMap.BBmerge56 is a tool in the BBTools suite that performs the merging of overlapping paired end reads (Bushnell, 2017). For assembling the metatranscriptome, the reads were first merged with the tool BBmerge, and then BBNorm was used to normalise the coverage so as to generate a flat coverage distribution. This type of operation can speed up assembly and can even result in an improved assembly quality.Rnnotator52 was employed for assembling the metatranscriptome samples 1–68. Rnnotator assembles the transcripts by using a de novo assembly approach of RNA-Seq data and it accomplishes this without a reference genome52. MEGAHIT57 was employed for assembling the metatranscriptome samples 69–82. The tool BBMap was used for reference mapping, the cleaned reads were mapped to metagenome/isolate reference(s) and the metatranscriptome assembly.Metatranscriptome analysisJGI performed the functional analysis on the metatranscriptomic dataset. JGI’s annotation system is called the Metagenome Annotation Pipeline (MAP) (v4.15.2)27. JGI used HMMER 3.1b258 and the Pfam v3059 database for the functional analysis of our metatranscriptomic dataset. This resulted in 11,205,641 genes assigned to one or more Pfam domain. This resulted in 8379 Pfam functional assignments and their gene counts across the 79 samples. The files were further normalised by applying hits per million.18S rDNA analysisA reference dataset of 18S rRNA gene sequences that represent algae taxa was compiled for the construction of the phylogenetic tree by retrieving sequences of algae and outgroups taxa from the SILVA database (SSUREF 115)60 and Marine Microbial Eukaryote Transcriptome Sequencing Project (MMETSP) database61. The algae reference database consists of 1636 species from the following groups: Opisthokonta, Cryptophyta, Glaucocystophyceae, Rhizaria, Stramenopiles, Haptophyceae, Viridiplantae, Alveolata, Amoebozoa and Rhodophyta. A diagram of the 18S classification pipeline can be found in Supplementary Fig. 1. In order to construct the algae 18S reference database, we first retrieved all eukaryotic species from the SILVA database with a sequence length of  > = 1500 base pairs (bp) and converted all base letters of U to T. Under each genus, we took the first species to represent that genus. Using a custom written script (https://github.com/SeaOfChange/SOC/blob/master/get_ref_seqs.pl), the species of interest (as stated above) were selected from the SILVA database, classified with NCBI taxa IDs and a sequence information file produced that describes each of the algae sequences by their sequence ID and NCBI species ID. Taxonomy from the NCBI database, eukaryote sequences from the SILVA database and a list of algal taxa including outgroups were used as input for the script. This information was combined with the MMETSP database excluding duplications.The algae reference database was clustered to remove closely related sequences with CD-HIT (4.6.1)62 using a similarity threshold of 97%. Using ClustalW (2.1)63 we aligned the reference sequences with the addition of the parameter iteration numbers set to 5. The alignment was examined by colour coding each species to their groups and visualising in iTOL64. It was observed that a few species were misaligning to other groups and these were then deleted using Jalview65. The resulting alignment was tidied up with TrimAL (1.1)66 by applying parameters to delete any positions in the alignment that have gaps in 10% or more of the sequence, except if this results in less than 60% of the sequence remaining. A maximum likelihood phylogenetic reference tree and statistics file based on our algae reference alignment was constructed by employing RaxML (8.0.20)67 with a general time reversible model of nucleotide substitution along with the GAMMA model of rate heterogeneity. For a description of the lineages of all species back to the root in the algae reference database, the taxa IDs were submitted for each species to extract a subset of the NCBI taxonomy with the NCBI taxtastic tool (0.8.4)68 Based on the algae reference multiple sequence alignment, with HMMER3 (3.1B1)69 a Profile HMM was created. A pplacer reference package using taxtastic was generated, which produced an organized collection of all the files and taxonomic information into one directory. With the reference package, a SQLite database was created using pplacer’s Reference Package PReparer (rppr). With hmmalign, the query sequences were aligned to the reference set and created a combined Stockholm format alignment. Pplacer (re-aligned to the reference set and created a combined Stockholm format alignment. Pplacer (1.1)70 was used to place the query sequences on the phylogenetic reference tree by means of the reference alignment according to a maximum likelihood model70 The place files were converted to CSV with pplacer’s guppy tool; in order to easily take those with a maximum likelihood score of  > = 0.5 and counted the number of reads assigned to each classification. This resulted in 6,053,291 reads that were taxonomically assigned being taken for analysis.Normalisation of 18S rDNA gene copy number18S rDNA gene copy number vary widely among eukaryotes. In order to create an estimate of abundances of the species in the samples the data had to be normalised. Previous work has explored the link between copy number and genome size71. However, there is not a single database of 18S rDNA gene copy numbers for eukaryote species. In order to address this, gene copy number and related genome sizes of 185 species across the eukaryote tree was investigated and plotted (Supplementary Fig. 2, Supplementary Table 4)68,71,72,73,74,75,76,77,78,79. Based on the log transformed data, a significant correlation with a R2 of 0.55 with a p-value  More

Macro- and micro-evolutionary integration between jaw complexesWe examined phenotypic associations between the lower oral and pharyngeal jaws (LOJ and LPJ, respectively) of 88 cichlid species from across Africa, primarily sampling from lakes in the East African Rift Valley: lakes Malawi, Tanganyika, and Victoria (Supplementary Data 1). We characterized jaw shapes based on 107 individuals using 3D geometric morphometrics by placing landmarks in positions that capture functionally (e.g., bony processes, sutures, etc.) and developmentally (e.g., distinct cellular origins) important aspects of morphology, including placing mirrored landmarks across midlines to gain symmetric configurations (Fig. 1e, Supplementary Fig. 1). We conducted a Procrustes superimposition, removed the effects of allometry to account for size differences, and then removed the effects of asymmetry to account for developmental noise. We performed a two-block partial least squares (PLS) analysis on the species mean landmark configurations and corrected for phylogenetic non-independence using a Bayesian time-calibrated tree31. We found the LOJ and LPJ were evolutionarily correlated (r-PLS = 0.482, P = 0.002, effect size (Z) = 2.585), but some taxa, particularly those with unique diets and/or modes of feeding, appeared to deviate from the best-fit line, indicating lower levels (or different patterns) of integration between jaws (Fig. 2a). Indeed, we found numerous taxa, typically from Lake Malawi, whereby covariation between the LOJ and LPJ appeared much different relative to other African cichlids. Taxa placed far from the best-fit line either (1) exhibited a specialized feeding morphology to better exploit an foraging niche shared with many competitors (i.e., Labeotropheus, algae; Copadichromis, zooplankton; Taeniolethrinops, insects), or (2) exhibited a specialized feeding morphology to take advantage of a more challenging food source (i.e., Trematocranus, snails). However, not all taxa that consume specialized prey were far from the best-fit line; Pungu, (primarily a sponge-feeder) and Perissodus (a scale-feeder), while exhibiting specialized feeding apparatuses to consume such prey, exhibited a relationship between their LOJ and LPJ that was in-line with other African cichlids (Supplementary Fig. 2). We also noted, that while Malawi cichlids exhibit a range of LOJ-LPJ relationships (from weak to strong), most Tanganyikan cichlids reside close to the best-fit line. However, when we examine the strength of integration in the Tanganyika group (n = 29, r-PLS = 0.698, P = 0.001, Z = 2.954) and Malawi group (n = 40, r-PLS = 0.541, P = 0.020, Z = 2.155), despite Tangyanika cichlids exhibiting higher Z-scores, consistent with stronger integration, a statistical comparison between groups finds no significant difference (Z pairwise = 1.188, P = 0.235). Comparisons between Tanganyikan and Malawi cichlids should not be influenced by sampling bias, as principal components analyses (PCA) on the LOJ and LPJ landmark data (Supplementary Data 2 and 3) showed that our sampling of Tanganyikan cichlids includes many species with extreme morphologies that reside at the outer edges of LOJ and LPJ morphospace (Supplementary Fig. 3). Indeed, cichlids from Lake Tanganyika exhibited similar LOJ morphological disparity (Malawi Procrustes variance (PV) = 0.074; Tanganyika PV = 0.057, P = 0.253) and greater LPJ morphological disparity (Malawi PV = 0.015; Tanganyika PV = 0.023, P = 0.012), relative to cichlids from Lake Malawi. Taken together, this indicates that while Tanganyikan cichlids exhibit comparable (i.e., LOJ), or greater (i.e., LPJ) morphological variation compared to Malawi cichlids, covariation between LOJ and LPJ shapes was generally similar between groups.Fig. 2: Phylogenetic two-block partial least squares analysis to assess macroevolutionary associations between lower oral and pharyngeal jaws.a Jaw shape associations across a broad sample of African cichlids (n = 88). Taxa from Lake Malawi are placed into two groups based on phylogenetic position: an mbuna ‘rock-dwellers’ group, and a non-mbuna group consisting of the utaka ‘sand-dwellers’ alongside other benthic species88. b Jaw shape associations across the Tropheops sp. species complex from across a depth gradient (n = 22). Oral and pharyngeal jaw wireframes denote morphologies at either end of the correlational axis. Source data are provided as a Source Data file.Full size imageWe next investigated the degree of integration at lower taxonomic levels. First, we analyzed the jaws of cichlids within the Tropheops species complex from Lake Malawi that is diverse and known to partition habitat by depth32,33. While Tropheops exhibited strong integration between jaws in on our macroevolutionary assessment, species within this genus occupy a broader niche space. Investigating integration within such a species complex provided an opportunity to understand whether habitat differences could lead to differences in integration between jaw complexes. Using the same landmarking procedure as described above, we characterized shape variation in the LOJs and LPJs of 22 wild-caught Tropheops taxa from 60 individuals, concentrating on members from localities across the southern portion of Lake Malawi (Supplementary Data 4). We again performed a two-block PLS analysis on the mean landmark configurations and accounted for phylogenetic non-independence using an amplified fragment length polymorphism tree33. Again, we found the LOJ and LPJ were evolutionarily correlated (r-PLS = 0.795, P = 0.006, Z = 2.521), indicating jaw integration does not appear to vary by habitat (Fig. 2b).Finally, we measured and compared integration between a species pair that exhibited relatively strong versus weak covariation between LOJ and LPJ shapes in our macroevolutionary assessments, Tropheops sp. ‘red cheek’ (TRC, relatively stronger integration) and Labeotropheus fuelleborni (LF, relatively weaker integration). Using the same landmarking protocol we performed separate two-block PLS analyses between LOJs and LPJs of LF and TRC (Supplementary Data 5). Notably, we found strong and significant integration between jaw complexes in TRC (n = 11, r-PLS = 0.957, P = 0.001, Z = 3.038; Fig. 3a) relative to LF (n = 17, r-PLS = 0.669, P = 0.22, Z = 0.794; Fig. 3b). Further, we found the effect sizes of jaw integration within TRC and LF to be statistically distinct (Z pairwise = 1.678, P = 0.047). Altogether, our data support the assertion that the LOJ and LPJ are evolutionarily integrated at multiple taxonomic levels, but they also appear to indicate that certain taxa, such as Labeotropheus, can more readily generate adaptive morphological variation in each jaw complex independently.Fig. 3: Two-block partial least squares analysis to assess microevolutionary associations between lower oral and pharyngeal jaws.a Shape associations among Tropheops sp. “red cheek” (TRC) individuals (n = 11). b Shape associations among Labeotropheus fuelleborni (LF) individuals (n = 17). c Shape associations among members of a hybrid cross between TRC and LF (n = 409). Oral and pharyngeal jaw wireframes denote morphologies at either end of the correlational axis. Source data are provided as a Source Data file.Full size imageGenetic basis for oral and pharyngeal jaw shape covariationTo understand whether phenotypic covariation between the LOJ and LPJ is genetically determined we performed a quantitative trait loci (QTL) analysis to identify prospective genomic regions involved in jaw shape variation for both the LOJ and LPJ. Specifically, we extended an existing genetic cross between the more strongly integrated TRC and the more weakly integrated LF to the F5 generation. Details of the pedigree may be found in34 and in the supplement. For this experiment, we genotyped 636 F5 hybrids and produced a genetic map containing 812 single-nucleotide polymorphisms (SNPs) spread across 24 linkage groups (Supplementary Data 6). With a total length of 1431 cM, our high-resolution linkage map contained a marker every 1.83 cM, on average, allowing us to leverage the increased number of recombination events that occurred to reach an F5 population. We then characterized LOJ and LPJ shape in 409 F5 hybrids using the same landmarking scheme described above, and performed a two-block PLS analysis. In concordance with our findings from natural populations, we documented an association between jaw complexes in this laboratory pedigree (r-PLS = 0.491, P = 0.001, effect size = 6.189; Fig. 3c).We next performed a PCA on the hybrid landmark configurations to distill the data down to a series of orthogonal axes that best explain shape variation among individuals. We extracted the first two PCs from the LOJ and LPJ as each axis represented more than 10% of the shape variation (Supplementary Data 7; Supplementary Figs. 4 and 5). The first axis of the LOJ reflected more general variation in depth, width, and length of the element (41.8% of variation), while the second axis reflected more specific variation in the length of the ascending arm of the articular––the process for which jaw closing muscles attach (12.7% of variation). The first axis of the LPJ reflected width, length, and wing process size (33.7% of variation), while the second axis reflected depth and the size of the anterior keel – the process for which the pharyngeal jaw pharyngohyoideus muscle attaches and controls jaw adduction (14.2% of variation). We then utilized these PC scores as traits to run in our QTL analyses to investigate the genetic basis for variation in these structures.QTL mapping implicates pleiotropic control of LOJ and LPJ shape variationIntegration between LOJ and LJP shapes in the F5 predicts that this pattern of covariation will be reflected in the genotype-phenotype map. Specifically, we predict that we will find overlapping QTL for both jaws. We used a multiple QTL mapping (MQM) approach to test this prediction. Specifically, we performed QTL scans for all four traits and quantitatively assessed the evidence for significant QTL marker(s) using a permutation procedure that reshuffles the phenotypic data relative to genotypic data 1000 times to generate a null distribution, disassociating any possible relationship between genotype and phenotype, to then compare the empirical distribution against35. Once candidate QTL markers were identified, we calculated an approximate Bayesian credible interval to determine the region in which a potential candidate locus would reside. We uncovered a total of five QTL for LOJ traits, and four QTL for LPJ traits (Fig. 4a; Supplementary Data 8). While most QTL localize to different linkage groups, we also identified some QTL that colocalized. Two traits (LOJ PC1, LPJ PC1) share a marker on LG4, while three traits (LOJ PC1, LOJ PC2, LPJ PC1) colocalized to the same markers on LG7. These data are consistent with pleiotropy on LG7 and possibly LG4.Fig. 4: Genetic analyses to identify regions of the genome responsible for major changes in jaw shape.All plots are based on 409 LFxTRC F5 hybrids. a QTL analysis to identify positions in the genome most associated with each trait. b Pleiotropy analysis on linkage group seven to determine whether the oral jaw PC1 trait colocalizes to the same region as the pharyngeal jaw PC1 trait. Significance was determined using a likelihood ratio test (LLRT). c Pleiotropy analysis on linkage group seven to determine whether the oral jaw PC2 trait colocalizes to the same region as the pharyngeal jaw PC1 trait. Significance was determined using a LLRT. d Fine mapping all traits across the entirety of LG7. Values furthest from 0 reflect the largest differences between hybrids with LF and TRC genotypes at a given marker. We find peak genotype-phenotype association at ~50 mb that coincides with our Bayes credible interval (grey bar). Intervals that surround the average phenotypic effect line denote standard error of the mean. e Fine mapping all traits across the Bayes credible interval. Population level genetic diversity (FST) data are applied to the map (black dots) with the opacity of each SNP dependent on the degree of segregation between LF and TRC, with those falling above an empirical Z-score threshold of 0.6 determined to be significant, and those above 0.9 deemed highly significant (green lines). Within the credible interval there are four SNPs with FST values of 1.0, but a single SNP that falls within a genotype-phenotype peak residing within an intron of dym (black circle). Source data are provided as a Source Data file.Full size imageWe then quantitatively assessed the evidence for pleiotropy using a likelihood ratio test (LLRT) to compare the null hypothesis of a common pleiotropic QTL to the alternative hypothesis that they are affected by separate QTL36,37. The overlap on LG4 at a single marker (43.57 cM) was deemed significant (LLRT = 1.85, P = 0.03), indicating that we can reject the null hypothesis and that these peaks likely represent separate QTL for each trait (Supplementary Fig. 6). The three traits that overlap on LG7 spanned two markers (19.12 cM–28.04 cM) and were all deemed non-significant (LOJpc1-LPJpc1: LLRT = 0.02, P = 0.66, Fig. 4b; LOJpc2-LPJpc1: LLRT = 0.20, P = 0.41, Fig. 4c), leading us to accept the null hypothesis and conclude that this interval likely contains a single pleiotropic QTL. Whether a single gene, or multiple closely linked genes drive this pleiotropic signal requires a fine-mapping approach.Notably, this locus on LG7 has been implicated previously in underlying LOJ and LPJ shape in another Lake Malawi cichlid cross between LF and Maylandia zebra38,39. Maylandia species, like Tropheops, were generally more integrated in our macroevolutionary analysis (Fig. 2a), and thus another cross between LF and a species with higher integration values point to the same locus. This suggests that the genetic mechanism of integration may be conserved.Fine mapping identifies two candidate genes critical for bone formationTo gain insights into which gene(s) may be pleiotropically regulating LOJ and LPJ jaw shape variation on LG7 we constructed a fine map with greater marker density to investigate genotype-phenotype associations with greater resolution. To that end, we anchored QTL intervals to particular stretches of physical sequence of the Maylandia zebra genome40. We then identified additional RAD-seq SNPs across the linkage group of interest and genotyped them in the F5. Based on this we created two fine maps: the first spanned the entirety of LG7 group with an average spacing of around one marker every 490 kb (Supplementary Data 9), the second matched the QTL marker range revealed by the Bayesian credible interval analysis with an average spacing of around one marker every 180 kb (Supplementary Data 10). We also calculated FST from a panel of wild-caught LF (n = 20) and TRC (n = 20), and primarily focused on FST values of 1.0 that would indicate complete segregation of a SNP between LF and TRC. At every marker on our LG7 fine maps, we calculated the difference in the values of our three colocalized traits between those hybrids homozygous for the LF allele and those homozygous for the TRC allele.We identified a small region on LG7 that exhibited large differences in the average phenotypic effect of those hybrids with either LF or TRC alleles. In our full LG7 map we identified a ~2 mb region (46.7 mb–48.7 mb) that peaked in all three traits (Fig. 4d; Supplementary Data 11). Notably, the traces of all three traits across our LG7 fine maps track together in an almost identical fashion. In our fine map that centered on the Bayes credible interval, we found evidence for both large phenotypic effects among all traits, and the presence of several FST markers approaching or equal to 1.0 (Fig. 4e; Supplementary Data 12). One marker combined an FST score of 1.0, indicating complete segregation of that allele between LF and TRC, with high average phenotypic effects across all traits (Supplementary Fig. 7). This SNP fell within an intron of dymeclin (dym), a gene that is necessary for correct organization of Golgi apparatus and controls endochondral bone formation during early development. Dym is critical for chondrocyte development and previous research using the zebrafish model found an expression pattern that spanned the presumptive mandibular and ceratobranchial regions at larval stages41. Mutations in this gene lead to profound effects on the size and shape of bones due to misregulated chondrocyte development42. Just downstream (8 kb, Supplementary Fig. 7) of dym is mothers against decapentaplegic homolog 7 (smad7), an antagonist of both TGF-β and BMP signaling and a suppressor of bone formation. As an inhibitory Smad, smad7 negatively regulates genes from the BMP and TGF-β signaling pathways (i.e., bmp-2, -4, -7, nodal, etc.) that are known to shape phenotypic differences in the craniofacial skeleton across a wide range of taxa including cichlids25,38,43, Geospiza finches44,45, and Anolis lizards46, primarily because these genes have the capacity to influence size in structures of trophic importance such as the mandible47. Both of these genes represent good candidates for controlling shape variation in the LOJ and the LPJ simultaneously. While two of the three traits peak at the dym SNP, when considering markers just outside the credible interval another peak is visible (especially for LOJ PC1) that sits close to notch1a, a gene involved in skeletal development and homeostasis. Notch1a is flanked by two fully segregated FST markers. The upstream marker is around ~60 kb from the promoter region, while the downstream marker resides around ~52 kb away from the gene within an intron of kcnt1, a gene involved in potassium channel development that appears to regulate brain function48. While kcnt1 reflects a poor candidate gene for our analysis, the intronic SNP could act as a distant enhancer of notch1a. Thus, given the combined results from QTL and fine-mapping, dym and smad7 represent strong candidates, but we cannot rule out notch1a.Correlated expression of key genes between LOJ and LPJWe used quantitative real-time PCR (qPCR) to assess and compare the expression levels of dym, smad7, and notch1a in the LOJ and LPJ of three mbuna genera from lake Malawi (Tropheops n = 6, Labeotropheus n = 8, Maylandia n = 8). We used Labeotropheus and Tropheops to complement our quantitative genetic analysis, and all three taxa were represented in our phenotypic assessments of integration, permitting a comparison between macroevolutionary associations of the LOJ and LPJ with the underlying genetic architecture and expression for jaw complex correlation. We collected tissue samples from young juveniles of these four taxa, taking the LOJ and LPJ, alongside the caudal fin to act as an internal control, and performed a phenol/chloroform RNA extraction. We designed primers with high amplification efficiency ( >90%) for our three genes (Supplementary Data 13), and used β-actin as our control gene. We calculated relative expression of the LOJ and LPJ using the 2-ΔΔCT method49, and compared expression across taxa and between tissues (Supplementary Data 14 and 15).We initially compared tissue level expression levels between Labeotropheus and Tropheops and found small differences in dym expression, with LF typically exhibiting slightly higher levels (t-test LOJ t = 2.863, P = 0.014; LPJ t = 1.212, P = 0.249; Fig. 5a). These results are consistent with previous expression studies that demonstrated how Labeotropheus typically has up-regulated bone and collagen markers and as a consequence has greater bone deposition and a more robust craniofacial skeleton50,51. Expression level differences were also noted for notch1a and smad7 (Fig. 5b-c); both showed reduced expression in LF, which is expected based on each genes role as negative regulators of bone formation52,53. While the differences between species were fairly small in smad7 between taxa (t-test LOJ t = −1.869, P = 0.086; LPJ t = −0.359, P = 0.726), they were more notable in notch1a (t-test LOJ t = −1.947, P = 0.080; LPJ t = −3.221, P = 0.009). Notch1a is involved in skeletal remodeling, previous research has shown LF exhibits a minimal plastic response to environmental stimuli51. Thus, the relatively low expression of notch1a in Labeotropheus compared to Tropheops is consistent with this observation. While only representing a single life-history stage, the expression differences between species suggest that all three genes may underlie the development of species-specific shapes for the LOJ and/or LPJ. However, visualizing the data this way cannot speak to whether one or more of these loci underlie the covariation of the jaws.Fig. 5: Comparing expression levels of dym, smad7, and notch1a via qPCR in the oral and pharyngeal jaws.a dym bar plot; (b) notch1a bar plot; (c), smad7 bar plot; (d), dym scatter plot; (e), notch1a scatter plot; (f), smad7 scatter plot. a–c bar plots depict mean relative expression levels, error bars denote standard error. d–f Scatterplots depict relative expression levels of the LOJ and LPJ, error bounds surrounding the linear regression line denote standard error. e inset, linear regression for each genus. Three cichlid taxa were included: Labeotropheus n = 8, Tropheops n = 6, Maylandia n = 8. Bar plot significance determined via t-tests: ●P  More

1.Gause, G. F. & Witt, A. A. Behavior of mixed populations and the problem of natural selection. Am. Nat. 69, 596–609 (1935).Article 

Google Scholar 
2.Hairston, N. G., Smith, F. E. & Slobodkin, L. B. Community structure, population control, and competition. Am. Nat. 94, 421–425 (1960).Article 

Google Scholar 
3.Ayala, F. J. Experimental invalidation of the principle of competitive exclusion. Nature 224, 1076–1079 (1969).ADS 
CAS 
PubMed 
Article 

Google Scholar 
4.Bengtsson, J. Interspecific competition increases local extinction rate in a metapopulation system. Nature 340, 713–715 (1989).ADS 
Article 

Google Scholar 
5.Bolnick, D. I. Intraspecific competition favours niche width expansion in Drosophila melanogaster. Nature 410, 463–466 (2001).ADS 
CAS 
PubMed 
Article 

Google Scholar 
6.Collins, S. Competition limits adaptation and productivity in a photosynthetic alga at elevated CO2. Proc. Biol. Sci. 278, 247–255 (2011).PubMed 

Google Scholar 
7.Osmond, M. M. & de Mazancourt, C. How competition affects evolutionary rescue. Philos. Trans. R. Soc. B 368, 20120085 (2013).Article 

Google Scholar 
8.Birch, L. C. Selection in Drosophila pseudoobscura in relation to crowding. Evolution 9, 389–399 (1955).Article 

Google Scholar 
9.Martin, M. J., Perez-Tome, J. M. & Toro, M. A. Competition and genotypic variability in Drosophila melanogaster. Heredity 60, 119–123 (1988).PubMed 
Article 

Google Scholar 
10.Harvey, J. A., Poelman, E. H. & Tanaka, T. Intrinsic inter- and intraspecific competition in parasitoid wasps. Annu. Rev. Entomol. 58, 333–351 (2013).CAS 
PubMed 
Article 

Google Scholar 
11.Pennacchio, F. & Strand, M. R. Evolution of developmental strategies in parasitic hymenoptera. Annu. Rev. Entomol. 51, 233–258 (2006).CAS 
PubMed 
Article 

Google Scholar 
12.Van Alphen, J. J. & Visser, M. E. Superparasitism as an adaptive strategy for insect parasitoids. Annu. Rev. Entomol. 35, 59–79 (1990).PubMed 
Article 

Google Scholar 
13.Varaldi, J. et al. Infectious behavior in a parasitoid. Science 302, 1930–1930 (2003).CAS 
PubMed 
Article 

Google Scholar 
14.Dorn, S. & Beckage, N. E. Superparasitism in gregarious hymenopteran parasitoids: ecological, behavioural and physiological perspectives. Physiol. Entomol. 32, 199–211 (2007).Article 

Google Scholar 
15.Gandon, S., Rivero, A. & Varaldi, J. Superparasitism evolution: adaptation or manipulation? Am. Nat. 167, E1–E22 (2006).PubMed 
Article 

Google Scholar 
16.Speirs, D. C., Sherratt, T. N. & Hubbard, S. F. Parasitoid diets: does superparasitism pay? Trends Ecol. Evol. 6, 22–25 (1991).CAS 
PubMed 
Article 

Google Scholar 
17.Tracy Reynolds, K. & Hardy, I. C. Superparasitism: a non-adaptive strategy? Trends Ecol. Evol. 19, 347–348 (2004).CAS 
PubMed 
Article 

Google Scholar 
18.Pan, M., Liu, T. & Nansen, C. Avoidance of parasitized host by female wasps of Aphidius gifuensis (Hymenoptera: Braconidae): the role of natal rearing effects and host availability? Insect Sci. 25, 1035–1044 (2018).PubMed 
Article 

Google Scholar 
19.Potting, R. P. J., Snellen, H. M. & Vet, L. E. M. Fitness consequences of superparasitism and mechanism of host discrimination in the stem borer parasitoid Cotesia flavipes. Entomol. Exp. Appl. 82, 341–348 (1997).Article 

Google Scholar 
20.Mackauer, B. B. Influence of superparasitism on development rate and adult size in a solitary parasitoid wasp, Aphidius ervi. Funct. Ecol. 6, 302–307 (1992).Article 

Google Scholar 
21.Keasar, T. et al. Costs and consequences of superparasitism in the polyembryonic parasitoid Copidosoma koehleri (Hymenoptera: Encyrtidae). Ecol. Entomol. 31, 277–283 (2006).Article 

Google Scholar 
22.Silva-Torres, C. S. A., Ramos, I. T., Torres, J. B. & Barros, R. Superparasitism and host size effects in Oomyzus sokolowskii, a parasitoid of diamondback moth. Entomol. Exp. Appl. 133, 65–73 (2009).Article 

Google Scholar 
23.Wylie, H. G. Delayed development of Microctonus vittatae (Hymenoptera: Braconidae) in superparasitized adults of Phyllotreta cruciferae (Coleoptera: Chrysomelidae). Can. Entomol. 115, 441–442 (1983).Article 

Google Scholar 
24.White, J. A. & Andow, D. A. Benefits of self-superparasitism in a polyembryonic parasitoid. Biol. Control 46, 133–139 (2008).PubMed 
PubMed Central 
Article 

Google Scholar 
25.Yamada, Y. Y. & Sugaura, K. Evidence for adaptive self-superparasitism in the dryinid parasitoid Haplogonatopus atratus when conspecifics are present. Oikos 103, 175–181 (2003).Article 

Google Scholar 
26.Varaldi, J., Fouillet, P., Bouletreau, M. & Fleury, F. Superparasitism acceptance and patch-leaving mechanisms in parasitoids: a comparison between two sympatric wasps. Anim. Behav. 69, 1227–1234 (2005).Article 

Google Scholar 
27.Varaldi, J., Patot, S., Nardin, M. & Gandon, S. A virus-shaping reproductive strategy in a Drosophila parasitoid. Adv. Parasitol. 70, 333–363 (2009).PubMed 
Article 

Google Scholar 
28.Carton, Y., Bouletreau, M., van Alphen, J. J. M. & van Lenteren, J. C. The Drosophila parasitic wasps. in The Genetics and Biology of Drosophila (eds Ashburner, M., Carson, H. L. & Thompson, J. N.) 347–394 (Academic Press, 1986).29.Kacsoh, B. Z., Lynch, Z. R., Mortimer, N. T. & Schlenke, T. A. Fruit flies medicate offspring after seeing parasites. Science 339, 947–950 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
30.Krzemien, J. et al. Control of blood cell homeostasis in Drosophila larvae by the posterior signalling centre. Nature 446, 325–328 (2007).ADS 
CAS 
PubMed 
Article 

Google Scholar 
31.Kraaijeveld, A. R. & Godfray, H. C. Trade-off between parasitoid resistance and larval competitive ability in Drosophila melanogaster. Nature 389, 278–280 (1997).ADS 
CAS 
PubMed 
Article 

Google Scholar 
32.Hwang, R. Y. et al. Nociceptive neurons protect Drosophila larvae from parasitoid wasps. Curr. Biol. 17, 2105–2116 (2007).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
33.Mortimer, N. T. et al. Parasitoid wasp venom SERCA regulates Drosophila calcium levels and inhibits cellular immunity. Proc. Natl Acad. Sci. USA 110, 9427–9432 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Huang, J. et al. Two novel venom proteins underlie divergent parasitic strategies between a generalist and a specialist parasite. Nat. Commun. 12, 234 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Martinson, E. O., Mrinalini, Kelkar, Y. D., Chang, C. H. & Werren, J. H. The evolution of venom by co-option of single-copy genes. Curr. Biol. 27, 2007–2013 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Jaffe, A. B. & Hall, A. Rho GTPases: biochemistry and biology. Annu. Rev. Cell. Dev. Biol. 21, 247–269 (2005).CAS 
PubMed 
Article 

Google Scholar 
37.Moon, S. Y. & Zheng, Y. Rho GTPase-activating proteins in cell regulation. Trends Cell Biol. 13, 13–22 (2003).CAS 
PubMed 
Article 

Google Scholar 
38.Xu, J. et al. RhoGAPs attenuate cell proliferation by direct interaction with p53 tetramerization domain. Cell Rep. 3, 1526–1538 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
39.Hinge, A. et al. p190-B RhoGAP and intracellular cytokine signals balance hematopoietic stem and progenitor cell self-renewal and differentiation. Nat. Commun. 8, 14382 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
40.Werner, E. GTPases and reactive oxygen species: switches for killing and signaling. J. Cell Sci. 117, 143–153 (2004).CAS 
PubMed 
Article 

Google Scholar 
41.Bailey, A. P. et al. Antioxidant role for lipid droplets in a stem cell niche of Drosophila. Cell 163, 340–353 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
42.Boguski, M. S. & McCormick, F. Proteins regulating Ras and its relatives. Nature 366, 643–654 (1993).ADS 
CAS 
PubMed 
Article 

Google Scholar 
43.Rittinger, K. et al. Crystal structure of a small G protein in complex with the GTPase-activating protein rhoGAP. Nature 388, 693–697 (1997).ADS 
CAS 
PubMed 
Article 

Google Scholar 
44.Simanshu, D. K., Nissley, D. V. & McCormick, F. RAS proteins and their regulators in human disease. Cell 170, 17–33 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
45.Outreman, Y., Le Ralec, A., Plantegenest, M., Chaubet, B. & Pierre, J. S. Superparasitism limitation in an aphid parasitoid: cornicle secretion avoidance and host discrimination ability. J. Insect Physiol. 47, 339–348 (2001).CAS 
PubMed 
Article 

Google Scholar 
46.Hofsvang, T. Discrimination between unparasitized and parasitized hosts in hymenopterous parasitoids. Acta Entomol. Bohemosl 87, 161–175 (1990).
Google Scholar 
47.van Lenteren, J. C. in Semiochemicals: Their Role in Pest Control (eds Nordlund, D. A., Jones, R. L. & Lewis, W. J.) 153–179 (Wiley and Sons, 1981).48.Ganesalingam, V. K. Mechanism of discrimination between parasitized and unparasitized hosts by Venturia canescens (hymenoptera: Ichneumonidae). Entomol. Exp. Appl. 17, 36–44 (2011).Article 

Google Scholar 
49.Hoffmeister, T. S. & Roitberg, B. D. To mark the host or the patch: decisions of a parasitoid searching for concealed host larvae. Evol. Ecol. 11, 145–168 (1997).Article 

Google Scholar 
50.Agboka, K. et al. Self-, intra-, and interspecific host discrimination in Telenomus busseolae Gahan and T. isis Polaszek (Hymenoptera: Scelionidae), sympatric egg parasitoids of the African cereal stem borer Sesamia calamistis Hampson (Lepidoptera: Noctuidae). J. Insect Behav. 15, 1–12 (2002).Article 

Google Scholar 
51.Liang, Q., Jia, Y. & Liu, T. Self- and conspecific discrimination between unparasitized and parasitized green peach aphid (Hemiptera: Aphididae), by Aphelinus asychis (Hymenoptera: Aphelinidae). J. Econ. Entomol. 110, 430–437 (2017).PubMed 

Google Scholar 
52.Gandon, S., Varaldi, J., Fleury, F. & Rivero, A. Evolution and manipulation of parasitoid egg load. Evolution 63, 2974–2984 (2009).PubMed 
Article 

Google Scholar 
53.Hughes, D. P. & Libersat, F. Neuroparasitology of parasite-insect associations. Annu. Rev. Entomol. 63, 471–487 (2018).CAS 
PubMed 
Article 

Google Scholar 
54.Sberro, H. et al. Large-scale analyses of human microbiomes reveal thousands of small, novel genes. Cell 178, 1245–1259 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
55.Zuzarte-Luis, V. & Mota, M. M. Parasite sensing of host nutrients and environmental cues. Cell Host Microbe 23, 749–758 (2018).CAS 
PubMed 
Article 

Google Scholar 
56.Cox, F. E. G. Parasites affect behavior of mice. Nature 294, 515–515 (1981).ADS 
CAS 
PubMed 
Article 

Google Scholar 
57.Elya, C. et al. Robust manipulation of the behavior of Drosophila melanogaster by a fungal pathogen in the laboratory. eLife 7, e34414 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
58.Hoover, K. et al. A gene for an extended phenotype. Science 333, 1401–140 (2011).ADS 
CAS 
PubMed 
Article 

Google Scholar 
59.Mcallister, M. K. & Roitberg, B. D. Adaptive suicidal-behavior in pea aphids. Nature 328, 797–799 (1987).ADS 
Article 

Google Scholar 
60.Maure, F., Brodeur, J., Droit, A., Doyon, J. & Thomas, F. Bodyguard manipulation in a multipredator context: different processes, same effect. Behav. Process. 99, 81–86 (2013).Article 

Google Scholar 
61.Mohan, P. & Sinu, P. A. Parasitoid wasp usurps its host to guard its pupa against hyperparasitoids and induces rapid behavioral changes in the parasitized host. PLoS ONE 12, e0178108 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
62.Muller, C. B. & Schmidhempel, P. Exploitation of cold temperature as defense against parasitoids in bumblebees. Nature 363, 65–67 (1993).ADS 
Article 

Google Scholar 
63.Noubade, R. et al. NRROS negatively regulates reactive oxygen species during host defence and autoimmunity. Nature 509, 235–239 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
64.Louradour, I. et al. Reactive oxygen species-dependent Toll/NF-κB activation in the Drosophila hematopoietic niche confers resistance to wasp parasitism. eLife 6, e25496 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
65.Sinenko, S. A., Shim, J. & Banerjee, U. Oxidative stress in the haematopoietic niche regulates the cellular immune response in. Drosoph. EMBO Rep. 13, 83–89 (2012).CAS 
Article 

Google Scholar 
66.Wang, Y. et al. Superoxide dismutases: dual roles in controlling ROS damage and regulating ROS signaling. J. Cell Biol. 217, 1915–1928 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
67.Colinet, D. et al. Extracellular superoxide dismutase in insects: characterization, function, and interspecific variation in parasitoid wasp venom. J. Biol. Chem. 286, 40110–40121 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
68.Colinet, D. et al. Extensive inter- and intraspecific venom variation in closely related parasites targeting the same host: the case of Leptopilina parasitoids of Drosophila. Insect Biochem. Mol. Biol. 43, 601–611 (2013).CAS 
PubMed 
Article 

Google Scholar 
69.Carton, Y., Frey, F. & Nappi, A. Genetic determinism of the cellular immune reaction in Drosophila melanogaster. Heredity 69, 393–399 (1992).PubMed 
Article 

Google Scholar 
70.Colinet, D., Schmitz, A., Depoix, D., Crochard, D. & Poirie, M. Convergent use of RhoGAP toxins by eukaryotic parasites and bacterial pathogens. PLoS Pathog. 3, 2029–2037 (2007).CAS 
Article 

Google Scholar 
71.Colinet, D. et al. The origin of intraspecific variation of virulence in an eukaryotic immune suppressive parasite. PLoS Pathog. 6, e1001206 (2010).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
72.Schlenke, T. A., Morales, J., Govind, S. & Clark, A. G. Contrasting infection strategies in generalist and specialist wasp parasitoids of Drosophila melanogaster. PLoS Pathog. 3, 1486–1501 (2007).PubMed 
Article 
CAS 

Google Scholar 
73.Anderl, I. et al. Transdifferentiation and proliferation in two distinct hemocyte lineages in Drosophila melanogaster larvae after wasp infection. PLoS Pathog. 12, e1005746 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
74.Forbes, A. A. et al. Revisiting the particular role of host shifts in initiating insect speciation. Evolution 71, 1126–1137 (2017).PubMed 
Article 

Google Scholar 
75.Allio, R. et al. Genome-wide macroevolutionary signatures of key innovations in butterflies colonizing new host plants. Nat. Commun. 12, 354 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
76.Araújo, M. S., Bolnick, D. I. & Layman, C. A. The ecological causes of individual specialisation. Ecol. Lett. 14, 948–958 (2011).PubMed 
Article 

Google Scholar 
77.Svanbäck, R. & Bolnick, D. I. Intraspecific competition drives increased resource use diversity within a natural population. P. Roy. Soc. B-Biol. Sci. 274, 839–844 (2007).
Google Scholar 
78.Laskowski, K. L. & Bell, A. M. Competition avoidance drives individual differences in response to a changing food resource in sticklebacks. Ecol. Lett. 16, 746–753 (2013).PubMed 
PubMed Central 
Article 

Google Scholar 
79.Huang, J., Reilein, A. & Kalderon, D. Yorkie and Hedgehog independently restrict BMP production in escort cells to permit germline differentiation in the Drosophila ovary. Development 144, 2584–2594 (2017).CAS 
PubMed 
PubMed Central 

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

Google Scholar 
81.Koren, S. et al. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res. 27, 722–736 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
82.Kolmogorov, M., Yuan, J., Lin, Y. & Pevzner, P. A. Assembly of long, error-prone reads using repeat graphs. Nat. Biotechnol. 37, 540–546 (2019).CAS 
PubMed 
Article 

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

Google Scholar 
84.Parra, G., Bradnam, K. & Korf, I. CEGMA: a pipeline to accurately annotate core genes in eukaryotic genomes. Bioinformatics 23, 1061–1067 (2007).CAS 
PubMed 
Article 

Google Scholar 
85.Waterhouse, R. M. et al. BUSCO applications from quality assessments to gene prediction and phylogenomics. Mol. Phylogenet. Evol. 35, 543–548 (2017).Article 
CAS 

Google Scholar 
86.Jurka, J. et al. Repbase Update, a database of eukaryotic repetitive elements. Cytogenet Genome Res. 110, 462–467 (2005).CAS 
PubMed 
Article 

Google Scholar 
87.Cantarel, B. L. et al. MAKER: an easy-to-use annotation pipeline designed for emerging model organism genomes. Genome Res. 18, 188–196 (2008).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
88.Li, H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34, 3094–3100 (2018).CAS 
PubMed 
PubMed Central 
Article 

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

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

Google Scholar 
91.Werren, J. H. et al. Functional and evolutionary insights from the genomes of three parasitoid Nasonia species. Science 327, 343–348 (2010).CAS 
PubMed 
Article 

Google Scholar 
92.Consortium, H. G. S. Insights into social insects from the genome of the honeybee Apis mellifera. Nature 443, 931–949 (2006).ADS 
Article 
CAS 

Google Scholar 
93.Geib, S. M., Liang, G. H., Murphy, T. D. & Sim, S. B. Whole genome sequencing of the braconid parasitoid wasp Fopius arisanus, an important biocontrol agent of pest tepritid fruit flies. G3-Genes Genom. Genet. 7, 2407–2411 (2017).CAS 

Google Scholar 
94.Standage, D. S. et al. Genome, transcriptome and methylome sequencing of a primitively eusocial wasp reveal a greatly reduced DNA methylation system in a social insect. Mol. Ecol. 25, 1769–1784 (2016).CAS 
PubMed 
Article 

Google Scholar 
95.Lindsey, A. R. et al. Comparative genomics of the miniature wasp and pest control agent Trichogramma pretiosum. BMC Biol. 16, 54 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
96.Stanke, M. et al. AUGUSTUS: ab initio prediction of alternative transcripts. Nucleic Acids Res. 34, W435–W439 (2006).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
97.Korf, I. Gene finding in novel genomes. BMC Biol. 5, 59 (2004).
Google Scholar 
98.Marchler-Bauer, A. et al. CDD: NCBI’s conserved domain database. Nucleic Acids Res. 43, D222–D226 (2015).CAS 
PubMed 
Article 

Google Scholar 
99.Hunter, S. et al. InterPro: the integrative protein signature database. Nucleic Acids Res. 37, D211–D215 (2008).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
100.Patro, R., Duggal, G., Love, M. I., Irizarry, R. A. & Kingsford, C. Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods 14, 417–419 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
101.Corpet, F. Multiple sequence alignment with hierarchical clustering. Nucleic Acids Res. 16, 10881–10890 (1988).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
102.Stamatakis, A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313 (2014).CAS 
PubMed 
PubMed Central 
Article 

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

Google Scholar 
104.Livak, K. J. & Schmittgen, T. D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C (T)) method. Methods 25, 402–408 (2001).CAS 
Article 

Google Scholar 
105.Zhang, X. S., Wang, T., Lin, X. W., Denlinger, D. L. & Xu, W. H. Reactive oxygen species extend insect life span using components of the insulin-signaling pathway. Proc. Natl Acad. Sci. USA 114, 7832–7840 (2017).Article 
CAS 

Google Scholar  More

1.Gaston, K. J. The Structure and Dynamics of Geographic Ranges (Oxford University Press, 2003).
Google Scholar 
2.Sexton, J. P., McIntyre, P., Angert, A. L. & Rice, K. J. Evolution and ecology of species range limits. Annu. Rev. Ecol. Evol. Syst. 40, 415–436 (2009).Article 

Google Scholar 
3.Holt, R. D., Keitt, T. H., Lewis, M. A., Maurer, B. A. & Taper, M. L. Theoretical models of species’ borders: single species approaches. Oikos 108, 18–27 (2005).Article 

Google Scholar 
4.Gaston, K. J. Geographic range limits: Achieving synthesis. Proc. R. Soc. Lond. B 276, 1395–1406 (2009).
Google Scholar 
5.Parmesan, C. et al. Empirical perspectives on species borders: From traditional biogeography to global change. Oikos 108, 58–75 (2005).Article 

Google Scholar 
6.Travis, J. M. J. & Dytham, C. In Dispersal Ecology and Evolution (eds Clobert, J. et al.) 337–348 (Oxford University Press, 2012).Chapter 

Google Scholar 
7.Kirkpatrick, M. & Barton, N. H. Evolution of a species’ range. Am. Nat. 150, 1–23 (1997).PubMed 
Article 
CAS 

Google Scholar 
8.Case, T. J. & Taper, M. L. Interspecific competition, environmental gradients, gene flow, and the coevolution of species’ borders. Am. Nat. 155, 583–605 (2000).PubMed 
Article 
CAS 

Google Scholar 
9.Rahbek, C. The role of spatial scale and the perception of large-scale species-richness patterns. Ecol. Lett. 8, 224–239 (2005).Article 

Google Scholar 
10.Hargreaves, A. L., Eckert, C. G. & Bailey, J. Evolution of dispersal and mating systems along geographic gradients. Implications for shifting ranges. Funct. Ecol. 28, 5–21 (2014).Article 

Google Scholar 
11.Hille, S. M. & Cooper, C. B. Elevational trends in life histories. Revising the pace-of-life framework. Biol. Rev. Camb. Philos. Soc. 90, 204–213 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
12.Boyle, W. A., Sandercock, B. K. & Martin, K. Patterns and drivers of intraspecific variation in avian life history along elevational gradients: A meta-analysis. Biol. Rev. 91, 469–482 (2016).Article 

Google Scholar 
13.Badyaev, A. V. & Ghalambor, C. K. Evolution of life histories along elevational gradients: Trade-off between parental care and fecundity. Ecology 82, 2948–2960 (2001).Article 

Google Scholar 
14.Bears, H., Martin, K. & White, G. C. Breeding in high-elevation habitat results in shift to slower life-history strategy within a single species. J. Anim. Ecol. 78, 365–375 (2009).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
15.Caro, L. M., Caycedo-Rosales, P. C., Bowie, R. C. K., Slabbekoorn, H. & Cadena, C. D. Ecological speciation along an elevational gradient in a tropical passerine bird?. J. Evol. Biol. 26, 357 (2013).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
16.Branch, C. L., Jahner, J. P., Kozlovsky, D. Y., Parchman, T. L. & Pravosudov, V. V. Absence of population structure across elevational gradients despite large phenotypic variation in mountain chickadees (Poecile gambeli). R. Soc. Open Sci. 4, 170057. https://doi.org/10.1098/rsos.170057 (2017).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
17.Chamberlain, D. E. et al. The altitudinal frontier in avian climate impact research. Ibis 154, 205–209 (2012).Article 

Google Scholar 
18.Hargreaves, A. L., Samis, K. E. & Eckert, C. G. Are species’ range limits simply niche limits writ large? A review of transplant experiments beyond the range. Am. Nat. 183, 157–173 (2014).PubMed 
Article 

Google Scholar 
19.Graham, C. H., Silva, N. & Velásquez-Tibatá, J. Evaluating the potential causes of range limits of birds of the Colombian Andes. J. Biogeogr. 37, 1863–1875 (2010).
Google Scholar 
20.Popy, S., Bordignon, L. & Prodon, R. A weak upward elevational shift in the distributions of breeding birds in the Italian Alps. J. Biogeogr. 37, 57–67 (2010).Article 

Google Scholar 
21.Chen, I.-C., Hill, J. K., Ohlemüller, R., Roy, D. B. & Thomas, C. D. Rapid range shifts of species associated with high levels of climate warming. Science 333, 1024–1026 (2011).ADS 
PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
22.Maggini, R. et al. Are Swiss birds tracking climate change?. Ecol. Model. 222, 21–32 (2011).Article 

Google Scholar 
23.Pearce-Higgins, J. W. & Green, R. E. Climate Change and Birds: Impacts and Conservation Responses (Cambridge University Press, 2014).Book 

Google Scholar 
24.Knaus, P. et al. Schweizer Brutvogelatlas 2013–2016. Verbreitung und Bestandsentwicklung der Vögel in der Schweiz und im Fürstentum Liechtenstein (Schweizerische Vogelwarte, 2018).
Google Scholar 
25.Chamberlain, D. E., Fuller, R. J., Bunce, R. G. H., Duckworth, J. C. & Shrubb, M. Changes in the abundance of farmland birds in relation to the timing of agricultural intensification in England and Wales. J. Appl. Ecol. 37, 771–788 (2000).Article 

Google Scholar 
26.Chamberlain, D. & Pearce-Higgins, J. Impacts of climate change on upland birds. Complex interactions, compensatory mechanisms and the need for long-term data. Ibis 155, 451–455 (2013).Article 

Google Scholar 
27.Sergio, F. & Newton, I. Occupancy as a measure of territory quality. J. Anim. Ecol. 72, 857–865 (2003).Article 

Google Scholar 
28.Fretwell, S. D. & Lucas, H. L. On territorial behavior and other factors influencing habitat distribution in birds. Acta Biotheor. 19, 16–36 (1969).Article 

Google Scholar 
29.Grüebler, M. U., Korner-Nievergelt, F. & von Hirschheydt, J. The reproductive benefits of livestock farming in barn swallows Hirundo rustica: Quality of nest site or foraging habitat?. J. Appl. Ecol. 47, 1340–1347 (2010).Article 

Google Scholar 
30.Schaub, M. & von Hirschheydt, J. Effects of current reproduction on apparent survival, breeding dispersal, and future reproduction in barn swallows assessed by multistate capture-recapture models. J. Anim. Ecol. 78, 625–635 (2009).PubMed 
Article 

Google Scholar 
31.Furrer, R. D. & Pasinelli, G. Empirical evidence for source-sink populations: A review on occurrence, assessments and implications. Biol. Rev. Camb. Philos. Soc. 91, 782–795 (2016).PubMed 
Article 

Google Scholar 
32.Plard, F., Turek, D., Grüebler, M. U. & Schaub, M. IPM2: Toward better understanding and forecasting of population dynamics. Ecol. Monogr. 8, e01364. https://doi.org/10.1002/ecm.1364 (2019).Article 

Google Scholar 
33.Grüebler, M. U. & Naef-Daenzer, B. Fitness consequences of pre- and post-fledging timing decisions in a double-brooded passerine. Ecology 89, 2736–2745 (2008).PubMed 
Article 

Google Scholar 
34.Grüebler, M. U., Morand, M. & Naef-Daenzer, B. A predictive model of the density of airborne insects in agricultural environments. Agric. Ecosyst. Environ. 123, 75–80 (2008).Article 

Google Scholar 
35.Jenni-Eiermann, S., Glaus, E., Grüebler, M. U., Schwabl, H. & Jenni, L. Glucocorticoid response to food availability in breeding barn swallows (Hirundo rustica). Gen. Comp. Endocrinol. 155, 558–565 (2008).PubMed 
Article 
CAS 

Google Scholar 
36.Schifferli, L., Grüebler, M. U., Meijer, H. A. J., Visser, G. H. & Naef-Daenzer, B. Barn Swallow Hirundo rustica parents work harder when foraging conditions are good. Ibis 156, 777–787 (2014).Article 

Google Scholar 
37.Shields, W. M. Factors Affecting nest and site fidelity in Adirondack barn swallows (Hirundo rustica). Auk 101, 780–789 (1984).Article 

Google Scholar 
38.Saino, N., Calza, S., Ninni, P. & Møller, A. P. Barn swallows trade survival against offspring condition and immunocompetence. J. Anim. Ecol. 68, 999–1009 (1999).Article 

Google Scholar 
39.Turner, A. The Barn Swallow (T & A D Poyser, 2006).
Google Scholar 
40.Newton, I. The Migration Ecology of Birds 1st edn. (Academic Press, 2007).
Google Scholar 
41.Ambrosini, R. & Saino, N. Environmental effects at two nested spatial scales on habitat choice and breeding performance of barn swallow. Evol. Ecol. 24, 491–508 (2010).Article 

Google Scholar 
42.Ambrosini, R. et al. The distribution and colony size of barn swallows in relation to agricultural land use. J. Appl. Ecol. 39, 524–534 (2002).Article 

Google Scholar 
43.Evans, K. L., Bradbury, R. B. & Wilson, J. D. Selection of hedgerows by Swallows Hirundo rustica foraging on farmland: the influence of local habitat and weather. Bird Study 50, 8–14 (2003).Article 

Google Scholar 
44.Newton, I. Population Limitation in Bird (Academic Press, 1998).
Google Scholar 
45.Paradis, E., Baillie, S. R., Sutherland, W. J. & Gregory, R. D. Patterns of natal and breeding dispersal in birds. J. Anim. Ecol. 67, 518–536 (1998).Article 

Google Scholar 
46.Scandolara, C. et al. Context-, phenotype-, and kin-dependent natal dispersal of barn swallows (Hirundo rustica). Behav. Ecol. 25, 180–190 (2014).Article 

Google Scholar 
47.Schaub, M., von Hirschheydt, J. & Grüebler, M. U. Differential contribution of demographic rate synchrony to population synchrony in barn swallows. J. Anim. Ecol. 84, 1530–1541 (2015).PubMed 
Article 

Google Scholar 
48.Camfield, A. F., Pearson, S. F. & Martin, K. Life history variation between high and low elevation subspecies of horned larks Eremophila spp. J. Avian Biol. 41, 273–281 (2010).Article 

Google Scholar 
49.Møller, A. P. Phenotype-dependent arrival time and its consequences in a migratory bird. Behav. Ecol. Sociobiol. 35, 115–122 (1994).Article 

Google Scholar 
50.Møller, A. P. Sexual Selection and the Barn Swallow (Oxford University Press, 1994).
Google Scholar 
51.Lerche-Jørgensen, M., Korner-Nievergelt, F., Tøttrup, A. P., Willemoes, M. & Thorup, K. Early returning long-distance migrant males do pay a survival cost. Ecol. Evol. 8, 11434–11449. https://doi.org/10.1002/ece3.4569 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
52.Pulliam, H. R. On the relationship between niche and distribution. Ecol. Lett. 3, 349–361 (2000).Article 

Google Scholar 
53.Pulliam, H. R. Sources, sinks, and population regulation. Am. Nat. 132, 652–661 (1988).Article 

Google Scholar 
54.Møller, A. P., de Lope, F. & Saino, N. Parasitism, immunity, and arrival date in a migratory bird, the barn swallow. Ecology 85, 206–219 (2004).Article 

Google Scholar 
55.Huntley, B., Green, R. E., Collingham, Y. C. & Willis, S. G. A Climatic Atlas of European Breeding Birds (Lynx Edicions, 2007).
Google Scholar 
56.Scridel, D. et al. A review and meta-analysis of the effects of climate change on Holarctic mountain and upland bird populations. Ibis 160, 489–515 (2018).Article 

Google Scholar 
57.Cormack, R. M. Estimates of survival from the sighting of marked animals. Biometrika 51, 429–438 (1964).MATH 
Article 

Google Scholar 
58.Jolly, G. Explicit estimates from capture-recapture data with both death and immigration-stochastic model. Biometrika 52, 225–247 (1965).MathSciNet 
PubMed 
MATH 
Article 
CAS 
PubMed Central 

Google Scholar 
59.Seber, G. A. F. A note on the multiple-recapture census. Biometrika 52, 249–259 (1965).MathSciNet 
PubMed 
MATH 
Article 
CAS 
PubMed Central 

Google Scholar 
60.Lebreton, J.-D., Burnham, K. P., Clobert, J. & Anderson, D. R. Modeling survival and testing biological hypotheses using marked animals: A unified approach with case studies. Ecol. Monogr. 62, 67–118 (1992).Article 

Google Scholar 
61.Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. https://doi.org/10.18637/jss.v067.i01 (2015).Article 

Google Scholar 
62.Saino, N., Martinelli, R. & Romano, M. Ecological and phenological covariates of offspring sex ratio in barn swallows. Evol. Ecol. 22, 659–674 (2008).Article 

Google Scholar 
63.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, 2017).64.Gelman, A. & Hill, J. Data Analysis Using Regression and Multilevel/Hierarchical Models (Cambridge University Press, 2007).
Google Scholar  More

1.Bowman, D. M. J. S. et al. Vegetation fires in the Anthropocene. Nat. Rev. Earth Environ. 1, 500–515 (2020).ADS 

Google Scholar 
2.Abatzoglou, J. T., Williams, A. P. & Barbero, R. Global emergence of anthropogenic climate change in fire weather indices. Geophys. Res. Lett. 46, 326–336 (2019).ADS 

Google Scholar 
3.Huang, Y., Wu, S. & Kaplan, J. O. Sensitivity of global wildfire occurrences to various factors in the context of global change. Atmos. Environ. 121, 86–92 (2015).ADS 
CAS 

Google Scholar 
4.van Oldenborgh, G. J. et al. Attribution of the Australian bushfire risk to anthropogenic climate change. Nat. Hazards Earth Syst. Sci. 21, 941–960 (2021).ADS 

Google Scholar 
5.Ward, M. et al. Impact of 2019–2020 mega-fires on Australian fauna habitat. Nat. Ecol. Evol. 4, 1321–1326 (2020).
Google Scholar 
6.Kablick III, G. P., Allen, D. R., Fromm, M. D. & Nedoluha, G. E. Australian PyroCb smoke generates synoptic-scale stratospheric anticyclones. Geophys. Res. Lett. 47, e2020GL088101 (2020).ADS 

Google Scholar 
7.Hirsch, E. & Koren, I. Record-breaking aerosol levels explained by smoke injection into the stratosphere. Science 371, 1269–1274 (2021).ADS 
CAS 

Google Scholar 
8.Schlosser, J. S. et al. Analysis of aerosol composition data for western United States wildfires between 2005 and 2015: Dust emissions, chloride depletion, and most enhanced aerosol constituents. J. Geophys. Res. Atmos. 122, 8951–8966 (2017).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
9.Barkley, A. E. et al. African biomass burning is a substantial source of phosphorus deposition to the Amazon, Tropical Atlantic Ocean, and Southern Ocean. Proc. Natl Acad. Sci. USA 116, 16216–16221 (2019).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
10.Guieu, C., Bonnet, S., Wagener, T. & Loÿe-Pilot, M.-D. Biomass burning as a source of dissolved iron to the open ocean? Geophys. Res. Lett. 32, L19608 (2005).ADS 

Google Scholar 
11.Ito, A. Mega fire emissions in Siberia: potential supply of bioavailable iron from forests to the ocean. Biogeosciences 8, 1679–1697 (2011).ADS 
CAS 

Google Scholar 
12.Abram, N. J., Gagan, M. K., McCulloch, M. T., Chappell, J. & Hantoro, W. S. Coral reef death during the 1997 Indian Ocean Dipole linked to Indonesian wildfires. Science 301, 952–955 (2003).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
13.Ito, A. et al. Pyrogenic iron: the missing link to high iron solubility in aerosols. Sci. Adv. 5, eaau7671 (2019).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
14.Jia, G. et al. in Climate Change and Land: an IPCC Special Report on Climate Change, Desertification, Land Degradation, Sustainable Land Management, Food Security, and Greenhouse Gas Fluxes in Terrestrial Ecosystems Ch. 2 (IPCC, in the press).15.Jiang, Y. et al. Impacts of wildfire aerosols on global energy budget and climate: the role of climate feedbacks. J. Clim. 33, 3351–3366 (2020).ADS 

Google Scholar 
16.Bowman, D. et al. Wildfires: Australia needs national monitoring agency. Nature 584, 188–191 (2020).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
17.New WWF report: 3 billion animals impacted by Australia’s bushfire crisis. WWF https://www.wwf.org.au/news/news/2020/3-billion-animals-impacted-by-australia-bushfire-crisis#gs.ebzve2 (2020).18.van der Velde, I. R. et al. Vast CO2 release from Australian fires in 2019–2020 constrained by satellite. Nature https://doi.org/10.1038/s41586-021-03712-y (2021).19.National Greenhouse Gas Inventory Report: 2018 (Australian Government, 2020); https://www.industry.gov.au/data-and-publications/national-greenhouse-gas-inventory-report-2018.20.Mahowald, N. M. et al. Aerosol impacts on climate and biogeochemistry. Annu. Rev. Environ. Res. 36, 45–74 (2011).
Google Scholar 
21.Boyd, P. W. et al. Mesoscale iron enrichment experiments 1993–2005: synthesis and future directions. Science 315, 612–617 (2007).ADS 
CAS 

Google Scholar 
22.Jickells, T. et al. Global iron connections between desert dust, ocean biogeochemistry, and climate. Science 308, 67–71 (2005).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
23.Martin, J. H. Glacial‐interglacial CO2 change: the iron hypothesis. Paleoceanography 5, 1–13 (1990).ADS 

Google Scholar 
24.Tagliabue, A. et al. Surface-water iron supplies in the Southern Ocean sustained by deep winter mixing. Nat. Geosci. 7, 314–320 (2014).ADS 
CAS 

Google Scholar 
25.Cassar, N. et al. The Southern Ocean biological response to aeolian iron deposition. Science 317, 1067–1070 (2007).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
26.Gabric, A. J., Cropp, R., Ayers, G. P., McTainsh, G. & Braddock, R. Coupling between cycles of phytoplankton biomass and aerosol optical depth as derived from SeaWiFS time series in the Subantarctic Southern Ocean. Geophys. Res. Lett. 29, 16-11–16-14 (2002).
Google Scholar 
27.Ardyna, M. et al. Hydrothermal vents trigger massive phytoplankton blooms in the Southern Ocean. Nat. Commun. 10, 2451 (2019).ADS 
PubMed 
PubMed Central 

Google Scholar 
28.Duprat, L. P. A. M., Bigg, G. R. & Wilton, D. J. Enhanced Southern Ocean marine productivity due to fertilization by giant icebergs. Nat. Geosci. 9, 219–221 (2016).ADS 
CAS 

Google Scholar 
29.Bixby, R. J. et al. Fire effects on aquatic ecosystems: an assessment of the current state of the science. Freshwater Sci. 34, 1340–1350 (2015).
Google Scholar 
30.Inness, A. et al. The CAMS reanalysis of atmospheric composition. Atmos. Chem. Phys. 19, 3515–3556 (2019).ADS 
CAS 

Google Scholar 
31.Shafeeque, M., Sathyendranath, S., George, G., Balchand, A. N. & Platt, T. Comparison of seasonal cycles of phytoplankton chlorophyll, aerosols, winds and sea-surface temperature off Somalia. Front. Marine Sci. 4, 384 (2017).
Google Scholar 
32.Cassar, N. et al. The influence of iron and light on net community production in the Subantarctic and Polar Frontal zones. Biogeosciences 8, 227–237 (2011).ADS 
CAS 

Google Scholar 
33.Mitchell, B. G. & Holm-Hansen, O. Observations of modeling of the Antartic phytoplankton crop in relation to mixing depth. Deep Sea Res. Part A 38, 981–1007 (1991).ADS 
CAS 

Google Scholar 
34.Longo, A. F. et al. Influence of atmospheric processes on the solubility and composition of iron in Saharan dust. Environ. Sci. Technol. 50, 6912–6920 (2016).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
35.Meskhidze, N., Nenes, A., Chameides, W. L., Luo, C. & Mahowald, N. Atlantic Southern Ocean productivity: fertilization from above or below? Global Biogeochem. Cycles 21, GB2006 (2007).ADS 

Google Scholar 
36.Sarmiento, J. L., Slater, R. D., Dunne, J., Gnanadesikan, A. & Hiscock, M. R. Efficiency of small scale carbon mitigation by patch iron fertilization. Biogeosciences 7, 3593–3624 (2010).ADS 
CAS 

Google Scholar 
37.Brzezinski, M. A., Jones, J. L. & Demarest, M. S. Control of silica production by iron and silicic acid during the Southern Ocean Iron Experiment (SOFeX). Limnol. Oceanogr. 50, 810–824 (2005).ADS 
CAS 

Google Scholar 
38.Lovenduski, N. S. & Gruber, N. Impact of the Southern Annular Mode on Southern Ocean circulation and biology. Geophys. Res. Lett. 32, L11603 (2005).ADS 

Google Scholar 
39.Cai, W., Cowan, T. & Raupach, M. Positive Indian Ocean Dipole events precondition southeast Australia bushfires. Geophys. Res. Lett. 36, L19710 (2009).ADS 

Google Scholar 
40.Chen, Y. et al. A pan-tropical cascade of fire driven by El Niño/Southern Oscillation. Nat. Climate Change 7, 906–911 (2017).ADS 
CAS 

Google Scholar 
41.Lim, E.-P. et al. Australian hot and dry extremes induced by weakenings of the stratospheric polar vortex. Nat. Geosci. 12, 896–901 (2019).ADS 
CAS 

Google Scholar 
42.Cai, W. et al. Increased frequency of extreme Indian Ocean Dipole events due to greenhouse warming. Nature 510, 254–258 (2014).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
43.Cropp, R. A. et al. The likelihood of observing dust-stimulated phytoplankton growth in waters proximal to the Australian continent. J. Mar. Syst. 117–118, 43–52 (2013).
Google Scholar 
44.Hamilton, D. S. et al. Impact of changes to the atmospheric soluble iron deposition flux on ocean biogeochemical cycles in the anthropocene. Global Biogeochem. Cycles 34, e2019GB006448 (2020).ADS 
CAS 

Google Scholar 
45.Duce, R. et al. Impacts of atmospheric anthropogenic nitrogen on the open ocean. Science 320, 893–897 (2008).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
46.Han, Y. et al. Asian inland wildfires driven by glacial-interglacial climate change. Proc. Natl Acad. Sci. USA 117, 5184–5189 (2020).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
47.van der Werf, G. R. et al. Global fire emissions estimates during 1997–2016. Earth Sys. Sci. Data 9, 697–720 (2017).ADS 

Google Scholar 
48.Orsi, A. H., Whitworth, T. & Nowlin, W. D. On the meridional extent and fronts of the Antarctic Circumpolar Current. Deep Sea Res. Part I 42, 641–673 (1995).
Google Scholar 
49.Sathyendranath, S. et al. An ocean-colour time series for use in climate studies: the experience of the Ocean-Colour Climate Change Initiative (OC-CCI). Sensors 19, 4285 (2019).ADS 
CAS 

Google Scholar 
50.Morcrette, J.-J. et al. Aerosol analysis and forecast in the European Centre for Medium-Range Weather Forecasts Integrated Forecast System: forward modeling. J. Geophys. Res. Atmospheres 114, D06206 (2009).ADS 

Google Scholar 
51.Levy, R. C. et al. Exploring systematic offsets between aerosol products from the two MODIS sensors. Atmos. Meas. Tech. 11, 4073–4092 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
52.Benedetti, A. et al. Aerosol analysis and forecast in the European Centre for Medium-Range Weather Forecasts Integrated Forecast System: 2. Data assimilation. J. Geophys. Res. 114, D13 (2009).
Google Scholar 
53.Kaiser, J. W. et al. Biomass burning emissions estimated with a global fire assimilation system based on observed fire radiative power. Biogeosciences 9, 527–554 (2012).ADS 
CAS 

Google Scholar 
54.Y. Bennouna et al. Validation Report of the CAMS Global Reanalysis of Aerosols and Reactive Gases, Years 2003–2019 (Copernicus Atmosphere Monitoring Service, 2020).55.Ito, A. et al. Evaluation of aerosol iron solubility over Australian coastal regions based on inverse modeling: implications of bushfires on bioaccessible iron concentrations in the Southern Hemisphere. Prog. Earth Planet. Sci. 7, 42 (2020).ADS 

Google Scholar 
56.Khaykin, S. et al. The 2019/20 Australian wildfires generated a persistent smoke-charged vortex rising up to 35 km altitude. Commun. Earth Environ. 1, 22 (2020).57.Haëntjens, N., Boss, E. & Talley, L. D. Revisiting Ocean Color algorithms for chlorophyll a and particulate organic carbon in the Southern Ocean using biogeochemical floats. J. Geophys. Res. Oceans 122, 6583–6593 (2017).ADS 

Google Scholar 
58.Boss, E. et al. The characteristics of particulate absorption, scattering and attenuation coefficients in the surface ocean; contribution of the Tara Oceans expedition. Methods Oceanogr. 7, 52–62 (2013).
Google Scholar 
59.de Boyer Montégut, C., Madec, G., Fischer, A. S., Lazar, A. & Iudicone, D. Mixed layer depth over the global ocean: an examination of profile data and a profile‐based climatology. J. Geophys. Res. Oceans 109, C12003 (2004).ADS 

Google Scholar 
60.Dong, S., Sprintall, J., Gille, S. T. & Talley, L. Southern Ocean mixed-layer depth from Argo float profiles. J. Geophys. Res. Oceans 113, C06013 (2008).ADS 

Google Scholar 
61.Cutter, G. A. et al. Sampling and Sample-handling Protocols for GEOTRACES Cruises, version 3.0 (2017).
Google Scholar 
62.Morton, P. L. et al. Methods for the sampling and analysis of marine aerosols: results from the 2008 GEOTRACES aerosol intercalibration experiment. Limnol. Oceanogr. Methods 11, 62–78 (2013).CAS 

Google Scholar 
63.Perron, M. M. G. et al. Assessment of leaching protocols to determine the solubility of trace metals in aerosols. Talanta 208, 120377 (2020).CAS 

Google Scholar 
64.Shelley, R. U., Landing, W. M., Ussher, S. J., Planquette, H. & Sarthou, G. Regional trends in the fractional solubility of Fe and other metals from North Atlantic aerosols (GEOTRACES cruises GA01 and GA03) following a two-stage leach. Biogeosciences 15, 2271–2288 (2018).ADS 
CAS 

Google Scholar 
65.Sanz Rodriguez, E. et al. Analysis of levoglucosan and its isomers in atmospheric samples by ion chromatography with electrospray lithium cationisation—triple quadrupole tandem mass spectrometry. J. Chromatogr. A 1610, 460557 (2020).CAS 

Google Scholar 
66.McLennan, S. M. Relationships between the trace element composition of sedimentary rocks and upper continental crust. Geochem. Geophys. Geosyst. 2, 1201 (2001).
Google Scholar 
67.Shelley, R. U. et al. Quantification of trace element atmospheric deposition fluxes to the Atlantic Ocean ( >40°N; GEOVIDE, GEOTRACES GA01) during spring 2014. Deep Sea Res. Part I 119, 34–49 (2017).CAS 

Google Scholar 
68.Sholkovitz, E. R., Sedwick, P. N., Church, T. M., Baker, A. R. & Powell, C. F. Fractional solubility of aerosol iron: synthesis of a global-scale data set. Geochim. Cosmochim. Acta 89, 173–189 (2012).ADS 
CAS 

Google Scholar 
69.Stein, A. F. et al. NOAA’s HYSPLIT atmospheric transport and dispersion modeling system. Bull. Am. Meteorol. Soc. 96, 2059–2077 (2016).ADS 

Google Scholar 
70.Kalnay, E. et al. The NCEP/NCAR 40-year reanalysis project. Bull. Am. Meteorol. Soc. 77, 437–471 (1996).ADS 

Google Scholar 
71.Tatlhego, M., Bhattachan, A., Okin, G. S. & D’Odorico, P. Mapping areas of the Southern Ocean where productivity likely depends on dust‐delivered Iron. J. Geophys. Res. Atmospheres 125, e2019JD030926 (2020).ADS 
CAS 

Google Scholar 
72.Stein, A. F., Rolph, G. D., Draxler, R. R., Stunder, B. & Ruminski, M. Verification of the NOAA smoke forecasting system: model sensitivity to the injection height. Weather Forecast. 24, 379–394 (2009).ADS 

Google Scholar 
73.Behrenfeld, M. J. & Falkowski, P. G. Photosynthetic rates derived from satellite‐based chlorophyll concentration. Limnol. Oceanogr. 42, 1–20 (1997).ADS 
CAS 

Google Scholar 
74.Behrenfeld, M. J., Boss, E., Siegel, D. A. & Shea, D. M. Carbon-based ocean productivity and phytoplankton physiology from space. Global Biogeochem. Cycles 19, GB1006 (2005).ADS 

Google Scholar 
75.Westberry, T., Behrenfeld, M. J., Siegel, D. A. & Boss, E. Carbon-based primary productivity modeling with vertically resolved photoacclimation. Global Biogeochem. Cycles 22, GB2024 (2008).ADS 

Google Scholar 
76.Silsbe, G. M., Behrenfeld, M. J., Halsey, K. H., Milligan, A. J. & Westberry, T. K. The CAFE model: a net production model for global ocean phytoplankton. Global Biogeochem. Cycles 30, 1756–1777 (2016).ADS 
CAS 

Google Scholar 
77.Laws, E. A., D’Sa, E. & Naik, P. Simple equations to estimate ratios of new or export production to total production from satellite‐derived estimates of sea surface temperature and primary production. Limnol. Oceanogr. Methods 9, 593–601 (2011).
Google Scholar 
78.Dunne, J. P., Armstrong, R. A., Gnanadesikan, A. & Sarmiento, J. L. Empirical and mechanistic models for the particle export ratio. Global Biogeochem. Cycles 19, GB4026 (2005).ADS 

Google Scholar 
79.Li, Z. & Cassar, N. Satellite estimates of net community production based on O2/Ar observations and comparison to other estimates. Global Biogeochem. Cycles 30, 735–752 (2016).ADS 
CAS 

Google Scholar 
80.Siegel, D. A. et al. Global assessment of ocean carbon export by combining satellite observations and food‐web models. Global Biogeochem. Cycles 28, 181–196 (2014).ADS 
CAS 

Google Scholar 
81.Marshall, G. J. Trends in the Southern Annular Mode from observations and reanalyses. J. Climate 16, 4134–4143 (2003).ADS 

Google Scholar 
82.Saji, N. H. & Yamagata, T. Possible impacts of Indian Ocean Dipole mode events on global climate. Climate Res. 25, 151–169 (2003).ADS 

Google Scholar  More

The Thau lagoon dataThe measurement campaign concerned four wastewater treatment plants (WWTP) in the Thau lagoon area in France, serving the cities of Sète, Pradel-Marseillan, Frontignan and Mèze. The measurements were obtained by using digital PCR20 (dPCR) to estimate the concentration of SARS-CoV-2 virus in samples taken weekly from 2020-05-12 to 2021-01-12. We provide further details about the measurement method in the “Methods” section.Figure 1Concentrations of SARS-CoV-2 (genome units per litre in logarithmic scale) from four WWTPs in Thau lagoon, measured weekly with dPCR technology from May 12th 2020 to January 12th, 2021. Note that there are some missing points.Full size imageFigure 1 shows the outcomes in a logarithmic scale over a nine months period. We summarise now their main features.

1.

An exponential phase starts simultaneously in Mèze and Frontignan WWTPs in early June.

2.

The genome units concentration curves in these two places reach, again simultaneously, a plateau. It has stayed essentially stable or slightly decreasing since then.

3.

The evolution at Sète and Pradel-Marseillan remarkably followed the previous two zones in a parallel way, with a two weeks lag. The measurements at Sète and Pradel-Marseillan continued to grow linearly (recall that this is in log scale, thus exponentially in linear scale), while the Mèze and Frontignan figures have stabilised ; then, after two weeks, they too stabilised at a plateau with roughly the same value as for the other two towns.

4.

The measurements seem to show a tendency to increase over the very last period.

The epidemiological model with heterogeneity and natural variability of population behaviourThe appearance of such plateaus and shoulders need not be explained either by psychological reactions or by public health policy effects. Indeed, the regulations were roughly constant during the measurement campaign and awareness or fatigue effects do not seem to have altered the dynamics over this long period of time. Witness to this is the fact that two groups of towns saw the same evolution, but two weeks apart one from the other. To understand this phenomena we propose a new model.Given the complexity and multiplicity of behavioural factors favouring the spread of the epidemic, we assume that the transmission rate involves a normalised variable (a in (0,1)) that defines an aggregated indicator of risky behaviour within the susceptible population. Thus, we represent the heterogeneity of individual behaviours with this variable. We take a as an implicit parameter that we do not seek to calculate. The classical SIR model is macroscopic and the type of model we discuss here can be viewed as intermediate between macroscopic and microscopic.The initial distribution of susceptible individuals (S_0(a)) in the framework of a SIR-type compartmental description of the epidemic can be reasonably taken as a decreasing function of a. We take the infection transmission rate (a mapsto beta (a)) to be an increasing function of a. In the Supplementary Information (SI) Appendix, the reader will find a more general version of this model involving a probability kernel of transition from one state to another. The model here can be derived as a limiting case of that more general version.Likewise, the behaviour of individuals usually changes from one day to another21. Many factors are at work in this variability: social imitation, public health campaigns, opportunities, outings, the normal variations of activity (e.g. work from home certain days and use of public transportation and work in office on others) etc. Therefore, the second key feature of our model is variability of such behaviours: variations of the population density for a given a do not only come from individuals becoming infected and leaving that compartment but also results from individuals moving from one state a to another21. In the simplest version of the model, variability is introduced as a diffusion term in the dynamics of susceptible individuals.The modelWe denote by S(t, a) the density of individuals at time t associated with risk parameter a, by I(t) the total number of infected, and by R(t) the number of removed individuals. We are then led to the following system:$$begin{aligned} frac{{partial S(t,a)}}{{partial t}} & = d{mkern 1mu} frac{{partial ^{2} S(t,a)}}{{partial a^{2} }} – beta (a)S(t,a)frac{{I(t)}}{N} \ frac{{{text{d}}I(t)}}{{{text{d}}t}} & = frac{{I(t)}}{N}{mkern 1mu} intlimits_{0}^{1} beta (a)S(t,a);da – gamma I(t), \ frac{{{text{d}}R(t)}}{{{text{d}}t}} = & gamma I(t), \ end{aligned}$$
(1)
where (gamma) denotes the inverse of typical duration (in days) of the disease and d a positive diffusion coefficient. System (1) is supplemented with initial conditions$$begin{aligned} S(0,a) = S_0(a), quad I(0) = I_0, quad hbox {and} quad R(0) = 0, end{aligned}$$
(2)
and with zero flux condition in a at (a=0, 1). In the “Methods” section below, we discuss the relation of this model with other current works.A more general modelIn a more general version of our model, we can consider the population of infected as also structured by the parameter a. The equations are as before but now we keep track of the class a in the infected population. The mechanism here is that a susceptible individual from class a can be infected by infectious from any class I(t, b) but then gives rise to an individual I(t, a) of the same parent class. We also assume that there is a diffusion of the infected behaviours. We denote by ({mathfrak {B}}(a,b)) the transmission rate of S(t, a) by I(t, b). For simplicity and because it is natural, we will assume that it is of the form$$begin{aligned} {mathfrak {B}}(a,b)= beta (a) beta (b) end{aligned}$$where (beta) is as before. For full generality, we can also envision multi-dimensional parameters (ain {mathbb {R}}^d), with (a_iin (0,1)). We are then led to the system:$$begin{aligned} frac{{partial S(t,a)}}{{partial t}} & = d;Delta _{a} S(t,a) – S(t,a)frac{{beta (a)}}{N}intlimits_{0}^{1} beta (b)I(t,b);db \ frac{{partial I(t,a)}}{{partial t}} & = d^{prime}Delta _{a} I(t,a) + S(t,a)frac{{beta (a)}}{N}intlimits_{0}^{1} beta (b)I(t,b)db – gamma I(t,a), \ frac{{{text{d}}R(t)}}{{{text{d}}t}} & = gamma intlimits_{0}^{1} I (t,b){mkern 1mu} db, \ end{aligned}$$
(3)
In the SI we write further, more general, forms of this model, with ({mathfrak {B}}(a,b)) and more general diffusion of behaviours, that can include jumps or non-local variations. The type of models we discuss here may also shed light on the initial phase of the epidemic. We plan to investigate these questions in future work.Patterns generated by the modelIn the next section, we will discuss how the model fits the data observed in the Thau lagoon measurements. But before that, we start by showing that the above model (1) can generate the different patterns we mentioned. For this we rely on numerical simulations without fitting real data. And indeed we obtain plateaus, shoulders, and oscillations. The latter can be interpreted as epidemic rebounds.The key parameter here is the diffusion coefficient d, which controls the amplitude of behavioural variability (see Fig. 2). Large values of d rapidly yield homogenised behaviours, leading to classical SIR-like dynamics of infectious individuals. For very small values of d, the system also has a simple dynamics, in the sense that I(t) has a unique maximum, and therefore has no rebounds. We derive this in the limit (d=0) for which we show in the SI that there are neither plateaus nor rebounds.For some intermediate range of the parameter d, plateaus may appear after an exponential growth, like in the initial phase of the SIR model. A small amplitude oscillation, called “shoulder”, precedes a temporary stabilisation on a plateau, followed by a large time convergence to zero of infectious population. We also show that for small enough d, time oscillations of the infectious population curve, i.e. epidemic rebounds, may be generated by Model (1). Such oscillations also appear after a plateau, in a similar way to what one can see in observations.Simulations in Fig. 2 illustrate the various patterns obtained on the dynamics of infected population as a function of the diffusion parameter. For small enough d, in the top left graph of Fig. 2, one can see oscillations of the fraction of infectious individuals. These oscillations cannot be achieved in the classical SIR model. In fact, the two lower graphs of that figure, for somewhat larger values of d, exhibit the SIR model outcomes. Indeed, for sufficiently large d, the system becomes rapidly homogeneous (i.e. constant with respect to a). Yet, such oscillations are standard in the dynamics of actual epidemics, like the current Covid-19 pandemic. The intermediate value of d, represented in the upper right corner of Fig. 2 shows the typical onset of a plateau at a rather high value of I. Note that this plateau is preceded by a first small dip and then a characteristic “shoulder-like” oscillation.Figure 2Model behaviour depending on diffusion parameter values: infected rate dynamics in logarithmic scale. From left to right and then top to bottom, graphs are associated with (d=10^{-5}), (d=5times 10^{-5}), (d=10^{-3}) and (d=5times 10^{-3}) (in (day^{-1}) unit).Full size imageSecondary epidemic peaks are of lower amplitude than the first one, as shown in the top graphs of Fig. 2. This empirical observation leads us to conjecture that, at least in many cases, it is a general property of this model (with (beta) independent of time). This property would then reflect a kind of dissipative nature of Model (1). It is natural to surmise that a change of behaviours in time may generate oscillations with higher secondary peaks. Such changes result for instance from lifting social distancing measures or from fatigue effects in the population.We illustrate this with numerical simulations in Fig. 3. We assume a collective time modulation of the (beta (a)) transmission profile. That is, we replace (beta (a)) by (beta (a)varphi (t)) for some time dependent function (varphi), the other parameters are the same as in the simulations shown in Fig. 2. We look at the effect of a “lockdown exit” type effect. Then, (varphi (t)) takes two constant values, 1 from (t=0) to (t={1000}) and 1.2 after (t={1100}). In between, that is, for (tin ({1000}, {1100})), (varphi (t)) changes linearly from the value 1 to 1.2.Figure 3Multiple epidemic rebounds: susceptible individuals are divided into 50 discrete groups in the case where relaxation of social distancing measures starts on Day (t=1000) and ends up on Day (t=1100). The fraction of infected individuals in the population is represented in the left graph in logarithmic scale and in linear scale in the right graph.Full size imageOne can clearly see a secondary peak with higher amplitude than the first one. Note that after this peak, a third one occurs, with a lower amplitude than the second one. This third peak happens in the regime when (beta) is again constant in time.The effect of variantsAnother important factor that yields secondary peaks with higher amplitudes is the appearance of variants. Consider the situation with two variants. We denote by (I_1(t)) and (I_2(t)) the corresponding infected individuals. The first variant, which we call the historical strain, is associated with (beta _1) and (I_1(0)) and starts at (t=0). The variant strain corresponds to (beta _2) and (I_2) and starts at Day (t=1000). In this situation, the system (1) is extended by the following system:$$begin{aligned} frac{{partial S(t,a)}}{{partial t}} & = d{mkern 1mu} frac{{partial ^{2} S(t,a)}}{{partial a^{2} }} – left( {beta _{1} (a)I_{1} (t) + beta _{2} (a)I_{2} (t)} right)frac{{S(t,a)}}{N}, \ frac{{{text{d}}I_{2} (t)}}{{{text{d}}t}} & = frac{{I_{2} (t)}}{N}{mkern 1mu} intlimits_{0}^{1} {beta _{2} } (a)S(t,a){mkern 1mu} da – gamma _{2} I_{2} (t), \ frac{{{text{d}}I_{1} (t)}}{{{text{d}}t}} & = frac{{I_{1} (t)}}{N}{mkern 1mu} intlimits_{0}^{1} {beta _{1} } (a)S(t,a){mkern 1mu} da – gamma _{1} I_{1} (t) \ frac{{{text{d}}R(t)}}{{{text{d}}t}} & = gamma _{1} I_{2} (t) + gamma _{1} I_{2} (t), \ end{aligned}$$
(4)
The total infected population is (I(t)=I_1(t)+I_2(t)). Figure 4 shows a simulation of this system. Before the onset of the second variant, i.e. for (t< 1000), we observe a peak, followed by a small shoulder and a downward tilted plateau. The second variant corresponds to a higher transmission coefficient: namely, we take here (beta _2(a) equiv frac{3}{2} beta _1(a)). When it appears at time (t=1000), initially there is no effect, because the initial number of infectious with variant 2 is very small. Then, there is an exponential growth caused by this second variant gaining strength. The secondary peak is then higher than the first one. A very small shoulder precedes another stabilisation on a downward plateau.Figure 4 also shows the dynamics of fractions of infected with each one of the variants. Note that the infectious with variant 1 very rapidly all but disappear at the onset of the second exponential growth phase. One might have expected that the historical strain would be gradually replaced by the new strain, merely tilting further downward the plateau. But that does not happen. Thus, it is remarkable that the historical strain gets nearly wiped out at the very beginning of the second exponential growth.Figure 4Multiple epidemic rebounds due to a variant virus: susceptible individuals are divided into 50 discrete groups in the case where a new variant appears at Day (t=1000). The transmission rate (beta _2) is taken as (beta _2(a) = 1.5 , beta _1(a)), (d=0.0002), (gamma _1=0.1) and (gamma _2= 0.05). The fraction of infected individuals in the population is represented in the left graph in logarithmic scale. The total infected population is represented in linear scale in the right graph (black curve), variant 1 in red and variant 2 in green.Full size imageApplication to the Thau lagoon measurementsModel (1) describes the dynamics of the fraction of infectious in the population, that is (t mapsto I(t)/N). Therefore, we need to derive this fraction from the wastewater measurements. To this end, we use an “effective proportionality coefficient” between the two quantities. This coefficient itself is derived from the measurements (compare Section “SARS-CoV-2 concentration measurement from wastewater with digital PCR” in the “Methods” part below). Calibration of model (1) also requires fitting the values of (gamma), the profiles (a mapsto beta (a)) and the initial distribution of susceptible individuals in terms of a.We carried this procedure and the resulting fitted curve is displayed in Fig. 5. Note that the outcome correctly captures the shoulder and plateau patterns.Figure 5Calibrated model on Sète area: blue dots are measures of SARS-CoV-2 genome units and black curve represents the total infected individuals as an output of the model discretized into (n_g=20) groups in a. Initial distribution of susceptible individuals and (beta) function are taken as described in supplementary information. Parameters d and (gamma) are taken as follows: (d=2.5 times 10^{-4}) (day^{-1}), and (gamma =0.1) (day^{-1}).Full size imageThe underlying dynamics of the rate of susceptible individuals is given in Fig. 6 below for (n_g=20) groups. The lower curve illustrates the competition phenomenon between diffusion and sink term due to new infections, depending on the level of risk a of each state.Figure 6Calibrated model on Sète WWTP: density of susceptible individuals of each group for (n_g=20). The densities of susceptible of each group is represented in colour curves as functions of time. The curves are ordered from top to bottom according to increasing a. The resulting average total susceptible population is represented in black. Susceptible individuals of highest a trait, which are represented in the bottom light blue curve exhibit a non monotonic behaviour.Full size image More