Ecology

Dispersal is assumed to contribute to microbiome composition and function; however, it is difficult to measure. Walters et al. now set out a 6-month experiment looking at different dispersal routes of environmental microorganisms to the surface soil layer. They set up different ‘traps’, either glass slides or freshly cut grass, to determine the number, identity and function of incoming microorganisms. The traps ‘recorded’ dispersal through air, from plants and their litter, or from below through the decomposing litter and bulk soil. This was achieved by placing the traps either on a pedestal, closing them off at the bottom or leaving them open, respectively. The authors found that the overall dispersal rate was low, with little influence of the route, with only 0.5% incoming bacterial cells per day compared with the number of resident cells. However, the dispersal routes did influence microbiome composition, at least if from above and close to the surface. Finally, without dispersal, the initial decomposition of the cut grass was slower.
This is a preview of subscription content More

Bailey, K. M. & Houde, E. D. Predation on eggs and larvae of marine fishes and the recruitment problem. Adv. Mar. Biol. 25, 1–83. https://doi.org/10.1016/S0065-2881(08)60187-X (1989).Article 

Google Scholar 
Houde, E. D. Fish early life dynamics and recruitment variability. Am. Fish. Soc. Symp. 2, 17–29 (1987).ADS 

Google Scholar 
Anderson, J. T. A review of size dependent survival during pre-recruit stages of fishes in relation to recruitment. J. Northw. Atl. Fish. Sci. 8, 55–66. https://doi.org/10.2960/J.v8.a6 (1988).Article 

Google Scholar 
McCarthy, I. D. Temporal repeatability of relative standard metabolic rate in juvenile Atlantic salmon and its relation to life history variation. J. Fish Biol. 57, 224–238. https://doi.org/10.1111/j.1095-8649.2000.tb00788.x (2000).Article 

Google Scholar 
Biro, P. A. & Stamps, J. A. Do consistent individual differences in metabolic rate promote consistent individual differences in behavior?. Trends Ecol. Evol. 25, 653–659. https://doi.org/10.1016/j.tree.2010.08.003,Pubmed:20832898 (2010).Article 
PubMed 

Google Scholar 
Endler, J. A. Natural Selection in the Wild (Princeton Univ. Pr., 1986).Meekan, M. G. & Fortier, L. Selection for fast growth during the larval life of Atlantic cod Gadus morhua on the Scotian Shelf. Mar. Ecol. Prog. Ser. 137, 25–37. https://doi.org/10.3354/meps137025 (1996).ADS 
Article 

Google Scholar 
Gilly, W. F. et al. Vertical and horizontal migrations by the jumbo squid Dosidicus gigas revealed by electronic tagging. Mar. Ecol. Prog. Ser. 326, 1–17 (2006).ADS 
Article 

Google Scholar 
Watanabe, H., Kubodera, T., Moku, M. & Kawaguchi, K. Diel vertical migration of squid in the warm core ring and cold water masses in the transition region of the western North Pacific. Mar. Ecol. Prog. Ser. 315, 187–197. https://doi.org/10.3354/meps315187 (2006).ADS 
Article 

Google Scholar 
Phillips, K. L., Jackson, G. D. & Nichols, P. D. Predation on myctophids by the squid Moroteuthis ingens around Macquarie and Heard Islands: stomach contents and fatty acid analyses. Mar. Ecol. Prog. Ser. 215, 179–189. https://doi.org/10.3354/meps215179 (2001).ADS 
CAS 
Article 

Google Scholar 
Field, J. C., Baltz, K., Phillips, A. J. & Walker, W. A. Range expansion and trophic interactions of the jumbo squid, Dosidicus gigas, in the California Current. CalCOFI Rep. 48, 131–146 (2007).
Google Scholar 
Ellis, T. & Gibson, R. N. Size-selective predation of 0-group flatfishes on a Scottish coastal nursery ground. Mar. Ecol. Prog. Ser. 127, 27–37. https://doi.org/10.3354/meps127027 (1995).ADS 
Article 

Google Scholar 
Takasuka, A., Aoki, I. & Oozeki, Y. Predator-specific growth-selective predation on larval Japanese anchovy Engraulis japonicus. Mar. Ecol. Prog. Ser. 350, 99–107. https://doi.org/10.3354/meps07158 (2007).ADS 
Article 

Google Scholar 
Tucker, S., Hipfner, J. M. & Trudel, M. Size- and condition-dependent predation: A seabird disproportionately targets substandard individual juvenile salmon. Ecology 97, 461–471. https://doi.org/10.1890/15-0564.1,Pubmed:27145620 (2016).Article 
PubMed 

Google Scholar 
Rodhouse, P. G. & Nigmatullin, C. M. Role as consumers. Phil. Trans. R. Soc. Lond. B 351, 1003–1022. https://doi.org/10.1098/rstb.1996.0090 (1996).ADS 
Article 

Google Scholar 
Hunsicker, M. E. & Essington, T. E. Size-structured patterns of piscivory of the longfin inshore squid (Loligo pealeii) in the mid-Atlantic continental shelf ecosystem. Can. J. Fish. Aquat. Sci. 63, 754–765. https://doi.org/10.1139/f05-258 (2006).Article 

Google Scholar 
Hunsicker, M. E. & Essington, T. E. Evaluating the potential for trophodynamic control of fish by the longfin inshore squid (Loligo pealeii) in the northwest Atlantic Ocean. Can. J. Fish. Aquat. Sci. 65, 2524–2535. https://doi.org/10.1139/F08-154 (2008).Article 

Google Scholar 
Wang, K. Y., Liao, C. H. & Lee, K. T. Population and maturation dynamics of the swordtip squid (Photololigo edulis) in the southern East China Sea. Fish. Res. 90, 178–186. https://doi.org/10.1016/j.fishres.2007.10.015 (2008).Article 

Google Scholar 
Sassa, C., Yamamoto, K., Tsukamoto, Y., Konishi, Y. & Tokimura, M. Distribution and migration of age-0 jack mackerel (Trachurus japonicus) in the East China and Yellow Seas, based on seasonal bottom trawl surveys. Fish. Oceanogr. 18, 255–267. https://doi.org/10.1111/j.1365-2419.2009.00510.x (2009).Article 

Google Scholar 
Tokai, T., Shiode, D., Sakai, T. & Yoda, M. Codend selectivity in the East China Sea of a trawl net with the legal minimum mesh size. Fish. Sci. 85, 19–32. https://doi.org/10.1007/s12562-018-1270-x (2019).CAS 
Article 

Google Scholar 
Sassa, C. & Konishi, Y. Vertical distribution of jack mackerel Trachurus japonicus larvae in the southern part of the East China Sea. Fish. Sci. 72, 612–619. https://doi.org/10.1111/j.1444-2906.2006.01191.x (2006).CAS 
Article 

Google Scholar 
Takahashi, M., Sassa, C. & Tsukamoto, Y. Growth-selective survival of young jack mackerel Trachurus japonicus during transition from pelagic to demersal habitats in the East China Sea. Mar. Biol. 159, 2675–2685. https://doi.org/10.1007/s00227-012-2025-3 (2012).Article 

Google Scholar 
Ishida, K. Feeding ecology of swordtip squid (Loligo edulis). Rep. Shimane Pref. Fish. Exp. Stan. 3, 31–35 (1981) (in Japanese).
Google Scholar 
Tashiro, M., Tokunaga, T., Machida, S. & Takata, J. Distribution of a squidfish, Loliogo edulis HOYLE, in the East China Sea. Bull. Nagasaki Pref. Inst. Fish. 7, 21–30 (1981) (in Japanese).
Google Scholar 
Jennings, S. & Warr, K. J. Smaller predator-prey body size ratios in longer food chains. Proc. Biol. Sci. 270, 1413–1417. https://doi.org/10.1098/rspb.2003.2392 (2003).Article 
PubMed 
PubMed Central 

Google Scholar 
Barnes, C., Maxwell, D., Reuman, D. C. & Jennings, S. Global patterns in predator-prey size relationships reveal size dependency of trophic transfer efficiency. Ecology 91, 222–232. https://doi.org/10.1890/08-2061.1 (2010).Article 
PubMed 

Google Scholar 
Cabana, G. & Rasmussen, J. B. Modelling food chain structure and contaminant bioaccumulation using stable nitrogen isotopes. Nature 372, 255–257. https://doi.org/10.1038/372255a0 (1994).ADS 
CAS 
Article 

Google Scholar 
Castilla, A. C., Hernández-Urcera, J., Gouranguine, A., Guerra, Á. & Cabanellas-Reboredo, M. Predation behaviour of the European squid Loligo vulgaris. J. Ethol. 38, 311–322. https://doi.org/10.1007/s10164-020-00652-4 (2020).Article 

Google Scholar 
Fiorito, G. et al. Guidelines for the Care and Welfare of Cephalopods in Research–A consensus based on an initiative by CephRes, FELASA and the Boyd Group. Lab. Anim. 49, 1–90. https://doi.org/10.1177/0023677215580006la.sagepub.com (2015).Article 
PubMed 

Google Scholar 
Percie du Sert, N. et al. Reporting animal research: Explanation and elaboration for the ARRIVE guidelines 2.0. PLoS Biol. 18, e3000411 (2020). https://doi.org/10.1371/journal.pbio.3000411Campana, S. E. How reliable are growth back-calculations based on otoliths?. Can. J. Fish. Aquat. Sci. 47, 2219–2227. https://doi.org/10.1139/f90-246 (1990).Article 

Google Scholar 
Xie, S. et al. Growth and morphological development of sagittal otoliths of larval and early juvenile Trachurus japonicus. J. Fish Biol. 66, 1704–1719. https://doi.org/10.1111/j.0022-1112.2005.00717.x (2005).Article 

Google Scholar 
Yasui, T. & Sakurai, Y. Gastric evacuation rate of Todarodes pacificus. Rep. Annu. Meet. Squid Res 32, 55–57 (2005) (in Japanese).
Google Scholar 
Šifner, S. K. & Vrgoč, N. Population structure, maturation and reproduction of the European squid, Loligo vulgaris, in the Central Adriatic Sea. Fish. Res. 69, 239–249. https://doi.org/10.1016/j.fishres.2004.04.011 (2004).Article 

Google Scholar 
Kono, N., Tsukamoto, Y. & Zenitani, H. RNA:DNA ratio for diagnosis of the nutritional condition of Japanese anchovy larvae Engraulis japonicus during the first-feeding stage. Fish. Sci. 69, 1096–1102. https://doi.org/10.1111/j.0919-9268.2003.00733.x (2003).CAS 
Article 

Google Scholar 
Booman, C., Folkvord, A. & Hunter, J. R. Responsiveness of starved northern anchovy Engraulis mordax larvae to predation attacks by adult anchovy. Fish. Bull. 89, 707–711 (1991).
Google Scholar 
Chick, J. H. & Van Den Avyle, M. J. Effects of feeding ration on larval swimming speed and responsiveness to predator attacks: Implications for cohort survival. Can. J. Fish. Aquat. Sci. 57, 106–115. https://doi.org/10.1139/f99-185 (2000).Article 

Google Scholar 
Hunsicker, M. E. et al. Functional responses and scaling in predator-prey interactions of marine fishes: Contemporary issues and emerging concepts. Ecol. Lett. 14, 1288–1299. https://doi.org/10.1111/j.1461-0248.2011.01696.x (2011).Article 
PubMed 

Google Scholar 
Chambers, R. C. & Miller, T. J. Evaluating fish growth by means of otolith increment analysis: spectral properties of individual-level longitudinal data in in Recent Developments in Fish Otolith Research (ed. Secor, D. H., Dean, J. M. & Campana, S. E.) 155–175 (University of South Carolina Press, 1995).Mizutani, T. et al. Diel variability in the catch composition of bottom trawl survey in East China Sea. Nippon Suisan Gakkaishi 71, 44–53 (2005). (in Japanese with English abstract). https://doi.org/10.2331/suisan.71.44.Sassa, C., Takahashi, M., Konishi, Y. & Tsukamoto, Y. Interannual variations in distribution and abundance of Japanese jack mackerel Trachurus japonicus larvae in the East China Sea. ICES J. Mar. Sci. 73, 1170–1185. https://doi.org/10.1093/icesjms/fsv269 (2016).Article 

