Ecology

Sampling
Killer whale biopsy samples of skin and blubber from 38 individuals were collected year-round from August 2017 to July 2018 in northern Norway. All whales were sampled according to relevant guidelines and regulations, and conducted under the permit FOTS-ID 10176 issued by Mattilsynet (the Norwegian Food Safety Authority, report nr. 2016/179856). Details of seasonal sampling locations, stable isotope dietary descriptors and classification of sampled individuals are described in a previous study14. In the current study, total Hg was analysed in skin from all individuals (n = 38), whereas organohalogen contaminants (OHC) was analysed in blubber of 31 individuals due to insufficient blubber for the remaining 7 individuals.
OHC analysis
OHC analysis was conducted at the Laboratory of Environmental Toxicology at the Norwegian University of Life Sciences, Oslo, Norway. We analysed a total of 83 OHCs: 49 organochlorines (OCs), including 34 PCBs and 15 organochlorine pesticides (OCPs), 18 brominated flame retardants (BFRs), including newer and unregulated compounds, and 16 hydroxylated metabolites (OH-metabolites) of PCBs and polybrominated diphenylethers (PBDEs). A full list of analysed compounds can be found in Supplementary Table S1.
We analysed OCs and BFRs using a multicomponent method, first described in 197842, and since modified for a range of compounds and biological matrices43,44,45,46. The analysis of the OH-metabolites was conducted according to previously published methods47,48. An outline of the method is described in the Supplementary Information. Reported concentrations were blank corrected based on the average concentration detected within blank samples. The limit of detection (LOD) was defined as three times the average noise in chromatograms, and ranged from 0.40 to 11.10 ng/g w.w. for OCs, 0.012 to 0.362 ng/g w.w. for BFRs and 0.013 to 0.040 ng/g w.w. for OH-metabolites (see Supplementary Table S2). Internal reference materials for OCs and BFRs (contaminated seal blubber, MTref01) and OH-metabolites (contaminated seal blood, MTref03) were also extracted in conjunction with sample material to assess method performance. Internal standard recoveries are listed in Supplementary Table S2.
Hg analysis
We analysed total Hg by atomic absorption spectrometry at the University of Oslo, using a Direct Mercury Analyser (DMA-80, Milestone Srl, Soirsole, Italy). Killer whale skin samples were freeze dried in a Leybold-Heraeus GT2 freeze dryer with a Leybold Vakuum GmbH vacuum pump (Leybold, Cologne, Germany) and then homogenised to a fine powder using an agate pestle and mortar. Approximately 0.002 g of killer whale skin were analysed in parallel with sample blanks and certified reference material (DORM-4, fish protein; DOLT-5, dogfish liver, National Research Council, Ottawa, Canada). If enough material, samples were analysed in duplicates to ensure precision of measurements and the arithmetic mean value used. Average recoveries of the certified reference materials were within 10% of the reported values. The detection limit of the instrument was 0.05 ng mercury.
Data treatment
We included OHC compounds found in levels above the instrument’s LOD in a minimum of 65% of the individual whale samples for statistical analysis (see Supplementary Table S1, Supporting Information for pollutants excluded). For individual concentrations below the LOD, we imputed left-censored data by replacing missing values with a random number between 0 and the LOD assuming a beta distribution (α = 5, β = 1) to retain the pattern of the dataset. In total, 95 values below the LOD were replaced, representing 6.52% of the OHC dataset. All total Hg samples were above the LOD.
We defined the ΣPCBs as the sum of all 28 PCB congeners detected in more than 65% of the whale samples (PCB-28, -66, -74, -87, -99, -101, 105, -110, -114, -118, -128, -137, -138, -141, -149, -151, -153, -156-, 157, 170, -180, -183, -187, -189, 194, -196, -206, -209). The definition for ΣPCBs varies within killer whale literature, with some studies analysing only a few core PCB congeners35, some all 209 of the possible congeners36, and others not providing a definition (e.g. for thresholds for possible health effects7). There will therefore inevitably be some errors in comparisons. However, since the ΣPCBs in killer whales is dominated by a few commonly reported congeners, typically PCB-153 and -13816,37, it is unlikely that the inclusion of other minor constituents will have a major influence on the total load. PCBs were further grouped according to the number of chlorine substitutions per molecule, i.e. homologue group to compare the pattern of PCBs. ΣDDTs was defined as the sum of p,p′-DDT, p,p′-DDD and p,p′-DDE, the ΣPBDEs as the sum of BDE-28, -47, -99, -100, -153 and -154 and the sum of chlordanes (ΣCHLs) as the sum of oxychlordane, trans-chlordane, cis-chlordane, trans-nonachlor and cis-nonachlor.
Statistical analyses
Statistical analyses were performed using R v. 3.4.149. The significance level was to set to α = 0.05, except in cases where the value was adjusted due to multiple testing, and was two-tailed. In addition to visual inspection, normality was tested using the Shapiro–Wilk’s test50 and homogeneity of variance by Levene’s test51 using the R package car52.
Whale dietary groups
The dietary groups used in this study are based on a previous study, which used stable isotope values inputted into a Gaussian mixture model to assign sampled individuals to two fish-eating groups: Herring-eaters and Lumpfish-eaters and one mammal-eating group Seal-eaters14. The three dietary groups were characterised by disparate, non-overlapping isotopic niches that were consistent with predatory field observations. The seal-eating group was defined by higher δ15N values than the two fish-eating groups.
We found that the herring and lumpfish-eating killer whales did not differ in either their OHC levels (Tukey’s HSD: p = 0.49) or total Hg levels (pairwise Welch’s t-test: p = 0.67). In this study, we thus combined the dietary groups Herring-eaters and Lumpfish-eaters into the group Fish-eaters, to enable easier comparison to the seal-eating killer whales.
We then used Welch’s t-test to compare the ΣPCB levels in the seal-eating and fish-eating dietary groups (using a log10 transformation), and to compare the total Hg levels in the skin between the two dietary groups.
OHC dataset
We used multivariate analysis to compare and visualise the differences in all the OHCs between the dietary groups, age and sex classes using the vegan package in R53. Principle Component Analysis (PCA) was used to visualise the main structure of the data: reducing the dimensions to two new, uncorrelated, latent variables termed principle components 1 and 2 (PC1 and PC2). We log-10 transformed contaminant levels to ensure normality and homogeneity of variance, and the presence of any influential outliers were checked by the Cook’s distance test. Redundancy Analysis (RDA) was used to extract and summarise the variation in the OHC levels constrained, and thereby explained, by a set of explanatory variables54. Significant associations between response variables and the explanatory variables were identified by an RDA based forward model selection, followed by a Monte Carlo forward permutation test (1,000 unrestricted permutations). The samples’ scores along PC1 were subject to one-way Analysis of Variance (ANOVA) followed by Tukey’s honestly significant difference post hoc test (Tukey’s HSD) to analyse differences between the three dietary groups. PC1 scores were also used to evaluate correlation to total Hg levels in the skin using a Spearman’s rank correlation test. Absolute concentrations were subject to PCA with lipid % as a covariate, after checking its significance using RDA, as lipid normalising data in inferential statistics can often lead to misleading conclusions55.
We lipid-normalised OHC values when comparing levels to threshold values for toxicity or other killer whale populations, and used the geometric mean as the average for each dietary group to reflect the log normal distribution of the data. In accordance with convention, efforts were made to only compare adult males with other worldwide populations, as reproductive female whales are known to transfer a substantial portion of their OHC burden to their calves35,36,38. In any case of comparison, similar metrics were compared (i.e. arithmetic mean, geometric mean, median) and all variables kept similar (i.e. sex, age, biopsy/stranded animals). We make the assumption in this study that the killer whales sampled in 2002 in Norway were fish-eaters for the following reasons: firstly, the whales were sampled on herring overwintering grounds, feeding on herring, and photographs were taken of five of the eight adults sampled and were identified as herring-eaters from previous field observations16. Secondly, the PCB pattern in the sampled whales showed 76% of ΣPCBs higher chlorinated congeners (hexaCBs or higher), which is more similar to the fish-eaters from our study (80% higher chlorinated congeners) than the seal-eaters (87% higher chlorinated congeners). Thirdly, the upper 95% confidence range of all pollutants reported in the 2002 killer whales falls below both the geometric and arithmetic mean values for seal-eaters from this study.
Total Hg dataset
The normal distribution of the data within each dietary group meant we used the arithmetic mean as an average. The three dietary groups (Herring-eaters, Lumpfish-eaters and Seal-eaters) were compared using a pairwise Welch’s t-test with a Benjamini–Hochberg False Discovery Rate correction to adjust for multiple testing. Because we found no difference between the Herring-eaters and Lumpfish-eaters (p = 0.67), we combined these two groups to a group called “Fish-eaters” for easier comparison with the seal-eaters. The total Hg levels in the skin of the two groups, Fish-eaters and Seal-eaters were compared using Welch’s t-test.
There is a strong positive correlation between Hg levels in the skin and liver in toothed whales, and this can be used to compare Hg levels measured in skin with hepatic toxicity threshold values56,57,58. To extrapolate to liver from skin in our samples, we chose an equation based on a model using concentrations in the liver (Hgliver μg/g w.w) and skin (Hgskin μg/g w.w) of bottlenose dolphins (Tursiops truncatus) (Eq. 1)58. We converted dry weight to wet weight using the water content for each individual whale measured during freeze drying.

$$ln left( {Hg_{liver} } right) = 1.6124 times ln left( {Hg_{skin} } right) + 2.0346$$
(1)

When comparing Hg concentrations to other worldwide populations, both male and female whales were included. This was due to a lack of information of sex in one of the populations for comparisons and because killer whales are unlikely to pass on Hg burdens to calves5,59. More

1.
Barbier, E. B. et al. The value of estuarine and coastal ecosystem services. Ecol. Monogr. 81, 169–193 (2011).
Article  Google Scholar 
2.
Douglas, E. J. et al. Macrofaunal functional diversity provides resilience to nutrient enrichment in coastal sediments. Ecosystems 20, 1324–1336 (2017).
CAS  Article  Google Scholar 

3.
Edgar, G. J. & Barrett, N. S. Effects of catchment activities on macrofaunal assemblages in Tasmanian estuaries. Estuar. Coast. Shelf Sci. 50, 639–654 (2000).
ADS  Article  Google Scholar 

4.
Murray, N. J. et al. The global distribution and trajectory of tidal flats. Nature 565, 222–225 (2019).
ADS  CAS  Article  Google Scholar 

5.
Hewitt, J. E., Thrush, S. F., Dayton, P. K. & Bonsdorff, E. The effect of spatial and temporal heterogeneity on the design and analysis of empirical studies of scale-dependent systems. Am. Nat. 169, 398–408 (2007).
Article  Google Scholar 

6.
Cadenasso, M. L., Pickett, S. T. A., Weathers, K. C. & Jones, C. G. A framework for a theory of ecological boundaries. Bioscience 53, 750 (2003).
Article  Google Scholar 

7.
Lohrer, A. M. et al. Biogenic habitat transitions influence facilitation in a marine soft-sediment ecosystem. Ecology 94, 136–145 (2013).
Article  Google Scholar 

8.
Schenone, S., O’Meara, T. A. & Thrush, S. F. Non-linear effects of macrofauna functional trait interactions on biogeochemical fluxes in marine sediments change with environmental stress. Mar. Ecol. Prog. Ser. 624, 13–21 (2019).
ADS  CAS  Article  Google Scholar 

9.
Thrush, S. F., Pridmore, R. D., Hewitt, J. E. & Cummings, V. J. Adult infauna as facilitators of colonization on intertidal sandflats. J. Exp. Mar. Biol. Ecol. 159, 253–265 (1992).
Article  Google Scholar 

10.
Borcard, D., Legendre, P. & Drapeau, P. Partialling out the spatial component of ecological variation. Ecology 73, 1045–1055 (1992).
Article  Google Scholar 

11.
Mermillod-Blondin, F., Rosenberg, R., François-Carcaillet, F., Norling, K. & Mauclaire, L. Influence of bioturbation by three benthic infaunal species on microbial communities and biogeochemical processes in marine sediment. Aquat. Microb. Ecol. 36, 271–284 (2004).
Article  Google Scholar 

12.
Dornhoffer, T., Waldbusser, G. & Meile, C. Modeling lugworm irrigation behavior effects on sediment nitrogen cycling. Mar. Ecol. Prog. Ser. 534, 121–134 (2015).
ADS  CAS  Article  Google Scholar 

13.
Braeckman, U. et al. Role of macrofauna functional traits and density in biogeochemical fluxes and bioturbation. Mar. Ecol. Prog. Ser. 399, 173–186 (2010).
ADS  CAS  Article  Google Scholar 

14.
Banta, G. T., Holmer, M., Jensen, M. H. & Kristensen, E. Effects of two polychaete worms, Nereis diversicolor and Arenicola marina, on aerobic and anaerobic decomposition in a sandy marine sediment. Aquat. Microb. Ecol. 19, 189–204 (1999).
Article  Google Scholar 

15.
Woodin, S. A. et al. Same pattern, different mechanism: Locking onto the role of key species in seafloor ecosystem process. Sci. Rep. https://doi.org/10.1038/srep26678 (2016).
Article  PubMed  PubMed Central  Google Scholar 

16.
Woodin, S. A., Wethey, D. S., Hewitt, J. E. & Thrush, S. F. Small scale terrestrial clay deposits on intertidal sandflats: Behavioral changes and productivity reduction. J. Exp. Mar. Biol. Ecol. 413, 184–191 (2012).
Article  Google Scholar 

17.
Thrush, S. F., Hewitt, J. E. & Pridmore, R. D. Patterns in the spatial arrangements of polychaetes and bivalves in intertidal sandflats. Mar. Biol. 102, 529–535 (1989).
Article  Google Scholar 

18.
Pridmore, R. D., Thrush, S. F., Hewitt, J. E. & Roper, D. S. Macrobenthic community composition of six intertidal sandflats in Manukau Harbour, New Zealand Macrobenthic community composition of six intertidal sandflats in Manukau Harbour, New Zealand. N. Z. J. Mar. Freshw. Res. 24, 81–96 (1990).
Article  Google Scholar 

19.
Turner, S. J. et al. Are soft-sediment communities stable? An example from a windy harbour. Mar. Ecol. Prog. Ser. 120, 219–230 (1995).
ADS  Article  Google Scholar 

20.
Zajac, R. N. et al. Responses of infaunal populations to benthoscape structure and the potential importance of transition zones. Limnol. Oceanogr. 48, 829–842 (2003).
ADS  Article  Google Scholar 

