Ecology

1.Howell, K. L. et al. A decade to study deep-sea life. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-020-01352-5 (2020).Article 

Google Scholar 
2.Howell, K. L. et al. A blueprint for an inclusive, global deep-sea ocean decade field program. Front. Mar. Sci. 7, 1–25. https://doi.org/10.3389/fmars.2020.584861 (2020).ADS 
Article 

Google Scholar 
3.Ramirez-Llodra, E. et al. Man and the last great wilderness: Human impact on the deep sea. PLoS ONE 6, 22588. https://doi.org/10.1371/journal.pone.0022588 (2011).ADS 
CAS 
Article 

Google Scholar 
4.Bernardino, A. F., Levin, L. A., Thurber, A. R. & Smith, C. R. Comparative composition, diversity and trophic ecology of sediment macrofauna at vents, seeps and organic falls. PLoS ONE 7, e33515. https://doi.org/10.1371/journal.pone.0033515 (2012).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
5.Thiel, M. & Gutow, L. The ecology of rafting in the marine environment. II. The rafting organisms and community. Ocean. Mar. Biol. 43, 279–418. https://doi.org/10.1201/9781420037449.ch7 (2005).Article 

Google Scholar 
6.McClain, C. & Barry, J. Beta-diversity on deep-sea wood falls reflects gradients in energy availability. Biol. Lett. 10, 20140129. https://doi.org/10.1098/rsbl.2014.0129 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
7.Knudsen, J. The Bathyal and Abyssal Xylophaga (Pholadidae, Bivalvia) (Danish Science Press Ltd., 1961).
Google Scholar 
8.Turner, R. Wood-boring bivalves, opportunistic species in the deep sea. Science 180, 1377–1379. https://doi.org/10.1126/science.180.4093.1377 (1973).ADS 
CAS 
Article 
PubMed 

Google Scholar 
9.Voight, J. R. Deep-sea wood-boring bivalves of Xylophaga (Myoida: Pholadidae) on the continental shelf: A new species described. J. Mar. Biol. Assoc. UK 88, 1459–1464. https://doi.org/10.1017/S0025315408002117 (2008).Article 

Google Scholar 
10.Turner, R. D. A survey and Illustrated Catalogue of the Teredinidae (Mollusca: Bivalvia) (Harvard University, 1966).Book 

Google Scholar 
11.Hoppe, K. N. Teredo Navalis—the Cryptogenic Shipworm. in Invasive Aquatic Species of Europe. Distribution, Impacts and Management. (ed. Leppäkoski E., Gollasch S., O. S.) 116–119, https://doi.org/10.1007/978-94-015-9956-6_12 (2002).12.Distel, D. L. & Roberts, S. J. Bacterial endosymbionts in the gills of the deep-sea wood-boring bivalves Xylophaga atlantica and X. washingtona. Biol. Bull. 192, 253–261. https://doi.org/10.2307/1542719 (1997).CAS 
Article 
PubMed 

Google Scholar 
13.Distel, D. L., Morrill, W., MacLaren-Toussaint, N., Franks, D. & Waterbury, J. Teredinibacter turnerae gen. nov., sp. Nov., a dinitrogen-fixing, cellulolytic, endosymbiotic gamma-proteobacterium isolated from the gills of wood-boring molluscs (Bivalvia: Teredinidae). Int. J. Syst. Evol. Microbiol. 52, 2261–2269 (2002).CAS 
PubMed 

Google Scholar 
14.O’Connor, R. M. et al. Gill bacteria enable a novel digestive strategy in a wood-feeding mollusk. Proc. Natl. Acad. Sci. U. S. A. 111, 5096–5104. https://doi.org/10.1073/pnas.1413110111 (2014).CAS 
Article 

Google Scholar 
15.Sabbadin, F. et al. Uncovering the molecular mechanisms of lignocellulose digestion in shipworms. Biotechnol. Biofuels 11, 1–14. https://doi.org/10.1186/s13068-018-1058-3 (2018).CAS 
Article 

Google Scholar 
16.Kooijman, S. A. L. M. Dynamic Energy Budget Theory for Metabolic Organisation (Cambridge University Press, 2010).
Google Scholar 
17.Sarà, G., Palmeri, V., Montalto, V., Rinaldi, A. & Widdows, J. Parameterisation of bivalve functional traits for mechanistic eco-physiological dynamic energy budget (DEB) models. Mar. Ecol. Prog. Ser. 480, 99–117. https://doi.org/10.3354/meps10195 (2013).ADS 
Article 

Google Scholar 
18.Sarà, G., Rinaldi, A. & Montalto, V. Thinking beyond organism energy use: A trait-based bioenergetic mechanistic approach for predictions of life-history traits in marine organisms. Mar. Ecol. 35, 506–515. https://doi.org/10.1111/maec.12106 (2014).ADS 
Article 

Google Scholar 
19.Mangano, M. C. et al. Moving toward a strategy for addressing climate displacement of marine resources: A proof-of-concept. Front. Mar. Sci. 7, 1–16. https://doi.org/10.3389/fmars.2020.00408 (2020).ADS 
Article 

Google Scholar 
20.Romano, C. et al. Wooden stepping stones: Diversity and biogeography of deep-sea wood-boring Xylophagaidae (Mollusca: Bivalvia) in the North-East Atlantic Ocean, with the description of a new genus. Front. Mar. Sci. https://doi.org/10.3389/fmars.2020.579959 (2020).Article 

Google Scholar 
21.Culliney, J. L. & Turner, R. D. Larval development of the deep-water wood boring bivalve, Xylophaga atlantica Richards (Mollusca, bivalvia, pholadidae). Ophelia 15, 149–161. https://doi.org/10.1080/00785326.1976.10425455 (1976).Article 

Google Scholar 
22.Romey, W., Bullock, R. & Dealteris, J. Rapid growth of a deep-sea wood-boring bivalve. Cont. Shelf Res. 14, 1349–1359. https://doi.org/10.1016/0278-4343(94)90052-3 (1994).ADS 
Article 

Google Scholar 
23.Gaudron, S. M. et al. Colonization of organic substrates deployed in deep-sea reducing habitats by symbiotic species and associated fauna. Mar. Environ. Res. 70, 1–12. https://doi.org/10.1016/j.marenvres.2010.02.002 (2010).CAS 
Article 
PubMed 

Google Scholar 
24.Gaudron, S. M., Haga, T., Wang, H., Laming, S. R. & Duperron, S. Plasticity in reproduction and nutrition in wood-boring bivalves (Xylophaga atlantica) from the Mid-Atlantic Ridge. Mar. Biol. 163, 1–12. https://doi.org/10.1007/s00227-016-2988-6 (2016).CAS 
Article 

Google Scholar 
25.Childress, J. J., Cowles, D. L., Favuzzi, J. A. & Mickel, T. J. Metabolic rates of benthic deep-sea decapod crustaceans decline with increasing depth primarily due to the decline in temperature. Deep Sea Res. Part A Oceanogr. Res. Pap. 37, 929–949. https://doi.org/10.1016/0198-0149(90)90104-4 (1990).ADS 
CAS 
Article 

Google Scholar 
26.Childress, J. J. Are there physiological and biochemical adaptations of metabolism in deep-sea animals?. Trends Ecol. Evol. 10, 30–36. https://doi.org/10.1016/S0169-5347(00)88957-0 (1995).CAS 
Article 
PubMed 

Google Scholar 
27.Tittensor, D. P., Rex, M. A., Stuart, C. T., Mcclain, C. R. & Smith, C. R. Species—energy relationships in deep-sea molluscs subject collections species—energy relationships in deep-sea molluscs. Biol. Lett. 7, 718–722 (2011).Article 

Google Scholar 
28.McClain, C. R., Allen, A. P., Tittensor, D. P. & Rex, M. A. Energetics of life on the deep seafloor. Proc. Natl. Acad. Sci. U. S. A. 109, 15366–15371. https://doi.org/10.1073/pnas.1208976109 (2012).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
29.Mickel, T. J. & Childress, J. J. Effects of pressure and temperature on the EKG and heart rate of the hydrothermal vent crab Bythograea Thermydron (Brachyura). Biol. Bull. 162, 70–82. https://doi.org/10.2307/1540971 (1982).Article 

Google Scholar 
30.Voight, J. R., Cooper, J. C. & Lee, R. W. Stable isotopic evidence of mixotrophy in Xylophagaids, deep-sea wood-boring bivalves. Front. Mar. Sci. 7, 50. https://doi.org/10.3389/fmars.2020.00050 (2020).Article 

Google Scholar 
31.Lika, K. et al. The ‘covariation method’ for estimating the parameters of the standard dynamic energy budget model I: Philosophy and approach. J. Sea Res. 66, 270–277. https://doi.org/10.1016/j.seares.2011.07.010 (2011).ADS 
Article 

Google Scholar 
32.Marques, G. M. et al. The AmP project: Comparing species on the basis of dynamic energy budget parameters. PLoS Comput. Biol. 14, 1–23. https://doi.org/10.1371/journal.pcbi.1006100 (2018).CAS 
Article 

Google Scholar 
33.Mariño, J., Augustine, S., Dufour, S. C. & Hurford, A. Dynamic Energy Budget theory predicts smaller energy reserves in thyasirid bivalves that harbour symbionts. J. Sea Res. 143, 119–127. https://doi.org/10.1016/j.seares.2018.07.015 (2019).ADS 
Article 

Google Scholar 
34.Brown, A. et al. Metabolic costs imposed by hydrostatic pressure constrain bathymetric range in the lithodid crab Lithodes maja. J. Exp. Biol. 220, 3916–3926. https://doi.org/10.1242/jeb.158543 (2017).Article 
PubMed 

Google Scholar 
35.Eisenmenger, M. J. & Reyes-De-Corcuera, J. I. High pressure enhancement of enzymes: A review. Enzyme Microb. Technol. 45, 331–347. https://doi.org/10.1016/j.enzmictec.2009.08.001 (2009).CAS 
Article 

