Ecology

1.
Hiernaux, J. Bantu expansion: the evidence from physical anthropology confronted with linguistic and archaeological evidence. J. Afr. Hist. 9, 505–515 (1968).
Article  Google Scholar 
2.
Vansina, J. New linguistic evidence and ‘the Bantu expansion’. J. Afr. Hist. 36, 173–195 (1995).
Article  Google Scholar 

3.
de Filippo, C., Bostoen, K., Stoneking, M. & Pakendorf, B. Bringing together linguistic and genetic evidence to test the Bantu expansion. Proc. Biol. Sci. 279, 3256–3263 (2012).
PubMed  PubMed Central  Google Scholar 

4.
Li, S., Schlebusch, C. & Jakobsson, M. Genetic variation reveals large-scale population expansion and migration during the expansion of Bantu-speaking peoples. Proc. R. Soc. B Biol. Sci. 281, 20141448 (2014).
Article  Google Scholar 

5.
Phillipson, D. W. The Later Prehistory of Eastern and Southern Africa (Heinemann, London, 1977).
Google Scholar 

6.
Holden, C. J. Bantu language trees reflect the spread of farming across sub-Saharan Africa: a maximum-parsimony analysis. Proc. Biol. Sci. 269, 793–799 (2002).
PubMed  PubMed Central  Article  Google Scholar 

7.
Bostoen, K. Pearl millet in early Bantu speech communities in Central Africa: a reconsideration of the lexical evidence. Afrika und Übersee 89, 183–213 (2007).

8.
Russell, T., Silva, F. & Steele, J. Modelling the spread of farming in the Bantu-speaking regions of Africa: an archaeology-based phylogeography. PLoS ONE 9, e87854 (2014).
PubMed  PubMed Central  Article  CAS  Google Scholar 

9.
Kahlheber, S., Bostoen, K. & Neumann, K. Early plant cultivation in the Central African rain forest: first millennium BC pearl millet from South Cameroon. J. Afr. Archaeol. 7, 253–272 (2009).
Article  Google Scholar 

10.
Grollemund, R. et al. Bantu expansion shows that habitat alters the route and pace of human dispersals. Proc. Natl Acad. Sci. USA 112, 13296–13301 (2015).
CAS  PubMed  Article  PubMed Central  Google Scholar 

11.
Bayon, G. et al. Intensifying weathering and land use in Iron Age Central Africa. Science 335, 1219–1222 (2012).
CAS  PubMed  Article  PubMed Central  Google Scholar 

12.
Bayon, G. et al. The roles of climate and human land-use in the late Holocene rainforest crisis of Central Africa. Earth Planet. Sci. Lett. 505, 30–41 (2019).
CAS  Article  Google Scholar 

13.
Garcin, Y. et al. Early anthropogenic impact on Western Central African rainforests 2,600 y ago. Proc. Natl Acad. Sci. USA 115, 3261–3266 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

14.
Ambrose, S. H. & DeNiro, M. J. The isotopic ecology of East African mammals. Oecologia 69, 395–406 (1986).
PubMed  Article  PubMed Central  Google Scholar 

15.
Ambrose, S. H. Stable carbon and nitrogen isotope analysis of human and animal diet in Africa. J. Hum. Evol. 15, 707–731 (1987).
Article  Google Scholar 

16.
Lee-Thorp, J. A., Sealy, J. C. & Morris, A. G. Iin Prehistoric Human Bone: Archaeology at the Molecular Level  (eds Lambert, J. & Grupe, G.) 99–120 (Springer Verlag, Berlin, 1993).

17.
Dlamini, N., Morris, A. G. & Sealy, J. Carbon isotopes and dental caries as evidence for regional variation in the diets of early farming communities from Katanga, Democratic Republic of the Congo. J. Afr. Archaeol. 14, 135–153 (2016).
Article  Google Scholar 

18.
Eggert, M. K. H. et al. Pits, graves and grains: archaeological and archaeobotanical research in southern Cameroon. J. Afr. Archaeol. 4, 273–298 (2006).
Article  Google Scholar 

19.
Eggert, M. K. H. In Europa im Geflecht der Welt. Mittelalterliche Migrationen in Globalen Bezügen (eds Borgolte, M., Dücker, J., Müllerburg, M., Predatsch, P. & Schneidmüller, B.) 193–216 (Akademie Verlag, Berlin, 2012).

20.
Eggert, M. K. H. In The Cambridge World Prehistory, Volume 1: Africa, South and Souheast Asia and the Pacific (eds Renfrew, C. & Bahn, P.) 183–203 (Cambridge University Press, Cambridge, 2014).

21.
Neumann, K., Bostoen, K., Höhn, A., Kahlheber, S., Ngomanda, A. & Tchiengue, B. First farmers in the Central African rainforest: a view from southern Cameroon. Quat. Int. 249, 53–62 (2012).

22.
Oslisly, R. et al. Climatic and cultural changes in the west Congo Basin forests over the past 5000 years. Philos. Trans. R. Soc. Lond. B Biol. Sci. 368, 20120304 (2013).
PubMed  PubMed Central  Article  Google Scholar 

23.
Oslisly, R. et al. West Central African peoples: survey of radiocarbon dates over the past 5000 years. Radiocarbon 55, 1377–1382 (2013).
CAS  Article  Google Scholar 

24.
Oslisly, R., White, L. & de Saulieu, G. In Pour une écologie historique en Afrique Centrale (eds de Saulieu, G., Elouga, M. & Sonké, B.) 121–139 (IRD/AUF, Yaoundé, 2016).

25.
Bostoen, K. et al. Middle to Late Holocene Paleoclimatic Change and the Early Bantu Expansion in the rain forests of Western Central Africa. Curr. Anthropol. 56, 354–384 (2015).
Article  Google Scholar 

26.
Wotzka, H.-P. In Trees, Grasses and Crops. People and Plants in Sub-Saharan Africa and Beyond (eds Eichhorn, B. & Höhn, A.) 407–429 (Dr. Rudolf Habelt GmbH, 2019).

27.
Wotzka, H.-P. In Auf dem Holzweg … Eine Würdigung für Ursula Tegtmeier. (eds Meurers-Balke, J., Zerl, T. & Gerlach, R.) 269–284 (Propylaeum, 2019).

28.
Kahlheber, S., Eggert, M. K. H., Seidensticker, D. & Wotzka, H.-P. Pearl millet and other plant remains from the Early Iron Age Site of Boso-Njafo (Inner Congo Basin, Democratic Republic of the Congo). Afr. Archaeol. Rev. 31, 479–512 (2014).
Article  Google Scholar 

29.
Schwartz, D. Assèchement climatique vers 3000 B.P. et expansion Bantu en Afrique centrale atlantique: quelques réflexions. Bull. Soc. Géol. Fr. 163, 353–361 (1992).
Google Scholar 

30.
Ngomanda, A., Neumann, K., Schweizer, A. & Maley, J. Seasonality change and the third millennium BP rainforest crisis in southern Cameroon (Central Africa). Quat. Res. 71, 307–318 (2009).
Article  Google Scholar 

31.
Mercader, J., Garralda, M. D., Pearson, O. M. & Bailey, R. C. Eight hundred-year-old human remains from the Ituri tropical forest, Democratic Republic of Congo: the rock shelter site of Matangai Turu Northwest. Am. J. Phys. Anthropol. 115, 24–37 (2001).
CAS  PubMed  Article  Google Scholar 

32.
Lee-Thorp, J. A., Sealy, J. C. & van der Merwe, N. J. Stable carbon isotope ratio differences between bone collagen and bone apatite, and their relationship to diet. J. Archaeol. Sci. 16, 585–599 (1989).
Article  Google Scholar 

33.
Kellner, C. M. & Schoeninger, M. J. A simple carbon isotope model for reconstructing prehistoric human diet. Am. J. Phys. Anthropol. 133, 1112–1127 (2007).
PubMed  Article  Google Scholar 

34.
Roberts, P. et al. Direct evidence for human reliance on rainforest resources in late Pleistocene Sri Lanka. Science 347, 1246–1249 (2015).
CAS  PubMed  Article  Google Scholar 

35.
Cerling, T. E., Hart, J. A. & Hart, T. B. Stable isotope ecology in the Ituri Forest. Oecologia 138, 5–12 (2004).
PubMed  Article  Google Scholar 

36.
van der Merwe, N. J. & Medina, E. The canopy effect, carbon isotope ratios and foodwebs in amazonia. J. Archaeol. Sci. 18, 249–259 (1991).
Article  Google Scholar 

37.
Verhegghen, A., Mayaux, P., Wasseige, Cde & Defourny, P. Mapping Congo Basin vegetation types from 300 m and 1 km multi-sensor time series for carbon stocks and forest areas estimation. Biogeosciences 9, 5061–5079 (2012).
CAS  Article  Google Scholar 

38.
Hart, T. B., Hart, J. A., Dechamps, R., Fournier, M. & Ataholo, M. C. Changes in forest composition over the last 4000 years in the Ituri basin, Zaire. In The Biodiversity of African Plants: Proceedings XIVth AETFAT Congress 22–27 August 1994, Wageningen, The Netherlands (eds van der Maesen, L. J. G., van der Burgt, X. M. & van Medenbach de Rooy, J. M.) 545–560 (Springer Netherlands, 1996).

39.
Mercader, J. et al. Phytoliths from archaeological sites in the tropical forest of Ituri, Democratic Republic of Congo. Quat. Res. 54, 102–112 (2000).
Article  Google Scholar 

40.
Schoeninger, M. J. & DeNiro, M. J. Nitrogen and carbon isotopic composition of bone collagen from marine and terrestrial animals. Geochim. Cosmochim. Acta 48, 625–639 (1984).
CAS  Article  Google Scholar 

41.
O’Connell, T. C., Kneale, C. J., Tasevska, N. & Kuhnle, G. G. C. The diet-body offset in human nitrogen isotopic values: a controlled dietary study. Am. J. Phys. Anthropol. 149, 426–434 (2012).
PubMed  PubMed Central  Article  Google Scholar 

42.
Bocherens, H., Koch, P. L., Mariotti, A., Geraads, D. & Jaeger, J.-J. Isotopic biogeochemistry (13 C, 18 O) of mammalian enamel from African Pleistocene hominid sites. Palaios 11, 306–318 (1996).

43.
Sponheimer, M. & Lee-Thorp, J. A. The oxygen isotope composition of mammalian enamel carbonate from Morea Estate, South Africa. Oecologia 126, 153–157 (2001).
PubMed  Article  PubMed Central  Google Scholar 

