Differential regulation of degradation and immune pathways underlies adaptation of the ectosymbiotic nematode Laxus oneistus to oxic-anoxic interfaces

Eukaryotes may experience oxygen deprivation under both physiological and pathological conditions. Because oxygen shortage leads to a reduction in cellular energy production, all eukaryotes studied so far conserve energy by suppressing their metabolism. However, the molecular physiology of animals that naturally and repeatedly experience anoxia is underexplored. One such animal is the marine nematode Laxus oneistus. It thrives, invariably coated by its sulfur-oxidizing symbiont Candidatus Thiosymbion oneisti, in anoxic sulfidic or hypoxic sand. Here, transcriptomics and proteomics showed that, whether in anoxia or not, L. oneistus mostly expressed genes involved in ubiquitination, energy generation, oxidative stress response, immune response, development, and translation. Importantly, ubiquitination genes were also highly expressed when the nematode was subjected to anoxic sulfidic conditions, together with genes involved in autophagy, detoxification and ribosome biogenesis. We hypothesize that these degradation pathways were induced to recycle damaged cellular components (mitochondria) and misfolded proteins into nutrients. Remarkably, when L. oneistus was subjected to anoxic sulfidic conditions, lectin and mucin genes were also upregulated, potentially to promote the attachment of its thiotrophic symbiont. Furthermore, the nematode appeared to survive oxygen deprivation by using an alternative electron carrier (rhodoquinone) and acceptor (fumarate), to rewire the electron transfer chain. On the other hand, under hypoxia, genes involved in costly processes (e.g., amino acid biosynthesis, development, feeding, mating) were upregulated, together with the worm’s Toll-like innate immunity pathway and several immune effectors (e.g., bactericidal/permeability-increasing proteins, fungicides). In conclusion, we hypothesize that, in anoxic sulfidic sand, L. oneistus upregulates degradation processes, rewires the oxidative phosphorylation and reinforces its coat of bacterial sulfur-oxidizers. In upper sand layers, instead, it appears to produce broad-range antimicrobials and to exploit oxygen for biosynthesis and development.