Google Scholar 
36.Kalenitchenko, D. et al. Bacteria alone establish the chemical basis of the wood-fall chemosynthetic ecosystem in the deep-sea. ISME J. 12, 367–379. https://doi.org/10.1038/ismej.2017.163 (2018).CAS 
Article 
PubMed 

Google Scholar 
37.Levesque, C., Limén, H. & Juniper, S. K. Origin, composition and nutritional quality of particulate matter at deep-sea hydrothermal vents on Axial Volcano NE pacific. Mar. Ecol. Prog. Ser. 289, 43–52. https://doi.org/10.3354/meps289043 (2005).ADS 
Article 

Google Scholar 
38.Limén, H., Levesque, C. & Kim Juniper, S. POM in macro-/meiofaunal food webs associated with three flow regimes at deep-sea hydrothermal vents on Axial Volcano, Juan de Fuca Ridge. Mar. Biol. 153, 129–139. https://doi.org/10.1007/s00227-007-0790-1 (2007).Article 

Google Scholar 
39.Culliney, J. L. Comparative larval development of the shipworms Bankia gouldi and Teredo navalis. Mar. Biol. 29, 245–251. https://doi.org/10.1007/BF00391850 (1975).Article 

Google Scholar 
40.Ramirez Llodra, E. Fecundity and life-history strategies in marine invertebrates. Adv. Mar. Biol. 43, 87–170. https://doi.org/10.1016/S0065-2881(02)43004-0 (2002).Article 
PubMed 

Google Scholar 
41.Fernandez-Arcaya, U. et al. Bathymetric gradients of fecundity and egg size in fishes: A Mediterranean case study. Deep Sea Res. Part A Oceanogr. Res. Pap. 116, 106–117. https://doi.org/10.1016/j.enzmictec.2009.08.001 (2016).ADS 
CAS 
Article 

Google Scholar 
42.Young, C. M., Emson, R. H., Rice, M. E. & Tyler, P. A. A paradoxical mismatch of fecundity and recruitment in deep-sea opportunists: cocculinid and pseudococculinid limpets colonizing vascular plant remains on the Bahamian Slope. Deep Sea Res. 92, 36–45. https://doi.org/10.1016/j.dsr2.2013.01.027 (2013).ADS 
Article 

Google Scholar 
43.Thorson, G. Reproductive and larval ecology of marine bottom invertebrates. Biol. Rev. 25, 1–45. https://doi.org/10.1111/j.1469-185X.1950.tb00585.x (1950).CAS 
Article 
PubMed 

Google Scholar 
44.Hitt, N. T. et al. Growth and longevity of New Zealand black corals. Deep. Res. Part I Oceanogr. Res. Pap. 162, e103298. https://doi.org/10.1016/j.dsr.2020.103298 (2020).Article 

Google Scholar 
45.McNichol, J. et al. Primary productivity below the seafloor at deep-sea hot springs. Proc. Natl. Acad. Sci. U. S. A. 115, 6756–6761. https://doi.org/10.1073/pnas.1804351115 (2018).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
46.Levin, L. A. et al. Hydrothermal vents and methane seeps: Rethinking the sphere of influence. Front. Mar. Sci. 3, 1–23. https://doi.org/10.3389/fmars.2016.00072 (2016).ADS 
Article 

Google Scholar 
47.Nedoncelle, K., Lartaud, F., de Rafelis, M., Boulila, S. & Le Bris, N. A new method for high-resolution bivalve growth rate studies in hydrothermal environments. Mar. Biol. 160, 1427–1439. https://doi.org/10.1007/s00227-013-2195-7 (2013).CAS 
Article 

Google Scholar 
48.Turekian, K. K., Cochran, J. K. & Bennett, J. T. Growth rate of a vesicomyid clam from the 21° N East Pacific Rise hydrothermal area. Nature 303, 55–56. https://doi.org/10.1038/303055a0 (1983).ADS 
Article 

Google Scholar 
49.Lutz, R. A. et al. Rapid growth at deep-sea vents. Nature 371, 663–664. https://doi.org/10.1038/371663a0 (1994).ADS 
Article 

Google Scholar 
50.Reed, A. J., Morris, J. P., Linse, K. & Thatje, S. Plasticity in shell morphology and growth among deep-sea protobranch bivalves of the genus Yoldiella (Yoldiidae) from contrasting Southern ocean regions. Deep. Res. Part I Oceanogr. Res. Pap. 81, 14–24. https://doi.org/10.1016/j.dsr.2013.07.006 (2013).ADS 
Article 

Google Scholar 
51.Oliver, G., Allen, J. A. & Yonge, M. The functional and adaptive morphology of the deep-sea species of the Arcacea (Mollusca: Bivalvia) from the Atlantic. Philos. Trans. R. Soc. London. B Biol. Sci. 291, 45–76. https://doi.org/10.1098/rstb.1980.0127 (1980).ADS 
Article 

Google Scholar 
52.Romano, C., Voight, J. R., Pérez-Portela, R. & Martin, D. Morphological and genetic diversity of the wood-boring Xylophaga (Mollusca, Bivalvia): New species and records from deep-sea Iberian canyons. PLoS ONE 9, 102887. https://doi.org/10.1371/journal.pone.0102887 (2014).ADS 
CAS 
Article 

Google Scholar 
53.Saulsbury, J. et al. Evaluating the influences of temperature, primary production, and evolutionary history on bivalve growth rates. Paleobiology 45, 405–420. https://doi.org/10.1017/pab.2019.20 (2019).Article 

Google Scholar 
54.Moss, D. K. et al. Lifespan, growth rate, and body size across latitude in marine bivalvia, with implications for phanerozoic evolution. Proc. R. Soc. B Biol. Sci. 283, 20161364. https://doi.org/10.1098/rspb.2016.1364 (2016).Article 

Google Scholar 
55.Tyler, P. A., Young, C. M. & Dove, F. Settlement, growth and reproduction in the deep-sea wood-boring bivalve mollusc Xylophaga depalmai. Mar. Ecol. Prog. Ser. 343, 151–159. https://doi.org/10.3354/meps06832 (2007).ADS 
Article 

Google Scholar 
56.Brown, J. H., Gillooly, J. F., Allen, A. P., Savage, V. M. & West, G. B. Toward a metabolic theory of ecology. Ecology 85, 1771–1789. https://doi.org/10.1890/03-9000 (2004).Article 

Google Scholar 
57.Maino, J. L., Kearney, M. R., Nisbet, R. M. & Kooijman, S. A. L. M. Reconciling theories for metabolic scaling. J. Anim. Ecol. 83, 20–29. https://doi.org/10.1111/1365-2656.12085 (2014).Article 
PubMed 

Google Scholar 
58.Gaudron, S. M., Demoyencourt, E. & Duperron, S. Reproductive traits of the cold-seep symbiotic mussel Idas modiolaeformis: gametogenesis and larval biology. Biol. Bull. 222, 6–16. https://doi.org/10.1086/bblv222n1p6 (2012).Article 

Google Scholar 
59.Hilário, A. et al. Estimating dispersal distance in the deep sea: Challenges and applications to marine reserves. Front. Mar. Sci. 2, 6. https://doi.org/10.3389/fmars.2015.00006 (2015).ADS 
Article 

Google Scholar 
60.Marsh, A. G., Mullineaux, L. S., Young, C. M. & Manahan, D. T. Larval dispersal potential of the tubeworm Riftia pachyptila at deep-sea hydrothermal vents. Nature 411, 77–80. https://doi.org/10.1038/35075063 (2001).ADS 
CAS 
Article 
PubMed 

Google Scholar 
61.Young, C. M. et al. Dispersal of deep-sea larvae from the intra-American seas: Simulations of trajectories using ocean models. Integr. Comp. Biol. 52, 483–496. https://doi.org/10.1093/icb/ics090 (2012).Article 
PubMed 

Google Scholar 
62.Yearsley, J. M., Salmanidou, D. M., Carlsson, J., Burns, D. & Van Dover, C. L. Biophysical models of persistent connectivity and barriers on the northern Mid-Atlantic Ridge. Deep. Res. Part II Top. Stud. Oceanogr. 180, 104819. https://doi.org/10.1016/j.dsr2.2020.104819 (2020).Article 

Google Scholar 
63.Levin, L. A. et al. Global observing needs in the deep ocean. Front. Mar. Sci. 6, 1–32. https://doi.org/10.3389/fmars.2019.00241 (2019).ADS 
Article 

Google Scholar 
64.McClain, C. R., Boyer, A. G. & Rosenberg, G. The island rule and the evolution of body size in the deep sea. J. Biogeogr. 33, 1578–1584. https://doi.org/10.1111/j.1365-2699.2006.01545.x (2006).Article 

Google Scholar 
65.Zonneveld, C. & Kooijman, S. A. L. M. Application of a dynamic energy budget model to Lymnaea stagnalis (L.). Funct. Ecol. 3, 269–278. https://doi.org/10.2307/2389365 (1989).Article 

Google Scholar 
66.Mueller, C. A., Augustine, S., Kooijman, S. A. L. M., Kearney, M. R. & Seymour, R. S. The trade-off between maturation and growth during accelerated development in frogs. Comp. Biochem. Physiol. A 163, 95–102. https://doi.org/10.1016/j.cbpa.2012.05.190 (2012).CAS 
Article 

Google Scholar 
67.MacArthur, R.H. & Wilson, E. The Theory of Island Biogeography (1967).68.Kooijman, S. A. L. M. Metabolic acceleration in animal ontogeny: An evolutionary perspective. J. Sea Res. 94, 128–137. https://doi.org/10.1016/j.seares.2014.06.005 (2014).ADS 
Article 

Google Scholar  More