44.
DeNiro, M. J. Postmortem preservation and alteration of in vivo bone collagen isotope ratios in relation to palaeodietary reconstruction. Nature 317, 806 (1985).
CAS  Article  Google Scholar 

45.
Ambrose, S. H. Preparation and characterization of bone and tooth collagen for isotopic analysis. J. Archaeol. Sci. 17, 431–451 (1990).
Article  Google Scholar 

46.
van Klinken, G. J. Bone collagen quality indicators for palaeodietary and radiocarbon measurements. J. Archaeol. Sci. 26, 687–695 (1999).
Article  Google Scholar 

47.
Sealy, J., Johnson, M., Richards, M. & Nehlich, O. Comparison of two methods of extracting bone collagen for stable carbon and nitrogen isotope analysis: comparing whole bone demineralization with gelatinization and ultrafiltration. J. Archaeol. Sci. 47, 64–69 (2014).
CAS  Article  Google Scholar 

48.
Harrington, N. R. The life habits of polypterus. Am. Nat. 33, 721–728 (1899).
Article  Google Scholar 

49.
Bocherens, H. & Drucker, D. Trophic level isotopic enrichment of carbon and nitrogen in bone collagen: case studies from recent and ancient terrestrial ecosystems. Int. J. Osteoarchaeol. 13, 46–53 (2003).
Article  Google Scholar 

50.
Cerling, T. E., Harris, J. M., Ambrose, S. H., Leakey, M. G. & Solounias, N. Dietary and environmental reconstruction with stable isotope analyses of herbivore tooth enamel from the Miocene locality of Fort Ternan, Kenya. J. Hum. Evol. 33, 635–650 (1997).
CAS  PubMed  Article  PubMed Central  Google Scholar 

51.
Roberts, P. et al. Fruits of the forest: Human stable isotope ecology and rainforest adaptations in Late Pleistocene and Holocene (~36 to 3 ka) Sri Lanka. J. Hum. Evol. 106, 102–118 (2017).
PubMed  Article  PubMed Central  Google Scholar 

52.
D’Andrea, A. C., Logan, A. L. & Watson, D. J. Oil palm and prehistoric subsistence in tropical West Africa. J. Afr. Archaeol. 4, 195–222 (2006).
Article  Google Scholar 

53.
Mbida, C. M., Van Neer, W., Doutrelepont, H. & Vrydaghs, L. Evidence for banana cultivation and animal husbandry during the first millennium bc in the forest of Southern Cameroon. J. Archaeol. Sci. 27, 151–162 (2000).
Article  Google Scholar 

54.
Soto, M. et al. Structural characterization and decontamination of dental calculus for ancient starch research. Archaeol. Anthropol. Sci. 11, 4847–4872 (2019).
Article  Google Scholar 

55.
Mercader, J. et al. Morphometrics of starch granules from sub-saharan plants and the taxonomic identification of ancient starch. Front. Earth Sci. 6, 146 https://doi.org/10.3389/feart.2018.00146 (2018).

56.
Höhn, A., Kahlheber, S., Neumann, K. & Schweizer, A. In Dynamics of Forest Ecosystems in Central Africa during the Holocene: Past-Present-Future (ed. Runge, J.) 29–41 (Taylor and Francis, London, 2007).

57.
Lupo, K. D. et al. An elusive record exposed: radiocarbon chronology of late Holocene human settlement in the northern Congo Basin, southern Central African Republic. Azania. Archaeol. Res. Afr. 53, 209–227 (2018).
Google Scholar 

58.
Schmitt, D. N. et al. An elusive record further exposed: additional excavations and chronometric data on human settlement in the northern Congo Basin rain forest, southern Central African Republic. Azania.: Archaeological Res. Afr. 54, 55–74 (2019).
Article  Google Scholar 

59.
Eichhorn, B. In Trees, Grasses and Crops. People and Plants in Sub-Saharan Africa and Beyond (eds Eichhorn, B. & Höhn, A.) 109–117 (Dr. Rudolf Habelt GmbH, 2019).

60.
Maley, J. & Brenac, P. Vegetation dynamics, palaeoenvironments and climatic changes in the forests of western Cameroon during the last 28,000 years B.P. Rev. Palaeobot. Palynol. 99, 157–187 (1998).
Article  Google Scholar 

61.
Vincens, A. et al. Forest response to climate changes in Atlantic Equatorial Africa during the last 4000 years BP and inheritance on the modern landscapes. J. Biogeogr. 26, 879–885 (1999).

62.
Mercader, J., Rovira, S. & Gómez-Ramos, P. Shared technologies: forager-farmer interaction and ancient iron metallurgy in the Ituri rainforest, Democratic Republic of Congo. Azania. Archaeol. Res. Afr. 35, 107–122 (2000).
Google Scholar 

63.
Wang, K. et al. Ancient genomes reveal complex patterns of population movement, interaction and replacement in sub-Saharan Africa. Sci. Adv. 6, eaaz0183 (2020).

64.
Mercader, J., Martı́, R., González, I. J., Sánchez, A. & Garcı́a, P. Archaeological site formation in rain forests: insights from the Ituri Rock Shelters, Congo. J. Archaeol. Sci. 30, 45–65 (2003).
Article  Google Scholar 

65.
Patin, E. et al. Inferring the demographic history of African farmers and pygmy hunter-gatherers using a multilocus resequencing data set. PLoS Genet. 5, e1000448 (2009).
PubMed  PubMed Central  Article  CAS  Google Scholar 

66.
Tishkoff, S. A. et al. The genetic structure and history of Africans and African Americans. Science 324, 1035–1044 (2009).
CAS  PubMed  PubMed Central  Article  Google Scholar 

67.
Roberts, P. et al. Late Pleistocene to Holocene human palaeoecology in the tropical environments of coastal eastern Africa. Palaeogeogr. Palaeoclimatol. Palaeoecol. 537, 109438 (2019).

68.
Malhi, Y. & Phillips, O. L. Spatial patterns and recent trends in the climate of tropical rainforest regions. Philos. Trans. R. Soc. Lond. B Biol. Sci. 359, 311–329 (2004).
PubMed  PubMed Central  Article  Google Scholar 

69.
Maslin, M. Atmosphere. Ecological versus climatic thresholds. Science 306, 2197–2198 (2004).
CAS  PubMed  Article  PubMed Central  Google Scholar 

70.
Zhou, L. et al. Widespread decline of Congo rainforest greenness in the past decade. Nature 509, 86–90 (2014).
CAS  PubMed  Article  Google Scholar 

71.
Dargie, G. C. et al. Age, extent and carbon storage of the central Congo Basin peatland complex. Nature 542, 86–90 (2017).
CAS  PubMed  Article  Google Scholar 

72.
Cairns, J. E. et al. Adapting maize production to climate change in sub-Saharan Africa. Food Security 5, 345–360 (2013).
Article  Google Scholar 

73.
Masuka, B. et al. Gains in maize genetic improvement in Eastern and Southern Africa: I. CIMMYT hybrid breeding pipeline. Crop Sci. 57, 168–179 (2016).
Article  Google Scholar 

74.
Bharucha, Z. & Pretty, J. The roles and values of wild foods in agricultural systems. Philos. Trans. R. Soc. Lond. B Biol. Sci. 365, 2913–2926 (2010).
PubMed  PubMed Central  Article  Google Scholar 

75.
Rasmussen, L. V., Watkins, C. & Agrawal, A. Forest contributions to livelihoods in changing agriculture-forest landscapes. Policy Econ. 84, 1–8 (2017).
Article  Google Scholar 

76.
Laurance, W. F., Sayer, J. & Cassman, K. G. Agricultural expansion and its impacts on tropical nature. Trends Ecol. Evol. 29, 107–116 (2014).
PubMed  Article  Google Scholar 

77.
Dargie, G. C. et al. Congo Basin peatlands: threats and conservation priorities. Mitig. Adapt. Strateg. Glob. Change 24, 669–686 (2019).
Article  Google Scholar 

78.
Nogherotto, R., Coppola, E., Giorgi, F. & Mariotti, L. Impact of Congo Basin deforestation on the African monsoon. Atmos. Sci. Lett. 14, 45–51 (2013).
Article  Google Scholar 

79.
Brooks, E. G. E. In The Status and Distribution of Freshwater Biodiversity in Central Africa (eds Brooks, E. G. E., Allen, D. J. & Darwall, W. R. T.) 110–121 (International Union for Conservation of Nature (IUCN), 2011).

80.
Wotzka, H.-P. Studien zur Archäologie des zentralafrikanischen Regenwaldes: die Keramik des inneren Zaïre-Beckens und ihre Stellung im Kontext der Bantu-Expansion (Heinrich-Barth-Institut, 1995).

81.
Eggert, M. K. H. Imbonga and Batalimo: ceramic evidence for early settlement of the equatorial rain forest. Afr. Archaeol. Rev. 5, 129–145 (1987).
Article  Google Scholar 

82.
Longin, R. New method of collagen extraction for radiocarbon dating. Nature 230, 241–242 (1971).
CAS  PubMed  Article  PubMed Central  Google Scholar 

83.
Brock, F., Higham, T. & Ramsey, C. B. Comments on the use of Ezee-FiltersTM and ultrafilters at Orau. Radiocarbon 55, 211–212 (2013).
CAS  Article  Google Scholar 

84.
Lee-Thorp, J. A. et al. Isotopic evidence for an early shift to C4 resources by Pliocene hominins in Chad. Proc. Natl Acad. Sci. USA 109, 20369–20372 (2012).
CAS  PubMed  Article  PubMed Central  Google Scholar 

85.
Sponheimer, M. et al. Hominins, sedges, and termites: new carbon isotope data from the Sterkfontein valley and Kruger National Park. J. Hum. Evol. 48, 301–312 (2005).
PubMed  Article  PubMed Central  Google Scholar 

86.
Dunbar, E., Cook, G. T., Naysmith, P., Tripney, B. G. & Xu, S. AMS 14C dating at the Scottish Universities Environmental Research Centre (SUERC). Radiocarb. Dating Lab. Radiocarb. 58, 9–23 (2016).
CAS  Google Scholar 

87.
Bronk Ramsey, C. Bayesian analysis of radiocarbon dates. Radiocarbon 51, 337–360 (2009).
Article  Google Scholar 

88.
Reimer, P. J. et al. IntCal13 and marine13 radiocarbon age calibration curves 0–50,000 years cal BP. Radiocarbon 55, 1869–1887 (2013).
CAS  Article  Google Scholar  More