Google Scholar 
Takahashi, M., Sassa, C., Nishiuchi, K. & Tsukamoto, Y. Interannual variations in rates of larval growth and development of jack mackerel (Trachurus japonicus) in the East China Sea: Implications for juvenile survival. Can. J. Fish. Aquat. Sci. 73, 155–162. https://doi.org/10.1139/cjfas-2015-0077 (2016).Article 

Google Scholar 
Takahashi, M., Sassa, C., Nishiuchi, K. & Tsukamoto, Y. Variability in growth rates of Japanese jack mackerel Trachurus japonicus larvae and juveniles in the East China Sea—effects of temperature and prey abundance in in Kuroshio Current, Physical, Biogeochemical and Ecosystem Dynamics (ed. Nagai, T., Saito, H., Suzuki, K. & Takahashi, M.) 295–307 (Wiley, 2019).Anraku, M. & Azeta, M. The feeding habits of larvae and juveniles of the yellowtail, Seriola quinqueradiata Temminck et Schlegel, associated with floating seaweeds. Bull. Seikai Reg. Fish. Res. Lab 33, 13–45 (1965) (in Japanese with English abstract).
Google Scholar 
Villanueva, R., Perricone, V. & Fiorito, G. Cephalopods as predators: a short journey among behavioral flexibilities, adaptations, and feeding habits. Front. Physiol. 8, 598. https://doi.org/10.3389/fphys.2017.00598,Pubmed:28861006 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Wang, R., Zuo, T. & Wang, K. The Yellow Sea Cold Bottom Water—an oversummering site for Calanus sinicus (Copepoda, Crustacea). J. Plankton Res. 25, 169–183. https://doi.org/10.1093/plankt/25.2.169 (2003).CAS 
Article 

Google Scholar 
Sassa, C., Kitajima, S., Nishiuchi, K. & Takahashi, M. Ontogenetic and inter-annual variation in the diet of Japanese jack mackerel (Trachurus japonicus) juveniles in the East China Sea. J. Mar. Biol. Assoc. U K 99, 525–538. https://doi.org/10.1017/S0025315418000206 (2019).Article 

Google Scholar 
Nakazawa, T., Ushio, M. & Kondoh, M. Scale dependence of predator–prey mass ratio: Determinants and applications. Adv. Ecol. Res. 45, 269–302. https://doi.org/10.1016/B978-0-12-386475-8.00007-1 (2011).Article 

Google Scholar 
Ohshimo, S., Tanaka, H., Nishiuchi, K. & Yasuda, T. Trophic positions and predator-prey mass ratio of the pelagic food web in the East China Sea and Sea of Japan. Mar. Freshw. Res. 67, 1692–1699. https://doi.org/10.1071/MF15115 (2016).Article 

Google Scholar 
Vidal, E. A. G. & Salvador, B. The tentacular strike behavior in squid: functional interdependency of morphology and predatory behaviors during ontogeny. Front. Physiol. 10, 1558. https://doi.org/10.3389/fphys.2019.01558 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Doubleday, Z. A. et al. Global proliferation of cephalopods. Curr. Biol. 26, R406–R407. https://doi.org/10.1016/j.cub.2016.04.002 (2016).CAS 
Article 
PubMed 

Google Scholar 
Overholtz, W. J., Link, J. S. & Suslowicz, L. E. Consumption of important pelagic fish and squid by predatory fish in the northeastern USA shelf with some fishery comparisons. ICES J. Mar. Sci. 57, 1147–1159. https://doi.org/10.1006/jmsc.2000.0802 (2000).Article 

Google Scholar 
Montevecchi, W. A. & Myers, R. A. Prey harvests of seabirds reflect pelagic fish and squid abundance on multiple spatial and temporal scales. Mar. Ecol. Prog. Ser. 117, 1–9. https://doi.org/10.3354/meps117001 (1995).ADS 
Article 

Google Scholar  More

The Sargasso Sea has long been considered as the spawning area for Atlantic eels, despite the absence of direct observations after more than a hundred years of the survey. We proposed a new insight on the location of Atlantic eels spawning areas eastward of the Sargasso Sea at the intersection between the Mid-Atlantic Ridge and the oceanic fronts1. Our hypothesis is based on a body of corroborating cues from literature. We suggested that European silver eels converge towards the Azores whatever their departure point from Europe and Northern Africa, then they follow the Mid-Atlantic Ridge south westerly until they reach oceanographic fronts where temperature and depths are favourable for reproduction. These orientation behaviours are potentially based on magnetic fields and odours that might be generated by the Mid-Atlantic Ridge volcanic activity and detected by eels during their diel vertical movements. The first favourable meeting point is then located at the crossing between the Mid-Atlantic Ridge and the oceanic thermic isotherms located around 45° W and 26° N. Our hypothesis is supported by (i) microchemical differences between the core of otoliths extracted from leptocephali collected in the Sargasso Sea and from glass eels collected across Europe suggesting that glass eels hatch in different chemical environments than leptocephali (ii) an asymmetric genetic introgression between American and European eels2 suggesting that the overlapping spawning areas favour transport of hybrids towards northern Europe rather than to America and to southern Europe. This supports the possible existence of several distinct spawning areas, where currents favour transport either westward (American eel), north eastward (hybrids and European eels) or eastward (European eels). To test this hypothesis, we developed a transport model and compared the dispersion dynamics of virtual leptocephali released from the Sargasso Sea and from above the Mid-Atlantic Ridge. The transport models showed that virtual eels released from the Mid-Atlantic Ridge reached Europe and America following similar patterns than those released from the Sargasso Sea thus supporting the Mid-Atlantic Ridge spawning hypothesis.Hanel et al.3 have raised several concerns, one of which being that “microchemical evidence was the only was the major argument supporting the Mid-Atlantic Ridge hypothesis”. This was their start point of a critical rebuttal of our findings to question our hypothesis. Instead, we consider that our regrettable error does not fundamentally contradict the possibility that eels do indeed successfully spawn outside the so-called Sargasso Sea.(Comment 1) The importance of seamounts as orientation and navigation cues towards a spawning area was hypothesized, no clear mechanism is proposed for how the migrating eels can detect the ridge.(Response 1) Our Hypothesis does not state that eels find a kind of shallow seamount where they spawn. Instead, we propose that orientation of silver eels during their spawning migration could be based on a combination of behavioural mechanisms including geomagnetism, odours, temperature and salinity gradients4,5,6,7,8. These environmental cues and related gradients are strongly controlled or influenced by the topography of the oceanic floor. The Mid-Atlantic Ridge and the Mariana areas have similarities with ridges and seamount chains oriented perpendicularly to temperature and salinity fronts surrounded by deep abyssal plains. Our Mid-Atlantic Ridge hypothesis proposed that Atlantic eels could use similar signposts as Japanese eel, which hatch near the Mariana Ridge9. Indeed, as for the Japanese eels, the orientation mechanism that lead Atlantic eels from the growth areas to the ridge are not understood, but the empirical observations from Righton et al.10 suggest that eels converge towards the Azores whatever their release point across Europe and that their diel vertical migration takes them down to 500–1000 m every day. The reasons for this behaviour are not elucidated, but since they cost energy, they are likely compensated by advantages such as orientation together with predator avoidance and sexual maturation11,12,13. Following our hypothesis, eels search for orientation cues during DVM. The geomagnetic fields are suggested to provide detectable information for silver eels on their oceanic spawning migration14. However, whether magnetic characteristics of the Mid-Atlantic ridge may provide detectable orientation cues still needs to be documented. Similarly, the existence of detectable odours that might be generated by the tectonic activity and hydrodynamics of the Mid-Atlantic ridge and serve as orientation cues for eels is still unknown. Hydrodynamic mesoscale turbulence and vertical flows have been shown to be generated along the Mid-Atlantic Ridge15, which we propose eels might be able to detect. There are no well supported spawning areas of freshwater eels other than A. japonica and one north Pacific population of A. marmorata. The spawning areas of the other species remain unknown. In the south west Indian Ocean, spawning areas of 3 species (A. mossambica, A. marmorata and A. bicolor) were proposed on the east of the Mascarene Ridge with a similar topography (although shallower) than along the Mid-Atlantic Ridge and the Mid-Pacific ridge and seamounts16,17. Inaccurate spawning areas were also proposed for the South Pacific A. diffenbachii between Fiji, New Caledonia and New Zealand; in the vicinity of a number of oceanic ridges and trenches18 that may also serve as landmarks. Because all eel species studied on their spawning migration show similar diel vertical migration behaviours, it is likely that common orientation mechanisms could lead to detection of oceanographic variability related to the topography of the sea floor and related geomagnetism, local hydrodynamic turbulence and odour caused by vertical currents. This kind of oceanic landscape (chains of seamounts) occurs on narrow areas which strongly increase the meeting probability of spawners searching for partners and favourable spawning places.(Comment 2) Drift simulation with departures from the Mid-Atlantic Ridge and from the Sargasso Sea showed similar results. This is not surprising since the modelling of larval drift seems essentially just to reflect the slow westward drift prevailing both in the Sargasso Sea and Mid-Atlantic Ridge areas. The assumption of using the intersection of the Mid-Atlantic Ridge by the two thermal fronts as presumed spawning places seems to have little basis. There is no indication neither of one nor two temperature fronts at depths where leptocephali are found along a 45  W latitudinal section in the middle of the Mid-Atlantic Ridge area.(Response 2) We agree with the comments that the similar distributions between the departure from the Sargasso Sea and the Mid-Atlantic Ridge are expected, as they mainly reflect the ocean circulation. This is also what we wanted to address, if different departures could lead to similar distributions, either Sargasso Sea or Mid-Atlantic Ridge could be candidates for the spawning area. We also agree that many eel larvae were collected at the two fronts in the Sargasso Sea, but not near the Mid-Atlantic Ridge. However, if the departure from the Sargasso Sea and the Mid-Atlantic Ridge led to similar distributions after 720 days, they were not the result of westward current, but the cause of a relatively quiet ocean in the Sargasso Sea and its surrounding area (i.e. Fig. 1). Without prevailing current, small larvae were mainly transported by ocean dispersion, and would later be transported by the major currents that lie in the north (Azores Current), south (North Equatorial Current), and west (Gulf Stream) of the Sargasso Sea. So, we compared departures at 100 km from west and east of the Mid-Atlantic Ridge. Subtle differences occurred (figure below). V-larvae departing from the east of the ridge dispersed relatively less northward compared to larvae released 100 km at the west of the ridge (this figure and original paper). Secondly v-larvae released at the south east of the study area (red dots on the figure, right panel) disperse relatively less towards the Caribbean Sea than when released at the west (red dots of the figure, left panel). This suggests that the dispersion of European eel larvae is optimum in an area comprised between the Mid-Atlantic Ridge and the Sargasso Sea (our previous simulation in the original paper), and declines eastward of the Mid-Atlantic Ridge (present simulation below).Figure 1Distribution of v-larvae released departure at the west (left) and east (right) of Mid-Atlantic Ridge. The tracking method is the same as described in the paper, v-larvae were release within 100 km west and east of the ridge.Full size imageHanel et al. also indicate that the convergence front weakens from West, in the Sargasso Sea, to East above the ridge. We consider that this constitutes an additional argument that the Mid-Atlantic Ridge is indeed at the edge of the convergence zone at the first area of the Atlantic Ocean where currents and temperatures are favourable for reproduction of eels.(Comment 3) Elevated manganese (Mn) concentrations in the otolith cores of glass eels as a hint for successful spawning only in areas with volcanic activity based on observations of Martin et al.18. However, the results from Martin et al.19 were entirely misread, resulting in a mis-interpretation of the data.(Response 3) Based on Martin et al.19, we stated that higher concentrations of Mn were found in glass eels’ otoliths collected across European estuaries than in otoliths of leptocephali larvae sampled in the Sargasso Sea. We suggested that this was the indication that glass eels were born in areas where volcanic activity produces high loads of Mn and other metals. This formed one of the arguments supporting our hypothesis that Atlantic eels could spawn in the proximity to the Mid-Atlantic Ridge. Thanks to Reinhold Hanel and colleagues, we realized that Martin et al.19 in fact showed that concentrations of Mn were higher in the center of otoliths of leptocephali larvae than in those of glass eels collected along the European coasts. Consequently, this argument is no longer valid. Nonetheless, otolith microchemical fingerprints significantly differ between young leptocephali sampled in the Sargasso Sea in 2008 and glass eels collected in Europe, hence suggesting that they have distinct spawning areas19. These authors indicated that the incorporation of elements from the environment to the otoliths needed to be better understood, namely as stated by Hanel et al., because of physiological and environmental control such as temperature and salinity. In addition, they outline that the dynamics of elements from the sea floor to the subsurface is not well understood and could be slow. We totally share these conclusions that are well known facts, and that simply confirm that environmental characteristics (trace element concentrations, salinity and temperature) are responsible for the elemental signature of the central part of otoliths. Hanel et al. also state that the composition of otoliths are also controlled by elemental maternal transfer from the egg to the otoliths. We are aware of this fact that has been shown is other fish species. However, the laser ablations were performed after the first feed check where maternal influence is reduced and is overruled by environment18. This supports the idea that glass eels collected in Europe do not originate from the same environments as leptocephali captured in the Sargasso Sea.(Comment 4) Insufficient sampling efforts and a limited area coverage of recent surveys as a possible reason for “false negative” observations along the Mid-Atlantic Ridge. This statement does not recognize the investigations by Johannes Schmidt as well as earlier and later surveys in the Mid-Atlantic Ridge area. The ICES “Eggs and Larvae database” records a total of 48 A anguilla leptocephali caught within the area 15–29 N and 43–48 W, at 10 stations between 1913 and 1970.Thanks for pointing out that larvae have been caught near the Mid-Atlantic Ridge, in which larvae were not newly hatched because of their relatively large size (23–45 mm). Ocean currents were weak and could flow either eastward or westward in this region, indicating that the spawning could occur from west to east of the ridge, without considering swimming. Note that ocean currents could change directions, so that it was also possible to spawn near the ridge after been transported eastward and westward.The observed distribution of small larvae  More