The crystalline nature and phase purity of the synthesized NiO, PdO and PdO–NiO nanocomposite were confirmed by XRD analysis and shown in Fig. 1a. All three XRD patterns show the highly crystalline nature, which confirms the purity of the samples. In PdO–NiO mixed metal oxide, The intense diffraction peaks at 2θ = 37.16, 43.24, 62.81, 75.32 and 79.34° indexed to (101), (012), (110), (113) and (202) planes of the cubic phase of NiO, which is highly consistent with standard JCPDS NO: 01-071-475120. On the other hand, the sharp peaks at 2θ = 34.54 and 55.79° indexed to (101) and (112) planes corresponding to the tetragonal crystalline phase of PdO with an average lattice parameter 3.043 Å, which is highly associated with standard JCPDS NO: 043-102421. Interestingly, metallic diffraction peaks and other impurity phases were not observed in the hybrid PdO–NiO nanocomposite. This confirming that metal source (Pd) completely oxidized to the metal oxide (PdO) and formed hybrid PdO–NiO nanocomposite. In addition, the single metal oxide (PdO, NiO) diffraction patterns were compared in Fig. 1a.Figure 1(a) X-ray diffraction pattern and (b) FTIR spectra of NiO, PdO and PdO–NiO (c) TGA curve of PdO–NiO nanocomposite.Full size imageThe sharp peaks of each metal oxide confirm the high crystallinity of the tetragonal phase of the PdO and cubic phase of NiO. Both diffraction pattern of the single metal oxides is well associated with standard JCPDS file: 01-071-4751 and 043-1024 respectively. The average crystallite size, micro-strains of the NiO, PdO and PdO doped NiO were calculated by using Scherrer analysis and W–H analysis and plotted in Fig. 1b. A small variation was observed in the crystallite size of the catalysts (PdO, NiO and PdO–NiO) in Fig. 1b, which is due to the difference in the distribution of the crystal in the catalysts. The average crystallite size of the PdO, NiO and PdO–NiO was found to be 10.8, 7.8 and 7.36 nm, respectively. The PdO doping reduces the crystallite size in the PdO–NiO composite, consistent with previous reports. The calculated crystallite sizes, d-spacing, micro-strains, and binding energies of the PdO, NiO, and PdO–NiO are shown in the Table 1. In addition, the experimentally calculated d-spacing value of the pure PdO and NiO was well correlated with theoretical values and shown in Table SI. 1. The synthesized catalyst showed the negative and positive slopes of ε are corresponds to the compressive and tensile stress, respectively.Table 1 Crystallite sizes, d-spacing, micro-strains, and binding energies of the PdO, NiO, and PdO–NiO.Full size tableFurthermore, the chemical bonding and functional groups were analyzed by FTIR spectrometer. Figure 1c show FTIR spectra of the synthesized NiO, PdO and PdO–NiO nanocomposite. Similar spectra were observed for the three catalyst, the absorption peaks at high-frequency region 3200–3400 cm−1 belongs to O–H stretching vibration of the water molecules, due to surface adsorption phenomenon. Furthermore, three absorption peaks appeared at 1398.2, 1237.8 and 1057.8 cm−1, which ascribed to the C–O, CH2 and C=O stretching vibrations, which is well associated with XPS analysis data. The metal oxide bonding peaks appeared in the frequency range of 480.2–702.6 cm−1. Hence the Pd–O and Ni–O stretching frequency in the PdO doped NiO sample confirmed the formation of hybrid PdO–NiO nanocomposite22. After that, the thermal stability of the PdO–NiO nanocomposite was studied by TGA analysis. Figure 1d shows the thermogram of the PdO–NiO nanocomposite. The first weight loss (7%) started in the temperature range of 65 to 180 °C. Due to the H2O molecules, evaporation and then the sustainable weight loss of around 10% was observed in the range of 187 to 574 °C. Beyond 600 °C a significant weight loss 20% was observed, which may be assigned to the unreacted CO3 combustion23. Hence, the XRD, FTIR and TGA spectral studies confirmed the formation of hybrid PdO–NiO nanocomposite.The electronic state and chemical bonding of the PdO doped NiO composite was analyzed by using XPS spectra. Figure 2 shows the XPS spectra of the PdO–NiO nanocomposite, the broad scan spectrum (Fig. SI. 1) of the PdO–NiO, which show the existence of the Pd (3d), Ni (2p), O (1s) and C(1s) elements. The deconvoluted Pd 3d XPS spectrum in Fig. 2a shows two major peaks at a binding energy of 336.4 and 342.2 eV corresponds to spin–orbit doublets of Pd 3d5/2 and Pd 3d3/2, respectively, which confirmed Pd2+ ions in the form of PdO in the PdO–NiO nanocomposite24,25. In addition, the satellite peak of Pd species appeared at a binding energy of 339.2 eV and 345.3 eV. In Ni 2p spectra (Fig. 2b), Ni 2p3/2 and Ni 2p3/2 spin–orbit doublets peaks were observed at 854.9 eV and 872.5 eV, which corresponded to Ni–O and Ni–OH, respectively and their corresponding satellite peaks are located at a binding energy of 867.2 eV and 879.1 eV26,27. Furthermore, O 1s spectra (Fig. 2c) show the two peaks at a binding energy of 529.3 eV and 534.2 eV, which ascribed M–O and M–OH species. The obtained XPS spectra of the PdO–NiO nanocomposite are well associated with the XRD and EDX analysis.Figure 2High resolution X-ray photoelectron spectroscopy (a) Pd 3d, (b) Ni 2p, (c) O 1s and (d) C 1s spectra of PdO–NiO nanocomposite.Full size imageThe morphology feature and elemental composition of the prepared NiO, PdO and PdO–NiO nanocomposite was scrutinised by FE-SEM. Figure 3 shows the SEM morphology images of NiO, PdO and PdO–NiO nanocomposite at low and high magnification. Pure metal oxides NiO and PdO samples in Fig. 3a–d shows the porous crystalline morphology with high purity of the respective elements. On the other hand, PdO doped NiO sample in Fig. 3e,f shows the uniform, monodisperse, spherical crystalline morphology. Which confirms that PdO uniformly distributed with NiO. Hence, the PdO doping enhances the surface area of the catalyst. In addition, the elemental composition and elemental mapping was analysed for PdO–NiO sample and shown in Fig. 3g,h. The EDX spectra and elemental mapping clearly confirms the presence of the Pd, Ni and O elements in the composite with high purity.Figure 3FESEM images (a,b) NiO (c,d) PdO, (e,f) PdO–NiO and (g,f) EDS spectra and elemental mapping of PdO–NiO nanocomposite.Full size imageThe detailed morphology and particle size distribution of the PdO–NiO NPs was measured by HR-TEM and the results are presented in Fig. 4. Figure 4a–c shows the typical HRTEM images of the as-synthesized PdO–NiO nanocomposite. The obtained TEM images confirmed the uniform distribution of the spherical PdO–NiO NPs, which agrees with FESEM results. In Fig. 4c (inset), the histogram reveals that formed PdO–NiO nanoparticles are uniformly distributed with an average particle size of about 9.64 ± 2.1 nm, which is well associated with XRD crystallite size. Furthermore, the SAED pattern was analyzed to understand the crystallinity and the crystal quality of the PdO–NiO nanoparticles are shown in Fig. 4d. The clear ring-like structure suggests the polycrystalline nature of PdO–NiO. The obtained diffraction rings d-spacing values are corresponding to the (101), (012), (110), (113) and (202) planes of the NiO nanoparticles. Figure 4e shows the lattice fringes of the PdO doped NiO nanoparticles. The fringes show the lattice planes for both metal oxides. The interplanar d-spacing value of 0.1997 nm to correspond to the (012) plane of the NiO phase and the d-spacing value of 0. 2145 nm to correspond to the (110) plane of the PdO in the composite. Which is well correlated with the XRD d-spacing values. Elemental mapping in Fig. 4f shows the presence of Ni, Pd and O elements with uniform distribution as similar as SEM mapping. The morphology results of the synthesized catalysts are well associated with XRD, FTIR and XPS analysis.Figure 4High resolution transmission electron microscopy HRTEM images (a–c), (d) SAED pattern, (e) interplanr spacing and (f) elemental mapping of PdO–NiO nanocomposite.Full size imageUV–Vis absorption spectra were analyzed for the as-synthesized catalysts NiO, PdO and PdO–NiO and the respective results are presented in Fig. 5a. The absorption spectra show the strongest absorption maxima at 234.8 nm for all three catalysts. In addition, the characteristic absorption band of NiO and PdO were observed at 338.2 nm and 422.1 nm respectively23, on the other hand, no characteristic absorption band was observed for PdO–NiO sample. Furthermore, the bandgap energy was calculated for three catalysts by using the Schuster-Kubelka–Munk function.$$(alpha {text{h}}nu ),{text{n}} = {text{A}}({text{h}}nu – {text{Eg}})$$
(1)
Figure 5(a) UV–Vis absorption spectra and (b) Plot of (αhν)2 Vs hν for NiO, PdO and PdO–NiO nanocomposite (Inset shows the bandgap energy of the catalyst).Full size imageThe bandgap energy (Eg) was achieved by extrapolating against the photon energy and the obtained results are shown in Fig. 5b. The calculated bandgap (Eg) of NiO, PdO and PdO–NiO are 4.05 eV, 3.84 eV and 4.24 eV, respectively25 (Fig. 5b, inset). The PdO doping with NiO increases its bandgap value. This suggests that the PdO interface and NiO interface are closely combined in the composite. The obtained band gap value of the catalysts is much higher than the reported bandgap energy. The bandgap energy is highly dependent on the particle’s size. The bandgap energy increases with decreasing particle size, which confirmed that the synthesized catalysts are in nanoscale. The bandgap energy (Eg) of the PdO–NiO catalyst is well associated with FESEM and TEM results.The photoluminescence (Pl) spectra of NiO, PdO and PdO–NiO materials were measured at 325 nm excitation wavelength and presented in Fig. 6. Figure 6a shows the Pl spectra of the pure NiO, PdO and PdO–NiO nanocomposite. The blue/violet emission was observed for all three samples at 364 nm due to the excitation of 3d8 electrons of Ni2+ ions from the conduction band to the valence band24. From Fig. 6a, it can be seen that the intensity of the PdO–NiO nanocomposite is lower than the pure NiO and PdO, which indicated the higher electron transfer between the NiO and PdO, which is well correlated with electrochemical results. The deconvoluted PL spectra of NiO, PdO and PdO–NiO materials are shown in Fig. 6b–d; four peaks have been fitted for each sample as shown in Fig. 6b–d. The UV emission at 364 nm (3.4 eV) corresponds to the near band edge (NBE) excitation of NiO25. The obtained PL spectra confirmed that the PdO–NiO nanocomposite has more conductivity than the pure metal oxides.Figure 6(a) PL spectra and (b–d) deconvoluted spectra of NiO, PdO and PdO–NiO nanocomposite respectively.Full size imageThe electrode kinetics of the NiO, PdO and PdO–NiO modified GC electrode was explored in 1 M KOH at different scan rate variation at room temperature. In addition, the resistance of the aforesaid electrodes was monitored in impedance analysis and shown in Fig. 7. Figure 7a, the NiO/GC show a pair of well-defined redox peaks at around 0.49 V and 0.44 V respectively, which corresponding to the reversible reaction between Ni2+ and Ni3+26. In addition, the redox peak currents linearly increase with increasing scan rate.Figure 7(a–c) Cyclic voltammetry and (d) EIS curves of NiO, PdO and PdO–NiO nanocomposite in 1 M KOH electrolyte solution.Full size imageFigure 7b shows CV pattern of PdO/GC electrode in 1 M KOH solution, which sows poor peaks palladium oxide and palladium reduction peak at 0.52 V and 0.53 V respectively due to the formation of oxyhydroxide on the electrode surface in basic medium. Whereas the PdO–NiO NPs modified GC electrode in Fig. 7c show well-defined Ni2+ and Ni3+ kinetics with eightfold higher peak current. Which confirms the Efficient electron transfer between NiO and PdO in the composite. Which is well associated with PL results. Furthermore, the impedance spectra were achieved for three electrodes at fixed over potential (500 mV s−1) in 1 M KOH electrolyte and presented in Fig. 7d. In Fig. 7d the Rct value of the pure metal oxides NiO/GC and PdO/GC electrode were obtained as 236 Ω and 2702 Ω respectively. Whereas the bimetal oxide PdO–NiO/GC show 425.7 Ω, which is lower than the pure PdO. Due to the superior electron transformation between each metal oxide.Catalytic reduction of Azo compoundsAzo compounds are highly toxic to the environment as well as human beings. Especially, nitrophenols are listed as the topmost hazardous chemical in the world. Hence the reduction of nitrophenols gains the most attention. Generally, the nitrophenol reduction reaction is thermodynamically favourable (E0 = − 0.76 V) at optimized conditions, whereas the NaBH4 acts as a reducing agent (E0 = − 1.33 V)27,28. However, the reduction rate is prolonged without the catalyst due to the kinetic barrier between the reducing agent and reactant. Hence, the catalytic reduction nitrophenols are a good way to convert to non-toxic aminophenol (AP) with the presence of NaBH4 as a reducing agent. The reduction reaction was easily monitored with a UV–Vis spectrometer. It is known that NaBH4 alone cannot reduce the nitrophenols into aminophenol’s. As shown in Fig. SI. 2 the fresh nitrophenols (NP, DNP and TNP) absorption peak appeared at 300–370 nm respectively. When the addition of reducing agent, the peak was shifted to 402–450 nm28. In addition, the solution color was turned light yellow to deep yellow, due to the formation of corresponding nitrophenolate ions in basic solution. However, no reduction was achieved over 2 h, indicating that nitrophenolate ions were very stable with NaBH4. Furthermore, the catalytic activity of the NiO was explored with three nitrophenols and shown in Fig. SI. 3. The pure NiO exhibits a poor catalytic reduction of nitro compounds. In contrast, the PdO–NiO catalyst show excellent activity on the nitrophenols, as shown in Fig. 8.Figure 8Catalytic activity and kinetic rate of PdO–NiO nanocomposite on reduction of nitrophenols with NaBH4 solution (a,b) NP, (c,d) DNP and (e,f) TNP.Full size imageIt can be seen that the nitrophenolate peak absorbance at 400 nm gradually decreases with reaction time, which confirmed that PdO–NiO promotes the electron and hydrogen transfer between the reactant. Due to the higher active sites of the PdO–NiO. The present PdO–NiO catalyst completes the reduction reaction of NP, DNP and TNP within 10, 13, 25 min respectively. In addition, the aminophenol absorption peak appeared around 300 nm for all three nitrophenols. On the other hand, the deep yellow solution turned colorless, indicating the formation of aminophenol. The rate constant κapp for each nitrophenol was calculated from the plot of ln (At/Ao) Vs. time. The proposed PdO–NiO catalyst exhibits excellent rate constant 0.1667, 0.0997, 0.0686 min−1 for NP, DNP and TNT, respectively, which is the higher rate constant than the previously reported catalyst (Table SI. 2). Generally, the reduction mechanism of nitrophenols to aminophenols follows many intermediate steps from nitro to nitroso and then to hydroxylamine and to final aminophenol. For these reaction required both electron transfer and active hydrogen atoms. Here, the BH4− ions produce the active hydrogen atoms on the surface of the catalyst and subsequently, the PdO–NiO catalyst enhances the electron transfer. As a result, the reduction of NP could be efficiently accelerated by the PdO–NiO catalyst. Furthermore, comparison of the catalytic reduction performance of nitrophenols with varies catalyst are shown in Table SI. 2. The reduction mechanism of the nitrophenols with PdO–NiO catalyst with NaBH4 as shown in Fig. SI. 6.Furthermore, the catalytic activity of PdO–NiO composite was explored by the reduction of organic azo dye compounds such as Methylene blue (MB), Rhodamine B (RhB) and Methyl orange (MO) with the addition of NaBH4 in the presence of PdO–NiO catalyst29. The reduction of each dye was monitored at different absorption peaks, as shown in Fig. 9. The intensity of each dye at respective wavelengths linearly reduced with time in the presence of the PdO–NiO. In addition, the reduction rate κapp for each dye was calculated from the plot of ln (At/Ao) Vs. time. PdO–NiO catalyst exhibited excellent reduction rate as 0.099, 0.0416 and 0.0896 min−1 for MB, RhB and MO, respectively, superior catalytic activity than previous reports. In addition, the azo dye solution turned into colourless, indicating the complete reduction occurs in the presence of PdO–NiO. The pure NiO catalyst exhibits poor reduction activity towards azo dyes with the presence of NaBH4 (Fig. SI. 4). Furthermore, comparison of the catalytic reduction performance azo dyes with varies catalyst are shown in Table. SI. 3. Additionally, the reduction of the mixture of nitrophenols (NP, DNP and TNP) and azo dyes (MB, RhB and MO) was tested with PdO–NiO nanocomposite and obtained results are shown in Fig. 10. Initially, the mixture of azo compounds is formed dark solution then rapidly turned into a colourless and became a transparent solution with the addition of PdO–NiO in the presence of NaBH430. The complete azo compounds reduction was achieved within 8 min. Hence, the proposed PdO–NiO is a promising catalyst for wastewater treatment. In addition, the effect of the catalyst loading on the reduction of mixture of azo compounds were studied with different loading amount of PdO–NiO catalyst (3–10 mg) and shown in the Table. SI. 4. Which show that the reduction rate increased with loading amount of the catalyst. In addition, the catalytic reduction performance of toxic azo compounds by various catalysts are shown in Table. 2. PdO–NiO catalyst exhibit superior reduction performance than the previously reported catalyst.Figure 9Catalytic activity and kinetic rate of PdO–NiO nanocomposite on reduction of azo dyes with NaBH4 solution (a,b) MB, (c,d) RhB and (e,f) MO.Full size imageFigure 10Reaction progress of an azo compounds mixture (4-NP, 2,4-DNP, 2,4,6-TNP, MB, RhB, MO) with PdO–NiO and NaBH4. Conditions: Dye: 100 ppm, 25 ml, Nitrophenol: 0.12 mM, 25 ml, NaBH4: 0.1 M, 5 ml and catalyst: 3 mg.Full size imageTable 2 Comparison of catalytic reduction performance of toxic azo compounds by various catalysts.Full size tableFurthermore, the reduction mechanism of azo dyes over PdO–NiO catalyst with reducing agent shown in Fig. SI. 7.After the complete reduction reaction, the catalyst property was analyzed to understand the stability of PdO–NiO. In Fig. 11a, FTIR spectra showed no noticeable changes before and after catalytic reduction of mixture reduction. Additionally, the SEM image (Fig. 11b) also showed no changes in the morphology of the PdO–NiO. In addition, HR-TEM image (Fig. 11c) was also analyzed to study the change in the particles size after catalytic reduction, show no considerable change in the particle size. Which proved that PdO–NiO is highly stable in the reduction conditions.Figure 11(a) FTIR spectra of PdO–NiO nanocomposite before and after reducing the mixture of azo compounds (b) FE-SEM image (c) HR-TEM image of PdO–NiO after reducing the mixture of azo compounds.Full size image More