Effective medium approximation
For a linearized version of Fisher , we first obtain the effect of disorder on invasion velocity using Effective Medium Approximation (EMA). Using the EMA approach, we can replace the spatially irregular diffusion constant with a uniform one in which the effective diffusivity is everywhere equal to (D_{e}) and can replace (D_0(1+xi f(x))) in Eq. (2)44. To obtain (D_{e}), one needs to discretize Eq. (2) using a finite-difference method, which leads to the following equation for the population density39:

$$begin{aligned} frac{partial C_i(t)}{partial t}=sum _{jin {i}}W_{ij}[C_j(t)-C_i(t)]+RC_i(t) ;, end{aligned}$$
(4)

where j belongs to nearest neighbors of i and (W_{ij}=D_{ij}/delta ^2) stands for the density flow rate between units i and j with distance of (delta). Due to spatially irregular diffusion constant we have (W_{ij}=W_0(1+xi f_{ij})). Following44, one has

$$begin{aligned} D_e=bigg ( int _{0}^{infty } frac{g(w) dw}{w} bigg )^{-1} end{aligned}$$
(5)

where g(w) is probability density function for (W_{ij}). Applying this approach to our case leads to

$$begin{aligned} D_{e}=D_0(1-|xi |^2/3) end{aligned}$$
(6)

where (|xi |) stands for dimensionless magnitude of (xi) ((xi ^2) has a physical dimension of length, meter).
Perturbation analysis for invasion front fluctuations
The first step towards the study of dynamics of a propagating front is linearizing (3) by neglecting the (C^2) term in an environment with effect diffusion constant, (D_e), as:

$$begin{aligned} frac{partial C}{partial t}=RC+frac{partial }{partial x} big ((D_e + xi f(x)) frac{partial }{partial x} Cbig ) end{aligned}$$
(7)

This is based on the fact that near the front, the cell density is (C ll 1). In other words, focusing on the dynamics of the front position, automatically grants us the possibility of linearizing (3).
The construction of the solution can be proceeded according to a valuable insight given in the classic paper43, where the particle density is written as follows

$$begin{aligned} C(zeta ,t) approx C_0(zeta +eta (t),t)+ delta C_1(zeta ,t) end{aligned}$$
(8)

where C is written in the comoving frame and (zeta = x-vt). It is assumed that (delta C_1 ll 1) and in the same order as the perturbing function. So that terms containing f(x) and (delta C_1) can be neglected. Furthermore, (C_0) is assumed to satisfy the linearized Eq. (3) with (xi =0), i.e.

$$begin{aligned} frac{partial C_0(zeta ,t)}{partial t}-hat{Gamma } C_0(zeta )=,& {} frac{partial C_0(zeta ,t)}{partial t} nonumber \&-bigg (D_e frac{d^2}{dzeta ^2} + vfrac{d}{dzeta } + Rbigg )C_0(zeta ,t) = 0 end{aligned}$$
(9)

Which has the following solution

$$begin{aligned} C_0(zeta ,t) = frac{1}{sqrt{4pi D_e t}}e^{-frac{1}{2}sqrt{frac{R}{D_e}}zeta }e^{-frac{zeta ^2}{4D_e t}} end{aligned}$$
(10)

The first term in (8) describes the effects of the perturbing function f(x) on the position of the propagating front, while the second term shows the change in the shape of the front. This approach has also been employed and well explained in a recent paper17. As shown in17,43, to determine the effective diffusion coefficient for the fluctuating front, it is sufficient to solve (3) using (8) for (eta (t)). Note also that since we are interested in the dynamics of the system in long times ((t gg frac{1}{R})), (v_e) can be assumed to be equal to (2sqrt{R D_e})45. Plugging (8) in Eq. (9) expressed in comoving coordinates and considering (xi =bar{xi }/D_e) yields

$$begin{aligned} frac{partial delta C_1}{partial t} – hat{Gamma }delta C_1 + dot{eta (t)} C_0(zeta ,t) = bar{xi } bigg (f(zeta ) C’_0(zeta ,t)bigg )’ end{aligned}$$
(11)

Noting that the operator (hat{Gamma }) is not self-adjoint (The adjoint of (hat{Gamma }) is: (hat{Gamma ^dagger }=D_e dfrac{d^2}{dzeta ^2} – v_e dfrac{d}{dzeta } + R)) and following43, we multiply Eq. (11) from the left in the eigenfunction of (hat{Gamma ^dagger }) with 0 eigenvalue (which is (e^{sqrt{frac{R}{D_e}}zeta })) and integrate. Thus,

$$begin{aligned} {,} & int _{-infty }^{infty } e^{sqrt{frac{R}{D_e}}zeta }frac{partial delta C_1(zeta ,t)}{partial t} dzeta + dot{eta (t)}int _{-infty }^{infty }e^{sqrt{frac{R}{D_e}}zeta }C’_0(zeta ,t) dzeta nonumber \&quad = bar{xi } int _{-infty }^{infty } e^{sqrt{frac{R}{D_e}}zeta }bigg (f(zeta ) C’_0(zeta ,t)bigg )’dzeta end{aligned}$$
(12)

Which yields

$$begin{aligned} dot{eta }(t) =bar{xi } dfrac{int _{-infty }^{infty } e^{sqrt{frac{R}{D_e}}zeta }bigg (f(zeta ) C’_0(zeta ,t)bigg )’ dzeta }{int _{-infty }^{infty }e^{sqrt{frac{R}{D_e}}zeta }C’_0(zeta ,t) dzeta } end{aligned}$$
(13)

Which can further be simplified into

$$begin{aligned} dot{eta }(t)=,& {} bar{xi }dfrac{int _{-infty }^{infty } e^{sqrt{frac{R}{D_e}}zeta }f(zeta ) C’_0(zeta ,t) dzeta }{int _{-infty }^{infty }e^{sqrt{frac{R}{D_e}}zeta }C_0(zeta ,t) dzeta } nonumber \=,& {} bar{xi } e^{-frac{R t}{4}}int _{-infty }^{infty } e^{sqrt{frac{R}{D_e}}zeta }f(zeta ) C’_0(zeta ,t) dzeta end{aligned}$$
(14)

Or equivalently,

$$begin{aligned} eta (t) =bar{xi } int _{0}^{t} dtau e^{-frac{Rtau }{4}} int _{-infty }^{infty } dzeta e^{sqrt{frac{R}{D_e}}zeta }f(zeta ) C’_0(zeta ,tau ) end{aligned}$$
(15)

According to17, the effective diffusion would be given by

$$begin{aligned} D_{C} = dfrac{langle eta ^2(t) rangle }{2t} end{aligned}$$
(16)

Or

$$begin{aligned} D_{C}=,& {} frac{bar{xi }^2}{2t}int _{0}^{t}d{T_1}int _{0}^{t}d{T_2}int _{-infty }^{infty } C’_0(zeta ,T_1)C’_0(zeta ,T_2) nonumber \&times e^{-frac{R T_1}{4}}e^{-frac{R T_2}{4}}e^{2sqrt{frac{R}{D_e}}zeta } dzeta end{aligned}$$
(17)

where we have performed an ensemble average over (eta ^2(t)) using the fact that (langle f(x)f(y) rangle = delta (x-y)). A numerical calculation of (16) can be readily computed using any mathematical software. However, valuable insight can still be obtained from (17), if we use dimensionless parameters (tau _i=frac{T_i}{t}) and (sigma =sqrt{frac{R}{D_0}}zeta). In other words

$$begin{aligned} D_{C}=,& {} bar{xi }^2 dfrac{sqrt{R}}{32pi D^{3/2}_0}int _{0}^{1}dtau _1int _{0}^{1}dtau _2int _{-infty }^{infty }dsigma nonumber \&times dfrac{left( 1+dfrac{sigma }{Rttau _1}right) left( 1+dfrac{sigma }{Rttau _2}right) }{sqrt{tau _1tau _2}}e^{-frac{sigma ^2}{4Rttau _1}}e^{-frac{sigma ^2}{4Rttau _2}}e^{sigma }e^{-R tfrac{(tau _1+tau _2)}{4}} end{aligned}$$
(18)

Equation (18) gives the effective diffusion coefficient for the stochastic behavior of the front. For a diffusive behavior, we would expect this effective diffusion coefficient to tend to a constant at large times. At large times, we can approximate the integral as follows,

$$begin{aligned} D_{C} approx bar{xi }^2 dfrac{sqrt{R}}{32pi D^{3/2}_e}int _{0}^{1}dtau _1int _{0}^{1}dtau _2int _{-infty }^{infty }dsigma dfrac{1}{sqrt{tau _1tau _2}}e^{-frac{sigma ^2}{4Rttau _1}} e^{-frac{sigma ^2}{4Rttau _2}}e^{sigma }e^{-R tfrac{(tau _1+tau _2)}{4}} end{aligned}$$
(19)

Luckily, Eq. (19) can be evaluated exactly to yield

$$begin{aligned} D_{C}&approx bar{xi }^2 dfrac{sqrt{R}}{32pi D^{3/2}_e} nonumber \&quad frac{2 pi (2 R t-1) text {Erf}left( frac{sqrt{R t}}{2}right) -2 pi e^{2 R t} text {Erf}left( frac{3 sqrt{R t}}{2}right) +2 pi e^{2 R t} text {Erf}left( sqrt{2} sqrt{R t}right) +8 sqrt{pi } e^{-frac{1}{4} (R t)} sqrt{R t}-4 sqrt{2 pi } sqrt{R t}}{R t} end{aligned}$$
(20)

Where Erf(x) is the error function. As (trightarrow infty) this gives the following simple relation for the diffusion constant for the wave front

$$begin{aligned} D_{C}(trightarrow infty ) = bar{xi }^2 dfrac{sqrt{R}}{8 D^{3/2}_e} end{aligned}$$
(21)

Substituting (xi =bar{xi }/D_e) and (D_e=D_0(1-|xi |^2/3)), we will get the following beautiful equation for the effective diffusion constant of the front at large times:

$$begin{aligned} D_{C}(trightarrow infty ) =dfrac{1}{8} xi ^2 sqrt{R D_0(1- |xi |^2/3)}. end{aligned}$$
(22) More

The early warnings of the outbreak
As early as Dec 31st, 2019, when Wuhan Municipal Health Commission first informed the public about the emerging pneumonia cases21, most of the cities (326 out of 346) exhibited at least some awareness of the emerging SARS-CoV-2 outbreak (Fig. 2b). However, awareness then decreased until Jan 19th, 2020, one day before the Chinese Centre for Disease Control and Prevention confirmed human-to-human transmissions of the novel coronavirus9. Since Jan 20th, 2020, overall awareness increased by a magnitude of at least five, demonstrating significant awareness across all cities (Fig. 2b). Awareness remained low as the epidemic spread, falling close to its lowest point on the starting day of Chunyun (Jan 10th, 2020). Considering cities that showed initial novel coronavirus awareness levels at least 1.5 times that of the search term “common cold”, we found a total of 166 alert cities as early as Dec 31st, 2019 (48 cities at a tighter threshold of (C = 3.0) times, illustrated in Fig. 2a). However, awareness decreased significantly during Chunyun.
Figure 2