21.
Kobayashi, G., Goto, R., Takano, T. & Kojima, S. Molecular phylogeny of Maldanidae (Annelida): Multiple losses of tube-capping plates and evolutionary shifts in habitat depth. Mol. Phylogenet. Evol. 127, 332–344 (2018).
Article  Google Scholar 

22.
Volkenborn, N. et al. Intermittent bioirrigation and oxygen dynamics in permeable sediments: An experimental and modeling study of three tellinid bivalves. J. Mar. Res. 70, 794–823 (2012).
CAS  Article  Google Scholar 

23.
Waldbusser, G. G., Marinelli, R. L., Whitlatch, R. B. & Visscher, P. T. The effects of infaunal biodiversity on biogeochemistry of coastal marine sediments. Limnol. Ocean. 49, 1482–1492 (2004).
CAS  Article  Google Scholar 

24.
Walker, B. H. Biodiversity and ecological redundancy. Conserv. Biol. 6, 18–23 (1992).
Article  Google Scholar 

25.
Volkenborn, N., Polerecky, L., Wethey, D. S. & Woodin, S. A. Oscillatory porewater bioadvection in marine sediments induced by hydraulic activities of Arenicola marina. Limnol. Oceanogr. 55, 1231–1247 (2010).
ADS  CAS  Article  Google Scholar 

26.
Thrush, S. F. et al. Changes in the location of biodiversity–ecosystem function hot spots across the seafloor landscape with increasing sediment nutrient loading. Proc. R. Soc. B. Biol. Sci. https://doi.org/10.1098/rspb.2016.2861 (2017).
Article  Google Scholar 

27.
O’Meara, T., Gibbs, E. & Thrush, S. F. Rapid organic matter assay of organic matter degradation across depth gradients within marine sediments. Methods Ecol. Evol. 9, 245–253 (2018).
Article  Google Scholar 

28.
Kana, T. M. et al. Membrane inlet mass spectrometer for rapid high-precision determination of N2, O2, and Ar in environmental water samples. Anal. Chem. 66, 4166–4170 (1994).
CAS  Article  Google Scholar 

29.
Thrush, S. F. et al. Changes in the location of biodiversity–ecosystem function hot spots across the seafloor landscape with increasing sediment nutrient loading. Proc. R. Soc. B Biol. Sci. https://doi.org/10.1098/rspb.2016.2861 (2017).
Article  Google Scholar 

30.
Legendre, P. & Legendre, L. F. J. Numerical Ecology (Elsevier, New York, 2012).
Google Scholar 

31.
Grömping, U. Relative importance for linear regression in R: The package relaimpo. J. Stat. Softw. 20, 17 (2006).
Google Scholar 

32.
Team RC. R: A language and environment for statistical computing. 2013. More

The backward and forward processes
Let us consider a branched tree, starting with (N_0) cells at time (t=0) and ending with N(t) cells at time t as shown on Fig. 1. We assume that all lineages survive up to time t, and therefore the final number N(t) of cells corresponds to the number of lineages in the tree.
The most natural way to sample the lineages is to put uniform weights on all of them. This sampling is called backward, (or retrospective) because at the end of the experiment one randomly chooses one lineage among the N(t) with a uniform probability and then one traces the history of the lineage backward in time from time t to 0, until reaching the ancestor population. The backward weight associated with a lineage l is defined as

$$begin{aligned} omega _{text {back}}(l)=N(t)^{-1} ,. end{aligned}$$
(1)

In a tree, some lineages divide more often than others, which results in an over-representation of lineages that have divided more often than the average. Therefore by choosing a lineage with uniform distribution, we are more likely to choose a lineage with more divisions than the average number of divisions in the tree.
The other way of sampling a tree is the forward (or chronological) one and consists in putting the weight

$$begin{aligned} omega _text {for}(l)= N_0^{-1} m^{-K(l)} ,, end{aligned}$$
(2)

on a lineage l with K(l) divisions, where m is the number of offspring at division. This choice of weights is called forward because one starts at time 0 by uniformly choosing one cell among the (N_0) initial cells, and one goes forward in time up to time t, by choosing one of the m offspring with equal weight 1/m at each division. The backward and forward weights are properly normalized probabilities, defined on the N(t) lineages in the tree at time t: (sum _{i=1}^{N(t)} omega _{text {back}}(l_i) = sum _{i=1}^{N(t)} omega _{text {for}}(l_i) =1).
Figure 1

Example of a tree with (N_0=1) and (N(t)=10) lineages at time t. Two lineages are highlighted, the first in blue with 2 divisions and the second in orange with 5 divisions. The forward sampling is represented with the green right arrows: it starts at time (t=0) and goes forward in time by choosing one of the two daughters lineages at each division with probability 1/2. The backward sampling is pictured by the left purple arrows: starting from time t with uniform weight on the 10 lineages it goes backward in time down to time (t=0).

Full size image

Single lineage experiments are precisely described by a forward process since experimentally, at each division, only one of the two daughter cells is conserved while the other is eliminated (for instance flushed away in a microfluidic channel9, 10). In these experiments, a tree is generated but at each division only one of the two lineages is conserved, with probability 1/2, while the rest of the tree is eliminated. This means that single lineage observables can be measured without single lineage experiments, provided population experiments are analyzed with the correct weights on lineages.
Link with the population growth rate
Since the backward weight put on a lineage depends on the number of cells at time t, it takes into account the reproductive performance of the colony but it is unaffected by the reproductive performance of the lineage considered. On the contrary, the forward weight put on a specific lineage depends on the number of divisions of that lineage but is unaffected by the reproductive performance of other lineages in the tree. Therefore, the difference between the values of the two weights for a particular lineage informs on the difference between the reproductive performance of the lineage with respect to the colony.
We now introduce the population growth rate:

$$begin{aligned} Lambda _t=frac{1}{t} ln frac{N(t)}{N_0} ,, end{aligned}$$
(3)

which is linked to forward weights by the relation

$$begin{aligned} frac{N(t)}{N_0}=sum _{i=1}^{N(t)} m^{K_i} omega _text {for}(l_i) = langle m^K rangle _text {for} ,, end{aligned}$$
(4)

where (langle cdot rangle _text {for}) is the average over the lineages weighted by (omega _text {for}), and (K_i=K(l_i)). Combining the two equations above, we obtain19:

$$begin{aligned} Lambda _t=frac{1}{t} ln langle m^K rangle _text {for} ,, end{aligned}$$
(5)

which allows an experimental estimation of the population growth rate from the knowledge of the forward statistics only.
Equation (4) can also be re-written to express the bias between the forward and backward weights of the same lineage

$$begin{aligned} frac{omega _{text {back}}(l)}{omega _text {for}(l)}=frac{m^{K(l)}}{langle m^K rangle _text {for}} ,, end{aligned}$$
(6)

which is the reproductive performance of the lineage divided by its average in the colony with respect to (omega _text {for}).
A similar relation is derived using the relation

$$begin{aligned} frac{N_0}{N(t)}=sum _{i=1}^{N(t)} m^{-K_i} omega _text {back}(l_i) = langle m^{-K} rangle _text {back} ,. end{aligned}$$
(7)

Combining Eqs. (5) and (7) we obtain:

$$begin{aligned} Lambda _t= – frac{1}{t} ln langle m^{-K} rangle _text {back} ,. end{aligned}$$
(8)

A similar equation as Eq. (6) can be obtained in terms of the backward sampling and reads: 

$$begin{aligned} frac{omega _{text {back}}(l)}{omega _text {for}(l)}=frac{langle m^{-K} rangle _text {back}}{m^{-K(l)}} ,. end{aligned}$$
(9)

Combining Eqs. (1) to (3), we obtain the fluctuation relation13,17:

$$begin{aligned} omega _{text {back}}(l)= omega _text {for}(l) e^{K(l) ln m – t Lambda _t} ,. end{aligned}$$
(10)

If we now introduce the probability distribution of the number of divisions for the forward sampling (p_text {for}(K)=sum _l delta (K-K(l)) omega _text {for}(l)) and similarly for the backward sampling, we can also recast the above relation as a fluctuation relation for the distribution of the number of divisions:

$$begin{aligned} p_{text {back}} (K,t)=p_{text {for}} (K,t) e^{K ln m – t Lambda _t} ,. end{aligned}$$
(11)

Let us now introduce the Kullback–Leibler divergence between two probability distributions p and q, which is the non-negative number:

$$begin{aligned} {{mathscr {D}}}_{text {KL}}(p||q)=int {mathrm {d}}x , p(x) ln frac{p(x)}{q(x)} ge 0 ,. end{aligned}$$
(12)

Using Eq. (10), we obtain

$$begin{aligned} {{mathscr {D}}}_{text {KL}}(omega _{text {back}}|| omega _text {for}) = langle K rangle _{text {back}} ln m – t Lambda _t ge 0 ,. end{aligned}$$
(13)

A similar inequality follows by considering ({{mathscr {D}}}_{text {KL}}(omega _{text {for}}|| omega _text {back})). Finally we obtain

$$begin{aligned} frac{t}{langle K rangle _{text {back}}} le frac{ln m}{Lambda _t} le frac{t}{langle K rangle _text {for}} ,. end{aligned}$$
(14)

In the long time limit, (lim nolimits _{t rightarrow + infty } t/langle K rangle _{text {back}} = langle tau rangle _{text {back}}), where (tau) is the inter-division time, or generation time, defined as the time between two consecutive divisions on a lineage. The same argument goes for the forward average. In the case of cell division where each cell only gives birth to two daughter cells ((m=2)), the center term in the inequality tends to the population doubling time (T_d). Therefore, this inequality reads in the long time limit:

$$begin{aligned} langle tau rangle _{text {back}} le T_d le langle tau rangle _text {for} ,. end{aligned}$$
(15)

Let us now mention a minor but subtle point related to this long time limit. For a lineage with K divisions up to time t, we can write (t=a + sum _{i=1}^{K} tau _i), where a is the age of the cell at time t and where (tau _i) is the generation time associated with the ith division. Then (t/ K= tau _m + a/K), where (tau _m) is the mean generation time along the lineage. For finite times, all we can deduce is (t/ K ge tau _m). Therefore the left inequality of Eq. (15) always holds

$$begin{aligned} langle tau rangle _{text {back}} le frac{t}{langle K rangle _{text {back}}} le frac{ln m}{Lambda _t} ,, end{aligned}$$
(16)

while the right inequality does not necessarily hold at finite time.
Inspired by work by Powell6, the inequalities of Eq. (15) have been theoretically derived in12 for age models. In our previous work17, we have replotted the experimental data of12 which confirm theses inequalities and we have shown theoretically that the same inequalities should also hold for size models. In fact, as the present derivation shows, the relation equation (14) is very general and only depends on the branching structure of the tree, while the relation equation (15) requires in addition the existence of a steady state. These inequalities and Eq. (11) express fundamental constraints between division and growth, which should hold for any model.
Stochastic thermodynamic interpretation
The results derived above have a form similar to that found in Stochastic Thermodynamics18. According to this framework, Eq. (5) is an integral fluctuation relation (similar to Jarzynski relation) while Eq. (11) is a detailed fluctuation relation (similar to Crooks fluctuation relation). Furthermore, the inequalities equation  (14) represent a constraint equivalent to the second law of thermodynamics, which classically follows from the Jarzynski or Crooks fluctuation relations. It is known that these inequalities take a slightly different form when expressed at finite time or at steady state, which is indeed the case here when comparing Eq. (14) with Eq. (15). A difference between work fluctuation relations like Crooks or Jarzynski and equations (5) and (11), is that Crooks or Jarzynski describe non-autonomous systems which are driven out of equilibrium by the application of a time-dependent protocol, whereas the relations for cell growth derived here concern autonomous systems, in the absence of any external protocol.
One of the main applications of Jarzynski or Crooks fluctuation relations concerns the thermodynamic inference of free energies from non-equilibrium fluctuations. Similarly, Eq. (5) or Eq. (11) can be used as estimators of the population growth rate. The specific advantage of Eq. (5) with respect to Eq. (11) is that it only requires single lineage statistics, which can be obtained from mother machine experiments. Let us now show how this can be done in practice. We use the data from20, where the growth of many independent lineages of E. coli have been recorded over 70 generations in a mother machine at three different temperatures (25 °C, 27 °C, and 37 °C), precisely 65 lineages for 25 °C, 54 for 27 °C, and 160 for 37 °C. For each temperature condition, we study the convergence of the estimator of the population growth rate based on Eq. (5), which we call (Lambda _{mathrm{lin}}) as a function of the length t of the lineages for a fixed number of independent lineages L, and as a function of the number of independent lineages for a fixed observation time.
Figure 2

Estimator of the population growth rate (Lambda _{mathrm{lin}}) based on Eq. (5), (a) as function of the the length t of the lineages and (b) as function of the number L of lineages used in the estimation. In (a), the curves for the three temperatures converge to a constant value. In (b), only the curve for 37 °C is shown and the horizontal dashed line represents the quantity (ln (2)/langle tau rangle _{text {for}}), which is smaller than the limit value of (Lambda _{mathrm{lin}}), as expected from the second law-like inequality, namely Eq. (15). In the inset, the purple histogram is the distribution of the number of divisions, while the green filled histogram is the histogram deduced from it by weighting it by a factor (2^K) and normalizing. All the 160 lineages were used to plot these histograms.

Full size image