Surface processes model and forcing mechanismsLandscape evolution over the last one million years interval is performed with the open-source modelling code Badlands34. It simulates the evolution of topography induced by sediment erosion, transport, and deposition (Fig. 1a). Amongst the different capabilities available in Badlands, we applied the fluvial incision and hillslope processes, which are described by geomorphic equations and explicitly solved using a finite volume discretisation. In this study, soil properties are assumed to be spatially and temporally uniform over the region, and we do not differentiate between regolith and bedrock. It is worth noting that the role of flexural responses induced by erosion and deposition is also not accounted for. Under these assumptions, the continuity of mass is governed by vertical land motion (U, uplift or subsidence in m/yr), long-term diffusive processes and detachment-limited fluvial runoff-based stream power law:$$frac{partial z}{partial t}=U+kappa {nabla }^{2}z+epsilon {(PA)}^{m}nabla {z}^{n}$$
(1)
where z is the surface elevation (m), t is the time (yr), κ is the diffusion coefficient for soil creep34 with different values for terrestrial and marine environments, ϵ is a dimensional coefficient of erodibility of the channel bed, m and n are dimensionless empirically constants, that are set to 0.5 and 1, respectively, and PA is a proxy for water discharge that numerically integrates the total area (A) and precipitation (P) from upstream regions34.Both κ and ϵ depend on lithology, precipitation, and channel hydraulics and are scale dependent34. All our landscape evolution simulations are running over a triangular irregular network of ~18. e6 km2 with a resolution of ~5 km, and outputs are saved every 1000 yr.The detachment-limited fluvial runoff-based stream power law is computed with a ({{{{{{{mathcal{O}}}}}}}}(n))-efficient ordering approach54 based on a single-flow-direction approximation where water is routed down the path of the steepest descent. The flow routing algorithm and associated sediment transport from source to sink depend on surface morphology, and sediment deposition occurs under three circumstances: (1) existence of depressions or endorheic basins, (2) if local slope is less than the aggregational slope in land areas and (3) when sediments enter the marine realm34. Submerged sediments are then transported by diffusion processes defined with a constant marine diffusion coefficient34.All landscape simulations are constrained with different forcing mechanisms, and five scenarios were tested (Supplementary Table 2).First, we impose precipitation estimates from the PaleoClim database38,39,40. These estimates are products from paleoclimate simulations (coupled atmosphere-ocean general circulation model) downscaled at approximately the same resolution as our landscape model (~5 km at the equator). Annual averages of precipitation rates are then used to provide rainfall trends in our simulations based on the ten specific snapshots available (from the mid-Pliocene warm period to late Holocene and present day). Between two consecutive snapshots, we assume that precipitation remains constant for the considered time interval. For exposed regions that are considered flooded in the PaleoClim database, we define offshore precipitations using a nearest neighbour algorithm where closest precipitation estimates are averaged from PaleoClim inland regions. To evaluate the role of precipitation variability on landscape dynamics, we also run a uniform rainfall scenario (2 m/yr obtained by averaging the annual precipitation rates from the PaleoClim database).Secondly, the models are forced with sea level fluctuations known to play a major role in the flooding history of the Sunda Shelf11,13,53. Two sea level curves are tested (Supplementary Fig. 1d). To account for the inherent uncertainty in reconstructed sea level variations, we chose a first curve37 obtained from a sea level stack constructed from five to seven individual reconstructions that agrees with isostatically adjusted coral-based sea level estimates at both 125 and 400 ka. The second one is taken from the global sea level curve reconstruction36 based on benthic oxygen isotope data and has been recently used to reconstruct the subsidence history of Sundaland17,18.The last forcing considered in our study is the tectonic regime. First, we chose to explore a non-tectonic model based on the default assumption of stability for the shelf17. Secondly, we assumed a uniform subsidence rate of −0.25 mm/yr recently derived from a combination of geomorphological observations, coral reef growth numerical simulations and shallow seismic stratigraphy interpretations17. Then, to represent the regional variations in the tectonic regime, we have compiled and digitised a number of calibration points (Supplementary Fig. 1b and Table 1) that were used to produce a subsidence and uplift map by geo-referencing calibration points and available tectonic polygons, and by Gaussian-smoothing and normalising the uplift and subsidence rates between the calibrated range to avoid sharp transitions in regions without observations. The resulting map does not account for fine spatial scale tectonic features such as fault systems43,55 or orogenic and sedimentary related isostatic responses. It rather represents a regional vertical tectonic trend with an overall uplift of Wallacea and NW Borneo regions and long-wavelength subsidence of Sunda Shelf and Singapore Strait17.Landscape evolution model calibrationThe landscape models start during the Calabrian in the Pleistocene Epoch, one million years before the present. At each time interval, the landscape evolves following Eq. (1) and the surface adjusts under the action of rivers and soil creep (Fig. 1a). In addition to surface changes, we extract morphometric characteristics such as drainage basins extents, river profiles lengths (Fig. 3 and Supplementary Fig. 2), distance between main rivers outlet (Supplementary Fig. 3) and tracks the cumulative erosion and deposition over time (Fig. 1b and Supplementary Fig. 1d).For model calibration, we perform a series of steps consisting in adjusting the initial elevation and the erosion–deposition parameters (i.e., κ and ϵ in Eq. (1)) to match with different observations.The initial paleo-surface is obtained by applying the uplift and subsidence rates backwards to calculate the total change in topography for the 1 Myr interval. Then, we test the simulated paleo-river drainages against results from a combination of phylogenic studies9,13 and paleo-river channels and valleys found from seismic and well surveys41,42,44. Iteratively, we modify our paleo-elevation to ensure those main river basins (e.g., Johore, Siam, Mekong, East Sunda) encapsulate the paleo-drainage maps reconstructed using lowland freshwater taxa described in13 (Supplementary Fig. 1a and Table 4) and that the major rivers follow paleo-rivers systems derived from both 2D and 3D seismic interpretations (Fig. 1b).For surface processes parametrisation, we tested different ranges of diffusion and erodibility coefficients and compared the final sediment accumulation across the Sunda Shelf (Fig. 1b) using estimated deposit thicknesses41,42,43,44. The Sunda Shelf is predominantly experiencing deposition over the past 500 kyr and increases in deposition are positively correlated with periods of sea level rise (i.e., Pearson’s coefficients for correlation with sea level above 80%, Supplementary Fig. 1d). After exploring a range of values, we set κ values to 1. e−2 and 8. e−2 m2/yr for terrestrial and marine environments and ϵ between 2.5 and 7.5 e−8yr−1 for the different scenarios to fit with chosen surveys dataset (Supplementary Table 2 and 3).Upon uniform subsidence case (−0.25mm/yr), flooding is limited, and the shelf only undergoes two full marine transgressions ( >60% of the shelf flooded) around 125 ka and during the last 10–20 kyr (Supplementary Fig. 1c). Upon spatially variable tectonics (non-uniform subsidence), partial flooding events are more pervasive, with higher magnitudes and greater temporal durations. Due to the shallow and flat physiography of Sundaland, we also note that even small increases in sea level amplitudes ( 0) and values higher than one and two standard deviations (zsc  > 1 and 2, respectively, Supplementary Fig. 1b). The approach provides a quantitative assessment of flow maps sensitivity to the chosen resistance maps.To gain additional insights into the distribution of connectivity regions across the shelf, we also employed a local spatial autocorrelation indicator, namely the Getis-Ord Gi⋆ index57. This hotspot analysis method assesses spatial clustering of the obtained current density maps, and the resultant z-scores provide spatially and statistically significant high or low clustered areas. The approach consists in looking at each local current value relative to its neighbouring one. From this spatial analysis, we extract both statistically significant hot and cold spots for each combination of resistance surfaces (Supplementary Fig. 5c). To extract statistically significant and persistent biogeographic connectivity areas across the exposed Sunda Shelf, we then combine all hotspots together and define preferential migration pathways as regions having a positive Gi⋆ z-scores for all resistance surfaces combination.We used the function zscore in the SciPy stats package to obtain the z-scores and the ESDA library for the Gi⋆ indicator computation.Modelling assumptions and limitationsThere are a number of important caveats for interpreting our modelling results.First, we made several assumptions related to our transient landscape evolution simulations. A single-flow direction algorithm54 was used to simulate temporal changes in river pathways. Recent work58 has shown that this algorithm might lead to numerical diffusion, fast degradation of knickpoints and underestimation of river captures particularly in flat regions. One way to address this would be to use a multiple flow direction method59 which allow for a better representation of flow distribution across the landscape. In this study, we also assumed a uniform and invariant soil erodibility coefficient for the entire domain and a detachment-limited erosion law. Even though the erodibility coefficient was calibrated independently for each simulation (Supplementary Table 3), soil cover and properties vary notably between Borneo, Sumatra, Java and the Malay Peninsula and soil conditions for the exposed sea floor would have changed significantly over successive flooding events12. Badlands software34 allows for multiple erodibility coefficients representing different soil compositions to be defined, and this functionality could be used to evaluate the impact of differential erosion on physiographic changes. Similarly, several transport-limited laws are also available and could be compared against our detachment-limited simulations.A second set of simplifications lies in the climatic conditions (i.e., rainfall regimes) used to force our simulations. We relied on the PaleoClim database40 which contains nine high-resolution paleoclimate dataset38,39,40 corresponding to specific time periods (4.2–0.3 ka, 8.326–4.2 ka, 11.7–8.326 ka, 12.9–11.7 ka, 14.7–12.9 ka, 17.0–14.7 ka, ca. 130 ka, ca. 787 ka and 3.205 Ma). The climate simulations from which these time periods are extracted do not consider emerged Sunda Shelf for the oldest inter-glacial events which can result in incorrect climatic pattern60. From 0.3–17 ka, precipitation fields in PaleoClim are obtained from the TRaCE21ka transient simulations of the last 21 kyr run with the CCSM3 model40. Although Fordham et al.39 show that precipitation errors range from 10–200% in their modern experiment, the paleoclim dataset provides a statistical downscaling method that includes a bias correction (namely the Change-Factor method, in which the anomaly between the modern simulation and observations is removed from the paleoclimate experiment) allowing the use of the model for paleoclimate studies40. The very same technique is applied for 130 ka and 787 ka fields that were obtained with different GCMs (namely HadCM3 and CCSM2). Given the absence of a million-year long transient climate simulation, we oversimplified the climatic conditions by considering that precipitation distribution and intensity remain constant between two consecutive intervals, generating an artificial stepwise evolution of rainfall through time. To evaluate the sensitivity of physiographic responses on the Sunda Shelf to precipitation, we ran a model with uniform rainfall over 1 Myr (scenario 4). Despite changes in the timing and extent of basins reorganisation (Supplementary Fig. 2 and Fig. 3b), we found limited differences in terms of flooding history and erosion/deposition patterns when compared with scenario 5 (Supplementary Fig. 1c, d and Supplementary Table 2). Recent work60 suggests clear regional responses induced by the emerged Sunda Shelf with seasonal enhancement of moisture convergence and continental precipitation induced by thermal properties of the land surface. This could significantly impact our simulation results. However, and at the time of writing, more continuous high-resolution paleoclimatic simulations considering the shelf as an emerged continental platform were still unavailable. Using high-resolution multi-model outputs would allow to target the uncertainty on climatic maps4 and will surely represent a significant advance for future studies. One approach would have used the orographic rainfall capability61 available in Badlands. The method is better suited to run generic simulations but falls short when applied to real cases as it relies on imposing paleo-environmental boundary conditions (e.g., temporal changes in wind direction and speed, moisture stability frequency or depth of moist layer) difficult to obtain for Earth-like model applied over geological time scales.Finally, our species-agnostic approach assumes an equally weighted cost between the three considered geomorphic features and does not account for additional factors (temperature, vegetation cover, solar radiation to cite a few), which are all important when assessing landscape connectivity for wildlife. Most importantly, we model connectivity at very large scales (5 km resolution). Often, species are highly influenced by microclimates and small-scale topography47. From our regional-scale simulations and hotspot analysis (Fig. 6), higher resolution models focusing on highly connected regions (across the Gulf of Thailand and Siam basin) could be applied to produce more detailed representations of species migration in the region. In addition, current flow field calculations from Circuitscape35 rely on randomly selecting nodes around the region of interest. For connectivity analysis, we used 33 terrestrial points located around the perimeter of the buffered Sundaland area (white contour line in Fig. 1b). Using a selection of nodes in a buffered region allows to reduce the bias in current density estimates46. However, bias might depend on the buffer size compared to the study area as well as the number of nodes selected46,47. Because of memory limitations and the great number of computed grids used to cover the past 500 kyr, we made a trade-off between buffer size and the number of selected points for pairwise calculations. Additional experiments could possibly be tested to evaluate bias in the proposed connectivity maps potentially using a tilling approach to reduce cell number45. More