Public awareness over time. (a) The frequency distributions of cities that exhibit the first significant signal of awareness over time. The number of cities for which searches for the combined term “Wuhan” and “pneumonia” exceed (user2{ C} = 3) times the search term “common cold” is reported every day. (b) Public awareness on the topic of “pneumonia” over time. All 346 cities exhibit at least some searches of the term “pneumonia” during the initial outbreak period. Of these, 326 cities recorded searches about it as early as Dec 31st, 2019. Cities are divided into two groups according to whether or not they had reported SARS cases in 2003–04. The mean values of awareness magnitude were computed on a daily basis for two groups of cities respectively. Accordingly, a paired t-test was performed on those two time-series, and we found the cities that had reported SARS cases had greater of awareness (t-statistic: 3.56; degrees of freedom: 23; p  More

1.
Waggoner, P. E. Agriculture and a climate changed by more carbon dioxide. In Changing Climate: Report of the Carbon Dioxide Assessment Committee (ed. National Research Council (NRC)) 383–418 (National Academy Press, Washington, DC, 1983). https://doi.org/10.1038/s41598-020-75350-9.
Google Scholar 
2.
Shinozaki, K., Yamaguchi-Shinozaki, K. & Seki, M. Regulatory network of gene expression in the drought and cold stress responses. J. Curr. Opin. Plant Biol. 6, 410–417 (2003).
CAS  Google Scholar 

3.
Deng, X. et al. Improving agricultural water use efficiency in arid and semiarid areas of China. J. Agric. Water Manag. 80, 23–40 (2006).
Google Scholar 

4.
Hussain, M. et al. Improving drought tolerance by exogenous application of glycinebetaine and salicylic acid in sunflower. J. Agron. Crop Sci. 194, 193–199 (2008).
CAS  Google Scholar 

5.
Wang, W., Vinocur, B. & Altman, A. Plant responses to drought, salinity and extreme temperatures: Towards genetic engineering for stress tolerance. J. Planta. 218, 1–14 (2003).
CAS  Google Scholar 

6.
Bartels, D. & Sunkar, R. Drought and salt tolerance in plants. J. Crit. Rev. Plant Sci. 24, 23–58 (2005).
CAS  Google Scholar 

7.
Rost, S. et al. Global potential to increase crop production through water management in rainfed agriculture. J. Environ. Res. Lett. 4, 1–9 (2009).
Google Scholar 

8.
Miflin, B. Crop improvement in the 21st century. J Exp Bot. 51, 1–8 (2000).
CAS  PubMed  Google Scholar 

9.
Islam, M. R. et al. Impact of water-saving superabsorbent polymer on oat (Avena spp.) yield and quality in an arid sandy soil. J. Sci. Res. Essays. 6, 720–728 (2011).
CAS  Google Scholar 

10.
Lu, H. et al. Earliest domestication of common millet (Panicum miliaceum) in East Asia extended to 10000 years ago. J. Proc. Natl. Acad. Sci. 106, 7367–7372 (2009).
ADS  CAS  Google Scholar 

11.
Anonymous. Plant area of millet in China in 2017. Available from: https://data.chinabaogao.com/nonglinmuyu/2019/0R43PA2019.html (2019).

12.
Tolk, J., Howell, T. & Evett, S. Effect of mulch, irrigation, and soil type on water use and yield of maize. J. Soil Tilll Res. 50, 137–147 (1999).
Google Scholar 

13.
Panda, R., Behera, S. & Kashyap, P. Effective management of irrigation water for maize under stressed conditions. J. Agric. Water Manag. 66, 181–203 (2004).
Google Scholar 

14.
Cattivelli, L. et al. Drought tolerance improvement in crop plants: An integrated view from breeding to genomics. J. Field Crop Res. 105, 1–14 (2008).
Google Scholar 

15.
Chaves, M. M. et al. How plants cope with water stress in the field? Photosynthesis and growth. J. Ann. Bot.-Lond. 89, 907–916 (2002).
CAS  Google Scholar 

16.
Xu, C. et al. Effect of biochar amendment on yield and photosynthesis of peanut on two types of soils. J. Environ. Sci. Pollut. R. 22, 6112–6125 (2015).
CAS  Google Scholar 

17.
Shah, N. & Paulsen, G. Interaction of drought and high temperature on photosynthesis and grain-filling of wheat. J. Plant Soil. 257, 219–226 (2003).
CAS  Google Scholar 

18.
Soda, W. et al. Co-composting of acid waste bentonites and their effects on soil properties and crop biomass. J. Environ. Qual. 35, 2293–2301 (2006).
CAS  PubMed  Google Scholar 

19.
Guiwei, Q., De Varennes, A. & Cunha-Queda, C. Remediation of a mine soil with insoluble polyacrylate polymers enhances soil quality and plant growth. J. Soil Use Manag. 24, 350–356 (2008).
Google Scholar 

20.
Hüttermann, A., Orikiriza, L. J. B. & Agaba, H. Application of superabsorbent polymers for improving the ecological chemistry of degraded or polluted lands. J. Clean Soil. Air. Water. 37, 517–526 (2009).
Google Scholar 

21.
Narjary, B. et al. Water availability in different soils in relation to hydrogel application. J. Geoderma 187, 94–101 (2012).
ADS  Google Scholar 

22.
Arbona, V. et al. Hydrogel substrate amendment alleviates drought effects on young citrus plants. J. Plant Soil. 270, 73–82 (2005).
CAS  Google Scholar 

23.
Wang, Y. et al. Influences of Millet–Peanut intercropping on photosynthetic characteristics and yield of Millet. J. Agric. Sci. Technol. 22, 153–165 (2020).
Google Scholar 

24.
Gao, H. et al. Diurnal change of photosynthetic characteristics and response tonight intensity of seven ornamental grasses. J. Acta Prataculturae Sinica 19, 87–93 (2010).
Google Scholar 

25.
Yordanov, I., Velikova, V. & Tsonev, T. Plant responses to drought, acclimation, and stress tolerance. J. Photosynthetica 38, 171–186 (2000).
CAS  Google Scholar 

26.
Yangwei, P. & Yan, S. Resources characteristics and market situation of bentonites at home and Abroad. J. Metal Mine. 95–99, 105 (2012).
Google Scholar 

27.
Fraenkel-Conrat, H., Singer, B. & Tsugita, A. Purification of viral RNA by means of bentonite. J. Virol. 14, 54–58 (1961).
CAS  Google Scholar 

28.
Lopez-Fernandez M, et al. Microbial community changes induced by uranyl nitrate in bentonite clay microcosms. In 16th International Clay Conference (ICC). Granada, Spain (2017).

29.
Bentahar, S. et al. Removal of a cationic dye from aqueous solution by natural clay. J. Groundw. Sustain. Devel. 6, 255–262 (2018).
Google Scholar 

30.
De Castro, M. L. F. A. et al. Adsorption of Methylene Blue dye and Cu(II) ions on EDTA-modified bentonite: Isotherm, kinetic and thermodynamic studies. J. Sustain. Environ. Res. 28, 197–205 (2018).
Google Scholar 

31.
Mi, J. et al. Effect of bentonite amendment on soil hydraulic parameters and millet crop performance in a semi-arid region. J. Field Crop Res. 212, 107–114 (2017).
Google Scholar 

32.
Hall, D. J. M. et al. Claying and deep ripping can increase crop yields and profits on water repellent sands with marginal fertility in southern Western Australia. J. Aust. J. Soil Res. 48, 178–187 (2010).
Google Scholar 

33.
Shi, Y., Chen, X. & Shen, S. Mechanisms of organic cementing soil aggregate formation and its theoretical models. J. Chin. J. Appl. Ecol. 13, 1495–1498 (2002).
Google Scholar 

34.
Mi, J. et al. Effect of sandy soil amendment on Dry-farmland water-conserving characteristic and millet seeding growth. J. Irrigat. Drain. 23, 92–96 (2015).
Google Scholar 

35.
Mi, J. et al. Effects of sandy soil amendment on soil moisture and growth status of millet with rainfed sandy soil in a semi-arid region. J. Adv. Mater. Res. 1092–1093, 1234–1242 (2015).
Google Scholar 

36.
Jiang, P. et al. Principles and experimental verification of capillary suction method for fast measurement of field capacity. J. Trans. CSAE. 22, 1–5 (2006).
CAS  Google Scholar 

37.
Andrenelli, M. et al. Field application of pelletized biochar: Short term effect on the hydrological properties of a silty clay loam soil. J. Agric. Water Manag. 163, 190–196 (2016).
Google Scholar 

38.
Jones, J. Jr. Kjeldahl Method for Nitrogen Determination (Micro-Macro Publishing, Athens, 1991).
Google Scholar 

39.
Nielsen, S. S. Food Analysis (Aspen Publishers Inc, New York, 1998).
Google Scholar 

40.
Ganzler, K., Salgo, A. & Valkó, K. Microwave extraction: A novel sample preparation method for chromatography. J. Chromatogr. A 371, 299–306 (1986).
CAS  Google Scholar 

41.
Min D B, Steenson D, Crude fat analysis. In Food Analysis. Vol. 2. New York: Kluwer Academic/Plenum Publishers. 113–131 (1998).

42.
Möller, J. Comparing Methods for Fibre Determination in Food and Feed, Vol. 1026712 (The Association of American Feed Control Officials, West Lafayette, 2014).
Google Scholar 

43.
Bakass, M., Mokhlisse, A. & Lallemant, M. Absorption and desorption of liquid water by a superabsorbent polyelectrolyte: Role of polymer on the capacity for absorption of a ground. J Appl Polym. Sci. 82, 1541–1548 (2001).
CAS  Google Scholar 

44.
Yu, J. et al. Superabsorbents and semiarid soil properties affecting water absorption. J. Soil. Sci. Soc. Am. J. 75, 2305–2313 (2011).
ADS  CAS  Google Scholar 

45.
Sun, Z. et al. Effect of ion-type and concentration on water-retention capacity of bentonite used in geosynthetic clay liner. J. Chin. Ceram. Soc. 38, 1826–1831 (2010).
Google Scholar 

46.
Suzuki, S. et al. Improvement in water-holding capacity and structural stability of a sandy soil in Northeast Thailand. J. Arid Land Res. Manag. 21, 37–49 (2007).
Google Scholar 

47.
Betti, G. et al. Size of subsoil clods affects soil–water availability in sand–clay mixtures. J. Soil Res. 54, 276–290 (2016).
Google Scholar 

48.
Al-Omran, A. et al. Impact of natural deposits of Saudi Arabia on selected physical properties of calcareous sandy soil. J. Dirasa Agric Sci. 29, 285–294 (2002).
Google Scholar 

49.
Tan, G. et al. Effect of super absorbent resin on the rate of maize emergence and soil moisture. J. Jilin Agric. Sci. 30, 26–27 (2004).
Google Scholar 

50.
Zhou, L. et al. Effect of bentonite-humic acid application on the improvement of soil structure and maize yield in a sandy soil of a semi-arid region. J. Geoderma. 338, 269–280 (2019).
ADS  CAS  Google Scholar 

51.
Fang, S. et al. Synthesis of chitosan derivative graft acrylic acid superabsorbent polymers and its application as water retaining agent. J. Int. J. Biol. Macromol. 115, 754–761 (2018).
CAS  Google Scholar 

52.
Zhang, X. et al. Dry matter, harvest index, grain yield and water use efficiency as affected by water supply in winter wheat. J. Irrigation Sci. 27, 1–10 (2008).
Google Scholar 

53.
Giri, G. S. & Schillinger, W. F. Seed priming winter wheat for germination, emergence, and yield. J. Crop Sci. 43, 2135–2141 (2003).
Google Scholar 

54.
Jamnická, G. et al. The soil hydrogel improved photosynthetic performance of beech seedlings treated under drought. J. Plant Soil Environ. 59, 446–451 (2013).
Google Scholar 

55.
Islam, M. R. et al. A lysimeter study of nitrate leaching, optimum fertilisation rate and growth responses of corn (Zea mays L.) following soil amendment with water-saving super-absorbent polymer. J. Sci/ Food Agric. 91, 1990–1997 (2011).
CAS  Google Scholar 

56.
Tahir, S. & Marschner, P. Clay amendment to sandy soil—Effect of clay concentration and ped size on nutrient dynamics after residue addition. J. Soils Sediments 16, 2072–2080 (2016).
CAS  Google Scholar 

57.
Croker, J. et al. Effects of recycled bentonite addition on soil properties, plant growth and nutrient uptake in a tropical sandy soil. J. Plant Soil. 267, 155–163 (2004).
CAS  Google Scholar 

58.
Brix, H. The effect of water stress on the rates of photosynthesis and respiration in tomato plants and loblolly pine seedlings. J. Physiol. Plantarum. 15, 10–20 (1962).
Google Scholar 

59.
Wang, L. et al. Progress in researches on effect of iron promoting accumulation of soil organic carbon. J. Acta Pedologica Sinica 55, 1041–1050 (2018).
Google Scholar 

60.
Taylor, J. R. N. & Duodu, K. G. Sorghum and millets: Grain-quality characteristics and management of quality requirements. In Cereal Grains (2nd edn) (eds Batey, I. & Miskelly, D.) 317–351 (Woodhead Publishing, Wrigley, Colin, 2017).
Google Scholar 

61.
Abbas, M. et al. Effect of some soil amendments on yield and quality traits of sugar beet (Beta vulgaris L.) under water stress in sandy soil. J. Egypt. J. Agron. 40, 75–88 (2018).
Google Scholar 

62.
Leila, K., Hassan, F. & Pooran, G. Effects of different irrigation and superabsorbent levels on physio-morphological traits and forage yield of millet (Pennisetum americanum L.). J. Am. Euras. J. Agric. Environ. Sci. 13, 1043–1049 (2013).
Google Scholar 

63.
Li, X. et al. Effects of super-absorbent polymers on a soil–wheat (Triticum aestivum L.) system in the field. J. Appl Soil Ecol. 73, 58–63 (2014).
ADS  Google Scholar 

64.
Sojka, R. et al. Polyacrylamide in agriculture and environmental land management. J. Adv Agron. 92, 75–162 (2007).
CAS  Google Scholar 

65.
Wang, A., Li, F. & Yang, S. Effect of polyacrylamide application on runoff, erosion, and soil nutrient loss under simulated rainfall. J. Pedosphere. 21, 628–638 (2011).
CAS  Google Scholar  More

Fish and vaccines
Two pedigree populations of Atlantic salmon (Salmo salar) called Fanad and Lochy were used in this study18 (Table 1). The two populations were managed separately and had different origins. Fish were provided in 2016 by the salmon fish farming company Salmones Camanchaca and pit tagged in April 2016 at an average weight of 26.4 ± 3.9 g and 30.2 ± 4.2 g, for populations Fanad and Lochy, respectively. During the freshwater growth period, salmon were immunized twice using commercial vaccines, following the strict Salmones Camanchaca protocols. First, fish were vaccinated by intraperitoneal (IP) injection with a pentavalent vaccine against P. salmonis, Vibrio ordalii, A. salmonicida, IPNV (infectious pancreatic necrosis virus) and ISAV (infectious salmon anemia virus). Second, fish were immunized by IP injection against P. salmonis using a monovalent live attenuated vaccine at the same time as the first vaccination. Since 2016, this double vaccination strategy has been a common practice in the Chilean salmon industry17. Fish were transferred as smolts to the Aquadvice experimental station in Puerto Montt, Chile. Unvaccinated fish were injected with PBS (phosphate-buffered saline) and used as control (Table 1). Prior to transferring the fish, a health check by RT-PCR was performed to verify that the fish were free of viral (ISAV and IPNV) and bacterial pathogens (Vibrio sp., Flavobacterium sp., P. salmonis, and Renibacterium salmoninarum). At the experimental station, all fish underwent a 15 days acclimatization period in seawater (salinity of 32% and a temperature of 15 ± 1 °C). Fish were fed daily ad libitum with a commercial diet.
Calculation of Piscirickettsia salmonis LD50
The median lethal dose (LD50) of P. salmonis (EM-90 type) was determined as previously described18. Briefly, animals from both populations were distributed in eight tanks of 350 L (n = 60 fish per tank). The LD50 was calculated in fish infected by IP injection with 200 μL of a P. salmonis suspension. Three dilutions were assessed from stock with concentrations of 1 × 106.63 TCID/mL (TCID = median tissue culture infective dose): 1 × 10–3 TCID/mL, 1 × 10–4 TCID/mL, and 1 × 10–5 TCID/mL. Controls were injected with 200 μL of PBS. Fish were monitored daily for 30 days, and mortalities were recorded. The presence of bacteria was assessed by qRT-PCR. In both infection scenarios, a single infection and coinfection, the highest dose of P. salmonis was used (1 × 10–3 TCID/mL) as a conservative measure because the fish grow about 100 g between LD50 and the main challenge (50 days).
Infection design, trait of resistance and protection added by vaccine
Fish were treated with two different types of infection, a single infection with P. salmonis (PS) or coinfection with both C. rogercresseyi and P. salmonis (CAL + PS) as previously described18. In short, infections against P. salmonis occurred at 822 ATU (accumulated thermal units) within the immunization period described by the vaccine manufacturer. Vaccinated and unvaccinated fish from populations Fanad and Lochy were equally distributed in four tanks of 6 m3, with two replicates per type of infection. For the single infection with P. salmonis, fish were IP injected. For the coinfection, fish were exposed first to sea lice and then to P. salmonis. A coinfection procedure was established based on our previous experience with this study model45,55. Infections with sea lice were performed by adding 60 copepodites per fish to each tank of coinfection. Copepodites were collected from egg-bearing females reared in the laboratory and confirmed as “pathogen-free” (P. salmonis, R. salmoninarum, IPNV, and ISAV) by RT-PCR diagnostic. After the addition of parasites, water flow was stopped for a period of 8 h, and tanks were covered to decrease light intensity, which favors a successful settlement of sea lice on fish55. A placebo procedure was applied to single infection tanks, keeping them in darkness and controlling the volume of water, temperature, oxygen levels, and fish density equivalent to those that were measured in coinfected tanks18. The secondary infection was performed with P. salmonis after seven days of sea lice infestation, and the establishment of the parasites was confirmed and quantified on all fish. Therefore, our experimental design had two types of treatments: (1) single infection (PS) or coinfection (CAL + PS); and (2) vaccinated or unvaccinated fish. Vaccinated and unvaccinated fish with a single infection were distributed in tanks 1 and 2, and vaccinated and unvaccinated fish with a coinfection were distributed in tanks 3 and 4. Further, fish were fasted for one day prior to each procedure to minimize the detrimental effects of stress on water quality parameters. Finally, fish were sedated with AQUI-S (50% Isoeugenol, 17 mL/100 L water) to reduce stress during handling. Fish were monitored daily for 30 days, and resistance to P. salmonis was measured individually as days to death. Protection added by vaccination was calculated as the difference of resistance between vaccinated fish and their unvaccinated full-sibs and represented under a single Genetic and Environment model (GxE, G = full-sib family; E = Vaccination treatment).
Comparison of moribund and survivor fish
Bacterial load, growth, and macroscopic lesions were evaluated in survivors and moribund fish. Moribund fish were obtained as dying fish when 50% of mortality was reached in both a single infection and coinfection treatments. Moribund fish were recognized and collected by three behavioral traits: lethargy, no response to stimuli, and slow swimming close to the tank wall. Resistance to P. salmonis was measured by days to death and mortality (alive versus dead) and monitored for 30 days15,45, survivors fish comprised those that lived at the end of experiment15. Forty fish were collected from each group of moribund and survivors, and from each treatment (PS and CAL + PS) and comparisons were performed between unvaccinated and vaccinated fish, twenty fish each group. However, due to the low number of unvaccinated survivors fish coinfected with P. salmonis and sea lice, it was not possible to compare with the vaccinated survivors fish.
Specific growth rate (SGR)
SGR was evaluated for moribund and survivors fish. The specific growth rate was calculated previous to infection, and post-infection as SGR = ((lnw2 − lnw1)*t−1)*100, where w2 corresponds to final weight, w1 to the initial weight, and t corresponds to the number of days between infection and death of the fish or the end of the trial if they survived56.
Piscirickettsia salmonis load
Piscirickettsia salmonis load was evaluated for moribund and survivors fish. P. salmonis load was estimated based on the amount of specific ribosomal RNA from the bacteria in the head kidneys of the infected fish, as measured by qRT-PCR. Dead fish were not used to evaluate bacterial load. Threshold cycle (CT) values from bacterial RNA was used as an indication of the bacterial load as previously described18. Head kidney samples were extracted from 20 moribund and survivors fish per group and preserved in RNAlater at − 80 °C until RNA extraction. RNA was extracted from tissue samples with the TRIzol reagent (Thermo Fisher Scientific, MA, USA) following the instructions provided by the manufacturer. DNA was removed through an additional step using a DNase incubation for 60 min at 37 °C. The quality of the RNA extraction was checked by visualizing the 28S and 18S rRNA bands resolved in 1% of agarose gels stained with SYBR Safe DNA gel stain (Invitrogen, CA, USA), and the total concentration of the RNA was measured spectrophotometrically in a MaestroNano device (MAESTROGEN, Hsinchu, Taiwan). One hundred nanograms of purified total RNA was used for the qRT-PCR reactions. The qRT-PCR reaction was prepared using the Brilliant III SYBR master mix (Agilent Technologies, CA, USA) by adding the template RNA, probes, and primers as described previously57. qRT-PCR was performed in the Eco Real-Time PCR system (Illumina, CA, USA), whose results were expressed in terms of CT. All samples were tested in triplicates and were calibrated to a plate standard that contained a combination of samples from all groups tested. Primers used for 23S gene of S. salar were forward primer TCTGGGAAGTGTGGCGATAGA and reverse primer TCCCGACCTACTCTTGTTTCATC.
Necropsy analysis
Macroscopic lesions from 20 fish per treatment were analyzed on moribund and survivors fish13; almost all survivors sampled fish were vaccinated, except one unvaccinated fish that survived to P. salmonis infection (data not shown). Fresh samples were analyzed by two veterinarians who were blinded to the treatments. Macroscopic lesions evaluated in the tissues were peeling or undergoing desquamation, congestion, and ecchymosis in the skin, paleness, and melanomacrophages in the gills, white hepatic nodules, hepatomegaly, spleen paleness, and splenomegaly. Macroscopic lesions were indicated as present or absent.
Statistical analysis
Significance levels of resistant to P. salmonis were obtained using a two-way ANOVA followed by a Tukey post-hoc test and unpaired t-test. The effects of populations and sex of fish on SGR and P. salmonis load were analyzed using a non-parametric Kruskal–Wallis test followed by a Dunn post-hoc test. Additionally, differences in the clinical signs of the P. salmonis infection between different treatments were analyzed using a non-parametric Chi-square proportion. All statistical analyses were performed using R Core Team (RStudio, Vienna, Austria). Graphs were designed with GraphPad Prism 8.0 software (GraphPad Software, CA, USA).
Quantitative genetic analysis
Each population in this study has a different genetic origin and has been managed as closed populations during the domestication process. Thus, (co) variance components of days to death were estimated independently for each population from the data of its genealogy (Table 1) using VCE 6.0 software by Groeneveld et al.58.
Heritability of days to death was estimated using the following univariate animal model:

$$y = upsilon 1 , + X_{1} t + X_{2} i + X_{3} v + Za + e,,,,{text{Model}},1$$

where y is the vector of the trait days to death, μ is the overall mean effect, t is the fixed effect of tank; i is the fixed effect of type of infection; v is the fixed effect of group of vaccination; a is the random effects vector of animal effects, with a ~ N(0, σa2A); and e is the random vector of errors, with e ~ N(0, σe2Ie). X1, X2, X3, and Z are incidence matrices, and A is the numerator relationship matrix obtained from pedigree information. The magnitude of estimated heritability was established following the classification of Cardellino and Rovira59: low (0.05–0.15), medium (0.20–0.40), and high (0.45–0.60) and very high ( > 0.65).
Genotype–environment interactions (GxE) were estimated by means of genetic correlations between the trait days to death measured in one environment (i.e., unvaccinated and single infection with P. salmonis) and the same trait measured in the other environment (i.e., vaccinated and coinfection).
Genetic correlations were estimated using the following bivariate animal model:

$$y_{1} ,;y_{2} = X_{1} d + X_{2} t(d) + X_{3} i(d) + X_{4} v(d) + Za(d) + e,,,,{text{Model}},1$$

where, y1 and y2 are the data vectors for the traits of interest (days to death in vaccinated and unvaccinated fish); d is the fixed vector of trait effects; t(d), i(d), v(d), are the fixed effects of tank, type of infection and group of vaccination effects within trait, respectively; a(d) is the random vector of animal effects within trait, with a(d) ~ N(0, A ⊗ G); and e is the random vector of errors, with e ~ N(0, I ⊗ R). The matrix G is a 2 × 2 variance–covariance matrix between traits defined by a genetic additive correlation term, rg, and a genetic variance (σgj2) for each trait. The matrix R is an unstructured 2 × 2 residual variance–covariance matrix with a different variance for each trait (σej2), and a covariance between traits (σeij). All other terms were previously defined. Correlations were classified as low (0–0.39), medium (0.40–0.59), high (0.60–0.79), and very high (0.80–1), regardless whether it was positive or negative. Significance testing of the estimates of heritability and genetic correlation were approximate as suggest by Åkesson et al.60. Thus, any genetic parameter value was considered significantly different from zero with P  More

1.
Thomas, C. D. et al. Extinction risk from climate change. Nature 427, 145–148 (2004).
CAS  Google Scholar 
2.
Alexander, J. M., Diez, J. M. & Levine, J. M. Novel competitors shape species’ responses to climate change. Nature 525, 515–518 (2015).
CAS  PubMed  Google Scholar 

3.
Seidl, R., Schelhaas, M.-J., Rammer, W. & Verkerk, P. J. Increasing forest disturbances in Europe and their impact on carbon storage. Nat. Clim. Change 4, 806–810 (2014).
CAS  Google Scholar 

4.
Post, E. & Forchhammer, M. C. Synchronization of animal population dynamics by large-scale climate. Nature 420, 168–171 (2002).
CAS  PubMed  Google Scholar 

5.
Koenig, W. D. & Liebhold, A. M. Temporally increasing spatial synchrony of North American temperature and bird populations. Nat. Clim. Change 6, 614–617 (2016).
Google Scholar 

6.
Shestakova, T. et al. Forests synchronize their growth in contrasting Eurasian regions in response to climate warming. Proc. Natl Acad. Sci. USA 113, 662–667 (2016).
CAS  PubMed  Google Scholar 

7.
Black, B. A. et al. Rising synchrony controls western North American ecosystems. Glob. Change Biol. 24, 2305–2314 (2018).
Google Scholar 

8.
Heino, M. Noise colour, synchrony and extinctions in spatially structured populations. Oikos 83, 368–375 (1998).
Google Scholar 

9.
Liebhold, A., Koenig, W. D. & Bjørnstad, O. N. Spatial synchrony in population dynamics. Annu. Rev. Ecol. Evol. Syst. 35, 467–490 (2004).
Google Scholar 

10.
Gouhier, T. C., Guichard, F. & González, A. Synchrony and stability of food webs in metacommunities. Am. Nat. 175, E16–E34 (2010).
PubMed  Google Scholar 

11.
Elton, C. S. Periodic fluctuations in the numbers of animals: their causes and effects. Br. J. Exp. Bot. 2, 119–163 (1924).
Google Scholar 

12.
Moran, P. A. P. The statistical analysis of the Canadian lynx cycle. Aust. J. Zool. 1, 291–298 (1953).
Google Scholar 

13.
Loreau, M. & de Mazancourt, C. Species synchrony and its drivers: neutral and nonneutral community dynamics in fluctuating environments. Am. Nat. 172, E48–E66 (2008).
PubMed  Google Scholar 

14.
Buma, B. et al. The value of linking paleoecological and neoecological perspectives to understand spatially-explicit ecosystem resilience. Landsc. Ecol. 34, 17–33 (2018).
Google Scholar 

15.
Hughes, B. B. et al. Long-term studies contribute disproportionately to ecology and policy. BioScience 67, 271–281 (2017).
Google Scholar 

16.
Gajewski, K., Viau, A. E., Sawada, M., Atkinson, D. E. & Fines, P. Synchronicity in climate and vegetation transitions between Europe and North America during the Holocene. Clim. Change 78, 341–361 (2006).
CAS  Google Scholar 

17.
Zhao, S. et al. The International Tree-Ring Data Bank (ITRDB) revisited: data availability and global ecological representativity. J. Biogeogr. https://doi.org/10.1111/jbi.13488 (2018).

18.
Babst, F., Poulter, B., Bodesheim, P., Mahecha, M. D. & Frank, D. C. Improved tree-ring archives will support earth-system science. Nat. Ecol. Evol. 1, 0008 (2017).
Google Scholar 

19.
Cook, E. R. & Peters, K. Calculating unbiased tree-ring indices for the study of climatic and environmental change. Holocene 7, 361–370 (1997).
Google Scholar 

20.
Gedalof, Z. E. & Berg, A. A. Tree ring evidence for limited direct CO2 fertilization of forests over the 20th century. Glob. Biogeochem. Cycles 24, GB3027 (2010).
Google Scholar 

21.
Gazol, A., Camarero, J. J., Anderegg, W. R. L. & Vicente-Serrano, S. M. Impacts of droughts on the growth resilience of Northern Hemisphere forests. Glob. Ecol. Biogeogr. 26, 166–176 (2016).
Google Scholar 

22.
Pomara, L. Y. & Zuckerberg, B. Climate variability drives population cycling and synchrony. Divers. Distrib. 23, 421–434 (2017).
Google Scholar 

23.
Briffa, M. et al. Trends in recent temperature and radial tree growth spanning 2000 years across northwest Eurasia. Phil. Trans. R. Soc. B 363, 2271–2284 (2008).
PubMed  Google Scholar 

24.
Ponocná, T. et al. Deviations of treeline Norway spruce radial growth from summer temperatures in East-Central Europe. Agric. Meteorol. 253, 62–70 (2018).
Google Scholar 

25.
Shestakova, T. A., Gutiérrez, E., Valeriano, C., Lapshina, E. & Voltas, J. Recent loss of sensitivity to summer temperature constrains tree growth synchrony among boreal Eurasian forests. Agric. Meteorol. 268, 318–330 (2019).
Google Scholar 

26.
Schurman, J. S. et al. Large-scale disturbance legacies and the climate sensitivity of primary Picea abies forests. Glob. Change Biol. 24, 2169–2181 (2018).
Google Scholar 

27.
Schweingruber, F. H. Tree Rings: Basics and Applications of Dendrochronology (Springer Science & Business Media, 1996).

28.
Speer, J. H. Fundamentals of Tree-Ring Research (Univ. of Arizona Press, 2010).

29.
Savolainen, O., Pyhäjärvi, T. & Knürr, T. Gene flow and local adaptation in trees. Annu. Rev. Ecol. Evol. Syst. 38, 595–619 (2007).
Google Scholar 

30.
Ripa, J. Analysing the Moran effect and dispersal: their significance and interaction in synchronous population dynamics. Oikos 89, 175–187 (2000).
Google Scholar 

31.
Hopson, J. & Fox, J. W. Occasional long distance dispersal increases spatial synchrony of population cycles. J. Anim. Ecol. 88, 154–163 (2019).
PubMed  Google Scholar 

32.
Johnson, C. A. et al. Effects of temperature and resource variation on insect population dynamics: the bordered plant bug as a case study. Funct. Ecol. 30, 1122–1131 (2016).
PubMed  Google Scholar 

33.
Parmesan, C. & Yohe, G. A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37–42 (2003).
CAS  PubMed  Google Scholar 

34.
St. George, S. The aberrant global synchrony of present-day warming. Nature 571, 483–484 (2019).
CAS  PubMed  Google Scholar 

35.
Neukom, R. et al. Consistent multidecadal variability in global temperature reconstructions and simulations over the common era. Nat. Geosci. 12, 643–649 (2019).
PubMed  PubMed Central  Google Scholar 

36.
Babst, F. et al. Twentieth century redistribution in climatic drivers of global tree growth. Sci. Adv. 5, eaat4313 (2019).
PubMed  PubMed Central  Google Scholar 

37.
Giguère-Croteau, C. et al. North America’s oldest boreal trees are more efficient water users due to increased [CO2], but do not grow faster. Proc. Natl Acad. Sci. USA 116, 2749–2754 (2019).
PubMed  Google Scholar 

38.
Duncan, A. B., Gonzalez, A. & Kaltz, O. Dispersal, environmental forcing, and parasites combine to affect metapopulation synchrony and stability. Ecology 96, 284–290 (2015).
PubMed  Google Scholar 

39.
Girardin, M. P. et al. No growth stimulation of Canada’s boreal forest under half-century of combined warming and CO2 fertilization. Proc. Natl Acad. Sci. USA 113, E8406–E8414 (2016).
CAS  PubMed  Google Scholar 

40.
Manzanedo, R. D. et al. Increase in CO2 concentration could alter the response of Hedera helix to climate change. Ecol. Evol. 8, 8598–8606 (2018).
PubMed  PubMed Central  Google Scholar 

41.
Venter, O. et al. Sixteen years of change in the global terrestrial human footprint and implications for biodiversity conservation. Nat. Commun. 7, 12558 (2016).
CAS  PubMed  PubMed Central  Google Scholar 

42.
Pederson, N. et al. Long-term drought sensitivity of trees in second-growth forests in a humid region. Can. J. Res. 42, 1837–1850 (2012).
Google Scholar 

43.
Kug, J. S., An, S. I., Ham, Y. G. & Kang, I. S. Changes in El Niño and La Niña teleconnections over North Pacific–America in the global warming simulations. Theor. Appl. Clim. 100, 275–282 (2010).
Google Scholar 

44.
Rahmstorf, S. & Coumou, D. Increase of extreme events in a warming world. Proc. Natl Acad. Sci. USA 108, 17905–17909 (2011).
CAS  PubMed  Google Scholar 

45.
Cai, W. et al. Increasing frequency of extreme El Niño events due to greenhouse warming. Nat. Clim. Change 4, 111–116 (2014).
CAS  Google Scholar 

46.
Fischer, E. M., Beyerle, U. & Knutti, R. Robust spatially aggregated projections of climate extremes. Nat. Clim. Change 3, 1033–1038 (2013).
Google Scholar 

47.
Trenberth, K. E., Fasullo, J. T. & Shepherd, T. G. Attribution of climate extreme events. Nat. Clim. Change 5, 725–730 (2015).
Google Scholar 

48.
Thurm, E. A. et al. Alternative tree species under climate warming in managed European forests. Ecol. Manage. 430, 485–497 (2018).
Google Scholar 

49.
Manzanedo, R. D. & Pederson, N. Towards a more ecological dendroecology. Tree Ring Res. 75, 152–159 (2019).
Google Scholar 

50.
Klesse, S. et al. Sampling bias overestimates climate change impacts on forest growth in the southwestern United States. Nat. Commun. 9, 5336 (2018).
CAS  PubMed  PubMed Central  Google Scholar 

51.
Ettinger, A. K., Kevin, R. F. & HilleRisLambers, J. Climate determines upper, but not lower, altitudinal range limits of Pacific Northwest conifers. Ecology 92, 1323–1331 (2011).
CAS  PubMed  Google Scholar 

52.
Grissino-Mayer, H. D. & Fritts, H. C. The International Tree-Ring Data Bank: an enhanced global database serving the global scientific community. Holocene 7, 235–238 (1997).
Google Scholar 

53.
Wilson, R. et al. Last millennium Northern Hemisphere summer temperatures from tree rings. Part I: the long term context. Quat. Sci. Rev. 134, 1–18 (2016).
Google Scholar 

54.
Cook, B. I., Anchukaitis, K. J., Touchan, R., Meko, D. M. & Cook, E. R. Spatiotemporal drought variability in the Mediterranean over the last 900 years. J. Geophys. Res. Atmos. 121, 2060–2074 (2016).
PubMed  PubMed Central  Google Scholar 

55.
Charney, N. D. et al. Observed forest sensitivity to climate implies large changes in 21st century North American forest growth. Ecol. Lett. 19, 1119–1128 (2016).
PubMed  Google Scholar 

56.
Bunn, A. G. A dendrochronology program library in R (dplR). Dendrochronologia 26, 115–124 (2008).
Google Scholar 

57.
R: A Language and Environment for Statistical Computing Version 3.5.0 (R Core Team, 2017).

58.
Sheffield, J., Goteti, G. & Wood, E. F. Development of a 50-yr high-resolution global dataset of meteorological forcings for land surface modelling. J. Clim. 19, 3088–3111 (2006).
Google Scholar 

59.
Lamarque, J.-F. et al. The Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP): overview and description of models, simulations and climate diagnostics. Geosci. Model Dev. 6, 179–206 (2013).
CAS  Google Scholar 