Firstly, for each temperature, we take into account all the lineages available and truncate them at an arbitrary time t smaller than the length of the shortest lineage of the set. On these portions of lineages of length t, we compute (Lambda _{mathrm{lin}}) versus the time t as shown in Fig. 2a. We see that the estimator (Lambda _{mathrm{lin}}) starts from zero, increases and eventually converges rather quickly towards a limiting value. The limit we found agree with the independent analysis carried out in19, with only one caveat, these authors reported that their estimator started at high values and then decreased towards the limit, while in our case, the estimator starts at zero and later increases towards the limit. In our case, the estimator needs to be zero at short times, before the first divisions occur.
Secondly, we truncate all the lineages at a fixed time equal to the length of the shortest lineage of the set, and compute (Lambda _{mathrm{lin}}) versus the number L of lineages considered for the estimation, which have been randomly selected from the ensemble of available lineages. As shown in Fig. 2b for the case at (37^{,circ } hbox {C}) (curves for the other temperatures look exactly the same), the convergence is also excellent in that case. Although the value of the population growth rate which is obtained in this way can not be measured independently from the evolution of the population in the mother machine setup, this convergence is indicative of the success of the method. The figure also confirms that the value of the population growth rate deduced from the estimator (Lambda _{mathrm{lin}}) is larger than (ln (2)/langle tau rangle _{text {for}}), as predicted by the right inequality of Eq. (15).
Here, the estimator is found to provide an excellent estimation, but this is not always so. For instance, for the inference of free energies from non-equilibrium work measurements, the exponential average of the estimator is often dominated by rare values, which are not accessible or not well sampled21. To understand why this problem does not arise here, we show in inset of Fig. 2b, the distribution P(K) of the number of divisions together with the same distribution weighted by the factor (2^K) and normalized. The peak of that modified distribution informs on the dominant values in the estimator21. Here, we observe that both distributions have a narrow support and are close to each other. The weighted distribution is peaked at (K=67) while P(K) is peaked at (K=66), therefore typical and dominating values are very close, which explains why the estimator is good.
Let us now further develop the Stochastic Thermodynamic interpretation of our results by analyzing the implications of the previous fluctuation relations when dynamical variables are introduced on the branched tree of the population. Let us introduce M variables labeled ((y_1,y_2, ldots ,y_M)) to describe a dynamical state of the system, then a path is fully determined by the values of these variables at division, and the times of each division. We call ({mathbf {y}}(t)=(y_1(t),y_2(t), ldots ,y_M(t))) a vector state at time t and ({{mathbf {y}}}={{mathbf {y}}(t_j)}_{j=1}^{K}) a path with K divisions. For cell growth models, the variables (y_i) can typically be the size and age of the cell, or the concentration of a key protein.
The probability ({{mathscr {P}}}) of path ({{mathbf {y}}}) is defined as the sum over all lineages of the weights of the lineages that follow the path ({{mathbf {y}}}):

$$begin{aligned} {{mathscr {P}}}({{mathbf {y}}},K,t)=sum _{i=1}^{N(t)} omega (l_i) , delta (K-K_i) delta ({{mathbf {y}}} – {{mathbf {y}}}_i) ,, end{aligned}$$
(17)

where ({{mathbf {y}}}_i) is the path followed by lineage (l_i). Using the normalization of the weights (omega) on the lineages, we show that ({{mathscr {P}}}) is properly normalized: (int mathrm {d}{{mathbf {y}}} sum _K {{mathscr {P}}}({{mathbf {y}}},K,t) = 1). We then define the number (n({{mathbf {y}}},K,t)) of lineages in the tree at time t that follow the path ({{mathbf {y}}}) with K divisions:

$$begin{aligned} n({{mathbf {y}}},K,t)=sum _{i=1}^{N(t)} delta (K-K_i) delta ({{mathbf {y}}} – {{mathbf {y}}}_i) ,. end{aligned}$$
(18)

This number of lineages is normalized as (int mathrm {d}{{mathbf {y}}} sum _K n({{mathbf {y}}},K,t) = N(t)). Then, the path probability can be re-written as

$$begin{aligned} {{mathscr {P}}}({{mathbf {y}}},K,t) = n({{mathbf {y}}},K,t) cdot omega (l) ,. end{aligned}$$
(19)

Since (n({{mathbf {y}}},K,t)) is independent of a particular choice of lineage weighting, we obtain

$$begin{aligned} frac{{{mathscr {P}}}_{text {back}}({{mathbf {y}}},K,t)}{{{mathscr {P}}}_text {for} ({{mathbf {y}}},K,t)}=frac{omega _{text {back}}(l)}{omega _text {for}(l)}= e^{K ln m – t Lambda _t} , , end{aligned}$$
(20)

which generalizes Eq. (11). In our previous work17, we have derived this relation for size models with individual growth rate fluctuations (i.e. ({mathbf {y}}=(x,nu ))) but we were not aware of the weighting method introduced by13, and for this reason, we used the term ‘tree’ to denote the backward sampling, and the term ‘lineage’ to denote the forward sampling.
This relation has a familiar form in Stochastic Thermodynamics. The central quantity called entropy production can indeed be expressed similarly as the relative entropy between probability distributions associated with a forward and a backward evolution. In this analogy, ({{mathbf {y}}}) is analog to the trajectory and (t Lambda _t – K ln m) is analog to the entropy production. Then, the equivalent of a reversible trajectory for which the entropy production is null is a lineage for which the number K of divisions is equal to (t Lambda _t / ln m), that is, a lineage having the same reproductive performance as that of the colony. When all the lineages in a tree have this property, there is no variability of the number of divisions among them. In that case, the forward and backward distributions are identical, and the cost function (t Lambda _t – K ln m) vanishes for all lineages. More