Author notesNina DombrowskiPresent address: Royal Netherlands Institute for Sea Research, Department of Marine Microbiology and Biogeochemistry, AB Den Burg, The NetherlandsKiley W. SeitzPresent address: EMBL Heidelberg, Meyerhofstraße 1, Heidelberg, GermanyThese authors contributed equally: Marguerite V. Langwig, Valerie De Anda.AffiliationsDepartment of Marine Science, University of Texas at Austin, Marine Science Institute, Port Aransas, TX, USAMarguerite V. Langwig, Valerie De Anda, Nina Dombrowski, Kiley W. Seitz, Ian M. Rambo & Brett J. BakerDepartment of Microbiology, Biomedicine Discovery Institute, Monash University, Clayton, VIC, AustraliaChris GreeningDepartment of Marine Sciences, University of North Carolina at Chapel Hill, Chapel Hill, NC, USAAndreas P. TeskeAuthorsMarguerite V. LangwigValerie De AndaNina DombrowskiKiley W. SeitzIan M. RamboChris GreeningAndreas P. TeskeBrett J. BakerCorresponding authorsCorrespondence to
Marguerite V. Langwig or Brett J. Baker. More

Literature search and data retrievalWe performed a systematic literature search on the internet (Google Scholar, Web of Science) using the following keywords: [gut microbiota], [animal microbiome], [gut microbiome 16S] and [captive AND wild AND microbiota]. This search yielded 222 articles on animal microbiomes published between 2014 and 2020. The materials and methods of these articles were analysed to ascertain whether the study met the following criteria: (i) all wild and captive samples were processed using identical procedures, (ii) compared wild and captive animals were phylogenetically closely related (members of the same species or species complex), (iii) captive individuals were born in captivity, or no information was provided about the origin of the captive animals; i.e., wild animals brought into captivity and sampled some time later were excluded, (iv) captive animals that underwent a deliberate selection process (e.g. inbred mice or domestic animals) were also excluded for considering them genetically not comparable to the wild counterparts, and (v) only datasets with sample sizes over 12 individuals were considered for analysis. Raw data were extracted from the databases and repositories indicated in the articles (accession numbers listed in the “Bioinformatic resources”).Bioinformatic sequencing data processingDatafiles from the different studies were (i) stored at the University of Copenhagen’s Electronic Research Data Repository (ERDA), (ii) assigned a unique study identifier and (iii) re-processed in the Danish National Supercomputer for Life Sciences ‘Computerome2’ using a new bioinformatic pipeline we developed for processing data with different characteristics, including sequencing mode, read length and 16S rRNA gene fragment. The entire code can be found in the “Bioinformatic resources”. In short, for each individual dataset, we quality-filtered (mean phred score of q = 25) and (if necessary) trimmed and merged the paired-end reads based on the sequence overlap using AdapterRemoval224. Primers (if present) were trimmed using Cutadapt25, and reads were dereplicated with USEARCH Derep26 using a relative minimum copy number threshold of 0.01% of the total sequencing depth. Reads were then converted into zero-ratio OTUs using the denoising algorithm UNOISE327, and USEARCH was used to map the reads back to the OTUs and create an OTU table. HS-Blast28 was used to assign taxonomy against the non-redundant Silva 132 database29, and taxonomic assignments were filtered using different identity thresholds for each taxonomic level: 97% for genus-level taxonomy, 95% for family-level taxonomy, 92% for order-level taxonomy and 90% for higher taxonomic levels30. To minimise the impact of incorrectly assigned taxa, taxonomic annotations below these identity thresholds were converted into unclassified, and not considered in downstream analyses. This pipeline yielded relative read abundances assigned to different taxa for each individual dataset analysed.Data quality filteringIndividual data files generated by the aforementioned pipeline were aggregated by study and host species into genus-level abundance tables. The two datasets of Sarcophilus harrisii retrieved from two different studies were processed independently. Taxonomic resolution was limited to the genus level to maximise taxonomic annotation rate and minimise biases introduced by the different 16S rRNA gene markers employed in the analysed studies. On the one hand, wild animals’ microbial communities often contain taxa that do not map to any catalogued species with enough molecular similarity to assign species-level annotation. On the other hand, the analysed datasets were generated based on the V4, V3–V4 and V1–V3 regions of the 16S rRNA gene (Supplementary Dataset), which hindered comparability at the ASV or zOTU level. We then proceeded to quality-filter the genus-level abundance tables of each species through filtering individuals by minimum sequencing depth, minimum diversity coverage and taxonomic annotation. Only individual datasets with more than 1000 reads and diversity coverage values over 99% were retained, and final genus-level abundance tables that contained at least five animals in each contrasting group were considered for analysis. Since the studied datasets contained traces of dietary items and host DNA, read counts assigned to taxonomic groups not assigned to Bacteria genera, or not present in the LTPs132_SSU release of the SILVA Living Tree (https://www.arb-silva.de/projects/living-tree) used for measuring the phylogenetic relationships among bacteria, were removed to ensure accurate measurements of phylogenetic diversities. In the cases where one group (either wild or captive) outnumbered the other, samples were randomly selected to ensure even sample sizes.Diversity and compositional analysesDiversity and compositional analyses were carried out in the R statistical environment v.3.6.331 and Python 3.8 based on the Hill numbers framework. The operations explained below were conducted using the R packages ape32, dendextend33, dmetar34, hilldiv35, meta36, metamicrobiomeR37, phylosignal38, phytools39, treedist40, vegan41, and the python package qdiv42. Hereafter functions and their respective packages are displayed as ‘package::function’. Statistical significance level was set at a FDR-adjusted p-value of 0.05. All charts and figures in the manuscript were originally generated either in R (full code of all original figures is included in “Bioinformatic resources”) and subsequently modified in Adobe Illustrator to achieve the desired layout without distorting the dimensions of the quantitative elements.Hill numbersThe Hill numbers framework encompasses the group of diversity measures that quantify diversity in units of equivalent numbers of equally abundant taxa43,44—in our context bacteria genera. Hill numbers provide a general statistical framework that is sufficiently robust and flexible to address a wide range of scientific questions that molecular ecologists regularly try to answer through measurement, estimation, partitioning and comparison of diversities45. To obtain a complete vision of the gut microbiome differences between wild and captive animals, we conducted all our diversity and compositional analyses based on three contrasting Hill numbers based metrics: the so-called dR, which only accounts for richness (i.e., order of diversity 0, whether bacteria taxa were present or not), dRE which considered Richness + Evenness of order of diversity 1 (i.e., the relative abundances of bacteria are proportionally weighed) and dRER, which considered Richness, + Evenness + Regularity (i.e., the phylogenetic relationships among bacteria are accounted for). Detailed explanations of these metrics can be found elsewhere17,46,47.Phylogenetic treesThe dRER metric required a Bacterial phylogenetic tree to compute the relatedness among bacterial taxa. As our datasets contained different fragments of the 16S rRNA gene, and thus we were unable to generate a phylogenetic tree directly from our DNA sequence data, we relied on the SILVA Living Tree, and used the LTPs132_SSU release as the reference phylogenetic tree. In addition, the time-calibrated host phylogeny required by the host phylogenetic signal and phylosymbiosis analyses was generated using Timetree48.Diversity metrics and meta-analysisWe computed individual-based diversity metrics using the function hilldiv::hill_div, and obtained average alpha diversity metrics per species, as well as wild and captive populations per species. We used a Kruskal–Wallis (KW) test as implemented in the function hilldiv::div_test to ascertain whether the mean diversity values varied across analysed host species, and a PERMANOVA (PMV) test using vegan::adonis function based on the pairwise dissimilarity matrix to test whether host species were compositionally distinct.Average alpha diversity metrics of wild and captive populations per species were used to conduct a random-effects-model (REM) meta-analysis with raw effect sizes using the function meta::metacont. We used the Sidik–Jonkman estimator for the between-study variance and the Knapp–Hartung–Sidik–Jonkman adjustment method. The overall effect was calculated using Hedge’s g (SMD) and its 95% confidence interval and p-value. An identical analysis was performed for the entire dataset and two representative subsets of five species, containing only datasets derived from primates and cetartiodactylans. Higgin’s and Thompson’s I2 test, Tau-squared T2 and Cochran’s Q were used for quantifying the heterogeneity between the included species. Due to the high heterogeneity found in the study, we evaluated whether the between-study heterogeneity was caused by outliers with extreme effect sizes, which could be distorting our overall effect. We defined an outlier if the species’s confidence interval did not overlap with the confidence interval of the pooled effect using dmetar::find.outliers function.The function detected three outliers in dR metric (GOGO, PEMA and TUTR), four in dRE (GOGO, PEMA, MOCH, EQKI) and seven in dRER (RHBR, PYNE, PEMA, TUTR, MOCH, CENI and AIME). Even when these outliers were excluded from the analysis the I2 heterogeneity value was substantial for dR (I2 from 79.3 to 70.3%) and moderate for dRE (I2 from 80.1 to 60.0%) and dRER (I2 from 86.9 to 54.2%) and significant for both (Cochran’s Q, p-value  More

Sample collectionSamples were collected from an evergreen broadleaf forest in central Taiwan (Lianhuachi Experimental Forest, Nantou County, 23°55′7″N 120°52′58″E) from January 2017 to March 2018. Permission to collect plants for the study was obtained from the Lianhuachi Research Center, Taiwan Forestry Research Institute, Council of Agriculture, Executive Yuan, Taiwan (Permission no.: 1062272538). The present study complies with the International Union for Conservation of Nature Policy Statement on Research Involving Species at Risk of Extinction and the Convention on the Trade in Endangered Species of Wild Fauna and Flora. Ant cadavers with fungal growth were collected from understory plants with a height of less than 3 m. Ant cadavers infected with O. unilateralis s. l. were removed carefully by cutting the leaf and placing it into a 50-mL conical centrifuge tube, which was then transported to the laboratory. Only cadavers in which the fungal growth stage preceded the development of perithecia, which theoretically has the highest biological activity, were collected (Fig. 1). In total, 24 infected P. moesta and 20 infected P. wolfi samples were collected.Figure 1Ophiocordyceps unilateralis sensu lato-infected (a) Polyrhachis moesta and (b) P. wolfi, with the stroma growing from the ant cadaver. The specimens were collected from the Lianhuachi Research Center, Taiwan and photographed in the laboratory by Wei-Jiun Lin.Full size imageIsolation and cultivation of bacteriaAnts on the leaves were first identified to species and then, using tweezers, each ant was placed carefully into a sterilized 1.5-mL microcentrifuge tube [see details in Lin et al. (2020)15. Samples were shaken one by one in 600 μL of sterilized water for a few seconds at 3000 revolutions/min (rpm) using a vortex mixer (AL-VTX3000L, CAE technology Co., Ltd., Québec, Canada), and were then soaked with 600 μL of 70% ethanol to sterilize the ant’s surface. The ethanol on the samples was washed twice with 600 μL of sterilized water, then vortexed in 400 μL of sterilized water. Next, 200 μL of the supernatant was spread homogeneously onto a Luria–Bertani (LB) agar plate (25 g Luria–Bertani broth and 15 g agar per liter) to confirm the absence of live bacteria.Bacteria from inside the ant host were released by homogenizing the ant host in 200 μL of water and culturing on LB agar plates at 28 °C for 2 days. Bacteria from each of the ant individuals were cultured independently and approximately equal numbers of the isolates were picked randomly with sterile toothpicks, and were suspended in the LB medium supplemented with 15% v/v glycerol and maintained at − 80 °C until the time of examination. In total, 247 bacterial isolates from P. moesta and 241 bacterial isolates from P. wolfi were collected.In addition to the bacterial isolates from the ant bodies, 60 bacterial isolates from soil, leaves, and air in the same forest were collected for the purpose of comparing their resistance to naphthoquinones (see below) by using the aforementioned procedure but without initial cleaning and sterilizing of the sample surface.Bacterial identificationBacteria collected from the ant hosts were identified by gene marker sequencing. Bacterial isolates were cultured in LB medium at 28 °C overnight to reach the log-phase, and genomic DNA was extracted following the methods described in Vingataramin and Frost (2015)20. Conspecies/strains of the bacterial isolates from the same host were determined using the randomly amplified polymorphic DNA (RAPD) method with the primer 5′-GAGGGTGGCGGTTCT-3′. PCR amplification was performed as follows: initial denaturation at 95 °C for 5 min, 40 cycles of amplification including denaturation at 95 °C for 1 min, annealing at 42 °C for 30 s, and extension at 72 °C for 1 min, followed by a final extension at 72 °C for 10 min. PCR products were run in 2% agarose gel and bacterial isolates were characterized by fragment patterns. For each of the ant hosts, bacterial isolates with the same RAPD pattern were considered to be the same strain. In total, 106 and 178 strains were found from P. moesta and P. wolfi, respectively. One of the bacterial isolates was selected at random to represent the strain and coded with “JYCB” followed by a series of numbers (e.g., JYCB191). Taxonomic status of each strain was determined to species by using the V3/V4 region of the 16S rDNA gene. PCR amplification with the primer set (8F: 5′-AGAGTTTGATCCTGGCTCAG-3′ and 1541R: 5′-AAGGAGGTGATCCAGCCGCA-3′)21,22 was performed under the following conditions: initial denaturation at 95 °C for 5 min, 40 cycles of amplification including denaturation at 95 °C for 1 min, annealing at 55 °C for 30 s, and extension at 72 °C for 1 min 45 s, followed by a final extension at 72 °C for 10 min. PCR products were first checked by running a gel, and were then sequenced at Genomics, Inc. (New Taipei City, Taiwan).The sequences of the bacterial strains from each of the ant hosts were first analyzed by the unweighted pair group method with arithmetic mean (UPGMA) analysis and clustered into clades according to the sequence dissimilarity ( More

Simulating competition for space using the “BacGo” modelTo investigate how spatial positioning of populations affects the outcome of microbial competition, we simulated two populations competing for space with a limiting size by building an individual-based model (named “BacGo”). The model was implemented in discrete grid boxes of a 20 × 20 array. As shown in Fig. 1a, our simulations were based on three basic assumptions. First, the two competing populations possess the same inherent growth rate and equal initial cell numbers, thus the only differences between them are their manners of colonizing free space. Second, the newly born daughter cell is located around its mother cell but with a random direction of spatial positioning [34], resulted in a microcolony with different spatial patterning. Lastly, if the selected box has been occupied, the newborn cell will compete for the box against the original occupants of the box and possesses a probability of 50% to survive [37].We first explored the outcome of spatial competition, which started by randomly distributing two populations on the grids with the same initial cell numbers of 10 for each. Based on our basic assumptions and the predictions of competitive exclusion theory [38], we hypothesized that only one population could win the competition and finally occupy all grids. As shown in 20,000 independent simulations with random initial distributions, we discovered that at the end of each simulation, only one population survived (Video S1 and Video S2). The Chi-square test showed no significant difference (P = 0.211) between the simulated winning times (10,177 of 20,000 simulations) and the random winning times (10,051 of 20,000 simulations) of the focus population. This result conformed with our initial assumption that cells possess a probability of 50% to survive in competing with original occupants. When we replicated simulations initiated with the same cell distribution, we found that the winning probabilities for each population changed in line with the initial distributions (Fig. S1). However, the winning probabilities never reached 100% no matter how the initial distribution changes. Together, these results suggested that unknown random factors may affect the final outcome of the competition.Next, we analyzed the dynamics of microbial colonization during our simulations. As summarized in Fig. 1b, we divided the competition process into two stages, the “occupation stage” and the “exclusion stage” (see Methods). To statistically characterize the competitive outcome at t3, we defined the winning asymmetry index, WinR, and the abundance asymmetry index, AbunR (see Methods). As shown in Fig. S2a, we found a strong positive correlation (R2 = 0.740, P  More

SamplingWater was sampled from Cep lake in Czechia, at a regular sampling site of 10 m depth (48.944 °N, 14.877 °E). The lake originates from sand mining in the 1970–80s. It is a permanent meso-oligotrophic (chlorophyll-a concentrations ranged from 1.4–16.4 µg L−1) seepage reservoir filled with groundwater penetrating from the nearby river Lužnice. The lake area is about 1.16 km2, with the maximum depth about 11–12 m. These characteristics are representative for most of temperate and boreal lakes [16].Samples were collected every four weeks from April till October in 2018, and from April till November in 2019. Ten liters of water were collected from 0.5 m depth using a Ruttner Water Sampler (model 11.003KC Denmark AS). Temperature and oxygen profiles were taken with an EXO1 multi-parameter probe (YSI Inc., Yellow Springs, USA). Water was transported to the laboratory within 2 h from the collection in a closed container made from high-density polyethylene, rinsed three times with the sampled water and stored in a cooled box.NutrientsSamples were filtered through glass fiber filters with 0.4 µm nominal porosity (GF-5, Macherey-Nagel, Düren, Germany). Concentrations of soluble reactive phosphorus (SRP) were determined spectrophotometrically [17, 18]. Concentrations of nitrate and ammonium were measured according to Procházková [19] and Kopáčkek and Procházková [20]. Dissolved organic carbon (DOC) and dissolved nitrogen (DN) were measured with a TOC 5000 A analyzer (Shimadzu, Kyoto, Japan).PigmentsSeston from 1.43 to 3.65 L of water was collected onto GF-5 glass fiber filters (diameter 47 mm, Macherey-Nagel). The filters were dried of excess water by gently pressing in a paper towel, and flush frozen in liquid nitrogen. Pigments were extracted and analyzed by HPLC as described in Piwosz et al. [21].Net primary production (NPP) and community respiration (oxygen measurements)Oxygen concentration was measured with the Winkler method [22]. It was chosen because it allows O2 concentration to be measured directly in the water without the need to consider carbonic equilibrium, which is the case when changes in CO2 concentration is measured [23]. Samples were unfiltered to avoid the removal of particle-associated bacteria and also of free-living AAP bacteria, which tend to be larger than average freshwater bacteria [12]. Glass stoppered Winkler type oxygen bottles (115 mL nominal capacity, VTR glass, Prague, Czechia) were filled with the sampled water directly from the sampler via a rubber tube. Each bottle was first rinsed three times and then filled without the formation of air bubbles. Water was allowed to overflow the neck of the bottle for about 1 min, and the bottle was closed with a glass stopper to avoid air bubbles. The closed bottles were kept in the dark in a cooled box. On the shore, three bottles were selected as T0, and 1.2 mL of manganese (II) chloride solution (concentration 3 mol L−1) was addded, followed by the addition of 1.2 mL of a mixture containing 4 mol L−1 of sodium iodide solution and 8 mol L−1 of sodium hydroxide solution. These samples were processed in the laboratory within 3 h. The remaining bottles were incubated for 24 h at in situ temperature in the IR-box prepared from the MAKROLON IR polycarbonate sheet (4 mm thickness, Professional Plastics, Inc. Fullerton, CA, USA). These panels have a maximum transmittance of 90% in the infrared region 850–2000 nm, 50% at 780 nm, and 0% 10 °C according to our measurements (180 days, Equation 3 in Supplementary File S1). Subsequently, the differences in the Cep Lake’s carbon budget for the surface layer (down to 0.5 m depth) were calculated by multiplying the integrated values by the volume of this layer (Equation 4 and 5 in Supplementary File S1).HCO3
− incorporationTriplicated water samples (32 mL) were incubated for 3.2–5.2 h in the IR light and the dark at in situ temperature, as described for respiration. Total activity added to each bottle was measured from 1 mL aliquot of the incubated sample that was transferred to a scintillation vial containing 20 μl of 5 mol L−1 NaOH (to prevent a loss of 14C-bicarbonate). Thirty mL of sample was filtered through 2.5 μm nitrate cellulose filters (Pragopor, Prague, Czechia, diameter 25 mm). Five mL of the filtered water was collected and subsequently filtered through a 0.17 μm nitrate cellulose filter. The resulting cell-free filtrate, which contained 14C-DOC was collected. The filtration was done at a low vacuum (0.02 MPa) to avoid cell breakage. The total CO2 assimilation rate was calculated as the sum of all these fractions.The filters were kept in an HCl-saturated atmosphere for 24 h at room temperature in a fume hood. They were placed in scintillation vials and dissolved in 1 mL of ethyl acetate (Penta, Prague, Czechia). Then, 5 mL of Ultima Golt LLT scintillation cocktail (PerkinElmer, Waltham, MA, USA) was added. Five mL of cell-free filtrates were acidified by adding 100 μL 5 mol L−1 HCl to volatilize non-incorporated H14CO3 and incubated 24 h at room temperature in a fume hood. Then, 10 mL of the scintillation cocktail was added. Finally, 5 mL of the scintillation cocktail was added to the total activity samples. Subsequently, the samples were gently mixed and left in the dark for 48 h. The radioactivity in the samples was measured using a Tri-Carb 2810 TR scintillation counter (PerkinElmer).To estimate carbon fluxes (μmol C L–1 h–1), a fraction of the added H14CO3 incorporated or released was corrected for the incubation time and multiplied by the concentration of total dissolved inorganic carbon (DIC). The DIC concentration was calculated based on temperature, pH, and alkalinity measurements (Inolab pH 720, WTW Xylem Inc. Rye Brook, NY, US) determined by Gran titration.Assimilation of organic monomersThe difference between microbial activity in the IR light and dark was also estimated based on assimilation rates of radiolabeled glucose, pyruvate, leucine and thymidine (American Radiolabeled Chemicals, St. Louis, MO, USA). Tritiated glucose (specific activity (SA): 2220 GBq mmol−1), leucine (SA: 4440 GBq mmol−1) and thymidine (SA: 2275.5 GBq mmol−1) were added to 5 mL samples to a final concentration of 5 nmol L−1. 14C-pyruvate (SA: 2.035 GBq mmol−1) was added to a final concentration of 10 nmol L−1. Trichloroacetic acid (TCA) was added to the killed controls to a final concentration of 1%. Samples were incubated for 1 h in the dark and IR light as described for respiration. The incubations were terminated as the killed controls and kept at 4 °C in the dark until processed within 65% except for the samples from 9th May and 29th Aug 2018 (10% each), 1st Aug 2018 (21%), 25th Nov 2018 (1%), and 14th Aug 2019 (22%), Supplementary Fig. S1A). Thus, we decided to concatenate the fastaq files and analyze both fractions together as the total community. This also facilitated statistical analysis, as the activity rates were measured for the whole community without fractionations.Reads quality was evaluated using FastQC v0.11.7 (Babraham Bioinformatics, Cambridge, UK). After primer sequences trimming using Cutadapt [29] (v1.16), the number of reads per sample ranged from 49,354 to 188,942. Subsequent analyses were done in the R/Bioconductor environment using the dada2 package (version 1.14.1) [30]. Forward and reverse reads were truncated to 225 bp and low quality sequences were filtered out with the filterAndTrim function (truncLen = c(225, 225), maxN = 0, maxEE = c(2, 2), truncQ = 2), which reduced the number of reads per sample to range from 30,190 to 143,552. After merging and chimera removal using the removeBimeraDenovo function, 4,893 amplicon sequence variants (ASV) were obtained. Rare ASVs (not seen >3 times in at least 20% of the samples) were removed, which reduced the number of ASVs to 658, and the number of reads to 14,613–69,046 per sample. Taxonomic assignment was done using SILVA 138.1 database [31, 32] released on August 27, 2020. ASVs identified as Chloroplast or Cyanobacteria were excluded from the analyses, giving the final number of 546 ASVs and from 10,819 to 54,799 reads per sample. The bacterial community composition graphs were done using phyloseq [33] and ggplot2 [34] packages.AAP community compositionThe composition of AAP community was analyzed by amplicon sequencing of pufM gene encoding the M subunit of bacterial type-2 reaction centers. This gene is routinely used for diversity studies of AAP bacteria [35].PufM gene amplicons (approx. 245 bp) were prepared using pufM_UniF (5′-GGN AAY YTN TWY TAY AAY CCN TTY CA-3′) and pufM_WAW (5′-AYN GCR AAC CAC CAN GCC CA-3′) primers [36]. PCR was performed in triplicate 20 μL reactions using Phusion™ High-Fidelity DNA Polymerase (Thermo Scientific, USA) with the following reaction conditions: 98 °C for 3 min, 27 cycles at 98 °C for 10 s, 58 °C for 30 s, 72 °C for 30 s, and a final extension at 72 °C for 5 min. The triplicate product reactions for each sample were pooled and gel purified using the kit Wizard SV Gel and PCR Clean-Up System (Promega, USA). The sequencing was performed on the Illumina MiSeq platform (2 × 150 bp) at Macrogen, South Korea.The fastq files were concatenated as described for bacteria communities. The Bray-Curtis similarity between two fractions for each sampling day was >70%, except for the samples from 1st Aug 2018 (47%) and 14th Aug 2019 (18%, Supplementary Fig. S1B).The samples were analyzed as described for bacterial communities. The number of reads per sample ranged from 192,360 to 239,418 after the cutadapt trimming. Forward and reverse reads were truncated to 130 bp, and the number of reads per sample after the quality filtering and denoising ranged from 189,432 to 235,311. Merging the forward and reverse reads with mergePairs function created 12,692 ASVs and reduced the number of reads to 183,136–221,281 per sample. The chimera removal lowered the number of ASVs to 1816, and the number of reads to 159,451–208,679. Rare ASVs (not seen >3 times in at least 20% of the samples) were removed, which resulted in the final 566 ASVs, and a number of reads ranging from 155,915 to 203,021 per sample. A manually curated taxonomic database was used for taxonomic assignment following the naïve Bayesian classifier method [37]. It contained 1580 unique pufM sequences, downloaded from the Fungene repository on May 16, 2019 (http://fungene.cme.msu.edu [38]), from metagenomes from the Římov Reservoir [39, 40] and from the Genome Taxonomy database accessed on September 16, 2020 [41].Statistical analysisLinear mixed-effects models were calculated in R (version 3.6.2) using lme function from the nlme package (version 3.1.143) on untransformed activity data and log10 transformed environmental variables [42]. Models’ parameters were estimated using maximum likelihood method and their significance was tested with ANOVA. Relationships between the activity measures, the environmental variables and the composition of AAP communities were investigated with distance-based linear models (DistLM) [43, 44] in Primer (version 7.0.13) with PERMANOVA + 1 add on (e-Primer, Plymouth, UK) [45]. The sequence reads were transformed with the varianceStabilizingTransformation function of the DESeq2 package [46] (version 1.14.1, blind = FALSE, fitType = “mean”).Data accessibilityThe sequences of 16S and pufM amplicons that support the findings of this study have been deposited in the EMBL database as the BioProject with the accession number PRJEB41596, together with most of the environmental metadata. The scripts and the remaining data supporting the results are included in the Supplementary Material. More