60.
GISTEMP Team GISS Surface Temperature Analysis (GISTEMP) Version 4 (NASA Goddard Institute for Space Studies, accessed 2 July 2018); https://data.giss.nasa.gov/gistemp/ More

1.
Manzanedo, R. D., HilleRisLambers, J., Rademacher, T. T. & Pederson, N. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-020-01306-x (2020).
2.
Shestakova, T. A. et al. Proc. Natl Acad. Sci. USA 113, 662–667 (2016).
CAS  Article  Google Scholar 

3.
Pugh, T. A. M. et al. Proc. Natl Acad. Sci. USA 116, 4382–4387 (2019).
CAS  Article  Google Scholar 

4.
Di Cecco, G. J. & Gouhier, T. C. Sci. Rep. 8, 14850 (2018).
Article  Google Scholar 

5.
Shestakova, T. A., Gutiérrez, E., Valeriano, C., Lapshina, E. & Voltas, J. Agric. For. Meteorol. 268, 318–330 (2019).
Article  Google Scholar 

6.
Hubau, W. et al. Nature 579, 80–87 (2020).
CAS  Article  Google Scholar 

7.
Adams, M. A., Buckley, T. N. & Turnbull, T. L. Nat. Clim. Change 10, 466–471 (2020).
CAS  Article  Google Scholar 