Ethics statement
This study did not involve endangered or protected species and was carried out in accordance with the guidelines of the French Polynesia code for animal ethics and scientific research (https://www.service-public.pf/diren/partager/code/). All protocols and experiments were approved by the CRIOBE-IRCP animal ethics committee (DL-20150214).
Model species
The convict surgeonfish Acanthurus triostegus (Linnaeus, 1758) is an abundant coral-reef-associated species found throughout the Indo-Pacific60, including around Moorea Island, French Polynesia11. It has a pelagic larval duration of ~53 days60 after which larvae move back to the reef61. A. triostegus is a well-studied species with regards to metamorphosis, with a well-defined developmental sequence11. The recruitment of A. triostegus larvae to reef habitats coincides with a true metamorphosis into juveniles, with this full process taking around 1 week11. Its metamorphosis is controlled by thyroid hormones (TH, the precursor thyroxine (T4), and the active hormone triiodothyronine (T3)), which is the same as other teleosts and other metamorphosing vertebrates such as amphibians11. TH signaling in larval A. triostegus is vulnerable to disruption by anthropogenic stressors, including the waterborne pesticide chlorpyrifos and artificial light at night11,17. Given the importance of TH during teleost metamorphosis10, and as metamorphosis coincides with recruitment in a taxonomically diverse range of coral-reef fishes (i.e., Acanthuridae, Apogonidae, Balistidae, Chaetodontidae, and Pomacentridae)11, we consider A. triostegus a representative teleost model for examining the effects of anthropogenic stressors on fish recruitment via their impacts on metamorphic processes.
Study period and site
This study was conducted from February 2015 to June 2018, at Moorea Island, French Polynesia (S17°32′16.4589″, W149°49′48.3018″). Sampling for the examination of sensory development under pharmacological treatments was conducted in 2015. Sampling for the investigation of behavioral preferences and survival to predation under pharmacological treatments, as well as sensory development under anthropogenic stressor exposure, were conducted in 2016. Sampling for the examination of survival to predation under anthropogenic stressors exposure, and TH levels under co-exposure to anthropogenic stressors were conducted in 2018. Sampling for the investigation of TH levels under single anthropogenic stressors was conducted in both 2016 and 2018.
Fish sampling
Settlement-stage A. triostegus (i.e., fully transparent individuals11, here define as day 0 (d0) individuals), were collected on the north–east coast of the island (S17°29′49.7362″, W149°45′13.899″) at night using a crest net11, as they transitioned from the ocean to the reef. Fish were then transferred to either in situ cages (for thyroid hormone signaling experiments) where they remained until d8, or to aquaria at the CRIOBE Marine Research Station (for increased temperature and chlorpyrifos exposure experiments) where they remained until d5.
Thyroid hormone signaling experiment: control, T3- and N3 treatments
To test the role of TH on sensory development, behavior (response to sensory cues), and survival (predation test) in metamorphosing A. triostegus, their TH pathway was pharmacologically manipulated. Fish were injected, at d0 immediately following capture, in their ventral cavity with 20 µl of a pharmacological treatment: (i) T3 + iopanoic acid (IOP) both at 10−6 M (T3 treatment), or (ii) NH3 at 10−6 M (N3 treatment). IOP was used as an inhibitor of deiodinase enzymes, following comparable work in mammals and amphibians62, and as routinely used in fish to prevent the immediate degradation of injected T348. The T3 treatment was therefore applied to promote TH signaling. NH3 is a known antagonist of TH receptors (TR) in vertebrates63 and in A. triostegus in particular11. NH3 prevents the binding of TH such as T3 to TR, therefore impairing the binding of transcriptional coactivators to TR, which therefore remain in an inactive and repressive conformation63,64. The N3 treatment was thus applied to repress TH signaling by disrupting the TH pathway. T3 and NH3 were initially suspended in dimethyl sulfoxide (DMSO), at 10−2 M, and then diluted in phosphate buffered saline (PBS) 1× to reach 10−6 M. Control fish were injected with 20 µL of DMSO diluted 10.000 times in PBS 1× to control for the effect of the solvent and injection. DMSO and NH3 non-toxicity has been previously determined11. Until subsequent sampling, all fish were re-treated each morning (i.e., at d2, d3, and d4 for fish sampled at d5) to maintain pharmacological activity11.
Thyroid hormone signaling experiment: fish husbandry
Following collection and subsequent treatment, larvae were transferred to a nursery area on the north coast of the island (S17°29′26.5378″, W149°53′29.2252″) where they were raised in in situ cages (cylindrical cages, diameter: 30 cm, height: 50 cm, 15 fish per cage). This allowed them to develop in in situ conditions11. As A. triostegus feeds on algal turf following settlement11,60, cages were stocked with a supply of turf-covered coral rubble that was replaced daily, ensuring both shelter and constant food availability. Fish were subsequently sampled on d2, d5, and d8, post collection to examine sensory development, behavior (response to sensory cues), and survival (predation test). Twelve in situ cages were used, and the use of any given cage was randomized prior each experiment.
Increased temperature and chlorpyrifos (CPF) exposures: fish husbandry and treatments
Following collection, d0 fish were maintained in groups of ten individuals at the CRIOBE Marine Research Station. Each group was held in a 30 L × 20 W × 20 H-cm aquaria containing 12 L of filtered (1-µm filter) seawater. All tanks were subject to a 12:12 h light–dark cycle (06:00–18:00 light period) and oxygenated with an air stone. Twelve aquaria were used (six for exposures to increased temperatures only, and six for exposures with CPF), and the use of any given aquarium was randomized prior each experiment.
For increased temperature treatments, seawater was in an open system, and water temperature was maintained at either 28.5 °C, 30.0 °C or 31.5 °C. 28.5 °C was chosen as the basal temperature as this was the mean temperature in the Moorea lagoon at the time of the study, and corresponds to the mean annual lagoon temperature in this region (http://observatoire.criobe.pf). Subsequent increases of +1.5 °C and +3.0 °C were selected, as these are in line with end-of-century projections for tropical Pacific sea surface temperatures37. Fresh coral rubble was added to the tanks and replaced each day to provide both food and shelter. Heaters controlled by thermostats were used to maintain the temperature treatments. Before the experiment, each tank was in open circuit with temperature maintained at 28.5 °C. At the beginning of the experiment, fish were introduced in the aquarium, and the thermostat temperature was then set up to 28.5 °C, 30.0 °C or 31.5 °C according to the treatment. The temperature of interest was reached within 2 h. Temperature in each tank was then visually checked (on the thermostat controller) at least five times per day, and never differed from the target temperature by more than 0.2 °C.
For CPF exposure, five different treatments were applied: unaltered seawater (control), seawater with acetone at a final concentration of 1:1.000.000 (CPF0, solvent control treatment, as CPF was made soluble using acetone), or seawater with CPF at a nominal concentration of either 1, 5, or 30 μg L−1 (CPF1, CPF5, and CPF30 treatments), based on the findings of recent studies of reef fishes exposed to CPF11,14,65. CPF was spiked in each tank from dilutions that were prepared in advance: 1 µg µL−1, 5 µg µL−1, and 30 µg µL−1. From these dilutions, 12 µL were pipetted and spiked in the 12-L exposure tanks, therefore reaching nominal concentrations of 1 µg L−1, 5 µg L−1, and 30 µg L−1. Similarly, 12 µL of acetone was spiked in the tank for the CPF0 condition. Spike was allowed to mix for 2 min (water mixing due to the air stone) before fish were introduced in the tank. At the end of the 32-h exposure, CPF concentrations in the water or in the fish tissues were not evaluated, as we were only interested in the effects of CPF spikes on fish metamorphic processes. Nevertheless, a previous study using similar methods and nominal concentrations of similar magnitude (i.e., ranging from 4 to 64 µg L−1) measured CPF levels corresponding to 80% of nominal concentrations after 24 h66, therefore suggesting a good stability of CPF levels in the condition of our study. Environmentally, these nominal concentrations represent high (CPF1) to extremely high (CPF30) exposures, as recorded contamination concentration in Australian and North American surface waters are generally below 1 µg L−1, with a few high outliers reaching up to 26 µg L−153. However, this shows that concentrations of up to 30 µg L−1 are possible on a short timescale (such as the 32-h exposure of this study), in particular in coastal and shallow areas such as fish nurseries. Aquaria used for CPF exposure treatments and associated controls were not equipped with coral rubble to prevent potential interaction with the pesticide11. This may explain the delay in trunk canal development observed in control fish from the CPF exposure treatments compared with control fish from the increased temperature treatments and control fish from the TH treatments (Figs. 1g and  3e). Water was replaced each day to ensure the maintenance of CPF concentrations11.
For combined increased temperature and CPF exposure treatments, fish were also maintained without coral rubble with water replaced daily.
For treatments where fish were exposed to an anthropogenic stressor and provided with supplemental T3, fish were maintained in either the elevated temperature or CPF exposure treatments as described above, but were also injected with either T3 (T3 treatment) or with DMSO (control, to control for solvent and injection). These injections were done as described above (thyroid hormone signaling experiment: control, T3- and N3-treatments).
Sample preparation and fish measurements
For thyroid hormone quantifications and histological analyses, fish were first euthanized in freshly prepared MS222 at 0.4 mg ml−1 in filtered seawater at 4 °C, and instantly placed on ice. Following euthanasia, all fish except those for histological analyses were weighed (W) then placed on a scale bar and photographed for standardized measurements (e.g., height (H) and standard length (SL)). These measurements were used to calculate Fulton’s condition factor K = 100*(W/SL3) (with SL, in cm, and W, in g) for each individual67. Weight measurements were performed using a precision balance (Ohaus Adventurer Precision), and length measurements were taken using the imageJ software.
Retinas
The retina is the light-sensitive tissue within the eye, and is composed of different cell layers organized ventrally and dorsally around the optic nerve (Fig. 1c). To analyze the retina’s cell layers, we removed the head region of euthanized fish in a way that allowed us to distinguish the right retina from the left retina, and the dorsal and ventral sides of the retina. This region was then fixed in Bouin solution before being embedded in paraffin for microtome sectioning (5 µm). Sections were then stained using hematoxylin and eosin. We chose cross-sections of the retina that were done at similar depth by selecting sections at the level of the optical nerve (Supplementary Fig. 5a). These sections enabled us to identify different cell segments, types, and layers, among which (i) photoreceptor external segments (perceiving light signals), (ii) photoreceptor nuclei, (iii) bipolar cells (which integrate the synaptic signals originating from the photoreceptors), and (iv) ganglion cells (which integrate signals from bipolar cells and create action potential toward the optic nerve) were easily distinguishable (Supplementary Fig. 5b). We only investigated the right eye of A. triostegus fish, as evidence suggests this species is visually lateralized at recruitment, and predominantly uses its right eye to examine visual predator stimuli14. Also, we only examined the dorsal side of the retina, as the ventral side (vs) was shown to not undergo change at metamorphosis in another coral-reef fish species, the goatfish Upeneus tragula13 (Supplementary Fig. 5a). To compare the developmental state of the retina between treatments, we looked at the peripheral area of the dorsal side (see the dotted square in Supplementary Fig. 5a, magnified in Supplementary Fig. 5b) as settlement-stage individuals in another acanthurid species (Naso brevirostris) showed weak if any spatial specialization of the retina, in particular on this axis68. We prioritized the examination of bipolar cell densities as important changes in bipolar cell density was observed during U. tragula metamorphosis13. However, we also examined the densities of photoreceptor external segments, photoreceptor nuclei, and ganglion cells along metamorphosis and across TH treatments, revealing for photoreceptor external segments (Supplementary Fig. 6) and nuclei (Supplementary Fig. 7) a TH-dependent density increase at metamorphosis, while ganglion cells showed no density variation during metamorphosis (Supplementary Fig. 8). Cell counting was performed in a 50-µm wide area perpendicular to the retina cell layers (Supplementary Fig. 5b). The measure of bipolar cell density corresponds to the number of bipolar cells per 0.001 mm² in the inner-nuclear cell layer (INL), after measuring the mean thickness of this INL (average of three measurements at random non-overlapping locations within the 50-µm wide area) (Supplementary Fig. 5b). The prn and ggc densities were obtained in a similar manner after measuring the mean width of their respective cell layer. Pes density corresponds to the exact number of pes in the 50-µm wide area as the pes cell layer is monocellular (Supplementary Fig. 5b). Cells that were located on the edges of the 50-µm wide area were included in the counts.
Olfactory organ lamellae
In fish, the olfactory organ is a fluid-filled blind sac that contains the ciliated sensory epithelium, which forms a rosette comprise of a number of folds, i.e., lamellae44. In A. triostegus, the left and right olfactory organs can be found in two cavities on the dorsal surface, between the eye and the snout edge, with water moving through each cavity from the anterior nostril to the posterior nostril (Supplementary Fig. 3). Pictures (e.g., Supplementary Fig. 3) were obtained using SEM microscopy following specialized tissue preparation involving the dissection of the olfactory organ (removal of the thin skin layer between the two nostrils) and fixation in a sodium cacodylate + glutaraldehyde solution (2.5% glutaraldehyde in 1 M sucrose and 0.1 M sodium cacodylate, pH 7.4) for a week at 4 °C. Samples were then washed (10% sucrose in 0.1 M sodium cacodylate solution, pH 7.4) 15 × 15 min at room temperature (RT). Samples were then dehydrated using successive ethanol (EtOH) baths (30%, 50%, 70%, 85%, 95%, 100%, 100%, 100%, each time during 15 min at RT). Samples were then placed in EtOH:HMDS (1:1) 15 min at RT, in 100% HMDS for 15 min at RT, and lastly in 100% HMDS for 30 min at RT. HMDS was used to replace critical point drying. Samples were not metalized. SEM pictures were obtained using a MiniMEB Hitachi TM3030 (University of French Polynesia). Given that lamellae are well formed and have a macroscopic size in A. triostegus at recruitment, lamellae counting was performed on tissues samples identically prepared, but observed under a simple binocular stereomicroscope. Lamellae counting was performed on both left and right olfactory organs to then have a measure of the average number of lamellae per olfactory organ. In the case where one of the left or right olfactory organ was damaged and thus accurate lamellae counts were not possible, we used the count from the other olfactory organ as the average number, as no significant difference was observed between the number of lamellae in the left and right olfactory organs in our preliminary observations.
Lateral line trunk canal pores
The lateral line system enables fish to detect water motions and pressure gradients, such as those caused by other fish (e.g., movement from other fish in the shoal or predator strikes). It is composed of superficial and canal neuromast receptor organs, which are the functional units of the lateral line system and are ciliary sensory organs, composed of hair cells, like those in the inner ear, located either on the skin or embedded in lateral line canals43. Canal neuromasts are found in the epithelium lining the bottom of the lateral line canals, and one canal neuromast is usually found within the short canal segments between two adjacent canal pores43. The development of neuromasts and the morphogenesis of lateral line canals and their pores initiate in late-stage larvae and continue through metamorphosis43. Counting the number of pores on the trunk canal is thus an appropriate way to rapidly characterize the development of the lateral line system when one cannot perform more advanced histological analyzes of the neuromasts. In d0 to d8 A. triostegus at recruitment, a fully formed trunk canal corresponding to a complete arched canal can be observed on each of the fish body flanks (Supplementary Fig. 4a). At this stage, A. triostegus does not exhibit scales69, and the trunk canal is therefore only composed of soft tissue (Supplementary Fig. 4a). Following the same preparation protocol as used for olfactory organs, we investigated the development of the lateral line system of A. triostegus by counting the number of pores on the trunk canal (Supplementary Fig. 4a-b). Counting was performed either on the left or the right side of the body, depending on tissue preservation. Indeed, as fish body frequently curved when dried to a critical point in HMDS, only one side of the body was generally accessible for counting trunk canal pores. We did not consider this an issue, as preliminary observations revealed no significant difference between the numbers of trunk canal pores present on the left or right flanks of the body. Variation in the number of pores cannot be attributed to variation in the length of the fully formed trunk canal (measured from the vicinity of the head to the base of the tail) as it does not change in A. triostegus during metamorphosis (Supplementary Fig. 9).
Thyroid hormone quantification
TH was extracted from frozen fish, following an extraction protocol adapted from previous studies70,71,72, including on coral-reef fishes11. Fish were individually crushed in 500 µl of methanol using a Minilys and glass beads (3*30 s, 5.000 rpm), centrifuged at 4 °C (10 min, 12.000 rpm), and supernatant reserved. These operations were performed twice, then fish were crushed one last time with 400 µl of methanol, 100 µl of chloroform, and 100 µl of barbital buffer (3*30 s, 5.000 rpm, RT), centrifuged at 4 °C (10 min, 12.000 rpm), and supernatant reserved. Pooled supernatants were then dried at 70 °C for 2 h. Dried pellets were re-suspended in 400 µl of methanol, 100 µl of chloroform, and 100 µl of barbital buffer, centrifuged at 4 °C (20 min, 12.000 rpm), and supernatant reserved. The same operation was again performed on the pellets. Pooled supernatants were once again dried out at 70 °C for 2 h with the final extract reconstituted in 2 ml of PBS 1× for quantification (Roche Elica kit on a Cobas analyzer, following the manufacturer’s standardized method). TH levels in pg g−1 of fish were then transformed into relative levels by selecting the respective control fish as standards in the increased temperature and CPF exposure treatment experiments.
Behavioral tests
A two-channel chemical choice flume73,74,75 was used to assess the responses of A. triostegus towards chemical predator cues (Supplementary Fig. 10). Each trial presented an individual fish with two water sources: control seawater (Ø = UV-sterilized and 1-µm filtered seawater from the collection site) vs “seawater containing chemical cues from predator” (P = UV-sterilized and 1-µm filtered seawater from the collection site in which five Lutjanus fulvus predators were soaked, into a 125-L tank, for 2 h prior to the experiment). Water was fed into the flume through two water inlets, with equivalent flow rates of 100 ml min−1 maintained using flow meters (MM Minimaster, Admi-France). This rate ensured laminar flow of each water source was maintained in parallel while allowing fish to swim naturally between the two water sources. Fine mesh and collimators also helped to ensure laminar flow of each water source was maintained75. Preliminary experiments were conducted to ensure the absence of unanticipated biased behavior within the two-channel flume (e.g., preference for one side of the flume over the other, irrespectively of the water sources, or preference for the drain area over the choice area). Prior to each trial, dye tests were conducted to confirm laminar flow, without eddies or areas of water mixing, within the choice area. After releasing the fish in the middle of the choice area, a 2-min acclimation period was observed, then fish position (left or right part of the choice area, or drain area) was recorded every 2 s for 5 min (Supplementary Fig. 10), using a camera (GoPro Hero 2) located above the edge of the flume tank. Water inlets were then switched (to account for any side preference due to fish’s immobility) followed by another 2-min acclimation period. A second 5-min test period followed, during which fish position was again assessed every 2 s. Preference or avoidance of water sources, and the absence of a side preference, were confirmed by comparing responses during the first and second 5-min test periods. These experiments were performed in the dark with red light, to limit potential visual perturbations such as the presence of an observer and to allow comparison between d2 and d0 fish, as d0 fish were tested immediately after collection (i.e., at night) as this is when they are actively moving from the ocean to the reef, and is thus the most biologically relevant time to do so. Fish that did not swim actively during the first acclimation period were removed from the analysis (n = 3) to prevent side preference bias. However, all remaining fish swam actively between the two water sources after water inlets were switched (either during the second acclimation period or during the second test period), ensuring no continuous side bias due to immobility. Fish that spent more than 50% of the test time in the drain area (i.e., where the two water sources mix) were also removed from the analysis (n = 4) as we considered that they did not show any particular preference or avoidance for any of the two water sources. Fish that did not make a clear choice between the two water sources but spent more than 50% of the time in the choice area were included in the analysis. This was done as we wanted to assess fish preference as well as the absence of preference.
A double-choice tank was used to assess the responses of A. triostegus in the presence of a visual predator cues. Each trial presented an individual fish with a choice of two visual stimuli, each contained in a separate aquarium placed at the end of the central rectangular choice tank. These were an aquarium, containing a 10-cm (standard length) L. fulvus (P condition) vs an empty aquarium, equipped with an air stone (Ø condition) (Supplementary Fig. 11). Fish were first placed into the central “no choice” area of the choice tank for a 2-min acclimation period. During this time, fish were not able to see the contents of either adjacent aquaria or access the “choice” areas of the choice tank as opaque panels were positioned at the edges of the “no choice” area (see the dotted lines in Supplementary Fig. 11). Following the acclimation period, the opaque panels were removed, allowing fish to observe the visual stimuli in the adjacent aquaria through the transparent walls of the choice tank and to access the choice areas. Fish position (i.e., choice area 1, no choice area, choice area 2; Supplementary Fig. 11) was then assessed every 2 s over a 10 min, using a camera (GoPro Hero 2) to limit any external visual disturbances such as an observer’s presence. The camera was located above the choice tank. Location of visual stimulus (left or right side of the double-choice tank) was switched between each fish to ensure the absence of a side preference. Fish that remained immobile and spent >50% of the test time in the “no choice” area were removed from the analysis (n = 11), as we considered that they did not show any particular preference or avoidance for any of the two visual stimuli.
Survival arena
The survival tests presented in Figs. 2c, 3f were conducted in an arena set up in situ in a nursery lagoon area of the north coast (S17°29’7.0272″, W149°49’51.1166″). The arena consisted of a 1-m3 cage with a hard bottom covered with sand and coral rubble and four lateral walls made from 5-mm fine mesh (Supplementary Fig. 12). Each trial consisted of 45 fishes, with 15 from each treatment group. Fish from each treatment group were tagged with a specific color using visible implant fluorescent filament (Northwest Marine Technology) 2 h prior being released into the arena. For each trial, all fishes were released simultaneously and allowed to acclimate for 30 min before the introduction of six L. fulvus (15–20 cm SL). After 2 h, predators were removed, and surviving A. triostegus was identified to treatment using their color tag. Color tags attributed to each fish group were randomly switched between each replicate to ensure no predation bias based on tag color. Differences in survival during this predation test cannot be attributed to variation in fish size or condition between treatments, as these did not differ significantly between groups (Supplementary Fig. 13).
Statistical analyses
All statistical analyses were conducted using R version 3.5.376. Conway–Maxwell–Poisson (COM-Poisson) generalized linear models (GLM) were used to assess if pharmacological and anthropogenic stressor treatments influenced the number of lamellae, the number of trunk canal pores, and the number of survivors to the predation experiment77. COM-Poisson GLM were also used to assess if pharmacological treatments influenced the number of photoreceptor external segments. Linear models (LM) were used to assess if pharmacological treatments influenced the bipolar cell, ganglion cell, and photoreceptor nuclei densities, and if pharmacological and anthropogenic stressor treatments influenced fish Fulton’s K condition factor. LM were also used to assess if the trunk canal length varied with age. Gamma generalized linear mixed-effect models (GLMEM) were used to assess if anthropogenic stressor exposures influenced TH levels and T3/T4 ratios78. TH level or T3/T4 ratios were used as the dependent variable, and replicate was included as a random factor to account for differences in TH levels only due to the two different Cobas analyzers that were used in the two different years. As preliminary experiments provided no evidence that season and lunar phase affected T3 levels in metamorphosing A. triostegus, we did not include them in our analyses (Gamma GLMEM, Supplementary Fig. 14). For each model, diagnostic plots were examined and outputs compared with raw data to confirm goodness-of-fit and residual homoscedasticity, and, when applicable, residual normality was assessed using Shapiro–Wilk normality test. Paired t tests or Wilcoxon signed-rank tests were used to assess whether fish spent more time in the no cue choice area vs predator-cue choice area, depending on residual normality (Shapiro–Wilk normality test).
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.
Ley RE, Peterson DA, Gordon JI. Ecological and evolutionary forces shaping microbial diversity in the human intestine. Cell. 2006;124:837–48.
CAS  PubMed  Google Scholar 
2.
Whitman WB, Coleman DC, Wiebe WJ. Prokaryotes: the unseen majority. Proc Natl Acad Sci USA. 1998;95 6578–83.

3.
Gentile CL, Weir TL. The gut microbiota at the intersection of diet and human health. Science. 2018;362:776–80.
CAS  PubMed  Google Scholar 

4.
Penesyan A, Kjelleberg S, Egan S. Development of novel drugs from marine surface associated microorganisms. Mar drugs. 2010;8:438–59.
CAS  PubMed  PubMed Central  Google Scholar 

5.
Chaparro JM, Sheflin AM, Manter DK, Vivanco JM. Manipulating the soil microbiome to increase soil health and plant fertility. Biol Fertil Soils. 2012;48:489–99.
Google Scholar 

6.
Mäder P, Fliessbach A, Dubois D, Gunst L, Fried P, Niggli U. Soil fertility and biodiversity in organic farming. Science. 2002;296:1694–7.
PubMed  Google Scholar 

7.
Classen AT, Sundqvist MK, Henning JA, Newman GS, Moore JA, Cregger MA, et al. Direct and indirect effects of climate change on soil microbial and soil microbial‐plant interactions: what lies ahead? Ecosphere. 2015;6:1–21.
Google Scholar 

8.
de Vries FT, Griffiths RI, Bailey M, Craig H, Girlanda M, Gweon HS, et al. Soil bacterial networks are less stable under drought than fungal networks. Nat Commun. 2018;9:3033.
PubMed  PubMed Central  Google Scholar 

9.
Van Der Heijden MG, Bardgett RD, Van Straalen NM. The unseen majority: soil microbes as drivers of plant diversity and productivity in terrestrial ecosystems. Ecol Lett. 2008;11:296–310.
PubMed  Google Scholar 

10.
Hawkes CV, Wren IF, Herman DJ, Firestone MK. Plant invasion alters nitrogen cycling by modifying the soil nitrifying community. Ecol Lett. 2005;8:976–85.
Google Scholar 

11.
Yuan J, Zhao J, Wen T, Zhao M, Li R, Goossens P, et al. Root exudates drive the soil-borne legacy of aboveground pathogen infection. Microbiome. 2018;6:156.
PubMed  PubMed Central  Google Scholar 

12.
De Corato U, Patruno L, Avella N, Lacolla G, Cucci G. Composts from green sources show an increased suppressiveness to soilborne plant pathogenic fungi: Relationships between physicochemical properties, disease suppression, and the microbiome. Crop Prot. 2019;124:104870.
Google Scholar 

13.
Finkel OM, Castrillo G, Paredes SH, González IS, Dangl JL. Understanding and exploiting plant beneficial microbes. Curr Opin Plant Biol. 2017;38:155–63.
PubMed  PubMed Central  Google Scholar 

14.
Berendsen RL, Pieterse CM, Bakker PA. The rhizosphere microbiome and plant health. Trends plant Sci. 2012;17:478–86.
CAS  PubMed  Google Scholar 

15.
Ploetz RC. Fusarium wilt of banana. Phytopathology. 2015;105:1512–21.
PubMed  Google Scholar 

16.
Cha J-Y, Han S, Hong H-J, Cho H, Kim D, Kwon Y, et al. Microbial and biochemical basis of a Fusarium wilt-suppressive soil. ISME J. 2016;10:119.
CAS  PubMed  Google Scholar 

17.
Gordon TR. Fusarium oxysporum and the Fusarium Wilt Syndrome. Annu Rev Phytopathol. 2017;55:23–39.
CAS  PubMed  Google Scholar 

18.
Laurence MH, Burgess LW, Summerell BA, Liew ECY. High levels of diversity in Fusarium oxysporum from non-cultivated ecosystems in Australia. Fungal Biol. 2012;116:289–97.
CAS  PubMed  Google Scholar 

19.
Nyvad B, Fejerskov O. An ultrastructural-study of bacterial invasion and tissue breakdown in human experimental root-surface caries. J Dent Res. 1990;69:1118–25.
CAS  PubMed  Google Scholar 

20.
Klein E, Ofek M, Katan J, Minz D, Gamliel A. Soil suppressiveness to Fusarium disease: shifts in root microbiome associated with reduction of pathogen root colonization. Phytopathology. 2012;103:23–33.
Google Scholar 

21.
Mendes LW, Mendes R, Raaijmakers JM, Tsai SM. Breeding for soil-borne pathogen resistance impacts active rhizosphere microbiome of common bean. ISME J. 2018;12:3038–42.
CAS  PubMed  PubMed Central  Google Scholar 

22.
Xiong W, Li R, Ren Y, Liu C, Zhao Q, Wu H, et al. Distinct roles for soil fungal and bacterial communities associated with the suppression of vanilla Fusarium wilt disease. Soil Biol Biochem. 2017;107:198–207.
CAS  Google Scholar 

23.
Ye XF, Li ZK, Luo X, Wang WH, Li YK, Li R, et al. A predatory myxobacterium controls cucumber Fusarium wilt by regulating the soil microbial community. Microbiome. 2020;8:17.
Google Scholar 

24.
Zhang S, Raza W, Yang X, Hu J, Huang Q, Xu Y, et al. Control of Fusarium wilt disease of cucumber plants with the application of a bioorganic fertilizer. Biol Fertil Soils. 2008;44:1073.
Google Scholar 

25.
Fu L, Penton CR, Ruan Y, Shen Z, Xue C, Li R, et al. Inducing the rhizosphere microbiome by biofertilizer application to suppress banana Fusarium wilt disease. Soil Biol Biochem. 2017;104:39–48.
CAS  Google Scholar 

26.
Shen Z, Ruan Y, Xue C, Zhong S, Li R, Shen Q, et al. Soils naturally suppressive to banana Fusarium wilt disease harbor unique bacterial communities. Plant. 2015;393:21–33.
CAS  Google Scholar 

27.
De Corato U, Patruno L, Avella N, Salimbeni R, Lacolla G, Cucci G, et al. Soil management under tomato-wheat rotation increases the suppressive response against Fusarium wilt and tomato shoot growth by changing the microbial composition and chemical parameters. Appl Soil Ecol. 2020;154:103601.
Google Scholar 

28.
Zhou D, Jing T, Chen Y, Wang F, Qi D, Feng R, et al. Deciphering microbial diversity associated with Fusarium wilt-diseased and disease-free banana rhizosphere soil. BMC Microbiol. 2019;19:161.
PubMed  PubMed Central  Google Scholar 

29.
da C Jesus E, Marsh TL, Tiedje JM, de S Moreira FM. Changes in land use alter the structure of bacterial communities in Western Amazon soils. ISME J. 2009;3:1004–11.
Google Scholar 

30.
Fierer N, Jackson RB. The diversity and biogeography of soil bacterial communities. Proc Natl Acad Sci USA. 2006;103:626–31.

31.
Mercado-Blanco J, Abrantes I, Barra Caracciolo A, Bevivino A, Ciancio A, Grenni P, et al. Belowground microbiota and the health of tree crops. Front Microbiol. 2018;9:1006.
PubMed  PubMed Central  Google Scholar 

32.
Ramirez KS, Knight CG, de Hollander M, Brearley FQ, Constantinides B, Cotton A, et al. Detecting macroecological patterns in bacterial communities across independent studies of global soils. Nat Microbiol. 2018;3:189–96.
CAS  PubMed  Google Scholar 

33.
Qiu M, Zhang R, Xue C, Zhang S, Li S, Zhang N, et al. Application of bio-organic fertilizer can control Fusarium wilt of cucumber plants by regulating microbial community of rhizosphere soil. Biol Fertil Soils. 2012;48:807–16.
CAS  Google Scholar 

34.
Wang B, Li R, Ruan Y, Ou Y, Zhao Y, Shen Q. Pineapple–banana rotation reduced the amount of Fusarium oxysporum more than maize–banana rotation mainly through modulating fungal communities. Soil Biol Biochem. 2015b;86:77–86.
CAS  Google Scholar 

35.
Alabouvette C. Fusarium wilt suppressive soils: an example of disease-suppressive soils. Australas Plant Pathol. 1999;28:57–64.
Google Scholar 

36.
Hornby D. Suppressive soils. Australas Plant Pathol. 1983;21:65–85.
Google Scholar 

37.
Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010. http://www.bioinformatics.babraham.ac.uk/projects/fastqc.

38.
Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7:335.
CAS  PubMed  PubMed Central  Google Scholar 

39.
Rognes T, Flouri T, Nichols B, Quince C, Mahé F. VSEARCH: a versatile open source tool for metagenomics. PeerJ. 2016;4:e2584.
PubMed  PubMed Central  Google Scholar 

40.
McDonald D, Clemente JC, Kuczynski J, Rideout JR, Stombaugh J, Wendel D, et al. The Biological Observation Matrix (BIOM) format or: how I learned to stop worrying and love the ome-ome. GigaScience. 2012;1:7.
PubMed  PubMed Central  Google Scholar 

41.
Liaw A, Wiener M. Classification and regression by randomForest. R N. 2002;2:18–22.
Google Scholar 

42.
Cortes C, Vapnik V. Support-vector networks. Mach Learn. 1995;20:273–97.
Google Scholar 

43.
Wright RE. Logistic regression. In Grimm LG, & Yarnold PR, (Eds), Reading and understanding multivariate statistics (p. 217–244). Washington, DC: American Psychological Association; 1995.

44.
Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc. 1996;58:267–88.
Google Scholar 

45.
Sing T, Sander O, Beerenwinkel N, Lengauer T. ROCR: visualizing classifier performance in R. Bioinformatics. 2005;21:3940–1.
CAS  PubMed  Google Scholar 

46.
Statnikov A, Henaff M, Narendra V, Konganti K, Li Z, Yang L, et al. A comprehensive evaluation of multicategory classification methods for microbiomic data. Microbiome. 2013;1:11.
PubMed  PubMed Central  Google Scholar 

47.
Zhang J, Liu Y-X, Zhang N, Hu B, Jin T, Xu H, et al. NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice. Nat Biotechnol. 2019;37:676–84.
CAS  PubMed  Google Scholar 

48.
Friedman J, Alm EJ. Inferring correlation networks from genomic survey data. PLoS Comput Biol. 2012;8:e1002687.
CAS  PubMed  PubMed Central  Google Scholar 

49.
Csardi G, Nepusz T. The igraph software package for complex network research. Int J Complex Syst. 2006;1695:1–9.
Google Scholar 

50.
Newman ME. The structure and function of complex networks. SIAM Rev. 2003;45:167–56.
Google Scholar 

51.
Walters W, Hyde ER, Berg-Lyons D, Ackermann G, Humphrey G, Parada A, et al. Improved bacterial 16S rRNA gene (V4 and V4-5) and fungal internal transcribed spacer marker gene primers for microbial community surveys. Msystems. 2016;1:e00009–15.
PubMed  Google Scholar 

52.
McKay G, Brown AE, Bjourson A, Mercer P. Molecular characterisation of Alternaria linicola and its detection in linseed. Eur J Plant Pathol. 1999;105:157–66.
CAS  Google Scholar 

53.
Adams RI, Bateman AC, Bik HM, Meadow JF. Microbiota of the indoor environment: a meta-analysis. Microbiome. 2015;3:49.
PubMed  PubMed Central  Google Scholar 

54.
Cornejo-Granados F, Gallardo-Becerra L, Leonardo-Reza M, Ochoa-Romo JP, Ochoa-Leyva A. A meta-analysis reveals the environmental and host factors shaping the structure and function of the shrimp microbiota. PeerJ. 2018;6:e5382.
PubMed  PubMed Central  Google Scholar 

55.
Koren O, Knights D, Gonzalez A, Waldron L, Segata N, Knight R, et al. A guide to enterotypes across the human body: meta-analysis of microbial community structures in human microbiome datasets. PLoS Comput Biol. 2013;9:e1002863.
CAS  PubMed  PubMed Central  Google Scholar 

56.
Duvallet C, Gibbons SM, Gurry T, Irizarry RA, Alm EJ. Meta-analysis of gut microbiome studies identifies disease-specific and shared responses. Nat Commun. 2017;8:1784.
PubMed  PubMed Central  Google Scholar 

57.
Rocca JD, Simonin M, Blaszczak JR, Ernakovich JG, Gibbons SM, Midani FS, et al. The Microbiome Stress Project: towards a global meta-analysis of environmental stressors and their effects on microbial communities. Front Microbiol. 2018;9:3272.
PubMed  Google Scholar 

58.
Gonzalez A, Navas-Molina JA, Kosciolek T, McDonald D, Vázquez-Baeza Y, Ackermann G, et al. Qiita: rapid, web-enabled microbiome meta-analysis. Nat Methods. 2018;15:796–8.
CAS  PubMed  PubMed Central  Google Scholar 

59.
Baxter NT, Ruffin MT, Rogers MA, Schloss PD. Microbiota-based model improves the sensitivity of fecal immunochemical test for detecting colonic lesions. Genome Med. 2016;8:37.
PubMed  PubMed Central  Google Scholar 

60.
Belk A, Xu ZZ, Carter DO, Lynne A, Bucheli S, Knight R, et al. Microbiome data accurately predicts the postmortem interval using random forest regression models. Genes. 2018;9:104.
PubMed Central  Google Scholar 

61.
Wilck N, Matus MG, Kearney SM, Olesen SW, Forslund K, Bartolomaeus H, et al. Salt-responsive gut commensal modulates TH17 axis and disease. Nature. 2017;551:585–9.
CAS  PubMed  PubMed Central  Google Scholar 

62.
Manici L, Caputo F. Fungal community diversity and soil health in intensive potato cropping systems of the east Po valley, northern Italy. Ann Appl Biol. 2009;155:245–58.
Google Scholar 

63.
Ploetz RC. Fusarium wilt of banana is caused by several pathogens referred to as Fusarium oxysporum f. sp. cubense. Phytopathology. 2006;96:653–6.
PubMed  Google Scholar 

64.
Shen Z, Xue C, Taylor PWJ, Ou Y, Wang B, Zhao Y, et al. Soil pre-fumigation could effectively improve the disease suppressiveness of biofertilizer to banana Fusarium wilt disease by reshaping the soil microbiome. Biol Fertil Soils. 2018;54:793–806.
CAS  Google Scholar 

65.
Wang B, Li R, Ruan Y, Ou Y, Zhao Y, Shen Q. Pineapple–banana rotation reduced the amount of Fusarium oxysporum more than maize–banana rotation mainly through modulating fungal communities. Soil Biol Biochem. 2015a;86:77–86.
CAS  Google Scholar 

66.
Wu X, Guo S, Jousset A, Zhao Q, Shen Q. Bio-fertilizer application induces soil suppressiveness against Fusarium wilt disease by reshaping the soil microbiome. Soil Biol Biochem. 2017;114:238–47.
Google Scholar 

67.
Forsyth LM, Smith LJ, Aitken EA. Identification and characterization of non-pathogenic Fusarium oxysporum capable of increasing and decreasing Fusarium wilt severity. Mycological Res. 2006;110:929–35.
Google Scholar 

68.
Lemanceau P, Alabouvette C. Biological control of fusarium diseases by fluorescent Pseudomonas and non-pathogenic Fusarium. Crop Prot. 1991;10:279–86.
Google Scholar 

69.
Liu L, Kong J, Cui H, Zhang J, Wang F, Cai Z, et al. Relationships of decomposability and C/N ratio in different types of organic matter with suppression of Fusarium oxysporum and microbial communities during reductive soil disinfestation. Biol Control. 2016;101:103–13.
CAS  Google Scholar 

70.
Pieretti I, Royer M, Barbe V, Carrere S, Koebnik R, Cociancich S, et al. The complete genome sequence of Xanthomonas albilineans provides new insights into the reductive genome evolution of the xylem-limited Xanthomonadaceae. BMC genomics. 2009;10:616.
PubMed  PubMed Central  Google Scholar 

71.
Li X, Zhang YN, Ding C, Jia Z, He Z, Zhang T, et al. Declined soil suppressiveness to Fusarium oxysporum by rhizosphere microflora of cotton in soil sickness. Biol Fertil soils. 2015;51:935–46.
Google Scholar 

72.
Wu L, Yang B, Li M, Chen J, Xiao Z, Wu H, et al. Modification of rhizosphere bacterial community structure and functional potentials to control pseudostellaria heterophylla replant disease. Plant Dis. 2019;104:25–34.
PubMed  Google Scholar 

73.
Shang Q, Yang G, Wang Y, Wu X, Zhao X, Hao H, et al. Illumina-based analysis of the rhizosphere microbial communities associated with healthy and wilted Lanzhou lily (Lilium davidii var. unicolor) plants grown in the field. World J Microbiol Biotechnol. 2016;32:95.
PubMed  Google Scholar 

74.
Abbasi S, Safaie N, Sadeghi A, Shamsbakhsh M. Streptomyces strains induce resistance to Fusarium oxysporum f. sp. lycopersici race 3 in tomato through different molecular mechanisms. Front Microbiol. 2019;10:1505.
PubMed  PubMed Central  Google Scholar 

75.
Liotti RG, da Silva Figueiredo MI, Soares MA. Streptomyces griseocarneus R132 controls phytopathogens and promotes growth of pepper (Capsicum annuum). Biol Control. 2019;138:104065.
CAS  Google Scholar 

76.
Tahvonen R. Microbial control of plant diseases with Streptomyces spp. 1. Eppo Bull. 2010;18:55–9.
Google Scholar 

77.
Yang Y, Zhang S-W, Li K-t. Antagonistic activity and mechanism of an isolated Streptomyces corchorusii stain AUH-1 against phytopathogenic fungi. World J Microbiol Biotechnol. 2019;35:145.
PubMed  Google Scholar 

78.
Cha J-Y, Han S, Hong H-J, Cho H, Kim D, Kwon Y, et al. Microbial and biochemical basis of a Fusarium wilt-suppressive soil. ISME J. 2015;10:119–29.
PubMed  PubMed Central  Google Scholar  More

Data collection and sequencing
In June 2016, diversity data were collected at 61 sites in a 740 km2 sub-catchment of the river Thur, northeastern Switzerland. No singular rain events took place during the sampling days, and thus all sampling was performed at base-flow conditions. Sampling sites were chosen in order to proportionally represent all stream orders in the river network (Supplementary Fig. 1) and to span the complete geographical extent of the catchment. At each site, a standardized35 3-min kicknet sampling applied to three microhabitats was performed to collect benthic macroinvertebrates. The subsamples from the three microhabitats were pooled and stored in ethanol. In the laboratory, debris was removed and all individuals belonging to may-, stone-, and caddisflies (Ephemeroptera, Plecoptera, and Trichoptera, abbreviated as EPT) were identified under a stereomicroscope to the genus or species level if applicable. One sample was lost due to handling error, and thus we subsequently only had EPT kicknet sample data from 60 sites.
At the same 61 sites at which kicknet data were collected, also eDNA samples were taken. We collected three independent samples of 250 mL of river water at each site (sampled below the surface and well above the river bottom). In small streams (up to about 1-m wide), the water was taken from the middle of the stream, while in larger streams the water was taken about 0.5 m away from the shore side. All three water samples per site were filtered on site on separate GF/F filters (pore size 0.7 µm Whatman International Ltd.), which were stored on ice immediately after filtering, and frozen at −20 °C within a few hours. Subsequently, these three samples were analyzed separately and as independent replicates. To prevent cross-contamination, the eDNA samples were collected a few meters upstream of the kicknet sampling, and kicknet sampling and eDNA sampling were performed by different people.
DNA was extracted with the DNeasy Blood and Tissue Kit (Qiagen GmbH) from the filter in a lab dedicated to low DNA-concentrated environmental samples (i.e., clean lab with positive air pressure, well-defined procedures and all work conducted under sterile benches, and no handling of high-DNA concentration samples or post-PCR products). A short barcoding region of the COI gene58,59 was targeted, and a dual-barcoded two-step PCR amplicon sequencing protocol for Illumina MiSeq was performed. In short, three PCR replicates were performed for each eDNA replicate with primers that contained an Illumina adapter-specific tail, a heterogeneity spacer and the amplicon target site. These three PCR replicates per sample replicate were pooled and indexed in a second PCR with the Nextera XT Index Kit v2 (lllumina). We measured the concentration of all indexed PCR reactions and pooled them in equimolar parts to a final library, which we ran twice on two consecutive Illumina MiSeq runs to increase sequencing depth. Raw data were demultiplexed and read quality was checked with FastQC60. Thereafter, end-trimming (usearch, version 10.0.240) and merging (Flash, version 1.2.11) of raw reads were performed, before primer sites were removed (cutadapt, version 1.12) and reads were quality‐filtered (prinseq‐lite, version 0.20.4). Next, we used UNOISE361, which has a built‐in error correction to reduce the influence of sequencing errors, to determine amplicon sequence variants (ZOTUs). To reduce sequence diversity, we implemented an additional clustering at 99% sequence identity. We checked ZOTUs for stop codons of the invertebrate mitochondrial code, to ensure an intact open reading frame. For the taxonomic assignment, all COI‐related sequences were downloaded from NCBI and ZOTUs were blasted against the NCBI COI collection. We extracted the top five best blast hits and then used the R packages “taxize”62 (version 0.9.7) and “rentrez”63 (version 1.2.2) to acquire taxonomic labels. We modified the selected COI sequences with the taxonomic labels in order to index the database. The ZOTUs were then assigned to taxa using Sintax (usearch) and the NCBI COI‐based reference. Details of the data collection and the bioinformatic pipeline can be further assessed in Mächler et al.36. The OTUs found in this dataset are generally saturating39. We acknowledge that the primers used are targeting eukaryotic diversity in general and were not only specific to EPT taxa, which might lead to an underestimation of detections of EPT by eDNA (i.e., false absences). This, however, is conservative with respect to our findings, as we would miss taxa with eDNA where the kicknet sample indicated their presence.
Network extraction and hydrological characterization
The river network was extracted from a 25-m Digital Elevation Model (DEM) of the region provided by the Swiss Federal Institute of Topography (Swisstopo) by applying the D8 method through a TauDEM routine in a GIS software64. A threshold of 800 pixels (0.5 km2) on contributing area was applied to identify the channelized (i.e., perennial, see ref. 64) portion of the drainage network, whose total length equaled 751 km. We then defined a reach as a sequence of pixels of the DEM matrix, originating from a source or a confluence, directed downstream, and ending at the following confluence or at the outlet. This resulted in a discretization of the river network comprising 760 reaches of median length 0.78 km (95% equal-tailed interval of the distribution: 0.07–3.18 km). Note that a lower threshold area would have resulted in a higher number of reaches, implying a more refined discretization of the river network but also an increased computational burden for the subsequent model runs. Hence, we chose the highest value of threshold area such that: i) the extracted river network retained all reaches where sampling sites were located; ii) a qualitative comparison with the vectorial hydrographic network provided by Swisstopo resulted adequate. For modeling purposes, reaches were considered as nodes of a graph with edges following flow direction. All variables referred to a node were considered homogeneously distributed within the corresponding reach.
In order to assess values of hydrological variables for all reaches, we made use of power-law scaling relationships, a well-established and universally applicable concept in hydrology42,43,44. In particular, river width w, river depth D and water velocity v are known to scale (both within a single cross-section and in the downstream direction) as a power-law of water discharge Q (refs. 42,43); along the flow direction, the relationships w~Q0.5, D~Q0.4, w~Q0.1 are valid over wide ranges of natural streams42. Moreover, water discharge scales linearly across a catchment with drainage area A (ref. 44). Strictly speaking, the relationship Q~A holds for mean annual values of Q; however, it can be reasonably extended to values of Q averaged over shorter time windows (say, at least 1 day), provided that the time scale of flow propagation is much shorter than 1 day, and that rainfall (and the resulting runoff generation) can be considered spatially homogeneous if aggregated at a daily scale. Both assumptions are reasonable for catchments up to 103 km2 (ref. 65), as the one here studied. Mean water discharges during the sampling days and stage-discharge relationships were available at four stations operated by the Swiss Federal Office for the Environment (FOEN) (see Fig. 1). River widths at these locations were estimated via aerial images. Power law relationships with drainage area for discharge and river width were then fitted on the four stations, yielding Q = 0.072 A1.056 and w = 1.586A0.526, where A is in km2, Qin m3s−1, and w in m. As for river depth, we discarded the station with lowest drainage area because we observed that the values of depth measured therein were highly overestimated with respect to the expected values, based on the other stations and the scaling exponent 0.4 (ref. 42). Hence, we limited the fit to the three other stations and obtained D = 0.073 A0.463, where D is in m. By assuming rectangular river’s cross-sections, we finally derived a power-law relationship linking water velocity v to drainage area: v = Q/(Dw) = 0.623 A0.067, where v is in ms−1. Notably, all scaling exponents obtained were very close to the literature values42,44. Details on the fit of these hydrological relationships are reported in Supplementary Fig. 2.
Choice of covariates
A first set of 18 covariates was chosen as representative of morphological, geological and land cover features of the catchment. These covariates can either reflect local or upstream characteristics. Land cover covariates were evaluated as local values, because they are assumed to potentially have a role in determining local taxon suitability (see also ref. 66). Geological covariates were calculated as upstream averaged values, because they are likely to affect the chemical composition of streamflow at a site; such process could affect local habitat suitability but is driven by the upstream catchment30. Morphological covariates were obtained by analyzing the DEM of the region, while geological and land cover raster maps were provided by Swisstopo67,68. The list of covariates is shown in Supplementary Table 1.
A second set of covariates was constituted by grouping all 760 reaches into 17 clusters according to geographical proximity and hydrological connectivity (see Fig. 1c). For each cluster, a corresponding covariate vector was defined with values equal to one for the reaches constituting that cluster, and zero otherwise. The addition of these ‘geographical’ covariates aimed at allowing eDITH to reproduce spatial patterns uncorrelated to any of the previous covariates (e.g., in the case of a taxon only inhabiting a single tributary of the catchment).
The 35 covariates were z-normalized and checked for possible multicollinearity. The 760-by-34 design matrix obtained by removing the covariate corresponding to cluster LUT of Fig. 1c (corresponding to the Luteren tributary, and arbitrarily chosen among the geographical covariates: each of them can indeed be expressed as a linear combination of the other 16 geographical covariates) yielded all pairwise correlation coefficients with absolute values lower than 0.80 and variance inflation factors lower than 10, which minimized the effects of multicollinearity69. The LUT covariate was however kept in the design matrix to facilitate model fitting.
The eDITH model
The eDNA transport component of the eDITH model is derived from Carraro et al.31