Stable transcriptional profile under hypoxic or anoxic sulfidic conditions. To understand the molecular mechanisms underlying L. oneistus response to oxygen, we subjected it to various oxygen concentrations. Namely, nematode batches were incubated under either normoxic (100% air saturation; O), hypoxic (30% air saturation; H) or anoxic (0% air saturation; A) conditions for 24 h. Additionally, given that L. oneistus thrives in reduced sand containing up to 25 µM sulfide 16,17 , we also incubated it in anoxic seawater supplemented with < 25 µM sulfide (anoxic sulfidic condition; AS; see 16 for the sulfide concentration experienced by L. oneistus during each incubation).
Whereas transcriptional differences of the symbiont Candidatus Thiosymbion oneisti incubated under normoxic (O) and hypoxic (H) conditions were negligible 16 , the expression profiles of nematode batches incubated under O conditions varied so much that they did not cluster (Fig. S1). Consequently, there was no detectable differential expression between the transcriptomes of O nematodes and any of the other transcriptomes (H, A or AS; Fig. S1B, C). We attribute the erratic transcriptional response of L. oneistus to normoxia to the fact that this concentration is not naturally experienced by L. oneistus 16,17 .
As for the expression profiles of nematodes subjected to the H, A or AS conditions, replicates of each condition behaved more congruently (Fig. S1B). However, we did not find any significant difference between the A and AS nematodes and only 0.05% of the genes (8 genes; Data S1) were differentially expressed between the H and A nematodes. Moreover, even at the proteome level, there was no significant difference between the H and A incubations (t-test, Benjamini-Hochberg correction, p < 0.05; Fig. S2A, Data S1). Finally, only the comparison of the AS and H transcriptomes resulted in 4.8% of the expressed genes (787 out of 16,526) being differentially expressed, with 434 genes being upregulated under the AS condition and 353 being upregulated under the H condition (Fig. S1C, Data S1).
Collectively, our data suggests that L. oneistus may be ill-equipped to handle normoxic sediment, but it maintains a largely stable physiological profile under all other conditions. Before discussing the subset of biological processes differentially upregulated in AS versus H nematodes and vice versa, we will present the physiological processes the worm appears to mostly engage with, irrespective of the environmental conditions we experimentally subjected it to.  19 and comprehensive literature search (Fig. 1, Data S2). Our manual classification was supported by automatic eggNOG classification (Data S2). Similarly, the H and A proteomes were pooled, and the 100 most abundant proteins out of 2,626 were detected (Fig. S2). Based on median gene expression values of the top 100 expressed genes, we found that some of the processes L. oneistus mostly engages with were ubiquitination (ubq-1 20 ), energy generation (e.g., globin glb-1-like 21  Lastly, 48 out of the top 100 most expressed genes, were also detected among the top 100 proteins (Fig. 1,  Fig. S2, and Data S2, Supplemental material). Despite the modest Spearman correlation between transcript and protein levels (ρ = 0.4, p-value < 0.01) (Fig. S3A), there was an overlap in the detected biological processes (e.g., energy generation, stress response or detoxification categories, carbohydrate metabolism, cytoskeleton, locomotion, nervous system) (Figs. S2, S3B). All in all, except for those encoding for immune effectors, top-transcribed L. oneistus genes could not be unambiguously ascribed to its symbiotic lifestyle. This differs to what has been observed for other chemosynthetic hosts, such as giant tubeworms and clams. It is perhaps because these animals acquire their endosymbionts horizontally and feed on them (as they are housed intracellularly) that they abundantly express genes involved in symbiont acquisition, proliferation control and digestion [24][25][26] . Notably, we did observe a partial overlap of the most expressed gene categories (e.g., oxidative stress, energy generation, immune response), when L. oneistus was compared to the marine gutless annelid Olavius algarvensis. We ascribe the overlap to the fact that, albeit endosymbiotic, O. algarvensis also inhabits shallow water sand (Fig. S4, Supplemental material) and, as hypothesized for L. oneistus, it may also acquire its symbionts vertically [27][28][29] .
To conclude, although both symbiont 16 and host transcriptomics do not suggest a high degree of inter-partner metabolic dependence in the L. oneistus ectosymbiosis, the nematode seems well-adapted to both anoxic sulfidic (AS) and hypoxic (H) sand (Fig. 2, Data S1). The transcriptional response of the worm to these two conditions is, however, significant (Fig. 2, Data S1), and it will be reported below.
Genes upregulated in anoxic sulfidic (AS) nematodes. Chaperones and detoxification. The expression of chaperone-encoding (e.g., hsp12.2, grpE, dnaJ/dnj-2, pfd-1, pfd-6 [30][31][32] ) and ROS-detoxification-related genes (e.g., superoxide dismutase sod-2 and a putative glutathione peroxidase, involved in the detoxification of superoxide byproducts and hydrogen peroxide, respectively 33 ) were higher in AS nematodes (Figs. 2 and 3). Notably, transcripts encoding for the heme-binding cytochrome P450 cyp-13B1 were also more abundant in AS (Fig. 3), perhaps to increase the worm's capacity to cope with putative ROS formation 34 . Indeed, as cells start being oxygen-depleted, mitochondrial ROS accumulate because of the inefficient transfer of electrons to molecular oxygen 11,12 . Alternatively, the upregulation of antioxidant-related genes in AS worms could represent 5 32 16 1 1 10 7 9 3 2 21 7 5 24 7 2 28 4 4 4 43 6 11 10 Figure 2. Median gene expression levels of selected L. oneistus metabolic processes among the differentially expressed genes between the hypoxic (H) and anoxic sulfidic (AS) conditions after 24 h. Individual processes among the differentially expressed genes are ordered according to their difference in median expression between the AS and H incubations. Namely, detoxification (far left) had the largest difference in median expression in the AS condition, whereas immune response (far right) had the largest median expression difference in the H condition. The absolute number of genes are indicated at the top of each process. Metabolic processes were manually assigned and confirmed with the automatic annotated eggNOG classification. For specific gene assignments see Data S1. Some genes are present in more than one functional category and processes comprising only one gene are not displayed in the figure but listed in Data S1. www.nature.com/scientificreports/ an anticipation response to an imminent reoxygenation. In animals alternating between anoxic and oxygenated habitats, the re-exposure to oxygen can be very dangerous, as it creates a sudden ROS overproduction that may overwhelm the oxidative defense mechanisms 1

Amino acids degradation
Heat schock   www.nature.com/scientificreports/ sion of ROS-counteracting genes is consistent with what has been reported for vertebrates and marine gastropods which, just like L. oneistus, alternate between oxygen-depletion and reoxygenation 1 .
Mitochondrial and cytoplasmic ribosome biogenesis. In the cellular stress imposed by oxygen deprivation, mitochondria are central to both death and survival (reviewed in 7 ). In this scenario, calcium regulation, the scavenging of ROS or the suppression of their production, and/or inhibition of the mitochondrial permeability transition pore (MPTP) opening, might help to preserve mitochondrial function and integrity 7,35 . In addition, removal of specific mitochondrial components (mitochondrial-associated protein degradation, MAD), might also arise to maintain the overall mitochondrial homeostasis 36 . Perhaps as a response to anoxia-induced stress (reviewed in 7 ), a gene involved in MAD (vms-1) 36 was upregulated in AS worms (Fig. 4). More abundant in this condition were also transcripts encoding for mitochondrial transmembrane transporters tin-44, slc-25A26 and C16C10.1 (UniProtKB O02161, Q18934, Q09461), putatively transporting peptide-containing proteins from the inner membrane into the mitochondrial matrix, such as S-Adenosyl methionine (Fig. 6). Surprisingly, although the translation elongation factor eef-1A.2 37 was downregulated in AS worms, not only various mitochondrial ribosome structural components (28S: mrps, 39S: mrpl) 38 , and mitochondrial translation-related genes (e.g., C24D10.6 and W03F8.3) 39 were upregulated in AS nematodes, but also several cytoplasmic ribosome biogenesis (40S: rps, 60S: rpl) 40 and subunit assembly genes (e.g., RRP7A−like) 41 (Fig. 4). Taken together, the maintenance of mitochondrial homeostasis, an anticipatory response to a potential upcoming ROS insult (see "Chaperones and detoxification" section) and/or their involvement in extra-ribosomal functions 42,43 might explain the upregulation of ribosomal biogenesis-related genes in AS nematodes. Although upregulation of ribosomal proteins has also been observed in gastropods exposed to anoxia 44 , increased ribosomal biogenesis (which oftentimes directly correlates with an increase in protein synthesis) is not expected in animals that must repress their metabolism to cope with oxygen deprivation 3 .
Energy generation. Equally surprising was the upregulation of all differentially expressed genes related to energy generation in AS nematodes (Fig. 4). Namely, besides putative oxygen-binding globulin-like genes (e.g., glb-1, glb-14) 21 , the following were upregulated in AS nematodes: key structural genes (e.g., atp-3, atp-5) 45 , assembly-related genes (H + -transport ATP synthase) 46 of the mitochondrial ATP synthase (complex V), genes related to complex I (lpd-5, nuo-2) 47,48 , a subunit of the succinate dehydrogenase involved in complex II (mev-1) 49 , a mitochondrial cytochrome c oxidase subunit II assembly gene related to complex IV (sco-1) 50 , and a mitochondrial gene (coq-5), involved in the synthesis of either ubiquinone (Q, aerobic) or rhodoquinone (RQ, anaerobic) electron carriers 51 (Fig. 4). This suggests that, under anoxia, the electron transfer chain (ETC) is rewired: the electrons still enter the ETC at complex I, but do not reach complex III and IV. Instead, complex II, which in aerobic conditions functions as a succinate dehydrogenase, is repurposed into a reductase, which shuttles electron from the electron donor rhodoquinone to the electron acceptor fumarate. This mechanism would maintain the flow of electrons through the ETC, and it would prevent mitochondrial ATP generation (complex V) from shutting down 51,52 .
In short, under AS, similarly to what has been observed in other free-living and parasitic nematodes, complex I appears to be the sole proton pump in this truncated form of ETC 51,52 . In accordance with this hypothesis, tryptophan (Trp) degradation-related genes (acsd-1, acsd-2) and the Trp RNA ligase (wars-1) 53 that might be required to synthesize RQ 51,52 were upregulated under AS. Intriguingly, an isocitrate dehydrogenase gene (idh-1) was also upregulated. This produces reducing equivalent (NADPH) carrying electrons that may fuel complex I 54 , but it might also add to the stimulation of the antioxidant capacity or to the maintenance of redox homeostasis by regenerating reduced glutathione 1,55 .
If glycolysis is a key process for ATP generation in anoxia 3 and if, consistently, hxk-2 was upregulated under this condition (Fig. 6), based on the expression levels of transcripts encoding for alpha-amylases (see "Carbohydrate metabolism" in Fig. 6), starch and/or glycogen 56 may be the prominent carbon sources under anoxic sulfidic conditions. Ubiquitin-proteasome system and proteases. Proteolysis supplies amino acids or polypeptides to the cells, while impeding the accumulation of damaged or misfolded proteins. The two main mechanisms of cellular proteolysis are the lysosome-mediated intracellular protein degradation (autophagy) and the proteasome-mediated protein degradation (ubiquitin-proteasome system, UPS). In the latter, ubiquitin-protein ligases covalently attach ubiquitin to proteins, allowing their recognition and further degradation by the proteasome 57 .
As shown in Fig. 1, transcripts encoding for polyubiquitin (ubq-1), had the highest median gene expression across all transcriptomes. However, all ubiquitination-related genes detected in the differential gene expression analysis between the AS and H conditions, were upregulated in AS worms ( Fig. 2 and 3, Data S1). For example, aos-1, encoding for a subunit of the ubiquitin-activating enzyme (E1) 58 , two ubiquitin-protein ligases (E3s without detected cullin domains) 57 , and kelch-like genes (e.g., kel-8-like and kel -20). The former are BTBdomain containing proteins known to interact with E3 enzymes, with kel-8 being involved in the degradation of glutamate neuroreceptors 59 . Additional ubiquitination-related genes upregulated in AS were csn-2, encoding for a component of the COP9 signalosome complex 60 , and proteasome genes (pas-2 and pas-3) 57 .
Among the proteases that were upregulated in AS worms, aspartyl proteases have been involved in neurodegeneration 61 , whereas plasminogen and the zinc matrix metalloproteinase ZMP-2 were both reported to mediate degradation of extracellular matrix (ECM) 62 (Fig. 3). C. elegans ZMP-2 was also shown to prevent the accumulation of oxidized lipoproteins 62 , and therefore, it may contribute to the enhanced antioxidant response observed in this condition. www.nature.com/scientificreports/ Autophagy and amino acid degradation. Besides acting coordinately to withstand stress, autophagy cooperates with apoptotic UPS for the recovery and supply of nutrients when these are scarce (reviewed in 63,64 ). Transcripts of two autophagy-related genes, bec-1 65 and the Ragulator complex protein LAMTOR4 (C7orf59-like) 66 were more abundant in AS nematodes (Fig. 3). While the former positively regulates autophagy 65 , the latter interacts with the mTOR Complex I (mTORC1), and tethers small GTPases (Rags and Rheb) to the lysosomal surface 66 . When amino acid levels are low, mTORC1 is not translocated to the lysosomal surface 66 , thereby favoring catabolic processes such as autophagy 67 . We propose that amino acid scarcity might result from the upregulation of genes involved in the degradation of lysin, glycin, tyrosin, cystein, leucin, isoleucin, valin or tryptophan (Fig. 3, Data S1). This would decrease mTORC1 activity and, in turn, stimulates nutrient recycling via autophagy in AS worms.  www.nature.com/scientificreports/ Conversely, we hypothesize that in H worms, active mTORC1 interacts with the ribosomal protein S6 kinase (S6K), encoded by the rsks-1 gene which is also up in H worms 68 (Fig. 3). This direct interaction, upon a cascade of phosphorylation events, would stimulate translation, and ultimately cell growth and proliferation 68,69 .

AS H A S H Translation Energy generation
All in all, although it is currently unclear whether increased autophagy is beneficial or detrimental, under AS conditions, the upregulation of genes involved in self-digestion might play a protective role and foster recovery from starvation 67 , pathogens 70 or from neuronal and muscular degeneration induced by oxygen deprivation 35 .
Lectins and mucins. Given that symbiont attachment may be mediated by Ca 2+ -dependent lectins [71][72][73] and that, under anoxia, the symbiont appeared to proliferate more 16 , we expected nematode C-type lectins to be upregulated under this condition. Indeed, nine C-type lectin domain (CTLD)-containing proteins were upregulated in AS L. oneistus adults and only two (clec-78 and clec-78-like-2) were upregulated in the presence of oxygen (Fig. 4). In addition to CTLD-containing proteins, mucins, a class of glycoproteins with more than 50% of its mass attributable to O-glycans, were also upregulated in AS nematodes. Considering that mucin glycans are used by vertebrate gut commensals for attachment, as well as being a source of nutrients 74 , it is conceivable that their upregulation in anoxia (Fig. 4), together with that of CTLD-containing proteins, would foster symbiont attachment.
We hypothesize that overexpression of two classes of putative symbiont-binding molecules, lectins and mucins, under conditions favoring symbiont proliferation (i.e., AS conditions) 16 , may mediate bacterial coat reinforcement.
Apoptosis. Mitochondria play an important role in apoptosis induction 54 . Indeed, MPTP opening due to ROS (or the severe ATP decline imposed by the absence of oxygen) may cause cytochrome c release from mitochondria and this, in turn, triggers caspase activation 7 . We observed that transcripts encoding for sco-1, a gene needed for the synthesis and assembly of mitochondrial cytochrome c 50 , were more abundant in AS worms (Fig. 4). Further, we observed upregulation of Caspase-3 (ced-3) which belongs to a family of cysteine proteases involved in apoptosis and which is activated upon mitochondrial cytochrome c release into the cytosol 54,75 . Additional apoptosis-related genes that appeared to be upregulated in AS worms were: bec-1 (Fig. 3), a gene that promotes autophagy and fine-tunes the Ced-3-mediated apoptosis 65 ; ttr-52, which mediates apoptotic cell recognition prior to engulfment 76 ; a BAG family molecular chaperone regulator 1 (BAG1-regulator); a celldeathrelated nuclease crn−2 77 and phagolysosome forming arl-8 78 , and a tyrosine-protein kinase (abl-1) that modulates apoptotic engulfment pathways 79 .
Lipid catabolism. Genes involved in lipid metabolism were similarly expressed between the AS and H conditions (Fig. 2, Data S1). In accordance, lipidomes of nematodes incubated in the presence or absence of oxygen were not significantly different (Fig. S5, Supplemental material). However, in line with the overall upregulation of degradation pathways, we observed upregulation of genes involved in fatty acid beta-oxidation (kat-1) 80 , in lipid digestion (the lipase lipl-6; UniProtKB E2S7J2), and lipid degradation (a peripilin-2-like protein) 81 . Moreover, a gene that might be involved in oxidative stress tolerance (a stearic acid desaturase fat-7 regulating the first step of the fatty acid desaturation pathway 82 was also upregulated in AS worms. Lipid degradation under anoxia might be a strategy to overcome starvation 83 . Notably, we also observed an upregulation of two genes involved in phosphatidylcholine (PC) synthesis (pmt-1, pmt-2) 84 (Fig. 5). Intriguingly, PC was more abundant in the anoxic symbiont 16 , although the latter cannot synthetize it. Thus, their upregulation in AS worms suggests worm-to-symbiont lipid transfer.
GABA-and glutamate-mediated neurotransmission. Upregulated genes related to GABA synthesis were unc-25, unc-104 and pdxk-1 (pyridoxal phosphate hexokinase) [85][86][87] (Fig. 5, Data S1). Consistent with an expected increase in glutamate requirement as a direct GABA precursor 88 , we observed downregulation of two glutamine synthetases and a delta-1-pyrroline-5-carboxylate synthase (gln-3 and alh-13 respectively; Fig. 6) 89,90 , known to convert glutamate to glutamine or to proline, respectively. Furthermore, an mgl-2-like gene encoding for a glutamate receptor, which is activated in the presence of glutamate 91 , was up in AS worms. Note that, when oxygen is limited, glutamate may act as a neurotoxic amino acid 92 . Therefore, increased GABA biosynthesis might, beneficially, prevent its accumulation 93 .
GABA-mediated neurotransmission has been documented for facultative anaerobic animals thriving in anoxic conditions 92,94 . Due to its inhibitory nature, it contributes to avoid membrane depolymerization 94 . Moreover, given that it relaxes muscles, the increment of GABA may impact the movement of the animal 85 . Therefore, upregulation of GABA-mediated neuronal activity might explain why L. oneistus did not form tight clusters in anoxia (Supplemental Movie 3).
Dopamine-mediated neurotransmission. A gene encoding for the tyrosine hydroxylase Cat-2 (cat-2), which is needed for dopamine biosynthesis 95 and two putative dopamine receptors (protein-D2-like and a G_PROTEIN_ RECEP_F1_2 domain-containing protein (dop-5)) 96 were upregulated in AS worms. Moreover, a dat-1-like gene mediating dopamine reuptake into the presynaptic terminals was downregulated 97 in AS worms (Fig. 5).
Calcium-binding and -sensing proteins. Finally, in AS worms several calcium-binding or -sensing proteins (e.g., ncs-2, cex-2, and a calbindin-like (CALB1 homolog); Fig. 5) 98,99 , as well as calcium transporters (cca-1; Transport category, Fig. 6) 100 were upregulated. On the one hand, we hypothesize their involvement in the inhibitory neural signaling described above (for example, Ncs-2 mediates the cholinergic and GABAergic expression of www.nature.com/scientificreports/ C. elegans 101 . On the other, they may protect cells against the stress inflicted by anoxia, which involves calcium overload and consequent cellular acidification 7 .

Genes upregulated in hypoxic (H) nematodes. Innate immune pathways and effectors. Animals rec-
ognize and respond to microbes by means of immunoreceptors including Toll-like receptors, conserved from sponges to humans 102 . We identified almost all genes belonging to this pathway, including the one encoding for the NF-kB transcription factor. This came as a surprise given that, up to now, the NF-kB has not been identified in any other nematode 103 . Equally surprising was the fact that not only two Toll-like receptors (tol-1 and tol-1like), but also genes encoding for antimicrobial proteins such as a peroxisome assembly factor involved in defense against Gram-(prx-11-like) 104 , a putatively antifungal endochitinase 105 and bactericidal/permeability-increasing proteins (BPIs) were also more abundant in H worms. BPIs may bind LPS and perforate Gram-membranes and have been shown to play a symbiostatic role in other invertebrates 106,107 . However, it is unclear whether activation of the L. oneistus Toll pathway leads to the nuclear NF-kB switching on the expression of antimicrobial genes or whether, as shown in C. elegans, the Toll pathway mediates behavioral avoidance of pathogens 108 . Overall, the apparent oxygen stimulation of a central innate immunity pathway and, directly or indirectly, of broad range antimicrobials could be adaptations to the fact that when crawling in superficial oxygenated sand, L. oneistus is not only exposed to predation by bigger animals, but also to pathogenic members of the bacterioplankton. Overexpression of broad-range antimicrobials in response to oxygen might therefore help L. oneistus to avoid colonization by potentially deleterious, fouling bacteria (e.g., vbrios, roseobacters and pseudoaltermonas/ alteromonadales) when crawling close to the water column 109  Development. Although development-related genes were some of the most expressed under all conditions (Fig. 1), many were upregulated in H nematodes (Figs. 2 and 5). Among the development-related genes upregu-  Genes are ordered by function in their respective metabolic pathways. For each process, the minority of genes that were upregulated in AS worms is shown in Data S1. Red denotes upregulation, and blue downregulation. FA, fatty acids; PC phosphatidylcholine; PL, phospholipids; metab: metabolism; synth: synthesis; assim: assimilation, oxid: oxidation; Transp: transporters. www.nature.com/scientificreports/ lated in H nematodes were those related to molting (e.g., nas-36, nas-38, chs-2, ptr-5, ptr-18, apl-1, myrf-1) [110][111][112][113][114] , germ line establishment (e.g., ccm-3, rsks-1) 115,116 , oogenesis/spermatogenesis (crt-1) 117 , embryonic development and yolk production (e.g., plt-1, mlc-5, crt-1, arrd-17, cpna-1, vit-6) [118][119][120][121][122][123] , and/or larval development (nmy-1, ifb-1) 124,125 , as well as male tail tip (cdt-1, plx-1, ver-3) [126][127][128] , and vulva morphogenesis (hda-1, unc-62) 129,130 (Fig. 5). Morever, transcripts encoding for a number of proteases shown to be involved in C. elegans molting (e.g., nas-36, nas-38) 110 , development (e.g., teneurin-a-like) 131 , neuronal regrowth or locomotion (tep-1) 132 and pharyngeal pumping (e.g., neprilysin nep-1) 133 were also more abundant in H worms. Remarkably, vav-1, which besides being involved in male tail tip and vulva morphogenesis 126 may also regulate the concentration of intracellular calcium 134 , was one of the few development-related genes to be downregulated in H nematodes (see previous section on Ca 2+ -binding proteins).
To sum up, and as expected, the host appears to exploit oxygen availability to undertake energetically costly processes, such as development and molting 135 .
Although glycolysis seems to generate ATP in both AS and H worms, it is not clear why the latter would prefer to respire cellulose or trehalose instead of starch. Given its role as a membrane stabilizer, we speculate that AS worms might prioritize the storage of trehalose over its degradation to preserve membrane integrity ( Fig. 6) 4,142 . Of note, based on its genome draft, the symbiont may synthetize and transport trehalose, but it may not use it 16 . Therefore, we hypothesize symbiont-to-host transfer of trehalose under hypoxia. Consistently, the symbiont's trehalose synthesis-related gene (otsB) 16 , and the host trehalase (tre-1; Fig. 6) were both upregulated under hypoxia and metabolomics could detect trehalose in both partners (Table S1). Metabolomics also detected sucrose in both the holobiont and the symbiont fraction (Table S1). Given that, based on transcriptomics and proteomics, the nematode can utilize sucrose but cannot synthesize it (Data S1), whereas the symbiont can 16 , we also hypothesize symbiont-to-host sucrose transfer.
Amino acid biosynthesis. Transcripts of genes involved in the synthesis of glutamine and proline (gln-3 and alh-13, respectively), aspartate (L-asparaginases) 157 and S-adenosyl-L-methionine (SAM) (sams-4) 158 were all upregulated in H worms (Fig. 6), as well as one encoding for the ornithine decarboxylase odc-1 which is involved in biosynthesis of the polyamine putrescin, and is essential for cell proliferation and tissue growth 159 . Moreover, polyamines, with their high charge-to-mass ratio may protect against superoxide radicals, which, as mentioned, harm cell membranes and organelles, oxidize proteins, and damage DNA 160 .
On the other hand, ceramides, which have anti-proliferative properties and which may mediate resistance to severe oxygen deprivation 165 , appeared to be mainly synthesized in AS worms, as indicated by the upregulation of genes involved in ceramide biosynthesis (asm-3, ttm-5, Fig. 6) 166 .
Sulfur metabolism. The mpst-7 gene, involved in organismal response to selenium and switched on in hypoxic C. elegans 175 , was upregulated in H nematodes (Fig. 6). Given that MPST-7 is thought to catalyze the conversion of sulfite and glutathione persulfide (GSSH) to thiosulfate and glutathione (GSH) 176 , hypoxia-experiencing L. oneistus might express this enzyme to recharge the cells with GSH and hence, help to cope with oxidative stress 177 . Also more abundant in H worms were transcripts encoding for the sulfatases 2 (sul-2) 178 and a PAPSproducing pps-1 (3′-phospho-adenosine-5′-phosphosulfate (PAPS), considered the universal sulfur donor) 141 , as well as for the chaperones pdi-6 and protein-disulfide-isomerase-A5-like which require oxygen to mediate correct disulfide bond formation in protein folding 6,179 (Fig. 6). Conversely, a putative sulfide-producing enzyme (mpst-1) which protects C. elegans from mitochondrial damage 180 was upregulated in AS nematodes, albeit together with two genes encoding for enzymes involved in its detoxification, a persulfide dioxygenase (ethe-1) and a cysteine dioxygenase (cdo-1) 179 . The first detoxifies sulfide by producing glutathione, while the latter catalyzes taurine synthesis via cysteine degradation. Of note, sulfide detoxification via taurine accumulation is a common strategy in chemosynthetic animals (reviewed in 181 ). All in all, L. oneistus appeared to limit excess accumulation of free sulfide in anoxia and to free sulfate when oxygen was available.

Conclusions
Overall, and irrespective of the conditions it was subjected to, L. oneistus mostly expressed genes involved in degradation processes, energy generation, stress response and immune defense. Astonishingly, L. oneistus did not enter suspended animation when subjected to anoxic sulfidic conditions for days. We hypothesize that in the absence of oxygen, ATP production is supported by symbiont-derived trehalose and cellulose catabolism, and by rewiring the ETC in such way as to use RQ as electron carrier, and fumarate as electron acceptor. Moreover, the nematode activates several degradation pathways (e.g., UPS, autophagy, and apoptosis) to gain nutrients from anoxia-damaged proteins and mitochondria. Further, AS worms also upregulated genes encoding for ribosomal proteins and putative symbiont-binding proteins (lectins). Finally, as proposed for other anoxia-tolerant animals, the worm seems to upregulate its antioxidant capacity in anticipation of reoxygenation. When in hypoxic conditions (Fig. 7, left), instead, we speculate that the worm uses starch for energy generation to engage in costly developmental processes such as molting, feeding, and mating, likely relying on excitatory neurotransmitters (e.g., acetylcholine), and it upregulates the Toll immune pathway and, directly or indirectly, the synthesis of broad range antimicrobials (e.g., fungicides, BPIs).
When looking at the Laxus-Thiosymbion symbiosis in light of what was recently published 16 , we could identify two signs of inter-partner metabolic dependence: under anoxia, worms might transfer lipids to their symbionts, and under hypoxia, the symbionts might transfer trehalose and sucrose to their hosts.
Furthermore, we may conclude that, wherever in the sand the consortium is, one of the two partners is bound to be stressed: in anoxia, the symbiont appears to proliferate more, while its animal host engages in degradation of damaged proteins and mitochondria and in detoxification. In the presence of oxygen, the situation is inverted: the symbiont seems massively stressed, while the host can afford energy-costly biosynthetic processes to develop and reproduce (Fig. 7). It is therefore fascinating that, in spite of the dramatically different needs a bacterium and animal must have, the Laxus-Thiosymbion symbiosis evolved.

Materials and methods
Sample collection. Laxus oneistus individuals were collected on multiple field trips (2016-2019) with cores of 60 cm length and 60 mm diameter (UWITEC, Mondsee, Austria) down to a depth of approximately 1 m depth in the sand bars off the Smithsonian Field Station, Carrie Bow Cay in Belize (16°48′11.01″N, 88°4′54.42″W). The collection of the nematodes, the incubations set up for RNA sequencing, lipidomics, proteomics and metabolomics, as well as the RNA extraction, and library preparation are in detailed described in 16 . Importantly, the nematodes had a bright white appearance and replicate incubations were started simultaneously. Note that the Supplemental material describes the metabolomics and sequencing data of Olavius algarvensis, as well as changes from 16 in the lipidomics and proteomics pipelines.
Host transcriptome de novo assembly. In preparation for the assembly, reads from each sample were first mapped to the symbiont as described before 16 , and remaining rRNA reads from all domains of life were removed from unmapped reads using sortmerna v2.1 in combination with the SSURef_NR99_119_ SILVA_14_07_14 and LSURef_119_SILVA_15_07_14 databases. Further, exact duplicate reads were removed using PRINSEQ lite's derep option. Read files free of symbiont reads, rRNA reads and exact duplicates were used as input for transcriptome sub-assemblies via Trinity v2.6.6 with the strand-specific option (--SS_lib_type F) 182 . Two sub-assemblies differing in the number and type of input read files were performed: (1) 9 input read files including biological triplicates from 3 incubation conditions (O, H, A) and (2) 4 input read files including a single replicate from 4 incubation conditions (O, H, A and hyper-O). Hyper-O refers to an incubation in which air was pumped directly into the exetainers for the entire incubation period to supersaturate the seawater (300% O 2 ). However, as this incubation condition yielded an incongruous transcriptional response by the symbiont (data not shown), these read data were only used to extend the host transcriptome's coding repertoire. The qualities of both sub-assemblies were assessed as described below.
We then performed an intra-assembly clustering step as described in 183 , during which identical transcripts were removed from the sub-assemblies using CD-HIT-EST 184  www.nature.com/scientificreports/ longest isoform for each 'gene' identified by Trinity was kept using Trinity's get_longest_isoform_seq_per_trin-ity_gene.pl utility. The remaining transcripts of each sub-assembly were then concatenated to produce a merged transcriptome assembly. The final assembly was created by applying another sequence clustering using CD-HIT-EST to avoid inter-assembly redundancy. Here, the identity parameter of 80% (-c 0.8) combined with a minimal coverage ratio of the shorter sequence of 80% (-aS 0.8) and minimal coverage ratio of the longest sequence of 0.005% (-aL 0.005) yielded the best-performing assembly in terms of number of transcripts (162,455) and contiguity (N50 value of 770) (data not shown). Assembly completeness was assessed by estimating completeness via BUSCO nematode single-copy orthologs (Simão et al. 2015). Importantly, the merged assembly yielded a higher BUSCO-based completeness compared with the two sub-assemblies; 79.2% of the BUSCO nematode single-copy orthologs were found to be present and complete in the final assembly (636 single-copy; 142 duplicated), whereas assembly (1) scored 77.8% (233 single-copy; 531 duplicated) and assembly (2) was 76.2% complete (314 single-copy; 434 duplicated). Further, assembled transcripts were filtered based on taxonomic classification. Transcripts were matched against the RefSeq protein database using blastx (E-value 1E−3), and the output was then used as input for taxonomic assignment via MEGAN v5 185 . Only transcripts classified as belonging to 'Eukarya' were kept (MEGAN parameters: Min Score: 50, Max Expected: 1E-2, Top Percent: 2), which reduced the number of putative L. oneistus transcripts to 30,562. Assembled transcripts were also functionally annotated using Trinotate 186 . Briefly, predicted protein coding regions were extracted using TransDecoder (https:// github. com/ Trans Decod er), both transcripts and predicted protein sequences were searched for protein homology via blastx and blastp, respectively, and predicted protein sequences were annotated for protein domains (hmmscan), signal peptides (signalP) and transmembrane domains (TMHMM). 85,859 transcripts exhibited at least one functional annotation. Finally, only taxonomyfiltered transcripts with at least one functional annotation were kept, thereby further reducing the number of putative host transcripts to 27,984, with 22,072 thereof predicted to contain protein coding regions. BUSCObased completeness for this filtered host transcriptome assembly was 78.8% (635 single-copy; 139 duplicated).

Gene expression analysis.
Raw sequencing reads quality assessment and preprocessing of data was followed as described in 16 . Trimmed reads were mapped to the de novo transcriptome assembly and transcript abundance was estimated using RSEM v1.3.1 187 Figure 7. Schematic representation of Laxus oneistus physiology in anoxic or hypoxic conditions. In anoxic sulfidic conditions (left), L. oneistus does not enter suspended animation. Instead, it upregulates the expression of genes mediating inhibitory neurotransmission, involved in symbiosis establishment (e.g., lectins, mucins) and in ribosome biogenesis. Metabolism may be supported by the degradation of starch and by rewiring the electron transfer chain: rhodoquinone (RQ) is used as electron carrier and fumarate as electron acceptor. Moreover, the worm activates degradation pathways (e.g., ubiquitin-proteasome system (UPS), autophagy, and apoptosis) and may anticipate reoxygenation by upregulating superoxide dismutase (SOD) and glutathione peroxidase (GP). In hypoxic conditions (right), instead, L. oneistus appears to use trehalose and cellulose for energy generation, while engaging in costly processes such as development, molting, feeding, and mating. Genes involved in excitatory neurotransmission are also upregulated, together with Toll-like receptors and immune effectors (e.g., fungicides, BPIs). www.nature.com/scientificreports/ the application of strandedness (--strandedness forward). Read counts per transcript were used for differential expression analysis, and TPM (transcripts per kilobase million) values were transformed to log 2 TPMs as described in 16 . Gene and differential expression analyses were conducted using R and the Bioconductor package edgeR v3.28.1 188,189 , and as shown in 16 . Here, we only describe the modifications that were made to the pipeline. Genes were considered expressed if at least ten reads in at least three replicates of one of the four conditions could be assigned. Excluding the replicates of the oxic condition, we found 74.9% of all predicted nematode proteinencoding genes to be expressed (16,526 genes out of 22,072). Log 2 TPM were used to assess sample similarities via multidimensional scaling based on Euclidean distances (R Stats package) 189 (Fig. S1B), and the average of replicate log 2 TPM values per expressed gene and condition was used to estimate expression strength. Median gene expression of entire metabolic processes and pathways per condition was determined from average log 2 TPM values.
Expression of genes was considered significantly different if their expression changed 1.5-fold between two treatments with a false-discovery rate (FDR) ≤ 0.05 190 . Throughout the paper, all genes meeting these thresholds are either termed differentially expressed or up-or downregulated. For the differential expression analyses between the AS, H and A conditions see Data S1. Heatmaps show mean-centered log 2 TPM expression values to highlight gene expression change.
All predicted L. oneistus proteins were automatically annotated using eggNOG-mapper v2 191 against eggNOG 5.0 192 using diamond v2.0.4 193 . All genes that are shown and involved in a particular process were manually curated by blasting (blastp) them against both the nr database 194 and the WormBase (https:// wormb ase. org/ tools/ blast_ blat) 195 . The Spearman's rank correlation was carried out with per-gene averaged log 2 TPM and orgNSAF% values over the oxic and anoxic incubations for the transcriptome and proteome, respectively, in R 189 .

Data availability
This Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession GJNO00000000. The version described in this paper is the first version, GJNO01000000. RNA-Seq data are available at the Gene Expression Omnibus (GEO) database and are accessible through accession number GSE188619. www.nature.com/scientificreports/