8.
Brienen R. J. W., Schöngart J. & Zuidema, P. A. in Tropical Tree Physiology Tree Physiology Vol. 6 (eds Goldstein, G. & Santiago L. S.) 439–461 (Springer, 2016).

9.
Hansen, B. B., Grøtan, V., Herfindal, I. & Lee, A. M. Ecography https://doi.org/10.1111/ecog.04962 (2020).

10.
Koenig, W. D. & Liebhold, A. M. Nat. Clim. Change 6, 614–617 (2016).
Article  Google Scholar  More

Experiments with robotic fish
We developed, and employed, a bio-mimetic robotic fish platform (Fig. 1, Supplementary Figs. 1–4 and Movie 1 and 2) in order to experimentally evaluate the costs and benefits of swimming together. We constructed two identical robotic fish, 45 cm in length and 800 g in mass. Each has three sequential servo-motors controlling corresponding joints, covered in a soft, waterproof, rubber skin. In addition, the stiffness of the rubber caudal fin decreases towards the tip33 (Supplementary Fig. 1). The motion of the servomotors is controlled using a bio-inspired controller called a central pattern generator (CPG)34,35 resulting in the kinematics that mimic normal real fish body undulations when swimming36 (see Supplementary Fig. 2 and Note 1). Here, due to the complexity of the problem (as discussed above) we consider hydrodynamic interactions between pairs of fish. We note that this is biologically meaningful as swimming in pairs is both the most common configuration found in natural fish populations7,10,37,38, and it has been found that even in schools fish tend to swim close to only a single neighbour7,37.
Fig. 1: The robotic fish platform employed to investigate hydrodynamic benefits of schooling.