$${hat{C}}_{j} = frac{1}{{Q_j}}sum_{i in gamma left( j right)} A_{S,i}exp left( { – frac{{L_{ij}}}{{v_{ij}tau }}} right)p_i,$$
(1)

where (hat C_j), the modeled eDNA concentration of a given taxon at site j, results from the summation of eDNA production rates pi across all network nodes upstream of j. These contributions are weighted by the relative source area AS,i and a first-order decay factor, where Lij is the along-stream length of the path joining i to j, vij the average water velocity along that path, and τ a characteristic decay time. Finally, Qj is a characteristic (mean during sampling days) value of water discharge across node j. Production rates are assumed to be proportional to the taxon’s density and expressed via an exponential link

$$p_i = p_0exp left( {{boldsymbol{beta }}^T{boldsymbol{X}}left( i right)} right),$$
(2)

where X(i) is a vector of (normalized) environmental covariates evaluated at node i, β a vector of covariate effects and p0 is a characteristic production rate31. Equation (2) represents a generalized linear model, a widely used class of species distribution models70. Note that Eq. (1) does not explicitly account for deposition of eDNA particles on the river bed. Decay and deposition are both relevant processes affecting eDNA removal, and thus detection, in stream water (see e.g., ref. 25). We chose not to model these processes independently because it would be hardly possible to disentangle the roles of degradation of genetic material and that of gravity-induced deposition. Decay time τ is therefore to be seen as a characteristic time (linked to half-life T1/2 by the relationship T1/2 = In(2)τ) for the presence of eDNA in water, regardless of the source of its depletion.
We here further hypothesize that the expected number of reads for a given taxon at a given site is proportional to the taxon’s eDNA concentration in the sample: (hat N_j = khat C_j). Equation (1) then becomes