Center for Microbial Communities, Department of Chemistry and Bioscience, Aalborg University, Aalborg, DenmarkMorten Kam Dahl Dueholm, Marta Nierychlo, Kasper Skytte Andersen, Vibeke Rudkjøbing, Simon Knutsson, Per H. Nielsen, Mads Albertsen & Per Halkjær NielsenEnvironmental Science Department, The Institute for Scientific and Technological Research of San Luis Potosi (IPICYT), San Luis Potosí, MexicoSonia ArriagaDepartment of Process, Energy and Environmental Technology, University College of Southeast Norway, Porsgrunn, NorwayRune BakkeCenter for Microbial Ecology and Technology, Ghent University, Ghent, BelgiumNico BoonInstitute for Water and Wastewater Technology, Durban University of Technology, Durban, South AfricaFaizal Bux & Sheena KumariVeolia Water Technologies AB, AnoxKaldnes, Lund, SwedenMagnus ChristenssonDepartment Of Chemical Engineering, Faculty of Engineering, University of Malaya, Kuala Lumpur, MalaysiaAdeline Seak May ChuaEnvironmental Engineering, Newcastle University, Newcastle, EnglandThomas P. CurtisThe Cytryn Lab, Microbial Agroecology, Volcani Center, Agricultural Research Organization, Rishon Lezion, IsraelEddie CytrynINGEBI-CONICET, University of Buenos Aires, Buenos Aires, ArgentinaLeonardo ErijmanDepartment of Biochemistry and Microbial Genetics, Biological Research Institute “Clemente Estable”, Montevideo, UruguayClaudia EtchebehereNIREAS-International Water Research Center, University of Cyprus, Nicosia, CyprusDespo Fatta-KassinosEnvironmental Engineering, McGill University, Montreal, QC, CanadaDominic FrigonSchool of Microbiology, Universidad de Antioquia, Medellín, ColombiaMaria Carolina Garcia-ChavesSchool of Civil and Environmental Engineering, Cornell University, Ithaca, NY, USAApril Z. GuWater Chemistry and Water Technology and DVGW Research Laboratories, Karlsruhe Institute of Technology (KIT), Karlsruhe, GermanyHarald HornDavid Jenkins & Associates, Inc, Kensington, CA, USADavid JenkinsInstitute for Water Quality and Resource Management, TU Wien, Vienna, AustriaNorbert KreuzingerWater Innovation and Research Centre, University of Bath, Bath, EnglandAna LanhamSingapore Centre of Environmental Life Sciences Engineering (SCELSE) Nanyang Technological University, Singapore, SingaporeYingyu LawWater Desalination and Reuse Center, King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi ArabiaTorOve LeiknesProcess Engineering in Urban Water Management, ETH Zürich, Zürich, SwitzerlandEberhard MorgenrothDepartment of Biology, Warsaw University of Technology, Warsaw, PolandAdam MuszyńskiEnvironmental Microbial Genetics Lab, La Trobe University, Melbourne, VIC, AustraliaSteve PetrovskiTechnologies and Evaluation Area, Catalan Institute for Water Research, ICRA, Girona, SpainMaite PijuanVA Tech Wabag Ltd, Chennai, IndiaSuraj Babu PillaiBiochemical Engineering Group, Universidade Nova de Lisboa, Lisboa, PortugalMaria A. M. ReisState Key Laboratory of Environmental Aquatic Chemistry, Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences, Beijing, ChinaQi RongWater Research Institute IRSA – National Research Council (CNR), Rome, ItalySimona RossettiLa Trobe University, Melbourne, VIC, AustraliaRobert SeviourDepartment of Civil and Environmental Engineering, University of Massachusetts Amherst, Amherst, MA, USANick TookerKemira Oyj, Espoo R&D Center, Espo, FinlandPirjo VainioEnvironmental Biotechnology, TU Delft, Delft, The NetherlandsMark van LoosdrechtVA Tech Wabag, Philippines Inc., Makati City, PhilippinesR. VikramanDepartment of Water Technology and Environmental Engineering, University of Chemistry and Technology, Prague, Czech RepublicJiří WannerEnvironmental Life Science Engineering, TU Delft, Delft, The NetherlandsDavid WeissbrodtSchool of Environment, Tsinghua University, Beijing, ChinaXianghua WenEnvironmental Biotechnology Lab, Department of Civil Engineering, The University of Hong Kong, Hong Kong, Hong KongTong Zhang More