a The reverse Karman vortices shedding by the robotic fish with dye flow visualisation. b Schematic view of the setup that allows setting various spatiotemporal differences between two robotic fish swimming in a flow tank (Front-back distance D ∈ [0.22, 1] BL (body length), Left-right distance G ∈ [0.27, 0.33] BL and Phase difference Φ ∈ [0, 2π]). A laser generator was used to visualise the hydrodynamic interactions (see Supplementary Note 1). c The phase difference Φ is evaluated by the difference between the undulation phase of the two robots. Undulation phase is evaluated based on the lateral position Lt of the tail tip. d Power cost (absolute value on the left y-axis, and relative value compared to the average power cost on the right y-axis), is shown as a function of the phase difference at D = 0.33 BL and G = 0.27 BL. Error bars are standard error of the mean.

Full size image

To evaluate the energetics of swimming together we conducted experiments on our pair of robotic fish in a flow tank (test area: 0.4-m-wide, 1-m-long and 0.45-m-deep; Fig. 1b and Supplementary Fig. 3). In order to conduct such an assessment we first measured the speed of our robots when freely swimming alone (we did so in a large tank 2-m-wide, 3-m-long and 0.4-m-deep). We then set the flow speed within our flow tank to this free-swimming speed (0.245 ms−1) allowing us to ensure the conditions in the flow tank are similar to those of the free-swimming robot. Unlike in the solitary free-swimming condition, to have precise control of spatial relationships in the flow tank we suspended each robotic fish by attaching a thin aluminium vertical bar to the back of each robot, which was then attached to a step motor above the flow tank (Supplementary Fig. 3 and Movie 2). To establish whether the robotic fish connected with a thin bar has similar hydrodynamics compared to when free swimming, we measured the net force (of the drag and thrust generated by the fish body in the front-back direction) acting on the robot in the flow tank. The measured net force over a full cycle (body undulation) was found to be zero; thus the bar is not measurably impacting the hydrodynamics of our robot fish in the front-back direction as they swim in the flow tank (Supplementary Fig. 5).
To further validate the utility of the platform, we also compared the power consumption of our robots swimming side-by-side, for different relative phase differences Φ, with equivalent measurements made with a simple 2D computational fluid dynamics (CFD) model of the same scenario (Supplementary Note 2). In both cases (see Supplementary Fig. 6a, c for robotic experiments and CFD simulations, respectively) we find that there exists an approximately sinusoidal relationship between power costs and phase difference which is defined as Φ = ϕleader − ϕfollower (Fig. 1c, d). Due to the 2D nature of the simulation, as well as many other inevitable differences between simulations and real world mechanics, the absolute power costs are different from those measured for the robots, but nevertheless the results from these two approaches are broadly comparable and produce qualitatively similar relative power distributions when varying the phase difference between the leader and follower. These results indicate that our robotic fish are both an efficient (making estimates of swimming costs is far quicker with our robotic platform than it is with CFD simulations) and effective (in that they capture the essential hydrodynamic interactions as well as naturally incorporate 3D factors) platform for generating testable hypotheses regarding hydrodynamic interactions in pairs of fish.
We subsequently utilise our robots to directly measure the energy costs associated with swimming together as a function of relative position (front-back distance D from 0.22 to 1 body length (BL) in increments of 0.022 BL and left-right distance G from 0.27 to 0.33 BL in increments of 0.022 BL) while also varying the phase relationships (phase difference, Φ) of the body undulations exhibited by the robots (the phase of the follower’s tailbeat ϕfollower relative to that of the leader’s ϕleader, Fig. 1c).
By conducting 10,080 trials (~120 h of data), we obtain a detailed mapping of the power costs relative to swimming alone associated with these factors (Fig. 2a). Such a mapping allows us to predict how real fish, that continuously change relative positions6,8, should correspondingly continuously adjust their phase relationship in order to maintain hydrodynamic benefits. To quantify the costs we determine the energy required to undulate the tail of each robot allowing us to define, and calculate, a dimensionless relative power coefficient as:

$$eta =frac{({P}_{1}^{{rm{Water}}}-{P}^{{rm{Air}}})-({P}_{2}^{{rm{Water}}}-{P}^{{rm{Air}}})}{{P}_{1}^{{rm{Water}}}-{P}^{{rm{Air}}}}=frac{{P}_{1}^{{rm{Water}}}-{P}_{2}^{{rm{Water}}}}{{P}_{1}^{{rm{Water}}}-{P}^{{rm{Air}}}},$$
(1)

where η is the relative power coefficient, PAir, ({P}_{1}^{{rm{Water}}}) and ({P}_{2}^{{rm{Water}}}) are the power costs of the robotic fish swimming in the air (an approximation of the dissipated power cost due to mechanical friction, resistance, etc. within the robot that are not related to interacting with the water), alone in water, or in a paired context in the water, respectively. ({P}_{1}^{{rm{Water}}}-{P}^{{rm{Air}}}) and ({P}_{2}^{{rm{Water}}}-{P}^{{rm{Air}}}) therefore represent the power costs due to hydrodynamics while swimming alone, and in a pair, respectively (see Methods section). Correspondingly, the coefficient η compares the energy cost of fish swimming in pairs to swimming alone. Positive values (blue in Fig. 2a) and negative values (red in Fig. 2a) respectively represent energy saving and energy cost relative to swimming alone. The difference between the maximum energy saving and maximum energy cost for the robots is ~13.4%.
Fig. 2: Robotic fish save energy by vortex phase matching (VPM).

a Relative power coefficient η shown as a function of the phase difference between the leader and the follower Φ and front-back distance D at left–right distance G = 0.31 BL. The dashed line (also in b) shows the functional relationship described in Eq. (2) that determines the theoretical phase relationship that maximally saves energy (Methods section). ({Phi }_{0}^{* }) is the optimal initial phase difference (fitted to the data points of maximum energy saving, as shown in b). The points marked by red square, blue circle and blue square indicate example cases depicted on panels c–e. b Location of maximal energy saving in the robotic trials. Point size and darkness denote the number of occurrences of each phase difference value at each front-back distance. c–e An illustration of important spatial configurations for vortex phase matching. Energy cost is related to how the follower moves its body relative to the direction of the induced flow of the vortices, in the opposite direction with Φ0 = ({Phi }_{0}^{* })+π (c) or in the same direction with Φ0 = ({Phi }_{0}^{* }) (d, e). Followers interact with the induced flow of vortices with the same body phase at any front-back distance (within the range of hydrodynamic interactions), termed vortex phase matching. (d, e; Φ0 = ({Phi }_{0}^{* }) describes the hydrodynamic interaction resulting in energy saving, see description in the text). As the front-back distance changes, the followers must dynamically adopt phase difference Φ, with respect to that of the leader.

Full size image

Our results indicate that there exists a relatively simple linear relationship between front-back distance and relative phase difference of the follower that minimises the power cost of swimming (as indicated by the dashed lines in Fig. 2a, b, the theoretical basis of which we will discuss below). This suggests that a follower can minimise energetic expenditure (and avoid substantial possible energetic costs) by continuously adopting a unique phase difference Φ that varies linearly as a function of front-back distance D (see Fig. 2b for example), even as that distance changes. We find that while left-right distance G does alter energy expenditure, this effect is minimal when compared to front-back positioning, and has little effect on the above relationship (Supplementary Figs. 7 and 8) in the range explored here.
Although we know fish generate reverse Kármán vortices at the Reynolds number (Re = Lu/ν ≈ 105, where L is the fish body length, u is the swimming or flow speed and ν is the kinematic viscosity) in our experiments39 (Supplementary Fig. 9 and Movie 1), turbulence will dominate over longer distances18. In accordance with this, we see a relatively fast decay in the benefits of swimming together as a function of D (e.g., D  > 0.7 BL, Supplementary Fig. 10), a feature we also expect to be apparent in natural fish schools (where it would likely be exacerbated by what would almost always be less-laminar flow conditions). Therefore, we expect, based on our results, that hydrodynamic interactions are dominated by short-distance vortex-body interactions (with D   2 BL), and it thus cannot benefit from neighbour-generated vortices. We also chose this method since isolating the fish would likely induce stress responses that could confound our results. To evaluate body kinematics in the presence of vortices we analysed the body undulations of the follower when in close proximity (within 0.4 BL), where hydrodynamic effects will be strongest (Supplementary Fig. 25). We find that in the vicinity of vortices, fish exhibit a higher tailbeat amplitude and lower tailbeat frequency (Supplementary Fig. 26), which indicates less power consumption48.
To further test if fish can save energy by adopting VPM with the typical vortex-body hydrodynamic interactions (Φ0 = −0.2π), we compared an estimation of the power consumption under different hydrodynamic interactions. Since the hydrodynamic interactions are mainly determined by the initial phase difference Φ0 (see above), we analysed performance in the full possible range from −π to π (see Supplementary Fig. 25 for the detailed method). We define relative energy saving when fish exhibit higher tailbeat amplitudes A (Fig. 4a) and lower tailbeat frequencies f (Fig. 4b) than average48, and find that the range is Φ0∈ [−0.5π, 0.5π] (the shaded area in Fig. 4). Figure 4 also shows that while fish adopting Φ0 ≈ 0 will save the most energy, those exhibiting Φ0 = −0.2π, as in our experiments, will save almost the same amount (thus they are very close to optimal in this respect).
Fig. 4: Relative energetic benefits to a follower in real fish pairs.

a, b Energy cost analysis was conducted by calculating the difference in amplitude A (a) and frequency f (b) at Φ0 and the same measurements with the opposite phase Φ0 + π (written as A+π and f+π respectively) as a function of initial phase difference Φ0 (Supplementary Fig. 25 and Note 4). Data are pooled from all pairs when the follower’s front-back positions are not >0.4 BL distance (where the hydrodynamic interactions are expected to be the strongest). The hatched areas show the energy saving zone of Φ0. The dashed line denotes ({Phi }_{0}^{* }), the most typically observed initial phase difference exhibited by our fish. (Average amplitude is 0.09 BL, average frequency is 2.3 Hz).

Full size image

Fish in our experiments (Fig. 3c at D = 0 BL) spent 59% of their time swimming with phase relationships (Φ0∈ [−0.5π, 0.5π]) that save energy, and the remaining 41% that imposes some (relative) energetic costs. However, because the energy cost has a sinusoidal relationship to the phase difference (Fig. 1d) simply calculating the percentage of time in each regime (in which there is either a benefit or a cost, regardless of the magnitude of each) is insufficient. By combining the frequencies (occurrences) of each phase difference Φ observed in Fig. 3c and the sinusoidal shape of the power cost as a function of Φ (Fig. 1d and Supplementary Fig. 6b, d), we can estimate that by behaving as they do, fish (in our flow conditions) save (by accumulating all benefits and extra costs; where a random behaviour would give 0) an overall 15% of the total possible (which would be achieved by perfectly adopting the optimal phase to the neighbour-generated vortices at all the time, Supplementary Fig. 27). It is possible that if fish are exposed to more challenging, stronger flow regimes (here we employed those of typical swimming), that this percentage will increase. However we would never expect fish motion to be completely dominated by a need to save energy as they must also move in ways as to obtain salient social and asocial information from their visual, olfactory, acoustic and hydrodynamic environment, such as to better detect food52, environmental gradients44 and threats16. Nevertheless, kinematic analysis suggests that they adopt VPM in a way that results in energy savings (dashed line in Fig. 4).
In summary, our bio-mimetic robots provided an effective platform with which we could explore the energetic consequences of swimming together in pairs and revealed that followers could benefit from neighbour-generated flows if they adjust their relative tailbeat phase difference linearly as a function of front-back distance, a strategy we term vortex phase matching. A model based on fundamental hydrodynamic principles, informed by our flow visualisations, was able to account for our results. Together, this suggests that the observed energetic benefit occurs when a follower’s tail movement coincides with the induced flow generated by the leader. Finally, experiments with real fish demonstrated that followers indeed employ vortex phase matching and kinematic analysis of their body undulations suggests that they do so, at least in part, to save energy. By providing evidence that fish do exploit hydrodynamic interactions, we gain an understanding of important costs and benefits (and thus the selection pressures) that impact social behaviour. In addition, our findings provide a simple, and robust, strategy that can enhance the collective swimming efficiency of fish-like underwater vehicles. More