$${hat{N}}_{j} = frac{1}{{Q_j}}sum_{i in gamma left( j right)} A_{S,i}exp left( { – frac{{L_{ij}}}{{v_{ij}tau }}} right)p_0^prime exp left( {{boldsymbol{beta }}^T{boldsymbol{X}}left( i right)} right),$$
(3)

where (p_0^prime = kp_0). Finally, by denoting with Li and wi the length and width of reach i, respectively, it is AS,i = Liwi, (v_{ij} = sum_{k in P_{ij}} L_k/sum_{k in P_{ij}} left( {L_k/v_k} right)), where Pij is the set of nodes constituting the path joining i to j. For the sake of clarity, all mathematical symbols used here and in the following are listed in Supplementary Table 2.
Model calibration
Measured numbers of reads at a given site are assumed to be distributed according to a geometric distribution with mean (hat N_j) obtained from Eq. (3). The geometric distribution is a special case of the negative binomial distribution, obtained by setting to one the parameter corresponding to the number of failures. Reasons for this choice are multifold: (i) it is a discrete distribution, as it is the case for read number data; (ii) it depends on a single parameter, thereby it reduces the complexity of the problem; (iii) it has null mode, which is in agreement with the data (in fact, the number of replicates with null read numbers ranged from 50 to 182 out of 183, depending on the genus); (iv) the distribution has fatter tails compared to the Poisson distribution, which is also a single-parameter, discrete distribution with null mode. Fat tails allow the interpretation of large observed numbers of reads.
Given these assumptions on the distribution of observed numbers of reads, the likelihood function reads