Study population and song recordingsAll animal procedures were carefully reviewed by the Williams College IACUC (WH-D), the Bowdoin College Research and Oversight Committee (2009–18), and the University of Guelph Animal Care Committee (08R601), and were carried out as specified by the Canadian Wildlife Service (banding permit 10789D).We studied Savannah sparrows (Passerculus sandwichensis) at the Bowdoin Scientific Station on Kent Island, New Brunswick, Canada (44.5818°N, 66.7547°W). Since 1988, individuals nesting within a 10 ha study area in the middle of the island (30–70 pairs each year; part of a larger population of 350–500 males breeding on Kent Island and two adjacent islands) have been colour-banded to facilitate visual identification, and complete demographic information is available for birds on the study site (though not for the entire population) for the years 1989–2004 and 2009–2013. Because of strong natal and breeding philopatry51, birds hatched on the study site itself represent 40–80% of adult breeders in that area, and because of the systematic banding program, ages are known. Each year adds a new generation to the population, with yearlings making up approximately half of the adult breeding males. The birds banded and recorded on the study site are estimated to make up 10–20% of the Savannah sparrow population on Kent Island and two nearby islands.Details of the recording methods used in this study (covering the years 1980, 1982, 1988-9, 1993-8, and 2003–13) can be found elsewhere36,49. Using digitally generated sound spectrograms (using SoundEdit Pro and Audacity), birds were scored as having either a) high note cluster=a final introductory segment interval including at least two different note types, or b) a click train=one or more introductory segment intervals including at least two clicks and no other note types, or c) both features36 (see Supplementary Fig. 1 for a full description of note types). Although a small proportion of birds (mean = 8.3%) did not include either feature in their songs (such birds either had no feature in the introductory segment intervals or one non-click note type in the final interval), we did not include this option in the model and omitted these birds from summaries of the data. We did not include data after the breeding year 2013 because of we began an experimental field tutoring study in the summer of 201364.ModellingWe used a dynamic, discrete time model which allowed us to focus our analysis to specific time points within the year that are related to song learning (the beginning and end of the breeding season). These were: (1) the return of older birds between breeding seasons, (2) the recruitment of young birds singing newly crystallized songs in the spring, and (3) reproduction, resulting in the addition of juveniles during the summer breeding season.Because survival data were not available for every year during the time span we studied, we captured the variation in survival rates observed in the field57 by using a binomial distribution centered on the average historical survival rate for each age class (addressing the possibility that cultural drift resulting from random differences in survival rates was responsible for the shift in song features). The model incorporates stochasticity to capture the variation in population dynamics and return rates by assigning parameter values for survival and return rates from empirically generated probability distributions.We did not include spatial distribution of song variants in the model; although spatial patterns can be important for the dynamics of language loss58, territories with birds singing click trains and high note clusters were intermixed and no spatial structure was apparent (Fig. 3).The model assumes that males choose which features to incorporate into the introductory sections of their songs during song development. Individuals fall into one of six mutually exclusive classes of male Savannah sparrows. The classes are defined by (1) the bird’s developmental stage in the song learning process: juvenile (J, the first year, when the song is plastic) or adult (A, after the first spring, when the song is crystallized), and (2) the variant or variants sung as part of the bird’s introduction (high note clusters, click trains, or both). Denoting note high note clusters with X and click trains with C, the adult classes are therefore AX, AC, and AXC, and the juvenile classes are JX, JC, and JXC. The sum of the individuals in these classes is the total male population.We used two times during each year – late spring and late summer – to correspond to stages in song development (Fig. 5). At a given time t, when breeding is underway in the late spring, the male population consists entirely of adults singing crystallized song, and therefore each juvenile class is empty. At the end of the summer, the population of males has been augmented by juveniles, which are initially assigned to the same variant class as their fathers. To capture these dynamics, we define an intermediate time step, denoted ti. Time t + 1 then corresponds to the following breeding season (late spring), when juvenile males hatched the previous year have completed song development, crystallized their songs, and joined the adult class.Fig. 5: Model of song development.We used two age classes (J = juvenile and A = adult) and three classes of introductions (C = click trains, X = high note clusters, and  XC = both). In the late spring of a given year (time = t), only adult males are present. In late summer, those adults have bred and both they and juvenile males are present; at this intermediate time (ti) each male is initially allocated the same introduction type as his father (solid lines). Then, as song development progresses and juvenile males can be influenced by other tutors, they may retain their initial introduction type or switch to either of the other two types (dashed lines) before they crystallize their songs late in the following spring (time = t+1), and join the breeding cohort, which also includes adult males from the previous year who returned to breed again.Full size imageIn the late summer the male population increases with the addition of juveniles hatched that year, some of which will return to join the singing population the following year; survivors will return to breed within a few hundred meters of where they hatched51. To fit the observed historical decline in the Kent Island population57, the total number of returning juveniles, r (including both those hatched on site and those immigrating from nearby populations at time), follows a Poisson distribution where m = 33.6 – .182x and x is the number of years since 1980 (this function results in a decline of 5 males per decade; the initial number on the study site used in the model, 70, was extrapolated from historical data). The size of each returning juvenile class at time ti then takes the form:$${{{{{{rm{JY}}}}}}}_{{{{{{{rm{t}}}}}}}^{{{{{{rm{i}}}}}}}} sim {{{{{rm{Poisson}}}}}}left(mright)frac{{{{{{rm{A}}}}}}{{{{{{rm{Y}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}}{{{{{{rm{A}}}}}}{{{{{{rm{X}}}}}}}_{{{{{{rm{t}}}}}}}+{{{{{rm{A}}}}}}{{{{{{rm{C}}}}}}}_{{{{{{rm{t}}}}}}}+{{{{{rm{AX}}}}}}{{{{{{rm{C}}}}}}}_{{{{{{rm{t}}}}}}}}$$
(1)
for each Y ∈ {X, C, XC}.After the following winter, the proportion of surviving adults at time t + 1 follows a binomial distribution where the mean survival rate s = 0.48 is derived from historical data. Therefore, each adult class takes the form:$${{{{{rm{A}}}}}}{{{{{{rm{Y}}}}}}}_{{{{{{rm{t}}}}}}+1} sim {{{{{rm{Binomial}}}}}}left({{{{{rm{AY}}}}}},{{{{{rm{s}}}}}}right)* {{{{{rm{A}}}}}}{{{{{{rm{Y}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}$$
(2)
At the beginning of the next breeding season, juveniles complete song learning64, choosing which variant to crystallize as part of the song, and enter an adult song class; thus all of the juvenile classes disappear at t + 1. Which adult class juveniles join depends on separate learning functions for each of the two variants, ϕX for the high note cluster and ϕC for the click train. The ϕ function takes values between 0 and 1 and gives the probability of crystallizing a song form during the transition from natal year to breeding, depending upon the frequency-dependent bias and selection parameters (see below). These functions define the proportion of features that appear in the next generation as compared to that of the previous generation. Therefore we have:$${{{{{rm{A}}}}}}{{{{{{rm{X}}}}}}}_{{{{{{rm{t}}}}}}+1}={left({{{upphi }}}_{{{{{{rm{X}}}}}}}right)}^{2}{{{{{rm{J}}}}}}{{{{{{rm{X}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}+{left(1-{{{upphi }}}_{{{{{{rm{C}}}}}}}right)}^{2}{{{{{rm{J}}}}}}{{{{{{rm{C}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}+{{{upphi }}}_{{{{{{rm{X}}}}}}}left(1-{{{upphi }}}_{{{{{{rm{C}}}}}}}right){{{{{rm{JX}}}}}}{{{{{{rm{C}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}+{{{{{rm{A}}}}}}{{{{{{rm{X}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}$$
(3)
$${{{{{rm{A}}}}}}{{{{{{rm{C}}}}}}}_{{{{{{rm{t}}}}}}+1}={left(1-{{{upphi }}}_{{{{{{rm{X}}}}}}}right)}^{2}{{{{{rm{J}}}}}}{{{{{{rm{X}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}+{left({{{upphi }}}_{{{{{{rm{C}}}}}}}right)}^{2}{{{{{rm{J}}}}}}{{{{{{rm{C}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}+left(1-{{{upphi }}}_{{{{{{rm{X}}}}}}}right){{{upphi }}}_{{{{{{rm{C}}}}}}}{{{{{rm{JX}}}}}}{{{{{{rm{C}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}+{{{{{rm{A}}}}}}{{{{{{rm{C}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}$$
(4)
$${{{{{rm{A}}}}}}{{{{{{rm{XC}}}}}}}_{{{{{{rm{t}}}}}}+1}=2{{{upphi }}}_{{{{{{rm{X}}}}}}}left(1-{{{upphi }}}_{{{{{{rm{X}}}}}}}right){{{{{rm{J}}}}}}{{{{{{rm{X}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}+2{{{upphi }}}_{{{{{{rm{C}}}}}}}left(1-{{{upphi }}}_{{{{{{rm{C}}}}}}}right){{{{{rm{J}}}}}}{{{{{{rm{C}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}+({{{upphi }}}_{{{{{{rm{X}}}}}}}{{{upphi }}}_{{{{{{rm{C}}}}}}}left(1-{{{upphi }}}_{{{{{{rm{X}}}}}}}right)left(1-{{{upphi }}}_{{{{{{rm{C}}}}}}}right){{{{{rm{JX}}}}}}{{{{{{rm{C}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}})+{{{{{rm{A}}}}}}{{{{{{rm{XC}}}}}}}_{{{{{{{rm{t}}}}}}}_{{{{{{rm{i}}}}}}}}$$
(5)
The sum of probabilities defining all of song crystallization outcomes for the songs of fathers with song type X is:$${left({{{upphi }}}_{{{{{{rm{X}}}}}}}right)}^{2}+{left(1-{{{upphi }}}_{{{{{{rm{X}}}}}}}right)}^{2}+2{{{upphi }}}_{{{{{{rm{X}}}}}}}left(1-{{{upphi }}}_{{{{{{rm{X}}}}}}}right)=1$$
(6)
Learning curvesTo define how young males’ song learning is influenced by the songs they hear, we used learning curves based on type III Holling response curves59 which provide a means to numerically capture functional responses. In our model, the type III curve models the response of juvenile to the song form of adults in the population based on two variables: (1) frequency-dependent bias that favors one form based on its prevalence within the adult population, and (2) selection that favors a particular form of the song.The learning curves, ϕx for the high note cluster and ϕc for the click train, are modified forms of the type III Holling response curve):$${{{upphi }}}_{{{{{{rm{x}}}}}}}=frac{{x}^{{{{{{rm{beta }}}}}}}/{{{{{rm{sigma }}}}}}}{{(1-x)}^{{{{{{rm{beta }}}}}}}+({x}^{{{{{{rm{beta }}}}}}}/{{{{{rm{sigma }}}}}})}$$
(7)
and$${{{upphi }}}_{{{{{{rm{c}}}}}}}=frac{{{{{{rm{sigma }}}}}},{c}^{{{{{{rm{beta }}}}}}}}{{(1-c)}^{{{{{{rm{beta }}}}}}}+{{{{{rm{sigma }}}}}}{{c}}^{{{{{{rm{beta }}}}}}}}$$
(8)
where x is the proportion of the high note cluster within the population, c is the proportion of the click train within the population, β is frequency-dependent bias (favoring learning the novel or retaining the common variant), and σ is selection on the novel variant (a preference for learning the variant that is not dependent on frequency of the variant and includes factors such as prestige bias, success bias, status, and content bias). Note that the two learning curves do not have identical equations, because selection is not frequency-dependent. In these equations, β  > 1 corresponds to conformist selection, and when β  1 correspond to selection for a novel variant and values of σ  More

Epidemiological description of exhibition farm outbreakThe index farm where highly pathogenic avian influenza (HPAI) H5N1 virus in captive birds occurred was an exhibition farm in St. John’s, Province of Newfoundland and Labrador, Canada. The farm housed 409 birds of different species, namely chickens, guineafowl, peafowl, emus, domestic ducks, domestic geese, and domestic turkeys. On 9th December 2021, the farm owner first noticed mortality. On 13th December, the farm owner reported the increased mortality to a local veterinarian. Autopsies on four chickens showed swollen heads and cutaneous haemorrhages. Oropharyngeal and cloacal swabs from these chickens tested positive for H5 avian influenza virus at the Atlantic Veterinary College, University of Prince Edward Island, and the Canadian Food Inspection Agency (CFIA) was notified. On 16th December, by which time 306 birds (mostly chickens, turkeys and guineafowl) had died, staff of the CFIA collected tissue samples from dead chickens, as well as oropharyngeal and cloacal swabs and sera from different species of captive birds still present (Table 1), after which all remaining captive birds were culled. All oropharyngeal and cloacal swabs tested positive for HPAI virus of the subtype H5N1 by real-time RT-PCR, and all sera tested positive for influenza nucleoprotein antibodies by ELISA. On 20th December, the CFIA confirmed the diagnosis of HPAI of the subtype H5N1.Table 1 List of samples for virological and serological analysis collected by CFIA on 17 December 2021 from different species of captive birds still present at the farm.Full size tableWild birds were frequently observed co-mingling with the captive population. Captive birds except emus were allowed to move freely in and out of the open pens in which they were housed. At an on-site pond, domestic ducks were reported to mingle with free-living mallards (scientific names of wild birds in Table 2), and a snowy egret had been observed around 2nd to 6th December. Other wild birds reported on the farm were common starlings, feral pigeons, and unspecified gulls.Table 2 Common and scientific species names of the birds mentioned in the text.Full size tableRetrospectively, HPAI H5N1 virus also was identified in tissues of a great black-backed gull found at a nearby pond in St. John’s. The gull had been found ill on 26th November 2021 and taken to a local wildlife rehabilitation centre, where it died the following day.Phylogenetic analysisPhylogenetic analyses were performed to compare the genome sequences of the Newfoundland viruses from the exhibition farm birds and a great black-backed gull found nearby to other influenza viruses in the database. Based on BLAST analysis all eight gene segments of the virus had a Eurasian origin, and the virus was identified as a clade 2.3.4.4b H5N1 virus. Based on maximum likelihood and time-resolved trees inferred by use of whole genome sequences, the Newfoundland viruses had a shared common ancestor with European viruses circulating in early 2021 (Figs. 1, 2). The dates for the most recent common ancestor (MRCA) of all gene segments ranged from December 2019 to April 2021 (Table 3). There was no evidence that the Newfoundland viruses were genetically closely related to other current or recent viruses circulating in Europe. In contrast to currently circulating European viruses, the sequences of the Newfoundland viruses had no evidence of reassortment with other avian influenza viruses after ancestral emergence (Fig. 3). The virus from the great black-backed gull was highly similar to those from the exhibition farm, except for a small number of nucleotide differences in the neuraminidase (N) gene segment.Figure 1Maximum likelihood phylogenetic tree of the H5 HA gene. Relationships among the European 2021 H5 2.3.4.4b HPAI strains (magenta) and the Newfoundland wild bird and outbreak strains (red) are shown. The tree is rooted by the outgroup and nodes placed in descending order. Clades are collapsed for clarity.Full size imageFigure 2Maximum likelihood phylogenetic tree of the H5 gene segments. Relationships among the European 2021 H5 2.3.4.4b HPAI strains (magenta) and the Newfoundland wild bird and outbreak strains (red) are shown. The tree is rooted by the outgroup and nodes placed in descending order; order: HA, NA, PA, PB1, PB2, NP, MP, NS.Full size imageTable 3 Dates for the most recent common ancestor (MRCA) of all gene segments.Full size tableFigure 3Phylogenetic incongruence analyses. Maximum likelihood trees for the H and N gene segments and internal gene segments from equivalent strains were connected across the trees. Tips and connecting lines are coloured according to the legend.Full size imageAnalysis of avian migration and potential routes for HPAI H5 virus to be carried across the Atlantic with migrating birdsThere are several pathways for HPAI H5N1 virus to be carried across the Atlantic with migrating birds, based on the multitude of migration routes of wild birds and their overlapping ranges at breeding, stop-over, and wintering sites. Ring-recovery data confirm the regular movements of wild birds from Europe to Iceland and other North Atlantic islands, and from there to North America, and also provide evidence for direct movements of for example seabirds and gulls (Supplementary Table 1). Ringed individuals with a European origin have been found on Newfoundland for 10 of the 24 species in Supplementary Table 1: barnacle goose (1 ringed individual), Eurasian wigeon (5), Eurasian teal (1), great skua (13), European herring gull (1), black-headed gull (1), black-legged kittiwake (102), purple sandpiper (1), Brunnich’s guillemot (15), and Atlantic puffin (50). Given that the most likely ancestor of the virus detected in Newfoundland was circulating in Northwest Europe between the beginning of the 2020/2021 outbreak in Europe in October 2020 and April 2021 (see above), likely routes include spring migration of bird species moving to Icelandic, Greenlandic or Canadian High Arctic breeding grounds, or migration directly across the Atlantic Ocean (Fig. 4).Figure 4Maps of transatlantic migration. Putative virus transmission pathways between Europe and Newfoundland via migratory waterfowl/shorebirds (a) and pelagic seabirds (b). Many Icelandic waterfowl and shorebirds (a) winter in Northwest Europe and return to Iceland to breed in spring (1), including whooper swans, greylag geese, pink-footed geese, Eurasian wigeons, Eurasian teals, northern pintails, common ringed plovers and purple sandpipers. Some bird populations use Iceland as a stopover site, and continue to breeding grounds in East Greenland (2; barnacle geese and pink-footed geese), the East Canadian Arctic (3; light-bellied brent geese, red knots, ruddy turnstones) and West Greenland (4; greater white-fronted geese). Migratory birds from Europe share these breeding areas with species that winter in North America, including Canada geese and snow geese from East Greenland and the East Canadian Arctic (5), and some Iceland-breeding species of duck, including small numbers of Eurasian wigeons, Eurasian teals, and tufted ducks (6). Several seabird species (b), such as gulls, skuas, fulmars and auks, have large breeding ranges in the Arctic. After the breeding season many species become fully pelagic and can roam large parts of the northern Atlantic. The mid-Atlantic ridge outside Newfoundland is an important non-breeding area for seabirds, and is frequented by auks from Iceland (7), Svalbard (8) and Norway (9), including large numbers of Atlantic puffins and Brünnich guillemots, and by black-legged kittiwakes and northern fulmars originating from Iceland, Norway and the United Kingdom (7–8, 10). There these birds are joined by seabirds from Canadian and Greenlandic waters (11). Direct migratory links to Newfoundland occurs through greater and lesser-black backed gulls as well as black-headed gulls from Iceland and Greenland (12, 13), and gulls also link the pelagic and the coastal zone around Newfoundland (14). Thickness of the lines highlights the relative approximate population sizes. Dashed lines show where small numbers of individuals, or vagrants, provide a potential pathway. For more details on species and population numbers see Table 2. This figure was prepared using the software R (version 4.0.5, https://www.r-project.org/) and the following packages: -ggplot2 (version 3.3.5, https://cran.r-project.org/web/packages/ggplot2/index.html), -sf (version 1.0.5, https://cran.r-project.org/web/packages/sf/index.html).Full size imageThe first possible route via Iceland seems to be the strongest link between Newfoundland and Europe14,15,16,17, because it is a meeting point of breeding bird populations which winter in Europe as well as along the East coast of North America. Numerous species, totaling almost two million individual birds, migrate annually from northwestern Europe to breeding grounds in Iceland and beyond. Several populations breed on Iceland, including swans (whooper swan) (Supplementary Table 1), geese (greylag goose, pink-footed goose), ducks (Eurasian wigeon, Eurasian teal, Northern pintail), gulls (great- and lesser black-backed gull, black-headed gull, black-legged kittiwake) and shorebirds (common ringed plover, purple sandpiper, Supplementary Table 1). In addition, several species (e.g. barnacle geese and pink-footed geese) migrating to breeding grounds further away (Greenland and/or Canadian High Arctic) make spring and autumn stopovers in Iceland18,19. This creates potential for the virus to have been spread northwards to Iceland (or further northward) in spring, where it could have circulated among breeding birds, or transmitted during autumn migration by species returning from the Arctic. Several Iceland-breeding species of ducks (Eurasian wigeon, Eurasian teal, tufted duck), gulls (lesser black-backed gull, black-legged kittiwake, black-headed gull) and alcids (Brunnich’s guillemot, Atlantic puffin) winter along the Atlantic coast of North America in variable numbers (Supplementary Table 1). If the virus was transmitted to one of these populations during their stay on Iceland, it could have been spread to Newfoundland during the subsequent autumn migration. Importantly, Iceland-breeding Eurasian wigeons or Eurasian teals could be responsible for both the journey to Iceland from European wintering grounds, as well as the journey from Iceland to Newfoundland, where these species are frequently encountered as vagrants (Supplementary Table 1)20,21.The second possible route is via species that migrate from northwestern Europe to the Canadian High Arctic and/or Northwest Greenland. These include shorebirds (e.g. ruddy turnstone, red knot) and some geese (light-bellied brent goose and greater white-fronted goose). If the virus circulated in these breeding populations and then moved to other coastal marine bird populations bordering Baffin Bay, which include huge numbers of colonial seabirds and marine waterfowl22,23, the virus could have followed a coastal or even pelagic route south with the large autumn migration of Arctic marine birds (various sea ducks, auks and larids)24,25 to emerge in Newfoundland. Alternatively, shorebirds and waterfowl may have played a role: several European-wintering populations have overlapping breeding grounds with populations wintering along the east coast of North America. Regarding geese, greater white-fronted geese share breeding grounds in western Greenland with Canada geese26,27, which migrate south along the Canadian Atlantic coast. Also, brent geese have overlapping breeding grounds with snow geese18. In addition, individual barnacle geese and pink-footed geese breeding in Greenland could also have travelled south to Newfoundland carrying the virus, as these birds are regular vagrants to the North American Atlantic coast28. While geese occur only in small numbers on Newfoundland, two barnacle geese and four pink-footed geese, probably originating from Greenland breeding grounds, were observed in the autumn of 2021. St. John’s is the first major population center on a coastal route south from Baffin Bay/Davis Strait and along the Labrador Shelf, so emergence in eastern Newfoundland is consistent with this route.Three wild bird species involved in the Iceland and/or Greenland/High Canadian Arctic routes deserve particular attention. Eurasian wigeon have been prominently involved in outbreaks in Eurasia, and are considered prime candidates for carrying HPAI virus over long distances29. Also, during the first stages of an outbreak they are one of the first species to be detected HPAI virus positive, often without clinical signs. Barnacle geese and greylag geese, which congregate in Iceland, were in the top three most abundant species detected H5-positive in Europe in late winter and early spring 20215. Given that both greylag and barnacle geese have populations breeding on Iceland/Greenland and wintering in Europe (particularly the UK), these two species are high on the list of probable vectors that transported the virus to Iceland/Greenland and finally to Newfoundland. The high involvement of infected geese in the HPAI dynamics, which was not seen before October 2020, together with the unusually high levels of HPAI H5 virus presence in wild birds in Northwest Europe in spring 2021, might also explain why HPAI H5 virus spread to Newfoundland this winter (2021/2022), and not in the previous winters (2020/2021, 2016/2017, 2014/2015, 2005/2006). It is, however, striking that no cases of HPAI H5 virus have been recorded on Iceland in 2021.A third possible, pelagic, route is directly across the Atlantic Ocean. Such a route could have started with coastal and pelagic seabirds in Northwest Europe, where the virus may have remained undetected for much of the summer of 2021. A subsequent migration of seabirds to Newfoundland waters in the autumn of 2021 could have brought the virus to North America. The large populations of black-legged kittiwakes and northern fulmars that breed in the U.K. have long been known to frequent Newfoundland waters30, and these movements have been corroborated by recent telemetry studies31. Further, recent telemetry information reveals that millions of pelagic seabirds breeding all across the Atlantic congregate over the Mid-Atlantic Ridge in the central North Atlantic at all times of year32, making a pelagic transmission route a possibility. From the pelagic wintering grounds off Newfoundland, a species that uses both pelagic and coastal habitats, possibly a gull, may have brought the virus to shore in St. John’s. Trans-Atlantic transmission via seabirds has been suggested for LPAI viruses, including detection of mosaic Eurasian-North American viruses in gulls and alcids12,33,34,35.For the time period and geographical frame considered, HPAI-H5-positive species included ducks (Eurasian wigeon, mallard, common eider), geese (barnacle, greylag, brent, pink-footed and greater white-fronted goose), swans (whooper), gulls (black-headed, herring, lesser black-backed, great black-backed), and shorebirds (red knot, ruddy turnstone) (Supplementary Table 2). Of these 15 species, ringed individuals with a European origin have been recorded on Newfoundland for barnacle goose (1 ringed individual), Eurasian wigeon (5), great skua (13), and black-headed gull (1) (Supplementary Table 1). Ringed individuals with a European origin have been found on Newfoundland for 5 species which were found to be HPAI-H5-positive between October 2020 and April 2021, such as Barnacle Goose (1), Eurasian Wigeon (5), Great Skua (13), Black-headed Gull (1). These species might be considered to be possible carriers of HPAI H5 virus from Europe in late winter 2020/2021 or early spring 2021 partly or all the way to Newfoundland. However, given the incompleteness of sampling and the possibility of wild birds carrying HPAI virus subclinically, the involvement of other wild bird species in transatlantic virus transport cannot be ruled out.Having reached the Avalon Peninsula of Newfoundland via one of above routes, the virus may have spread further within the abundant local populations of ducks and gulls wintering in the city of St. John’s. The peridomestic populations of some of these species may be candidates for incursion of the virus into the farm in St John’s. More

Discovery of a novel archaeal lineage in wetland sedimentsTo examine archaeal community composition and function in a mangrove ecosystem, we analyzed metagenomic data from 13 sediment samples taken from mangrove wetlands in Techeng Island of Zhanjiang and Dongzhai Harbour of Haikou, China (Supplementary Fig. 1). De novo assembly of these sequencing data (60–120 Gbp for each sample) and genome binning resulted in 242 archaeal MAGs ( >70% complete; 85.4% AAI) detected in two metagenomes in IMG database generated from sediments of Lake Towuti, Indonesia (Supplementary Fig. 1 and Supplementary Table 2). Two additional related MAGs (TDP8 and TDP10, Table 1) encoding the complete Mcr complex were subsequently recovered from these metagenomes. For these MAGs (with exception of HK01M), the mcrABG operon and other genes related to methane metabolism were located on long contigs (≥11,476 bp) whose sequence composition features were consistent with their corresponding genomes (Supplementary Fig. 3), supporting the accurate assignment of these contigs to each MAG. The estimated genome size range for the seven MAGs recovered was 1.06–2.55 Mbp with total number of coding sequences ranging from 1151 to 3291. We examined vertical distribution of these MAGs in sediment cores of two sampling sites and found that their relative abundance increased gradually as depth increased from 15 to 100 cm (Supplementary text; Supplementary Fig. 4). Subsequent searches of public sequencing databases using the 16S rRNA and mcrA gene sequences annotated in these MAGs identified related species in freshwater lake sediments, hot springs, mangrove wetlands, rice paddy soils, hydrothermal vents, and deep-sea sediments distributed in different regions of the world (Supplementary text; Supplementary Table 3 and Supplementary Fig. 5).Phylogenomic analysis using 122 concatenated archaeal-specific marker proteins revealed that the seven MAGs and “Ca. M. tengchongensis” formed a distinct lineage that is sister to the order Nitrososphaerales (Fig. 1a and Supplementary Fig. 2b). Phylogenetic analyses of the 16S and 23S rRNA genes recovered from these MAGs supported the novelty of this lineage (Supplementary Table 4 and Supplementary Fig. 2a), with pairwise nucleotide comparisons of 16S rRNA genes revealing an identity of 79.1–87.3% to publicly available Nitrososphaeria genomes (Supplementary Table 5). The seven MAGs belonging to the novel lineage had an AAI of 44.0–52.3% to all other genomes of the Nitrososphaeria (Supplementary Table 6), further supporting their classification as a separate order [21, 22]. Collectively, these phylogenetic analyses indicate that these MAGs represent four different genera of the recently described family “Ca. Methylarchaceae” within a novel order—designated here as “Ca. Methylarchaeales” (Fig. 1a and Supplementary Fig. 2 and Supplementary Tables 5 and 6). H03B1, HK01M, HK01B, and HK02M1 represents one genus (69.7–80% AAI to other MAGs), HK02M2 represents the second (68.9–80% AAI to other MAGs), TDP8 and TDP10 represent the third (70.2–82.5% AAI), and “Ca. M. tengchongensis” represents the fourth (68.9–82.5% AAI); the former three genera are named here “Ca. Methanoinsularis”, “Ca. Methanoporticola”, and “Ca. Methanotowutia”, respectively.The “Ca. Methylarchaeales” are potentially methyl-reducing methanogens with b-type cytochromesAnnotation of the eight “Ca. Methylarchaeales” MAGs confirmed genes involved in archaeal methane metabolism (Supplementary Table 7 and Fig. 2), including those encoding the Mcr complex (mcrABG and auxiliary genes mcrCD), and the ATP-binding protein AtwA (component A2) required for Mcr activation [23]. The “Ca. Methylarchaeales” harbor genes for methane production from methanol and methylamines (mtaA, mtbA, mttB, mtbB, and mtmB) (Supplementary Table 7 and Fig. 2), suggesting that the “Ca. Methylarchaeales” have potential to perform methyl-reducing methanogenesis, as previously suggested for “Ca. M. tengchongensis” [7], and members of the orders Methanomassiliicoccales [15], “Ca. Methanofastidiosales” [19] and “Ca. Methanomethylicales” [4]. All of the “Ca. Methylarchaeales” MAGs encoded a tetrahydromethanopterin (H4MPT) S-methyltransferase subunit H (MtrH), and either a MtrX or MtrA, that are homologous to those of Methanosarcina barkeri (Supplementary Table 7). Phylogenetic analysis revealed that the “Ca. Methylarchaeales” MtrH subunits are more closely related to a MtrH (BP07_RS03240) of Methermicoccus shengliensis than to the MtrH subunits of Methanosarcina (Supplementary Fig. 6). It is likely that the “Ca. Methylarchaeales” MtrH may be involved in methyl transfer directly to H4MPT, as previously shown in M. shengliensis for utilization of methoxylated aromatic compounds [24]. The absence of a complete gene operon for Mtr complex suggests that the “Ca. Methylarchaeales” cannot use the CO2 reduction or aceticlastic pathway for methanogenesis.Fig. 2: Proposed metabolic pathways in the “Ca. Methylarchaeales”.Genes found in H03B1/HK01M/HK01B/ HK02M1 (blue-green dots), HK02M2 (pink dots), TDP8/TDP10 (green dots), and JZ-2-bin_220 (brown dots) or missing from all bins (gray) are indicated. Genes associated with these pathways and their full name are provided in Supplementary Table 7. MP Methanophenazine, Fd ferredoxin, b-cyt b-type cytochrome.Full size imageIn contrast to the “Ca. Methanomethylicales”, all genes for the Wood-Ljundahl pathway (WLP) and acetyl-CoA decarbonylase/synthase: CO dehydrogenases (ACDS/CODH) are also present in all the genomes (Supplementary Table 7 and Fig. 2). However, we did not identify the energy-converting hydrogenase complex and F420-reducing hydrogenase complex, both of which are required for the oxidation of the methyl groups to CO2 via the WLP [12]. This suggests that the “Ca. Methylarchaeales” cannot utilize the methylotrophic pathway for methanogenesis. Similar to methyl-reducing methanogens of the Methanonatronarchaeales [17], function of the defective WLP remains a mystery.The “Ca. Methylarchaeales” MAGs contain one or two copies of a gene encoding heterodisulfide reductase subunit D (HdrD) (Supplementary Fig. 7 and Supplementary Table 7), one of which was co-located with a b-type cytochrome gene (Fig. 3a and Supplementary Fig. 7), which is similar to the hdrDE operon of Methanosarcina barkeri [25]. The b-type cytochromes in the HdrDE-like complex of the “Ca. Methylarchaeales” are integral membrane proteins with five transmembrane helical segments that harbor a nitrate reductase gamma subunit domain (PF02665) (Fig. 3c and Supplementary Figs. 7 and 8). Sequence analysis of these b-type cytochromes revealed two histidine residues located in Helix 2 of these proteins in all the “Ca. Methylarchaeales” genomes, two histidine residues located in Helix 5 for H03B1, and single histidine and methionine residues located in Helix 5 for “Ca. Methanotowutia” and “Ca. Methanoinsularis” (Supplementary Fig. 7b and Fig. 3c). These residues are suggested to be involved in the binding of two heme groups [26], similar to the NarI of E. coli [27] and HdrE of M. barkeri [25]. It is assumed that the two heme groups ligated to histidine or methionine residues of Helix 1 and Helix 5 are on the periplasmic and cytoplasmic side of the membrane bilayer respectively, and are responsible for electron transfer. In addition, the hdrDE operon is adjacent to the mcrABDG operon in all the “Ca. Methylarchaeales” MAGS (Fig. 3a), supporting their role in methanogenesis for these microorganisms. Collectively, these findings strongly indicate that members of the “Ca. Methylarchaeales” are b-type cytochrome-containing methanogens that use the HdrDE complex to reduce the heterodisulfide CoM-S-S-CoB of Coenzymes M and B generated in the final step of methanogenesis [28] (Fig. 2).Fig. 3: Gene composition and structural model of HdrDE and VhtAGC complexes in the “Ca. Methylarchaeales”.a Gene composition of contigs/scaffolds containing the gene cluster of heterodisulfide reductase (HdrDE) complex. Genes related to methane metabolism are highlighted with red, blue, yellow, and cyan. The hdrDE complex gene cluster is always adjacent to mcrABDG operon. b Gene composition of methanophenazine-reducing hydrogenase (VhtAGC) complex. Genes for VhtAGC were collocated on the same contig/scaffold, forming a transcriptional unit. c Structural model of b-type cytochromes in HdrDE and VhtAGC complexes showing the proposed heme ligation.Full size imageWe identified a homolog of a 11-subunit NADH-quinone oxidoreductase complex in each “Ca. Methylarchaeales” genome (Supplementary Table 7) whose gene cluster resembles to the F420H2 dehydrogenase (Fpo) found in Methanosarcina [29] (Supplementary Fig. 9b). Phylogenetic analysis of the large subunit revealed that the “Ca. Methylarchaeales” complex is more closely related to the Fpo and Fpo-like complexes of Methanosarcinales and Methanomassiliicoccales than to group 4 [NiFe] hydrogenases (Supplementary Fig. 10). The absence of the typical [NiFe]-binding motifs in the catalytic subunit excludes the possibility that the complex is a group 4 [NiFe] hydrogenase (Supplementary Fig. 9a). In addition, the complex also lack the FpoF subunit required for binding and oxidation of F420H2 [15]. This suggests that this Fpo-like complex is unable to interact with F420H2, and instead may use reduced ferredoxin as an electron donor, similar to its proposed role for the Methanomassiliicoccales [15] and Methanosaeta thermophila [30]. In six MAGs from “Ca. Methanoinsularis”, “Ca. Methanoporticola”, and “Ca. M. tengchongensis”, genes for soluble methyl viologen-reducing hydrogenase/heterodisulfide reductase complex (MvhADG/HdrABC) and methanophenazine-reducing hydrogenase complex (VhtAGC) are missing. It is extremely unlikely that genes encoding all MvhADG/HdrABC and VhtAGC complex subunits are present in these near-complete genomes but were missed by sequencing. Thus, it is proposed that these microorganisms may use the Fpo-like complex directly to accept electrons from reduced ferredoxin, and subsequently channel these electrons to the HdrDE complex coupled to the reduction of CoM-S-S-CoB (Fig. 2), as shown previously for Methanosaeta thermophila [30]. The reduced ferredoxin may be produced by some unidentified hydrogenases or an unknown pathway. The H03B1 MAG also encodes a formate dehydrogenase subunit A gene (fdhA) co-located with a fdhB gene (Supplementary Table 7) and a putative b-type cytochrome with five transmembrane helices and a prokaryotic b561 domain (PF01292) binding two heme groups (Supplementary Fig. 11c) that is similar to FdhC of E. coli. “Ca. M. tengchongensis” contained fdhAB genes, with the fdhB gene collocated with a gene for a cytochrome b561 with four transmembrane helices and two heme groups (Supplementary Fig. 11b). It is likely that these microorganisms may be able to use formate dehydrogenase to reduce methanophenazine pool which could then transfer electrons to the membrane-bound HdrDE complex (Fig. 2). We identified a geranylfarnesyl diphosphate synthase homolog in each “Ca. Methylarchaeales” genome. Phylogenetic analysis revealed that these enzymes cluster together with the geranylfarnesyl diphosphate synthase of M. mazei, likely suggesting that the “Ca. Methylarchaeales” may be able to synthesize methanophenazine, as previously shown in M. mazei [31] (Supplementary Fig. 12).The “Ca. Methanotowutia” (TDP8 and TDP10) MAGs encode the small and large subunits for a [NiFe] active site-containing hydrogenase co-located with a gene for membrane-spanning b561 domain (PF01292) cytochrome b (Fig. 3b), which is similar to the operon of VhtAGC complex found in Methanosarcina with cytochromes [12]. The b-type cytochrome harbors five transmembrane helices, with histidine or methionine residues located in Helix 1, 2, 5 for the ligation of two heme groups (Supplementary Fig. 11a). It has been proposed that the VhtA is guided to the cell membrane with the help of twin-arginine signal peptide of VhtG and its [NiFe] active site faces periplasmic side [32, 33]. As a result, two H+ ions generated by H2 oxidation are released into the periplasm while two electrons are transferred to heme groups of VhtC through Fe-S clusters of VhtG [12, 34]. Furthermore, the electron carrier methanophenazine connects VhtAGC with HdrDE, and its reduction and reoxidation results in the release of two additional H+ ions into the periplasm (Fig. 2) [34, 35]. Altogether, four electrogenic protons are generated in the system, which can be used to drive the synthesis of one ATP via an archaeal A-type ATP synthase. The HdrDE complex that receives electrons from the methanophenazine can be used to reduce CoM-S-S-CoB (Fig. 2), enabling the coupling of methane production with energy conservation. This is the first report of a VhtAGC complex and an HdrDE complex found in an mcr-containing archaeal lineage outside the Euryarchaeota superphylum (Fig. 1) and indicates that “Ca. Methanotowutia” may be capable of performing H2-dependent methyl-reducing methanogenesis. The membrane-bound electron transport chain is more efficient than electron bifurcation that is used by methanogens without cytochromes [12].Sequence analysis revealed that key conserved residues of the McrA sequences of the “Ca. Methylarchaeales”, including the binding sites for F430 cofactors, coenzyme M, and coenzyme B [36], are the same as those in McrA sequences of members of the Euryarchaeota superphylum, with exception that the cysteine at site α452 is replaced with an alanine or serine (Supplementary Fig. 13 and Supplementary Table 8). Phylogenetic trees of concatenated and individual McrABG were reconstructed, showing that the “Ca. Methylarchaeales” encode canonical Mcr complexes that cluster with those of putative methane-metabolizing archaea and are divergent from those of short-chain alkane-oxidizing archaea (Fig. 4 and Supplementary Fig. 14). These results support the view that the “Ca. Methylarchaeales” metabolize methane.Fig. 4: Phylogeny of the Mcr/Mcr-like complex showing the relationship with their species tree.a Maximum-likelihood tree (IQ-TREE, LG + C60 + F + G) based on an alignment of concatenated McrABG/McrABG-like subunits from 167 archaeal genomes. The Mcr-like complex is found in short-chain alkane-oxidizing archaea. b Maximum-likelihood tree (IQTREE, LG + C60 + F + G) based on concatenated 122 archaeal-specific marker proteins using the same genomes with those of Mcr/Mcr-like tree. Ultrafast bootstraps values ≥95 are indicated with green filled squares.Full size imageWe also explored the possibility that the “Ca. Methylarchaeales” may be able to oxidize methane. In reported anaerobic methanotrophic archaea (ANME), methane oxidation is coupled to the reduction of several electron acceptors (nitrate, sulfate or metal oxides). Known ANME are predicted to utilize canonical terminal respiratory reductases or multi-heme c-type cytochromes (MHCs) to transfer electrons to a syntrophic partner microorganism [37], metal oxides [38, 39] or humics [40]. We could not identify any terminal reductases or MHCs in the “Ca. Methylarchaeales” genomes. Previous studies have hypothesized that formate or acetate might act as potential syntrophic electron carriers between methane-oxidizing archaea and their partners [41, 42], and members of the “Ca. Methylarchaeales” possesses the genetic potential for the production of these electron carriers. However, to our knowledge, these electron-transferring mechanisms have never been experimentally verified for ANME. Collectively, these analyses suggest that these “Ca. Methylarchaeales” are more likely methanogens, although empirical studies are required to confirm this.Similar to all described methanogens [15], the “Ca. Methylarchaeales” do not encode a complete tricarboxylic acid cycle, with citrate synthase, fumarase and succinate dehydrogenase absent from these MAGs. The “Ca. Methylarchaeales” lack a canonical pyruvate kinase for glycolysis (Supplementary Fig. 15 and Supplementary Table 7). However, pyruvate-water dikinase or pyruvate phosphate dikinase in gluconeogenesis may replace pyruvate kinase to catalyze the reversible interconversion of phosphoenolpyruvate and pyruvate, as shown in cultivated methanogens Methanomassiliicoccales [15]. The identification of sugar transport proteins and a variety of extracellular and intracellular carbohydrate-active enzymes (CAZymes) including glycoside hydrolases (EC 3.2.1.1 and 5.4.99.16) and glycosyltransferases (EC 2.4.1, 2.4.1.83, and 2.4.99.18, etc.) in the “Ca. Methylarchaeales” (Supplementary Fig. 15) suggests that they may be able to utilize sugars as an alternative carbon and energy source, as previously hypothesized for the “Ca. Methanomethylicales” and “Ca. Bathyarchaeia” [4, 8]. However, comparative genomics revealed that cultured methanogens that do not utilize sugars also encode similar proteins (Supplementary Fig. 15) [12, 13], and they may instead be involved in biosynthetic pathways. In addition, peptide and amino acid transporters, and enzymes related to peptide fermentation including extracellular peptidases, endopeptidases, 2-oxoglutarate ferredoxin oxidoreductase (kor), 2-ketoisovalerate ferredoxin oxidoreductase (vor), indolepyruvate ferredoxin oxidoreductase (ior), and pyruvate ferredoxin oxidoreductase (por) are present in both the “Ca. Methylarchaeales” and cultured methanogens (Supplementary Fig. 15). Nevertheless, to our best knowledge, peptide fermentation has never been reported in these isolated methanogens to date. Thus, the genes may be involved in assimilation and metabolism of amino acids in the “Ca. Methylarchaeales” and other newly discovered uncultured methanogens [4, 8, 12].Evolution of the b-type cytochrome-containing methanogensThe rapid increase in the number and diversity of MAGs has greatly expanded the known diversity and distribution of Mcr genes in archaea. To investigate the evolutionary history of the Mcr complexes in methanogens, we inferred the phylogeny of concatenated McrABG subunits based on all mcr-containing archaeal genomes available in public databases. In accordance with previous studies [43, 44], lineages in Class I and Class II methanogens within the Euryarchaeota superphylum appear congruent between McrABG and species trees while H2-dependent methylotrophic methanogens Methanomassiliicoccales and Methanonatronarchaeia, and methanotroph “Ca. Methanophagales” (ANME-1) are not (Fig. 4). The results were further supported by the phylogeny of the six conserved markers (m4–m9) in this (Supplementary Fig. 16) and previous studies [44]. These markers are solely present in archaea containing Mcr or Mcr-like complexes and suggested to be involved in activation, folding and assembly of Mcr subunits [44]. The Mcr genes of “Ca. Methanomethylicales” and “Ca. Korarchaeia” within the phylum Thermoproteota were previously suggested to be acquired via HGTs, since they are closely related with those of methylotrophic methanogens of the Euryarchaeota superphylum in McrABG tree [44]. However, analyses including our “Ca. Methylarchaeales” MAGs and several others with an Mcr complex revealed good congruence between the concatenated McrABG, m4-m9 genes, and the genome-based trees for the lineages within the Thermoproteota (including the “Ca. Methanomethylicales”, “Ca. Korarchaeia”, “Ca. Nezhaarchaeales”, and our “Ca. Methylarchaeales”; Fig. 4 and Supplementary Fig. 16) support vertical inheritance and evolution independent of the Euryarchaeota superphylum. Wide distribution of mcr genes in archaea (Supplementary Fig. 17 and Supplementary Table 9) and their congruence with the genome-based tree for many lineages within the Euryarchaeota superphylum and the Thermoproteota suggest that these genes likely have originated before the divergence of these two major archaeal lineages.Recently, amalgamated likelihood estimation (ALE) has been used to estimate presence probability of McrA in each internal node in a rooted archaeal species tree, supporting the presence of McrA with high confidence in the common ancestor of Class I and Class II methanogens, “Ca. Methanofastidiosales”/“Ca. Nuwarchaeales” in Euryarchaeota superphylum, as well as “Ca. Methanomethylicales”, “Ca. Korarchaeia”, and “Ca. Nezhaarchaeales” in the Thermoproteota [45]. Compared to the previous study [45], our ALE results support the presence of McrA with high confidence [presence probability (pp) >0.9] at the basal node of “Ca. Methanomethylicales”, “Ca. Nezhaarchaeales”, “Ca. Korarchaeia”, and “Ca. Methylarchaeales” in the Thermoproteota (Supplementary Fig. 17), suggesting an earlier origin of Mcr complex in Thermoproteota. The difference is likely attributed to the addition of “Ca. Methylarchaeales”. Confidence in evolutionary inferences from ALE analyses will require expansion of genome coverage of some of the poorly represented or yet-to-be-discovered Mcr-containing lineages. A previous study showed that an ancestral McrA sequence were more closely related to McrA from “Ca. Methanodesulfokores washburnensis” in the “Ca. Korarchaeia” compared to any other lineages [6], possibly supporting our inference that methane metabolism may have evolved relatively early in Thermoproteota.The b-type cytochrome in HdrDE complex belongs to the protein family of nitrate reductase gamma subunit (PF02665, NarI). Using all publicly available archaeal genomes, we found that the NarI domain-containing cytochromes (NarI-Cyt) are primarily used in three electron transfer complexes: HdrDE, dissimilatory nitrate reductase (NarGHI) [46], and sulfite reductase (DsrABCJKMOP). For the HdrDE and NarGHI complexes, the genes encoding the subunits are co-localized in archaeal genomes, each forming a transcriptional unit. However, in the Dsr complex, only a DsrK is co-localized with a DsrM (b-type cytochrome) while other subunits are usually not adjacent to the DsrKM but separated by few genes [6]. We examined distribution of the three complexes in archaea. A total of 101 genomes were found to encode these complexes (66 for HdrDE, 16 for Nar, 23 for Dsr), and they are distributed across the Euryarchaeota superphylum, Thermoproteota, and Asgardarchaeota (Supplementary Fig. 17 and Supplementary Table 9). Among these genomes, the HdrDE is found in methanogens and methanotrophs belonging to the class “Ca. Methanosarcinia”, the orders Methanomicrobiales and Methanonatronarchaeales, and in alkane-oxidizing archaea belonging to the orders Archaeoglobales, “Ca. Syntropharchaeales”, and Methanosarcinales (GoM-Arc1) (Supplementary Fig. 17). In Mcr-containing archaea outside of the Euryarchaeota superphylum, the complex is exclusively found in the “Ca. Methylarchaeales” (Fig. 1 and Supplementary Fig. 17).Phylogenetic analyses of the NarI-Cyt were conducted to investigate the evolution of these genes in archaea (Fig. 5a). The results showed that these cytochromes have experienced frequent horizontal gene transfer, especially DsrM. The DsrM sequences annotated in members of the Thermoproteota form a distinct cluster. In the cluster, Archaeoglobi and “Ca. Hydrothermarchaeota” DsrM branch far from their Euryarchaeota superphylum relatives, and have potentially gained their cytochromes from a member of the “Ca. Korarchaeia”. Similarly, the “Ca. Methanoperedenaceae” and Archaeoglobi might have acquired their NarI genes from a member of Thermoproteia. Congruence between the cytochrome and genome-based trees for members of the Thermoproteota suggest that these cytochromes might have evolved before the diversification of this phylum. We further inferred a gene tree using concatenated HdrDE complex (Fig. 5b). The topological structure of this tree exhibits high congruence with the genome-based tree for all lineages except the Methanonatronarchaeia, supporting an early presence of the complex in archaea. This suggestion is supported by ALE analyses which indicate the presence of NarI-Cyt with high confidence in the common ancestor of Thermoproteota (pp = 0.69) and in the common ancestor of “Ca. Halobacteriota” (pp = 0.70) (Supplementary Fig. 17).Fig. 5: Phylogeny of NarI-domain-containing b-type cytochromes and concatenated HdrDE complexes in archaea.a Phylogeny of NarI-domain-containing b-type cytochromes in archaea. b Phylogeny of concatenated HdrDE complexes in archaea. The maximum-likelihood trees of NarI-domain-containing b-type cytochromes (a (a’)) and concatenated HdrDE subunits (b (a’)) from representative archaea are inferred with IQ-TREE (LG + C60 + F + G, -bb: 10,000 for NarI-domain, 1000 for HdrDE). The b-type cytochromes comprising different enzyme complexes are indicated by different color dots (light red for HdrE, yellow for DsrM, and blue for NarI). HdrE heterodisulfide reductase E subunit, DsrM sulfite reductase M subunit, NarI dissimilatory nitrate reductase I subunit. The maximum-likelihood trees of a concatenated set of 122 archaeal-specific marker proteins using the same genomes as those of NarI-domain-containing b-type cytochrome tree (a (b’)) and HdrDE complexe tree (b (b’)), respectively. The trees were computed with IQ-TREE using LG + C60 + F + G model. These genomes or clades with Mcr complexes are marked by pink dots. The bootstrap support values ≥95 are indicated with green filled squares.Full size imageAs mentioned above, b-type cytochromes are classified into different protein families, and form part of many membrane-bound electron transfer complexes in bioenergetic pathways [47, 48]. Aside from HdrDE, Nar, and Dsr, such complexes also include Vht, Fdh, b6f complex, bc1 complex, and succinate dehydrogenase (Sdh). We examined the distribution of different families of b-type cytochromes in 416 representative archaea covering 41 orders or phyla of the Euryarchaeota superphylum, Thermoproteota, and Asgardarchaeota (Supplementary Fig. 17 and Supplementary Table 9). A total of 246 genomes contained these b-type cytochromes that were distributed across 23 archaeal lineages. In total, 11 of the 13 lineages of the Thermoproteota, and 11 of the 24 orders in Euryarchaeota superphylum, had b-type cytochrome, suggesting its pervasiveness in archaea. We conducted phylogenetic analyses of the b-type cytochromes from different families (Fig. 6a). The result indicates that cytochromes from Fdh and Sdh complexes form two large clusters. Within each cluster, lineages from Thermoproteota or the Euryarchaeota superphylum were essentially grouped together, suggesting that these cytochromes may have evolved before the divergence of these major archaeal lineages. The cluster of cytochromes of the b6f complex is close to those of the bc1 complex, consistent with the suggestion that bacterial cytochromes in bc1 complex might originate from cytochromes in b6f complex [48]. A phylogenetic analysis of concatenated VhtAGC showed clustering of lineages from Thermoproteota with Archaeoglobi (Fig. 6b), suggesting ancient exchanges of the Vht complex among these lineages. Taken together, these results support an early origin of b-type cytochromes in archaea. Previous studies also imply that some core enzymes for bioenergetic pathways, including membrane-integral b-type cytochrome, formate dehydrogenase, [NiFe]-hydrogenase, the Rieske/cytb complexes, and NO-reductases, were present in the Last Universal Common Ancestor of Bacteria and Archaea [48, 49].Fig. 6: Phylogeny of b-type cytochromes and concatenated VhtAGC complexes in archaea.a Maximum-likelihood tree of b-type cytochromes of representative archaea (NarI-domain-containing b-type cytochromes not included) inferred using IQ-TREE (the best model: cpREV + F + G4). Different families of b-type cytochromes are shown. vht methanophenazine-reducing hydrogenase complex, fdh formate dehydrogenase, sdh succinate dehydrogenase. b Maximum-likelihood tree of concatenated VhtAGC subunits retrieved from representative archaea, inferred with IQ-TREE using LG + C60 + F + G model. The bootstrap support values ≥95 are indicated with filled squares. Genomes or clades with Mcr complexes are marked by green filled dots. The number of sequences for branches is given in parenthesis. The pink branches represent members of Thermoproteota phylum while the black branches represent members of Euryarchaeota superphylum.Full size imageAs the heme is indispensable to b-type cytochrome [47], we also investigated distribution of its biosynthetic pathway in archaea. Although there are 11 genes involving in the heme biosynthesis, the three genes (Ahb-NirDH, Ahb-NirJ1, and Ahb-NirJ2), responsible for conversion from precorrin-2 to heme, are the key to this pathway. Thus, these three genes were used as markers denoting the presence of heme biosynthetic pathway. Among 41 archaeal lineages, 32 had this pathway including the “Ca. Methylarchaeales” (Supplemental text, Fig. 2, Supplementary Fig. 17 and Supplementary Table 9). Phylogenetic analyses reveal that these lineages from Thermoproteota largely cluster together for Ahb-NirDH (Supplementary Fig. 18). However, for Ahb-NirJ1 and Ahb-NirJ2, lineages from the Euryarchaeota superphylum, the Thermoproteota, and Asgardarchaeota are tangled up, suggesting frequent HGTs of these genes between these lineages. The wide distribution of this pathway across the Euryarchaeota superphylum, the Thermoproteota, and Asgardarchaeota (Supplementary Fig. 17 and Supplementary Table 9) suggests that a common ancestor may have been able to synthesize heme. This observation further supports the possibility of the early presence of b-type cytochromes in archaea.Here we described the discovery of the novel archaeal order “Ca. Methylarchaeales”, expanding known methanogen and archaeal diversity. Members of the lineage are methyl-reducing methanogens that can conserve energy via membrane-bound electron transport chains. The “Ca. Methylarchaeales” are globally distributed in anoxic lake and marine sediments, suggesting that they make an important contribution to global methane emissions. Our broader analyses suggest that methanogens who use b-type cytochrome-containing complexes to transfer electrons may have originated before diversification of Thermoproteota or “Ca. Halobacteriota” phyla based on a conservative estimation for the origin of McrA and NarI-Cyt genes in the ALE analysis. A previous study using molecular clock analyses to indicate that the diversification of Thermoproteota likely occurred in the early Archean Eon [45]. Archean oceans are thought to have been anoxic and contain abundant ferrous iron from hydrothermal volcanics [50, 51], which would have provided sufficient raw materials for heme synthesis by methanogens. In addition, CO2, H2, and organic compounds produced by volcanic activity are transported to the early oceans [52], which provides adequate carbon and energy sources for methanogenic growth. Compared to hydrogenotrophic methanogens using electron bifurcation, methanogens using the membrane-bound electron chain have a higher energy production efficiency and growth yield, providing an advantage for members of the “Ca. Methylarchaeales” described here.Taxonomic proposals“Ca. Methanotowutia igneaquae” (gen. nov., sp. nov.)Methanotowutia (Me.tha.no.to.wu’ti.a. N.L. pref. methano-, pertaining to methane; N.L. fem. n. Methanotowutia methanogenic organism named after the lake Towuti in Indonesia where members of the genus were first discovered).Methanotowutia igneaquae (ig.ne.a’quae. L. masc. adj. igneus, of fire; L. fem. n. aqua, freshwater, pertaining to freshwater habitats; N.L. gen. n. igneaquae from/of water of fire, referring to the volcanic lake environment). This organism is deduced to be able to use methylated compounds for methanogenesis. Representative genomes are near-complete bins TDP8 (Accession No. SAMN15658089) and TDP10 (Accession No. SAMN15658091) recovered from freshwater sediments in Lake Towuti in Indonesia with the latter the type genome for the species.“Ca. Methanoinsularis halodrymi” (gen. nov., sp. nov.)Methanoinsularis (Me.tha.no.in.su.la’ris. N.L. pref. methano-, pertaining to methane; L. fem. adj. insularis, from an island; N.L. fem. n. Methanoinsularis methanogenic organism from an island, specifically referring to Techeng Island in China where these microorganisms were discovered).Methanoinsularis halodrymi (ha.lo.dry’mi. Gr. masc. n. hals (gen. halos) salt; Gr. masc. n. drymos coppice; N.L. gen. n. halodrymi of salty woodland, referring to the mangrove wetland environment). This uncultivated microorganism is assumed to be able to perform methylotrophic methanogenesis. The type genome for the species is the bin H03B1 (Accession No. SAMN15658086) recovered from mangrove wetlands in Techeng Island in China.“Ca. Methanoinsularis haikouensis” (gen. nov., sp. nov.)Methanoinsularis haikouensis (hai.kou.en’sis. N.L. fem. adj. haikouensis, pertaining to Haikou). This uncultivated microorganism is assumed to be able to perform methylotrophic methanogenesis. Representative genomes are the bins HK01M, HK01B, HK02M1 (Accession No. SAMN25131447, SAMN25131448, SAMN25131449) recovered from mangrove wetlands in Dongzhai Harbour in Haikou, China.“Ca. Methanoporticola haikouensis” (gen. nov., sp. nov.)Methanoporticola (Me.tha.no.por.ti’co.la. N.L. pref. methano-, pertaining to methane; L. masc. n. portus, harbour; L. suff. -cola (from L. masc. or fem. n. incola), inhabitant, dweller; N.L. masc. n. Methanoporticol, a methane-forming dweller of a harbor, specifically referring to Dongzhai Harbour in China where these microorganisms were discovered).Methanoporticola haikouensis (hai.kou.en’sis. N.L. masc. adj. haikouensis, pertaining to Haikou). This uncultivated microorganism is assumed to be able to perform methylotrophic methanogenesis. The type genome for the species is the bin HK02M2 (Accession No. SAMN25131450) recovered from mangrove wetlands in Dongzhai Harbour in Haikou, China.“Ca. Methylarchaeales” (ord. nov.)Methylarchaeales (Me.thyl.ar.cha.ea’les. N.L. neut. n. Methylarchaeum (Candidatus) type genus of the order; -ales, ending denoting an order; N.L. fem. pl. n. Methylarchaeales, the order of the genus “Ca. Methylarchaeum”); Methylarchaeaceae (Me.thyl.ar.chae.a.ce’ae. N.L. neut. n. Methylarchaeum (Candidatus) type genus of the family); -aceae, ending denoting a family; N.L. fem. pl. n. Methylarchaeaceae, the family of the genus “Ca. Methylarchaeum”). More

Experimental designData collection took place in different sites disseminated along the mainland French coastline in sectors dedicated to Pacific oyster farming. Over the years, the number of sites monitored varied from 43 sites until 2009, to 13 between 2009 and 2013, and finally to 8 sites since 2015. Here, we focus on 13 sites (Fig. 1 & Table 1) that were almost continuously monitored since 1993. All these sites stand in tidal areas except Marseillan, located in the Mediterranean Thau lagoon, for which tidal variations are only tenuous and Men-er-Roué which is in subtidal deep-water oyster culture area in the Bay of Quiberon. Sentinel oysters were reared in plastic meshed bags fixed on iron tables, mimicking the oyster farmers practices. In Marseillan, half-grown oysters were cemented onto vertical ropes (from 1993 to 2007 and from 2015 to 2018), reared in Australian baskets (from 2008 to 2011), or put in bags fixed on iron tables (2012, 2013, 2014). As for spat oysters, they were reared in pearl-nets between 2008 and 2011 or put in bags since 2012.Fig. 1Site locations (coordinates in WGS84) along the French coastline. The site numbers refer to Table 1.Full size imageTable 1 Site identification and coordinates in WGS84.Full size tableDuring the 1993–2013 period, at the beginning of each annual campaign, one batch of diploid spat (three in 2012 and 2013) and one batch of diploid half-grown oysters were bought from an oyster farmer (i.e., wild-caught individuals) and then deployed simultaneously on all sites of the monitoring network. Here, the term “batch” designates a group of oysters born from the same reproductive event (spatfall or hatchery cohort), having experienced strictly the same zootechnical route. One batch could eventually be reared in several different bags (up to 3) deployed in the same site. Different batches were never mixed in the same bag.During the 2009–2013 period, up to three additional batches of triploid spat were bought in commercial hatcheries and included in the survey strategy (for a maximum of 6 batches of spat per site in 2012 and 2013). In 2009, the batches that were bought had already been exposed to a first wave of mortality before being followed by the network. Thus, the data collected this year should be interpreted with caution. Since 2014, the origin of spat and half-grown oysters has changed notably to better control the initial health status of oysters (no contact with the natural environment before deployment in all sites). The hatchery facility of Ifremer-Argenton now produces the sentinel diploid spat used in the monitoring network (one batch for all sites per campaign), whereas, the half-grown oysters was composed of spat reared on the same location the previous year but not monitored.Data collectionAfter the deployment of the different batches at the beginning of the campaign (seeding dates from February to April depending on the year), growth and mortality were longitudinally monitored yearly. Until 1999, annual campaigns usually ended in the winter of the year the monitoring began (i.e. in December), whereas, during the period 2000–2018, all sites frequently extended the campaign to end in the winter (February to March) of the following year.Observations were collected on each site quarterly until 2008 but then monthly to bimonthly depending on the season. At each sampling date, local operators carefully emptied each bag in separate baskets, counted the dead individuals (those with open or empty shells) and alive ones, and removed the dead individuals. Then local operators weighed all alive individuals in each basket (mass taken at the bag level, protocol mainly used between 1993 and 1998 and since 2004) and/or collected 30 individuals to individually weigh them in the laboratory (mass taken at the individual level, protocol used between 1995 and 2010 for spat and since 1996 for half-grown oysters).Data cleaningDuring the 2009–2013 period, several batches of spat were monitored per site and campaign. Some had a similar background to the batches monitored before 2009 (i.e., wild-caught spat from natural spatfall collected in the bay of Arcachon). To ensure the continuity of the time-series, we thus decided to remove all mass and mortality data of spat that did not originate from natural spatfall in the Bay of Arcachon, as well as triploid spat bought in hatcheries (see Table 2 for the origin and number of batches kept per site and campaign). To ensure that the life-cycle indicators are as comparable as possible between campaign and site (i.e. estimated in a common restricted time window), we removed data collected after December 31 of the year the monitoring began, as well as the site × campaign combinations when monitoring ended before October because the growth or mortality could still be in the exponential phase during this end-of-follow-up periods26. As the protocol of mass data collection changed over the years, we could not only use the mass data taken at the bag level or that at the individual level without greatly breaking the continuity of the time-series. We thus kept data taken at the individual level until 2008 and those taken at the bag level since 2009. We then checked for nonsense or missing data (e.g., the mass of a bag was equal to 0 or missing although they were still alive oysters in the bag), duplicated values and removed data for bags not part of the protocols or incorrectly identified. Finally, we removed site × campaign combinations for which we had fewer than four mass or mortality data because more data is necessary to study the temporal pattern of growth and mortality.Table 2 Origins of the different oyster batches retained after data cleaning.Full size tableData processing and analysisAt this point, the available data were, therefore, the number of living individuals per bag, the number of dead individuals since the last visit, the individual mass (g) of oysters (until 2008) and the total mass (g) of the living individuals per bag (since 2009).For mass data collected until 2008, we calculated the mean of the individual mass per date × site × age class combination by averaging the mass of the individuals. In other cases (mass data collected since 2009), we calculated the mean mass of individuals for each bag × date × site × age class combination by dividing the total mass of living oysters by the number of living individuals and then averaged data by date × site × age class combination. Our mass data, hereafter called mean mass data, is thus composed of the mean of the individual mass until 2008 and the mean mass of individuals since 2009.For mortality data, we could not calculate a cumulative mortality per bag × date × site × age class combination as (1-frac{number;of;alive;oysters;at;sampling;date}{number;of;oysters;at;previous;sampling;date}) because the total number of oysters (dead and alive) on a specific date often differed from the number of alive oysters at the previous date (e.g., because oysters were lost from the bags, or were sampled for complementary analyses such as pathogen detection). We thus took into account changes in oyster numbers between visits and calculated cumulative mortality using the following formula: CMt = 1 − ((1 − CMt-1) × (1 − IMt)). CMt = Cumulative mortality at time t; CMt-1 = Cumulative mortality at time t-1; IMt = Mortality rate at time t. IMt was obtained by dividing the number of dead oysters by the sum of alive and dead oysters at time t. When several bags were followed, we then averaged the cumulative mortality per date × site × age class combinations.We modeled the evolution of the mean mass and cumulative mortality data as a function of time to cope with changes in data frequency acquisition during annual monitoring campaigns. According to previous studies, annual mortality and growth curves in C. gigas follow a sigmoid curve11,26. Therefore, we fitted a logistic model, Eq. (1), and a Gompertz model, Eq. (2), which correspond to the most commonly used sigmoid models for growth and other data27, to describe Yt = mean mass (in grams) and cumulative mortality at time t.$${Y}_{t}=frac{a}{left(1+{e}^{left(-btimes left(t-cright)right)}right)}$$
(1)
$${Y}_{t}=atimes {e}^{left(-eleft(-btimes left(t-cright)right)right)}$$
(2)
These equations estimate three parameters: the upper asymptote (a), the slope at inflection (b), and the time of inflection (c).As the mean mass of half-grown individuals at the beginning of the campaign was higher than 0, we also fitted a four-parameter version of the logistic model, Eq. (3), and Gompertz model, Eq. (4), which is commonplace in the growth-curve analysis of bacterial counts27, and estimated (d) which represents the lowest asymptote of the curve. This parameter also moves the model curve vertically without changing its shape. The upper asymptote thus becomes equal to d + a.$${Y}_{t}=d+frac{a}{left(1+{e}^{left(-btimes left(t-cright)right)}right)}$$
(3)
$${Y}_{t}=d+atimes {e}^{left(-eleft(-btimes left(t-cright)right)right)}$$
(4)
Model fitting was carried out using non-linear least squares regressions (R package nls.multstrat28). This method allows running 5000 iterations of the fitting process with start parameters drawn from a uniform distribution and retaining the fit with the lowest score of Akaike Information Criterion (AIC). The sigmoid curve (i.e. logistic or Gompertz) with the lowest mean AIC of all models was selected as the best curve describing the data (see technical validation section). More