$$fleft( {{boldsymbol{N}}{mathrm{|}}{boldsymbol{beta }},p_{0}^{prime} ,tau } right) = mathop {prod }limits_{j in S} left[ {mathop {prod }limits_{o = 1}^3 frac{1}{{1 + hat N_j}}left( {frac{{{hat{N}}_j}}{{1 + {hat{N}}_j}}} right)^{N_{jo}}} right],$$
(4)

where S is the set of sampling sites used to calibrate the model and N is a |S|-by-3 matrix whose entries are observed numbers of reads Njo at site j and replicate o (|S| is the cardinality of S, equal to 61 in this case).
The posterior distribution of parameters β, (p_0^prime), τ was sampled by means of an Adaptive Metropolis algorithm71. The prior distribution adopted for β was a multivariate normal with null mean and covariance matrix equal to 9 ∗ I, where I is the identity matrix of order equal to the number of components of β. We then adopted a uniform prior distribution for (p_0^prime) and a lognormal prior distribution for τ with mean equal to 2.55 h and standard deviation equal to 1.36 h, equivalent to a normal distribution for ({mathrm{ln}}left( tau right)) (with τ expressed in seconds) with mean equal to 9 and standard deviation equal to 0.5; such prior distribution was derived from a previous study31. For each of the 50 model runs (corresponding to the 50 EPT genera detected in the eDNA samples), Markov chains were randomly initialized; the first 5000 elements of the chains built by the Adaptive Metropolis algorithm were then discarded (burn-in phase), while the following 10,000 were retained. With the above-mentioned settings, eDITH was used to generate spatial patterns of relative density for all of the 50 genera.
Assessing detection probability and presence maps
In order to enable the comparison among maps of relative density for the different genera, we resorted to the evaluation of detection probability. This was defined as the probability that, given the relative density of a genus at a given reach predicted by eDITH, an eDNA sample taken from that reach would yield a nonzero number of reads, if that reach were unconnected from the river network. From Eq. (3), the expected number of reads of a sample taken from an unconnected reach reads

$$hat N_{U,i} = p_i^prime frac{{A_{S,i}}}{{q_i}}exp left( { – frac{{L_i}}{{v_itau }}} right),$$
(5)

where qi is the discharge directly contributing to reach i (such that (Q_i = sum_{k in gamma left( i right)} q_k)), Li and vi are length and water velocity relative to reach i, respectively. Thus, according to the assumption of geometric distribution for observed read numbers, the probability of a non-zero read number for a sample with expected read number (hat N_{U,i}) is

$$P_{D,i} = frac{{hat N_{U,i}}}{{1 + hat N_{U,i}}}.$$
(6)

Detection probabilities are evaluated by using median posterior values for parameters β, (p_0^prime), τ. Finally, maps of detection probability are converted into presence maps by applying the threshold36PD,i ≥ 2/3. The summation of presence maps for all genera yields the genus richness map displayed in Fig. 5a.
A number of methods for threshold selection in species distribution models exist in the literature72; however, preliminary analyses (not reported) showed that methods such as the prevalence approach and the average probability/suitability approach (see ref. 72) resulted in lower accuracy estimates (see below) as opposed to the fixed threshold approach. Moreover, the fixed threshold approach PD,i ≥ 2/3 was used by Mächler et al.36; since the eDITH model is to be seen as a way to translate “upstream-averaged” eDNA measurements into local equivalent eDNA measurements, it appears appropriate to keep the same criterion already used with the raw eDNA data.
For a given genus that was found both in the eDNA and kicknet samples, the accuracy of a presence map obtained by eDITH with respect to the kicknet data was determined as the fraction of true positives (i.e., sites where both model and kicknet assessed presence) and true negatives (i.e., sites where both model and kicknet assessed absence) over all (60) sites where kicknet was performed. On the contrary, false positives are sites where eDITH predicts presence while kicknet indicates absence; the vice versa holds for false negatives.
We underline that the eDITH model was run for all genera that were found in at least one replicate at one site. Such an inclusive choice was operated in a bid to maximize the available information and avoid false absences. However, the estimates of occurrence of taxa performed by our model remain conservative: indeed, for the 8 genera that were never found at any site with at least 2 out of 3 nonzero read numbers (see Supplementary Data 1), the fraction of river reaches where eDITH predicted presence was always lower than 3%.
Note that calculating the accuracy between model predictions and eDNA data with the method used for the model vs. kicknet comparison would be formally wrong: indeed, eDNA data are an aggregate measure of the upstream taxon distribution, while eDITH-based presence estimates refer to a local variable, because upstream contributions have been disentangled by the model. Conversely, an adequate measure of consistency between model predictions and eDNA data must compare quantities of the same type (i.e., both referred to the upstream catchment). To this end, we utilized the goodness-of-fit test detailed below.
Goodness-of-fit test
An ad hoc goodness-of-fit test was required owing to the use of a geometric distribution to represent the observed numbers of reads. The test relies on a bootstrapping technique derived from that proposed by Mi et al.73, to which the reader is referred for details on the theoretical foundation of the method.
For a given site j, let (hat N_j) be the read number predicted by eDITH (see Eq. (3)) and Njo, o = 1,2,3 the triplet of observed read numbers at that site. The sample standard deviation of the data with respect to the predicted mean reads

$$s_j = sqrt {mathop {sum }limits_{o = 1}^3 left( {N_{jo} – hat N_j} right)^2} ,$$
(7)

while Pearson’s residuals are given by (r_{jo} = left( {N_{jo} – hat N_j} right)/s_j). Now, let (tilde N_jsim {mathrm{Geom}}left( {hat N_j} right)) be a random variable that is geometrically distributed with mean equal to (hat N_j). Let us generate a large (h = 1, …, 105) number of triplets (tilde N_{jho}), o = 1, 2, 3 and compute their sample standard deviation (tilde s_{jh}) and residuals (tilde r_{jho}):

$$tilde s_{jh} = sqrt {mathop {sum }limits_{o = 1}^3 left( {tilde N_{jho} – hat N_j} right)^2} ;,tilde r_{jho} = frac{{tilde N_{jho} – hat N_j}}{{tilde s_{jh}}}.$$
(8)

Let now dj and (tilde d_{jh}) be the sum of squared deviations of ordered residuals (of the data and of the sampled distributions, respectively) from the medians of their sampled distributions

$$d_j = mathop {sum }limits_{o = 1}^3 left( {r_{jo} – tilde r_{jo}^{left( {50} right)}} right)^2;,tilde d_{jh} = mathop {sum }limits_{o = 1}^3 left( {tilde r_{jho} – tilde r_{jo}^{left( {50} right)}} right)^2,$$
(9)

where (tilde r_{jo}^{left( {50} right)}) is the median of the residuals of the o-th components of the generated triplets (tilde N_{jho}). The p-value can thus be computed as

$$p_j^{{mathrm{GOF}}} = frac{{1 + mathop {sum }nolimits_{h = 1}^{10^5} {mathbf{1}}_{tilde d_{jh}!ge! d_j}left( h right)}}{{1 + 10^5}},$$
(10)

where ({mathbf{1}}_{tilde d_{jh}! ge! d_j}left( h right)) is an indicator function equal to one if (tilde d_{jh} ge d_j) and null otherwise. We finally assumed that the null hypothesis H0 that the triplet Njo, o = 1, 2, 3 is geometrically distributed with mean (hat N_j) (namely, that the model correctly reproduces the data at site j) cannot be rejected if (p_j^{GOF}) > 0.05.
Cross-validation analysis of the effect of sample size
In order to evaluate how model predictions vary as a function of data availability, we performed 9 additional simulations (termed AS), in which only subsets of the sampling sites were used to calibrate eDITH: 3 of them (simulation group termed AS1) were based on 49 out of 61 sites (~80%), 3 (AS2) were based on 37 sites (~60%) and 3 (AS3) were based on 25 sites (~40%). For each simulation, a quasi-random choice of the excluded subset of sites was operated. In particular, to ensure spatial coverage of the catchment, we first imposed that, for each AS group, the distribution of Strahler stream order values of the calibration subset was proportional to that of the complete set of sites (shown in Supplementary Fig. 2). The allocation of number of excluded sites per stream order value was determined via the D’Hondt method74. The excluded sites were then randomly sampled while respecting the constraint on stream order-based allocation. Note that, for a given thus-obtained calibration set, eDITH was run for all 50 genera.
For these additional simulations, calibration was performed as described with regards to the complete model (CM) (see “Model calibration” section), except that Markov Chains were here shorter (burn-in phase: 1000 elements; length of the retained chain: 5000 elements) to speed up the computational process.
Finally, performance indices were calculated in analogy to CM. Specifically, these were goodness of fit and accuracy. Goodness of fit is defined as the ability of the modeled spatial patterns of expected read numbers to reproduce the observed read numbers, and is evaluated as the fraction of sites (either all sampling sites, calibration or validation subsets) where the null hypothesis that observed read numbers come from the hypothesized distribution cannot be rejected (see “Goodness-of-fit test” section). Accuracy is defined as the agreement between model-based genus presence predictions and kicknet observations, and is evaluated as the fraction of sites marked as true positives or true negatives (see “Assessing detection probability and presence maps” section). We then calculated loss of goodness of fit as the difference between the goodness of fit of the AS simulations and that of CM. Loss of accuracy is analogously calculated.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

Media and culture conditions
All bacterial culturing steps, with or without phages, were undertaken at room temperature (~ 22 °C) using R2A agar and R2B broth (Alpha Biosciences Inc., Baltimore, Maryland, USA) as growth media. Liquid cultures were grown with agitation at 200 rpm. Between uses, bacterial isolates were stored on R2A agar at 4 °C or as cryostocks at − 80 °C (~ 20% glycerol). Purified phage suspensions were stored at 4 °C between uses. The buffer solution (SM buffer) for phage resuspension contained 100 mM NaCl, 8 mM MgSO4 and 50 mM Tris–Cl (pH 7.5) in autoclaved Milli-Q water.
Bacterial isolation from sampled groundwater
Groundwater from 11 wells at the Skelstofte test sites (Vedby, Denmark, 54°52′34.1ʺN 11°16′18.7ʺE) was sampled at depths of between 7 and 12 m below the surface. These test sites are used to study bioremediation and monitoring strategies of groundwater contamination24. Briefly, to obtain water samples representative of the reservoir, wells were purged of ~ 13–22 times their volume (well ø = 4 mm) using a peristaltic pump with a 0.2 L min−1 flow rate to avoid disruption of the natural water movement, and stable values were observed for O2, pH, conductivity, temperature and other redox parameters before sampling (individual parameters varied between wells) (data not shown). Subsequently, 1-L samples were extracted into sterile glass bottles and stored at 4 °C until use. From the samples, heterotrophic bacteria were isolated on R2A agar (Alpha Biosciences Inc., Baltimore, Maryland, USA), which gave a collection of 84 isolates representing diverse colony morphotypes.
In brief, sample bottles with groundwater were shaken thoroughly and 100 µL of each well sample was plated and left to grow for 4–14 days. Colony growth from samples varied greatly and library building was based on acquiring all the colony morphotypes throughout the incubation period. Where possible, several colonies of a distinct morphotype (up to a maximum of eight) were picked. To obtain pure isolates, colonies were re-streaked three times.
Phage screening, isolation and amplification
Phage activity was detected in a pre-screening of phage enrichment cultures containing sample water, concentrated media and the relevant isolate. Screening hits, indicated by clear or turbid zones in the top agar layer as a result of plausible phage activity, were then followed by a new enrichment to verify phage activity and subsequently obtain pure phage isolates.
The initial enrichment and pre-screening were performed by mixing 2.5 mL 2 × R2B, 2.5 mL sample water and then inoculating with 50 µL culture of the relevant isolate (isolated from that sample water). Following five days of growth, enrichment cultures were centrifuged at 10,000 g for three minutes and supernatant was filtered through a 0.22-µm sterile syringe PVDF membrane filter (Millex Durapore, Burlington, MA, USA). Phage activity in the filtrate from each enrichment was screened in a standard double agar overlay assay. Briefly, 10 µL of each filtered enrichment was drop-plated in 10 replicates onto a semisolid agar, R2B + 50 mM CaCl2/MgCl2 + 0.6% agarose, inoculated with 100 µL culture of relevant bacterial isolate, using R2A as the bottom agar layer. Plates were incubated at room temperature and inspected daily for one week, with any potential plaque forming noted as ‘hits’ for groundwater sample and isolate.
The follow-up enrichment cultures of phage-host ‘hits’ contained 4 mL 10 × R2B broth, 35 mL groundwater sample, 1 mL culture of relevant host and 10 mM CaCl2/MgCl2. Enrichment cultures were incubated for 24 h. Subsequently, the culture was supplemented with 1 M NaCl and incubated for 30 min. The culture was then centrifuged at 5,000 g for five minutes and the supernatant was filtered through a 0.45-µm PVDF syringe filter (Millex Durapore, Burlington, MA, USA). Then 100 µL of the filtered supernatant was used in a double agar overlay assay as described above, and individual plaques were picked and re-plated three times to obtain pure phage samples.
A high titre of both phages was obtained by polyethylene glycol (PEG) precipitation as described in9 with modifications. Briefly, 200 mL R2B was inoculated with 100 µL host culture, infected with 100 µL of pure phage lysate and left overnight. The following day, 1 M NaCl was added to cultures and left on an orbital shaker for one hour to burst bacterial cells. Cultures were then spun at 12,000 g for 10 min, supernatant was collected and PEG was added to reach 10% w/v. The mixture was left for two hours on an orbital shaker for phage-PEG adsorption. Finally, the mixture was centrifuged at 12,000 g for 10 min and the pellets with phage virions were resuspended and collected in 3–5 mL SM buffer. Following PEG precipitation, titre was determined by PFU counts by drop-plating dilution series of collected particles on an agar overlayer, as described above. A titre range between 1010–1011 PFU mL−1 was obtained. Single plaques or high-titre phage samples (1010–1011 PFU mL−1) were used in downstream analysis to characterise and determine the phage’s (i) morphology, (ii) burst size, (iii) general genome features (iv) homology with other related phages, and (v) virion proteins. For TEM imaging and peptide sequencing, phage samples were further purified by a caesium chloride gradient according to the protocol of Clokie & Kropinski25.
Host DNA extraction, sequencing and identification
DNA was extracted from 5 mL of each host culture (grown for five days) using an Ultraclean Microbial DNA Isolation Kit (Mo bio Laboratories, Carlsbad, California, USA), following the manufacturer’s protocol. DNA libraries were prepared with the Nextera XT DNA kit (Illumina, San Diego, USA) according to the manufacturer’s protocol. The prepared libraries were then sequenced in a 2 × 251 paired-end sequencing run, as part of a flow cell, using the Illumina MiSeq v2 kit (Illumina, San Diego, CA, USA). In assembling host draft genomes, CutAdapt (v1.8.3) was used to quality-trim sequence reads (bases with  More

1.
Cerling, T. E. et al. Stable isotope-based diet reconstructions of Turkana Basin hominins. Proc. Natl. Acad. Sci. USA 110, 10501–10506 (2013).
ADS  CAS  PubMed  Google Scholar 
2.
Wynn, J. G. et al. Geological and palaeontological context of a Pliocene juvenile hominin at Dikika, Ethiopia. Nature 443, 332–336 (2006).
ADS  CAS  PubMed  Google Scholar 

3.
van der Merwe, N., Masao, F. & Bamford, M. Isotopic evidence for contrasting diets of early hominins Homo habilis and Australopithecus boisei of Tanzania. South Afr. J. Sci. 104, 153–155 (2008).
Google Scholar 

4.
White, T. D. et al. Ardipithecus ramidus and the paleobiology of early hominids. Science 326, 64–86 (2009).
ADS  Google Scholar 

5.
Cerling, T. E. et al. Diet of Paranthropus boisei in the early Pleistocene of East Africa. Proc. Natl. Acad. Sci. USA 108, 9337–9341 (2011).
ADS  CAS  PubMed  Google Scholar 

6.
Henry, A. G. et al. The diet of Australopithecus sediba. Nature 487, 90–93 (2012).
ADS  CAS  PubMed  Google Scholar 

7.
Levin, N., Haile-Selassie, Y., Frost, S. R. & Saylor, B. Z. Dietary change among hominins and cercopithecids in Ethiopia during the early Pliocene. Proc. Natl. Acad. Sci. USA 112, 12304–12309 (2015).
ADS  CAS  PubMed  Google Scholar 

8.
Lee-Thorp, J. 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).
ADS  CAS  PubMed  Google Scholar 

9.
Sponheimer, M. et al. Isotopic evidence for dietary variability in the early hominin Paranthropus robustus. Science 314, 980–982 (2006).
ADS  CAS  PubMed  Google Scholar 

10.
Scott, R. S. et al. Dental microwear texture analysis reflects diets of living primates and fossil hominins. Nature 436, 693–695 (2005).
ADS  CAS  PubMed  Google Scholar 

11.
Martin, J. E., Vance, D. & Balter, V. Magnesium stable isotope ecology using mammal tooth enamel. Proc. Natl. Acad. Sci. USA 112, 430–435 (2015).
ADS  CAS  PubMed  Google Scholar 

12.
Martin, J. E., Tacail, T. & Balter, V. Non-traditional isotope perspectives in vertebrate palaeobiology. Palaeontology 60, 485–502 (2017).
Google Scholar 

13.
Jaouen, K. & Pons, M. L. Potential of non-traditional isotope studies for bioarchaeology. Archaeol. Anthropol. Sci. 9, 1389–1404 (2017).
Google Scholar 

14.
Tacail, T., Le Houedec, S. & Skulan, J. L. New frontiers in calcium stable isotope geochemistry: perspectives in present and past vertebrate biology. Chem. Geol. https://doi.org/10.1016/j.chemgeo.2020.119471 (2020).

15.
Bourgon, N. et al. Zinc isotopes in Late Pleistocene fossil teeth from a Southeast Asian cave setting preserve paleodietary information. Proc. Natl. Acad. Sci. USA 117, 4675–4681 (2020).
CAS  PubMed  Google Scholar 

16.
Balter, V. et al. Were Neandertalians essentially carnivores? Sr and Ba preliminary results of the mammalian palaeobiocoenosis of Saint-Césaire. C. R. Acad. Sci. IIA 332, 59–65 (2001).
Google Scholar 

17.
Sponheimer, M. & Lee-Thorp, J. A. Enamel diagenesis at South African Australopith sites: implications for paleoecological reconstruction with trace elements. Geochim. Cosmochim. Acta 70, 1644–1654 (2006).
ADS  CAS  Google Scholar 

18.
Balter, V., Braga, J., Télouk, P. & Thackeray, J. F. Evidence for dietary change but not landscape use in South African early hominins. Nature 489, 558–560 (2012).
ADS  CAS  PubMed  Google Scholar 

19.
Joannes-Boyau, R. et al. Elemental signatures of Australopithecus africanus teeth reveal seasonal dietary stress. Nature 572, 112–115 (2019).
CAS  PubMed  Google Scholar 

20.
Cerling, T. E. et al. Global vegetation change through the Miocene/Pliocene boundary. Nature 389, 153–158 (1997).
ADS  CAS  Google Scholar 

21.
Ungar, P. S., Grine, F. E. & Teaford, M. F. Dental microwear and diet of the Plio-Pleistocene hominin Paranthropus boisei. PLoS ONE 3, e2044 (2008).
ADS  PubMed  PubMed Central  Google Scholar 

22.
Skulan, J. & DePaolo, D. J. Calcium isotope fractionation between soft and mineralized tissues as a monitor of calcium use in vertebrates. Proc. Natl. Acad. Sci. USA 96, 13709–13713 (1999).
ADS  CAS  PubMed  Google Scholar 

23.
Clementz, M. T., Holden, P. & Koch, P. L. Are calcium isotopes a reliable monitor of trophic level in marine settings? Intl. J. Osteoarchaeol. 13, 29–36 (2003).
Google Scholar 

24.
Martin, J. E., Tacail, T., Adnet, S., Girard, C. & Balter, V. Calcium isotopes reveal the trophic position of extant and fossil elasmobranchs. Chem. Geol. 415, 118–125 (2015).
ADS  CAS  Google Scholar 

25.
Reynard, L. M., Henderson, G. M. & Hedges, R. E. M. Calcium isotope ratios in animal and human bone. Geochim. Cosmochim. Acta. 74, 3735–3750 (2010).
ADS  CAS  Google Scholar 

26.
Li, Q., Thirwall, M. & Müller, W. Ca isotopic analysis of laser-cut microsamples of (bio)apatite without chemical purification. Chem. Geol. 422, 1–12 (2016).
ADS  CAS  Google Scholar 

27.
Melin, A. D. et al. Calcium and carbon stable isotope ratios as paleodietary indicators. Am. J. Phys. Anthropol. 154, 633–643 (2014).
PubMed  Google Scholar 

28.
Chu, N. C., Henderson, G. M., Belshaw, N. S. & Hedges, R. E. Establishing the potential of Ca isotopes as proxy for consumption of dairy products. Appl. Geochem. 21, 1656–1667 (2006).
CAS  Google Scholar 

29.
Martin, J. E., Tacail, T., Cerling, T. E. & Balter, V. Calcium isotopes in enamel of modern and Plio-Pleistocene East African mammals. Earth Planet. Sci. Lett. 503, 227–235 (2018).
ADS  CAS  Google Scholar 

30.
Hassler, A. et al. Calcium isotopes offer clues on resource partitioning among Cretaceous predatory dinosaurs. Proc. Roy. Soc. B 285, 20180197 (2018).
Google Scholar 

31.
Heuser, A., Tütken, T., Gussone, N. & Galer, S. J. Calcium isotopes in fossil bones and teeth—Diagenetic versus biogenic origin. Geochim. Cosmochim. Acta. 75, 3419–3433 (2011).
ADS  CAS  Google Scholar 

32.
Tacail, T. et al. Assessing human weaning practices with calcium isotopes in tooth enamel. Proc. Natl. Acad. Sci. USA 114, 6268–6273 (2017).
CAS  PubMed  Google Scholar 

33.
Tacail, T. et al. Calcium isotopic patterns in enamel reflect different nursing behavior among South African early hominins. Sci. Adv. 5, eeax3250 (2019). https://doi.org/10.1126/sciadv.aax3250.
ADS  Google Scholar 

34.
Beynon, A. D., Dean, M. C. & Reid, D. J. Histological study on the chronology of the developing dentition in Gorilla and Orangutan. Am. J. Phys. Anthropol. 86, 189–203 (1991).
Google Scholar 

35.
Kelley, J. & Schwartz, G. T. Dental development and life history in living African and Asian apes. Proc. Natl. Acad. Sci. 107, 1035–1040 (2010).
ADS  CAS  PubMed  Google Scholar 

36.
Fletcher, A. & Nowell, A. The development of feeding behaviour in wild western lowland gorillas (Gorilla gorilla gorilla). Behaviour 145, 171–193 (2008).
Google Scholar 

37.
Prado-Martinez, J. et al. Great ape genetic diversity and population history. Nature 499, 471–475 (2014).
ADS  Google Scholar 

38.
Li, Q. et al. Spatially-resolved Ca isotopic and trace element variations in human deciduous teeth record diet and physiological change. Env. Archaeol. In Press, 1–10, https://doi.org/10.1080/14614103.2020.1758988 (2020).

39.
Scott, R. S., Teaford, M. F. & Ungar, P. S. Dental microwear texture and anthropoid diets. Am. J. Phys. Anthropol. 147, 551–579 (2012).
PubMed  Google Scholar 

40.
Martin, F. et al. Dietary niches of terrestrial cercopithecines from the Plio-Pleistocene Shungura Formation, Ethiopia: evidence from Dental Microwear Texture Analysis. Sci. Rep. 8, 14052 (2018).
ADS  PubMed  PubMed Central  Google Scholar 

41.
Fashing, P. J., Nguyen, N., Venkataraman, V. V. & Kerby, J. T. Gelada feeding ecology in an intact ecosystem at Guassa, Ethiopia: variability over time and implications for theropith and hominin dietary evolution. Am. J. Phys. Anthropol. 155, 1–16 (2014).
PubMed  Google Scholar 

42.
Souron, A. Morphology, diet, and stable carbon isotopes: on the diet of Theropithecus and some limits of uniformitarianism in paleoecology. Am. J. Phys. Anthropol. 166, 261–267 (2018).
PubMed  Google Scholar 

43.
Cerling, T. E., Chritz, K. L., Jablonski, N. G., Leakey, M. G. & Manthi, F. K. Diet of Theropithecus from 4 to 1 Ma in Kenya. Proc. Natl. Acad. Sci. USA 110, 10507–10512 (2013).
ADS  CAS  PubMed  Google Scholar 

44.
Grine, F. E., Sponheimer, M., Ungar, P. S., Lee‐Thorp, J. & Teaford, M. F. Dental microwear and stable isotopes inform the paleoecology of extinct hominins. Am. J. Phys. Anthropol. 148, 285–317 (2012).
PubMed  Google Scholar 

45.
Leakey, M. G. et al. New fossils from Koobi Fora in northern Kenya confirm taxonomic diversity in early Homo. Nature 488, 201–204 (2012).
ADS  CAS  PubMed  Google Scholar 

46.
Peters, C. R. & Vogel, J. C. Africa’s wild C4 plant foods and possible early hominid diets. J. Hum. Evol. 48, 219–236 (2005).
PubMed  Google Scholar 

47.
Schmitt, A.-D. 2016. Earth-Surface Ca Isotopic Fractionations. In Calcium Stable Isotope Geochemistry (eds Tipper, E. T., Schmitt, A. D. and Gussone, N.) 145–172. (Springer, Berlin Heidelberg, 2016).

48.
Yeakel, J. D., Bhat, U., Elliott Smith, E. A. & Newsome, S. Exploring the isotopic niche: isotopic variance, physiological incorporation, and the temporal dynamics of foraging. Front. Ecol. Evol. 4, https://doi.org/10.3389/fevo.2016.00001 (2016).

49.
Wood, B. & Strait, D. Patterns of resource use in early Homo and Paranthropus. J. Hum. Evol. 46, 119–162 (2004).
PubMed  Google Scholar 

50.
Wood, B. & Schroer, K. Paranthropus: where do things stand? in Human paleontology and prehistory (eds Marom, A. and Hovers, E.) 95–107 (Springer, Cham, 2017).

51.
Cerling, T. E. et al. Dietary changes of large herbivores in the Turkana Basin, Kenya from 4 to 1 million years ago. Proc. Natl. Acad. Sci. USA 112, 11467–11472 (2015).
ADS  CAS  PubMed  Google Scholar 

52.
Tacail, T., Albalat, E., Télouk, P. & Balter, V. A simplified protocol for measurement of Ca isotopes in biological samples. J. Anal. Atom. Spectrom. 29, 529–535 (2014).
CAS  Google Scholar 

53.
Martin, J. E. et al. Calcium isotopic evidence for vulnerable marine ecosystem structure prior to the K/Pg extinction. Curr. Biol. 27, 1641–1644 (2017).
CAS  PubMed  Google Scholar 

54.
Balter, V. et al. Calcium stable isotopes place Devonian conodonts as first level consumers. Geochem. Persp. Lett. 10, 36–39 (2019).
Google Scholar 

55.
Tacail, T., Télouk, P. & Balter, V. Precise analysis of calcium stable isotope variations in biological apatites using laser ablation MC-ICPMS. J. Anal. Atom. Spectrom. 31, 152–162 (2016).
CAS  Google Scholar 

56.
R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/ (2018).

57.
Heuser, A., Eisenhauer, A., Scholz-Ahrens, K. E. & Schrezenmeir, J. Biological fractionation of stable Ca isotopes in Göttingen minipigs as a physiological model for Ca homeostasis in humans. Isotopes Environ. Health Stud. 52, 633–648 (2016).
CAS  PubMed  Google Scholar 

58.
Heuser, A. & Eisenhauer, A. The calcium isotope composition (δ44/40Ca) of NIST SRM 915b and NIST SRM 1486. Geostand. Geoanal. Res. 32, 311–315 (2008).
CAS  Google Scholar 

59.
Heuser, A., Schmitt, A., Gussone, N. & Wombacher, F. Analytical Methods. In Calcium Stable Isotope Geochemistry. 23–73 (Springer, Berlin, Heidelberg, 2016).

60.
Tipper, E.T., Schmitt, A. & Gussone, N. Global Ca Cycles: Coupling of Continental and Oceanic Processes. In Calcium Stable Isotope Geochemistry. 173–222 (Springer, Berlin, Heidelberg, 2016). https://doi.org/10.1007/978-3-540-68953-9_6.

61.
Sponheimer, M. et al. Diets of southern African Bovidae: stable isotope evidence. J. Mammal. 84, 471–479 (2003).
Google Scholar  More