Brian W Strehlow

Astrid Schuster

Warren R Francis

Lisa Eckford-Soper

Beate Kraft

Rob Mcallen

Ronni Nielsen

Susanne Mandrup

Donald E Canfield

Transcriptomic responses of sponge holobionts to in situ, seasonal anoxia and hypoxia

Deoxygenation can be fatal for many marine animals; however, some sponge species are tolerant of hypoxia and anoxia. Indeed, two sponge species, Eurypon sp. 2 and Hymeraphia stellifera, survive seasonal anoxia for months at a time. To understand their tolerance mechanisms, we performed differential gene expression analyses on the sponges, their mitochondria and their microbial symbionts under in situ conditions of normoxia, hypoxia and anoxia. Each species possessed a unique microbiome, but the microbiomes of each species were dominated by a species-specific Thaumarchaeon and a Gammaproteobacterium. Holobiont gene expression was species-and oxygen-level dependent, though there were some shared interspecific responses to deoxygenation. In general, few changes occurred in the expression of sponge metabolic genes as a function of oxygenation level, indicating that they may remain metabolically active under anoxia. However, ATP synthesis genes were significantly upregulated under hypoxia when compared to normoxia, and genes for DNA replication were downregulated. Mitochondrial gene expression was effectively unchanged under both hypoxia and anoxia. Nevertheless, both anoxia and hypoxia caused upregulation of heat shock proteins (HSPs), indicating cellular level adaptations to deoxygenation stress. A meta-analysis demonstrated that sponge transcriptional responses to anoxia were distinct from those displayed by other invertebrates while dormant, and the hypothesis of sponge dormancy under anoxia was not supported. Thaumarchaeota symbionts also upregulated stress response genes in hypoxia, while maintaining expression of oxygen-dependent metabolic pathways under hypoxia and anoxia. Gammaproteobacteria symbionts showed relatively few noteworthy changes in gene expression in response to anoxia but decreased metabolic gene expression in hypoxia. There was no clear evidence of upregulated anaerobic respiration in the transcriptomes of the sponge holobionts under anoxia or hypoxia. The tolerance of some sponges to prolonged anoxia warrants further investigation and could give them an advantage in future oceans following climate change as well as in ancient oceans when oxygen concentrations were lower than at present.

Introduction

Ocean oxygen levels are predicted to decline by 1-7% over the next century [START_REF] Keeling | Ocean Deoxygenation in a Warming World[END_REF]. Climate change is expanding marine oxygen-minimum zones, while enhanced nutrient runoff is causing coastal eutrophication and subsequent deoxygenation [START_REF] Kirchman | Dead Zones: The Loss of Oxygen from Rivers, Lakes, Seas, and the Ocean[END_REF]. Deoxygenation is expected to negatively affect the physiology, metabolism and survival of many marine organisms (Diaz and Rosenberg 2008). As a result of deoxygenation and other anthropogenic stressors, we are facing large scale habitat and biodiversity losses, but the full effects of these stressors are still unknown. Sessile organisms that cannot move or migrate in response to deoxygenation are particularly vulnerable. Sponges (Phylum: Porifera Grant, 1836), sessile, ecologically important filterfeeding animals, are hypothesized to suffer from increasing hypoxic/anoxic areas worldwide [START_REF] Breitburg | Declining Oxygen in the Global Ocean and Coastal Waters[END_REF][START_REF] Levin | Manifestation, Drivers, and Emergence of Open Ocean Deoxygenation[END_REF][START_REF] Stramma | Expanding Oxygen-Minimum Zones in the Tropical Oceans[END_REF], and some in situ mortality of sponges (along with other taxa) following deoxygenation has been reported [START_REF] Altieri | Tropical Dead Zones and Mass Mortalities on Coral Reefs[END_REF][START_REF] Stachowitsch | Mass Mortality in the Gulf of Trieste: The Course of Community Destruction[END_REF]). However, species-specific responses to deoxygenation in Porifera are variable, and some species are tolerant of hypoxia. Lab experiments have shown that the demosponge Halichondria (Halichondria) panicea Pallas, 1766 can survive low oxygen levels down to 0.5-4% of air saturation (a.s.) [START_REF] Mills | Oxygen Requirements of the Earliest Animals[END_REF], and Tethya wilhelma can survive down to 0.25-2% a.s. for 4 days and <0.05% a.s. for at least one hour [START_REF] Mills | The Last Common Ancestor of Animals Lacked the HIF Pathway and Respired in Low-Oxygen Environments[END_REF]. Additional demosponge species that tolerate similar levels of hypoxia include Geodia barretti Bowerbank, 1858 with 7% a.s. (20 µM O2, Leys and[START_REF] Leys | Oxygen and the Energetic Requirements of the First Multicellular Animals[END_REF] and Haliclona pigmentifera Dendy, 1905 (0.50-0.68% a.s., Gunda and[START_REF] Gunda | Effects of Dissolved Oxygen Levels on Survival and Growth in Vitro of Haliclona pigmentifera (Demospongiae)[END_REF]. In a recent experimental study, four taxonomically distinct sponge species (Cliona celata Grant, 1826, Suberites carnosus Johnston, 1842, Polymastia crocea Kelly-Borges and [START_REF] Kelly-Borges | Revision of Southwest Pacific Polymastiidae (Porifera: Demospongiae: Hadromerida) with descriptions of new species of Polymastia Bowerbank, Tylexocladus Topsent, and Acanthopolymastia gen. nov. from New Zealand and the Norfolk Ridge, New Caledonia[END_REF] and Suberites australiensis [START_REF] Bergquist | The Marine Fauna of New Zealand: Porifera, Demospongiae, Part 1 (Tetractinomorpha and Lithistida)[END_REF] were exposed to various levels of hypoxia (1.5-40% a.s., 0.13-3.3 mg O2 L -1 ) in aquaria. These sponges exhibited species-specific physiological, morphological and behavioral responses to deoxygenation. Sub-lethal oxygen thresholds for these sponges were in the range of 6-20% a.s. (0.48-1.56 mg O2 L -1 ) while mortality was only observed in one of the two species (P. crocea) exposed to severe hypoxia lower than 5% a.s. (0.13 mg O2 L -1 ) (Micaroni et al. 2021a). The sponge species most sensitive to deoxygenation in controlled experiments, P. crocea, had a median lethal time (time until death, LT50) of 286 h at ~1.5% a.s. (0.13 mg O2 L -1 , (Micaroni et al. 2021a), which was still longer than most marine organisms under comparable conditions (reviewed in Micaroni et al. 2021a and[START_REF] Hughes | Coral Reef Survival under Accelerating Ocean Deoxygenation[END_REF]. Few studies have investigated the deoxygenation tolerance of sponges in situ. [START_REF] Bayer | Microbial Strategies for Survival in the Glass Sponge Vazella Pourtalesii[END_REF] found that the glass sponge Vazella pourtalesii Schmidt, 1870 lives in reduced oxygen conditions (< 4 mL/L) on the Scotian Shelf, Canada. [START_REF] Schuster | Effects of Seasonal Anoxia on the Microbial Community Structure in Demosponges in a Marine Lake in Lough Hyne, Ireland[END_REF] reported at least nine encrusting demosponge species of the orders Axinellida [START_REF] Lévi | Sur une nouvelle classification des Démosponges[END_REF]Poecilosclerida Topsent, 1928 that tolerate seasonal anoxia in a Marine Lake in Lough Hyne, Ireland for ~three months at a time in the summer, and there are more sponge species at this site [START_REF] Bell | The Distribution and Prevalence of Sponges in Relation to Environmental Gradients within a Temperate Sea Lough: Vertical Cliff Surfaces[END_REF]. The anoxia-tolerant sponges of Lough Hyne therefore have LT50s that exceed their ca. three month exposure to anoxia, and, compared to other animals, their tolerance is rivaled only by the anemone Bunodosoma cavernatum Bosc, 1802 (LT50 in anoxia >1.5 months, [START_REF] Ellington | Metabolic Responses of the Sea Anemone Bunodosoma Cavernata (Bosc) to Declining Oxygen Tensions and Anoxia[END_REF]. Moreover, the abundance (percent cover) of encrusting sponges in the seasonally anoxic site of Lough Hyne only slightly decreased over 30 years despite widespread sponge losses in the well oxygenated parts of the lough (Micaroni et al. 2021b), indicating that these sponges are uniquely suited to their seasonally anoxic niche. Many invertebrates exposed to harsh conditions like anoxia or hypoxia have evolved molecular response pathways to survive in these environments. For instance, under anoxia, marine mollusks shift from lipid, carbohydrate and amino acid metabolism to carbohydrate catabolism as the primary fuel for ATP generation (Larade and Storey 2009). Some mollusks also maximize anoxic survival time (mollusc LT50s at 0 mg L -1 : 360-672 h, see Micaroni, Strano, et al. 2021b) by building up a fixed reserve of internal fuels (primarily glycogen) to increase ATP output and to reduce acidity. A small suite of upregulated genes also plays a critical role in cell preservation, facilitating long-term viability under hypoxia (Larade and Storey 2009). Some anemones may even have elements of a facultative anaerobic metabolism [START_REF] Sassaman | Adaptations To Environmental Oxygen Levels In Infaunal And Epifaunal Sea Anemones[END_REF].

Only one study has investigated the transcriptomic response of a sponge to hypoxia. After four days of exposure to hypoxia (0.25-2% a.s.), 128 genes within the sponge T. wilhelma were differentially expressed compared to controls, but none of these genes were clearly related to metabolic or stress responses [START_REF] Mills | The Last Common Ancestor of Animals Lacked the HIF Pathway and Respired in Low-Oxygen Environments[END_REF]). However, one hour exposure to near anoxia (<0.05% a.s., 105 nM O2) elicited stress responses in T. wilhelma's transcriptome [START_REF] Mills | The Last Common Ancestor of Animals Lacked the HIF Pathway and Respired in Low-Oxygen Environments[END_REF], demonstrating that transcriptomic responses were mostly observed under very low oxygen levels. In most metazoans, the hypoxia-inducible factor (HIF) pathway is essential for maintaining oxygen homeostasis in response to hypoxia [START_REF] Samenza | Life with Oxygen[END_REF]. However, in sponges and ctenophores, the likely sister-lineages to other metazoans [START_REF] Feuda | Improved Modeling of Compositional Heterogeneity Supports Sponges as Sister to All Other Animals[END_REF], key components of the HIF pathway, including orthologs for HIFa and HIFb, are absent [START_REF] Mills | The Last Common Ancestor of Animals Lacked the HIF Pathway and Respired in Low-Oxygen Environments[END_REF]. HIFa is a transcription factor responsible for triggering metabolic adaptation to deoxygenation in most metazoans [START_REF] Samenza | Life with Oxygen[END_REF]. Due to the absence of HIFa in sponges, the mechanisms of tolerance to hypoxia and anoxia are largely unknown [START_REF] Mills | The Last Common Ancestor of Animals Lacked the HIF Pathway and Respired in Low-Oxygen Environments[END_REF]). The sponge microbiome may play a role in sponge deoxygenation tolerance. Sponges can host diverse microbiomes that perform vital functions and can constitute large portions of their biomass [START_REF] Wilkinson | Microbial Associations in Sponges. II. Numerical Analysis of Sponge and Water Bacterial Populations[END_REF]. The host and microbiome together are referred to as a 'holobiont' [START_REF] Webster | Marine Sponges and Their Microbial Symbionts: Love and Other Relationships[END_REF]. Some sponge symbionts have anaerobic capabilities, including sulfate reduction [START_REF] Hoffmann | An Anaerobic World in Sponges[END_REF], denitrification, and anaerobic ammonium oxidation (anammox) [START_REF] Hoffmann | Complex Nitrogen Cycling in the Sponge Geodia barretti[END_REF]). In V. pourtalesii, the microbial community is thought to be specialized for the low-oxygen conditions in which it lives [START_REF] Bayer | Microbial Strategies for Survival in the Glass Sponge Vazella Pourtalesii[END_REF]. A previous study has shown that the microbiomes of anoxia-tolerant sponges from Lough Hyne were unchanged under anoxia, hypoxia and normoxia despite drastic changes in the microbial communities of the water and sediment, indicating that the holobiont as a whole is tolerant to deoxygenation [START_REF] Schuster | Effects of Seasonal Anoxia on the Microbial Community Structure in Demosponges in a Marine Lake in Lough Hyne, Ireland[END_REF]. Specifically, the microbiomes of the anoxia-tolerant sponge species Eurypon sp. 2 and Hymeraphia stellifera Bowerbank, 1864 were both dominated by two microbesone Thaumarchaeon and one Gammaproteobacterium, though each sponge hosted a different species from each clade. Eurypon sp. 2 harbors the Thaumarchean OTU1 (closest related isolate: Nitrosopumilus adriaticus strain NF5), and H. stellifera harbors OTU3 (closest related isolate: Nitrosopumilus oxyclinae strain HCE1). For the Gammaproteobacteria, Eurypon sp. 2 contains OTU2 (closest related isolate: Arenicella chitinivorans strain ALEN 1); whereas H. stellifera contains OTU7 (closest related isolate: A. chitinivorans strain ALEN 1, [START_REF] Schuster | Effects of Seasonal Anoxia on the Microbial Community Structure in Demosponges in a Marine Lake in Lough Hyne, Ireland[END_REF]. Although their location within Lough Hyne sponges has not yet been determined, Gammaproteobacteria and Thaumarchaeota symbionts are abundant and widely distributed throughout the sponge tissue in other species [START_REF] Moeller | Characterization of a Thaumarchaeal Symbiont That Drives Incomplete Nitrification in the Tropical Sponge Ianthella basta[END_REF]. These symbionts also participate in integrated metabolic processes with their host species [START_REF] Moitinho-Silva | Integrated metabolism in sponge-microbe symbiosis revealed by genome-centered metatranscriptomics[END_REF]. [START_REF] Schuster | Effects of Seasonal Anoxia on the Microbial Community Structure in Demosponges in a Marine Lake in Lough Hyne, Ireland[END_REF] hypothesized potential metabolic pathways within these two symbiont taxa that might function under deoxygenation. Based on the genomes of their close relatives, the Thaumarchaeota in these sponges were hypothesized to have components of a facultatively anaerobic metabolism, including fumarate, nitrite, and sulfite respiration and the ability to ferment [START_REF] Bayer | Microbial Strategies for Survival in the Glass Sponge Vazella Pourtalesii[END_REF], which could aid in deoxygenation tolerance. A recent study demonstrated that another close Thaumarchaeota relative, Nitrosopumilus maritimus, can produce trace levels of oxygen (100-200 nM) under environmental anoxia [START_REF] Kraft | Oxygen and Nitrogen Production by an Ammonia-Oxidizing Archaeon[END_REF]. The pathway of oxygen production and the genes involved are not currently known. However, comparative proteomics suggested that proteins upregulated during oxygen production included some AsnC family proteins and nitroreductases [START_REF] Kraft | Oxygen Production by an Ammonia-Oxidizing Archaeon[END_REF]. We therefore hypothesized that nitroreductases and AsnC family proteins would be upregulated under hypoxia and anoxia in sponge Thaumarchaeota. Gammaproteobacteria symbionts were hypothesized to use sulfide or thiosulfate as an electron donor anaerobically coupled to nitrate reduction, removing exogenous, toxic sulfide from the holobiont under anoxia [START_REF] Schuster | Effects of Seasonal Anoxia on the Microbial Community Structure in Demosponges in a Marine Lake in Lough Hyne, Ireland[END_REF]. To understand the mechanisms of sponge holobiont tolerance and survival under extended anoxia, we performed differential gene expression analysis of the sponges Eurypon sp. 2 and H. stellifera as well as their respective mitochondria and microbial symbionts including Thaumarchaeota and Gammaproteobacteria. Sponges were collected from Lough Hyne under in situ normoxia, hypoxia and anoxia. Since neither of these holobionts are established model species, we generated reference transcriptomes for the sponges and assembled the genomes of their symbionts and mitochondria from metagenomic sequences. Additionally, we performed a meta-analysis of the sponge differential expression data to determine if the responses of Lough Hyne sponges to deoxygenation are similar to those of T. wilhelma under deoxygenation, other invertebrates under heat stress [START_REF] Dixon | Coral Reefs. Genomic Determinants of Coral Heat Tolerance across Latitudes[END_REF], in a dauer state [START_REF] Sinha | Divergent Gene Expression in the Conserved Dauer Stage of the Nematodes Pristionchus Pacificus and Caenorhabditis Elegans[END_REF], and during diapause [START_REF] Gong | Global Transcriptome Analysis of Orange Wheat Blossom Midge, Sitodiplosis Mosellana (Gehin) (Diptera: Cecidomyiidae) to Identify Candidate Transcripts Regulating Diapause[END_REF]. Finally, we considered the implications of sponge deoxygenation tolerance for future oceans as well as early animal evolution under low atmospheric oxygen concentrations.

Materials and Methods

Sampling

Sponges were sampled between July 2018 and August 2019 at Labhra Cliff (Lough Hyne Nature Reserve, Ireland; N51°30.05309 W9°18.17679) as described in [START_REF] Schuster | Effects of Seasonal Anoxia on the Microbial Community Structure in Demosponges in a Marine Lake in Lough Hyne, Ireland[END_REF] under permit no. R23-27/2018 issued by the Irish Department of Environment, Heritage, and Local Governments. Samples were taken under various in situ oxygen conditions that were categorized as follows: normoxic (5.3-12 mg L -1 , 49.3-111% a.s.), hypoxic (1.30-3.56 mg L -1 , 12.1-33.1% a.s.), and anoxic (0.00-0.01 mg L -1 , 0.00-0.93% a.s., i.e. instrument detection limit). All samples were frozen in liquid nitrogen within ~15 minutes of collection and kept at liquid nitrogen temperatures in a dry shipper until returned to the laboratory where samples were stored in a -80°C freezer until extraction. The sponges were identified using molecular barcodes (see [START_REF] Schuster | Effects of Seasonal Anoxia on the Microbial Community Structure in Demosponges in a Marine Lake in Lough Hyne, Ireland[END_REF] as Eurypon sp. 2 and H. stellifera. DNA/RNA extraction Total RNA was isolated using TRIzol TM Reagent from Invitrogen (ThermoFisher Scientific, Denmark). For RNA sequencing (RNAseq), eight samples of H. stellifera (anoxic = 2, hypoxic = 5, normoxic = 1) were taken, and sixteen samples of Eurypon sp. 2 (anoxic = 2, hypoxic = 7, normoxic = 7) were collected (n = 24 in total). Metadata for individual samples, including oxygen concentration, season, depth and collection date are included in Supplemental Table 1. Since these two species are indistinguishable in situ, the samples taken for each species were unequal in number. Remaining DNA in the RNA extracts were removed using the DNase MAX kit (Qiagen). Additionally, high molecular weight (HMW) DNA was extracted from one individual of each species using the protocol from [START_REF] Vargas | Slime Away: A Simple CTAB-Based High Molecular Weight DNA and RNA Extraction Protocol for 'Difficult' Invertebrates[END_REF] for metagenomic sequencing. The Eurypon sp. 2 sample used for metagenomic sequencing was taken under anoxia (DC24), and RNA was also sequenced from this sample. Due to limited availability of samples and DNA, a sample of H. stellifera that was not included in the RNAseq (DC297, sampled under normoxia) was used for metagenome sequencing. DNA and RNA yields were determined using a Qubit TM (Thermo Fisher Scientific) and a NanoDrop TM spectrophotometer (Thermo Fisher Scientific). RNA quality was evaluated on an Agilent TapeStation using an RNA ScreenTape (Agilent Technologies).

Reference transcriptomes and metatranscriptomes

All RNA samples were sequenced twice with two different library preps -one for polyadenylated (poly A) messenger RNA (mRNA, from eukaryotes) and one for total RNA (to also capture prokaryotic RNA). Library preparation for the eukaryotic (sponge) RNA was done using NEBNext® Ultra™ II Directional RNA Library Prep Kit for Illumina (New England Biolabs, USA). For the prokaryotic RNAseq, the supernatant from the poly A enrichment was treated with NEBNext ® rRNA Depletion Kit (Bacteria) and then NEBNext ® rRNA Depletion Kit (Human/Mouse/Rat, New England Biolabs, USA) to eliminate ribosomal RNA (rRNA) that causes sequence biases. This combination of poly A exclusion and rRNA depletion eliminated the most rRNA compared to other combinations of methods (or single methods in isolation) in preliminary tests, yielding sequences with only ~25% rRNA. For the eukaryotic (n = 24) and prokaryotic (n = 24) RNA samples, library preparation was completed using the NEBNext Ultra II Directional RNA library preparation kit (New England Biolabs, Denmark). Libraries pooled in equimolar concentrations were sequenced (5M reads per sample, 2 x 50 bp) at the Center for Functional Genomics and Tissue Plasticity (ATLAS) at SDU, using a NovaSeq 6000 (Illumina, USA).

To ensure good coverage for the final reference transcriptomes, two samples from each species were sequenced deeper in additional runs. For H. stellifera, the sample for deeper sequencing (DC129, sampled under normoxia) was also present in the sample set above. In contrast, for Eurypon sp. 2 another sample (DC128, sampled under normoxia) was used due to the limited availability of sample tissue. As in the sequencing above, two libraries were created and sequenced for each sample -one with eukaryotic RNA and one with prokaryotic RNA. For this deeper sequencing, mRNA in the sponge transcriptomes was purified using the NucleoTrap mRNA kit (Macherey-Nagel, Germany). Prokaryotic RNA was rRNA depleted using the Ribozero Magnetic kit (Illumina, USA). An additional purification step was performed for all RNA samples using the standard protocol for the CleanNGS kit (CleanNA, Netherlands) before sequencing preparation using the NEBNext Ultra II Directional RNA library preparation kit (New England Biolabs, Denmark). Samples were diluted in equimolar concentrations, pooled and sequenced (2 x 150 bp) on a HiSeq X platform (Illumina, USA) by DNASense ApS (Aalborg, Denmark). For the eukaryotic RNA, shallow (2 x 50 bp) sequencing generated 9.20 x 10 8 and 4.38 x 10 8 reads for Eurypon sp. 2 and H. stellifera, respectively, across all sampled sponges, and deeper sequences (2 x 150 bp) generated 2.05 x 10 8 and reads 3.45 x 10 8 for each sponge of the two species, respectively. For the prokaryotic RNA, shallow sequencing yielded 1.27 x 10 9 and 6.68 x 10 8 reads for Eurypon sp. 2 and H. stellifera, respectively. The deeper prokaryotic sequencing generated 1.52 x 10 8 and 1.20 x 10 8 for Eurypon sp. 2 and H. stellifera, respectively (totals are shown in Supplemental Table 2). To obtain reference transcriptomes for each sponge species without prokaryotic symbionts, all eukaryotic RNA sequences were combined into a single file. These sequences were trimmed for quality control using public scripts [START_REF] Matz | Assembly and Annotation of Transcriptomes de Novo[END_REF][START_REF] Meyer | A Collection of Scripts for Analyzing Transcriptome Assemblies[END_REF]) and assembled with Trinity v 2.8.5 [START_REF] Grabherr | Full-Length Transcriptome Assembly from RNA-Seq Data without a Reference Genome[END_REF]. This assembly was filtered for contaminants and annotated using egg-NOG mapper [START_REF] Huerta-Cepas | eggNOG 5.0: A Hierarchical, Functionally and Phylogenetically Annotated Orthology Resource Based on 5090 Organisms and 2502 Viruses[END_REF]. The final sponge transcriptomes included only genes with metazoan annotations or no annotations. For a detailed description of all filtering and assembly steps, see [START_REF] Strehlow | Novel Reference Transcriptomes for the Sponges Carteriospongia Foliascens and Cliona Orientalis and Associated Algal Symbiont Gerakladium Endoclionum[END_REF] and (https://zenodo.org/badge/latestdoi/573040118) for all scripts used in the current paper. Sponge reference transcriptomes fit within the expected range of size, gene number, and N50 (Supplemental Table 2) of published sponge transcriptomes [START_REF] Strehlow | Novel Reference Transcriptomes for the Sponges Carteriospongia Foliascens and Cliona Orientalis and Associated Algal Symbiont Gerakladium Endoclionum[END_REF][START_REF] Riesgo | The Analysis of Eight Transcriptomes from All Poriferan Classes Reveals Surprising Genetic Complexity in Sponges[END_REF]. Additionally, sponge reference transcriptomes were effectively complete based on Benchmarking Universal Single-Copy Orthologs (BUSCO) analysis (Supplemental Table 2, [START_REF] Simão | BUSCO: Assessing Genome Assembly and Annotation Completeness with Single-Copy Orthologs[END_REF] and were therefore suitable for use for subsequent differential expression analysis (see below). Sponge reference transcriptome assemblies are available through NCBI under BioProject number PRJNA89319.

Reference Metagenomes

Two reference metagenomes were sequenced for Eurypon sp. 2 and H. stellifera at DNASense ApS (Aalborg, Denmark). DNA was fragmented to ~550 bp (M220 Focused-ultrasonicator, Covaris, USA) and used for library preparation following the NEBNext Ultra II DNA library preparation kit. The DNA library was paired-end sequenced (2 x 150 bp) on a NovaSeq platform (Ilumina, USA).

Sequences were assembled using Megahit v.1.2.9 [START_REF] Li | MEGAHIT v1.0: A Fast and Scalable Metagenome Assembler Driven by Advanced Methodologies and Community Practices[END_REF]. The final contig fasta files were used to plot GC (%) vs. mean coverage of mapped reads to identify individual bins. For this, we used the custom python scripts spadescontigstocovgc.py (Francis 2020a) as well as the R script contig_gc_coverage.R [START_REF] Francis | Make Interactive Plot of Contigs in Genomes or Metagenomes[END_REF] to extract bins for the Thaumarchaeota and Gammaproteobacteria symbiont of each sponge species. Final statistics for each genome are shown in Supplemental Table 2. Raw sequences are available from NCBI under BioProject number PRJNA893197. Metagenome assemblies are available on https://github.com/bstrehlow/Sponge-transcriptomic-responses-to-seasonal-anoxia.

Mitogenomes

Using only the deep sequences of total ('prokaryotic') RNA, the mitochondrial genomes of single individuals of Eurypon sp. 2 (DC128) and H. stellifera (DC129) were assembled following the pipeline in [START_REF] Plese | Trimitomics: An Efficient Pipeline for Mitochondrial Assembly from Transcriptomic Reads in Nonmodel Species[END_REF]. For future utilization, single contigs that likely represent the mitochondrial genomes of Eurypon sp. 2 and Hymeraphia stellifera were retrieved from the reference metagenome in addition to the coding regions used for differential expression below. The assembled contigs for each metagenome were blasted against the mitochondrial genome of DC128 (Eurypon sp. 2) using blastn, with an e-value cutoff of 1x10 -50 for Eurypon sp. 2 (DC97) and 1x10 -30 for Hymeraphia stellifera (DC297). In general, there was only a single hit for each protein coding gene, confirming a single contig as the mitochondrial genome for each sample. Based on the blast, these two contigs contained all protein coding genes retrieved from the transcriptome that were the focus of this analysis. These 'mitogenomes' can be accessed here: https://zenodo.org/badge/latestdoi/573040118.

Differential Expression Analysis

Trimmed reads from the shallow sequencing for each individual were mapped back to their respective reference transcriptome or genome using bowtie2 (version 2.3.5.1, Langmead and Salzberg 2012) and mapped reads were counted for each contig using samcount.pl (see [START_REF] Matz | Assembly and Annotation of Transcriptomes de Novo[END_REF]. Eukaryotic reads were mapped to the sponge reference transcriptome for each species. Prokaryotic reads were mapped to their respective metagenomic bins -metagenome assembled genomes (MAGs), i.e. Thaumarchaeota or Gammaproteobacteria, within each sponge species. The mitochondrial reference was broken into its 14 genes and prokaryotic (total RNA) reads from both species were mapped to their matching mitochondrial reference and counted as above.

Analyses were performed to compare the gene expression of each target (sponge species, symbiont species, and mitochondria; n=8) under the three different oxygen conditions (anoxia, hypoxia and normoxia). Differential gene expression analysis was performed with the package DESeq2 v. 1.22.2 [START_REF] Love | Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2[END_REF]) in R (R Core Team 2016) with (RStudio Team 2020). Three pairwise comparisons of oxygen condition (normoxia vs. hypoxia, normoxia vs. anoxia, and anoxia vs. hypoxia) were computed from DESeq2 results. Significantly differentially expressed genes between each of these pairwise comparisons were identified using a false discovery rate (FDR) adjusted p-value threshold of 0.1. The ReportingTools software package [START_REF] Huntley | ReportingTools: An Automated Result Processing and Presentation Toolkit for High-Throughput Genomic Analyses[END_REF]) was used to create PDF pages displaying 4 plots of significantly differentially expressed genes and their annotations in various oxygen conditions, as well as box plots, p-values, and LFCs (https://zenodo.org/badge/latestdoi/573040118). These files are searchable in case particular genes are of interest to future studies. We note that comparisons to H. stellifera under normoxia were made with only one normoxic sample due to sampling restrictions. Despite this sampling limitation, the normoxic data for H. stellifera are still shown as a baseline and reference. For each sponge species, log fold changes (LFC, from DESeq2) in all genes were calculated so that positive LFCs were upregulated and negative values were downregulated in hypoxia or anoxia relative to normoxia. For the comparisons between hypoxia and anoxia, positive LFCs were upregulated in anoxia and vice versa.

KOG analysis

EuKaryotic Orthologous Groups (KOG class) assignments were used to compare large scale trends and potential correlations in differential gene expression between Lough Hyne sponges under varied oxygen conditions and external gene expression data from other species and conditions. Delta-ranks, a measure of each KOG class enrichment with up or downregulated genes, were calculated from LFC data using the KOGMWU package in R [START_REF] Dixon | Coral Reefs. Genomic Determinants of Coral Heat Tolerance across Latitudes[END_REF], and delta-ranks were plotted using the package pheatmap [START_REF] Kolde | Pheatmap: Pretty Heatmaps[END_REF]). While it is also possible to calculate delta-ranks with signed p-values (see KOGMWU supporting documentation), LFC was used here to include changes between normoxia and other conditions for H. stellifera. Since there was only one normoxic replicant for this species, LFC was a more appropriate metric than signed p values. Within each pairwise comparison (i.e. hypoxic, anoxic, and anoxic vs hypoxic), significantly enriched KOGs were also calculated using KOGMWU (FDR-adjusted p < 0.05). These KOG enrichments were also compared to those calculated from the sponge Tethya wilhelma after four days of exposure to hypoxia and one hour of exposure to anoxia [START_REF] Mills | The Last Common Ancestor of Animals Lacked the HIF Pathway and Respired in Low-Oxygen Environments[END_REF]. To test whether KOG enrichment in sponges exposed to variable oxygen conditions was correlated with gene expression in other invertebrate stress responses, delta rank data was included from adult Acropora millepora coral exposed to heat stress (31.5 °C for three days, [START_REF] Dixon | Coral Reefs. Genomic Determinants of Coral Heat Tolerance across Latitudes[END_REF]. Similarly, to test the similarity of the sponge data to noted decreases in metabolic activity in invertebrates, KOGs enriched in the dauer (dormant) state in Caenorhabditis elegans [START_REF] Sinha | Divergent Gene Expression in the Conserved Dauer Stage of the Nematodes Pristionchus Pacificus and Caenorhabditis Elegans[END_REF]) and the diapause state of the midge Sitodiplosis mosellana [START_REF] Gong | Global Transcriptome Analysis of Orange Wheat Blossom Midge, Sitodiplosis Mosellana (Gehin) (Diptera: Cecidomyiidae) to Identify Candidate Transcripts Regulating Diapause[END_REF] were added using delta rank tables calculated by [START_REF] Strader | Red Fluorescence in Coral Larvae Is Associated with a Diapause-like State[END_REF]. Correlations between the KOG delta ranks in different datasets were performed using Pearson's Correlation in KOGMWU. Within the significantly enriched KOGs 'replication, recombination and repair' and 'energy production and conversion,' genes that were significantly differentially expressed were plotted in heatmaps to further investigate sponge transcriptomic responses to deoxygenation at the gene level in Lough Hyne sponges. KOG enrichment analysis was performed as above for each of the two major symbionts within the two sponge species; however, Clusters of Orthologous groups (COGs) were used since they are designed for prokaryotes rather than eukaryotic KOGs and are fewer in number.

Gene ontology (GO) term enrichment analysis

Log fold changes calculated above were used for rank-based gene ontology (GO) analysis with adaptive clustering (GO_MWU, [START_REF] Wright | Gene Expression Associated with White Syndromes in a Reef Building Coral, Acropora hyacinthus[END_REF] to determine what GO terms were significantly enriched under various oxygen conditions. These analyses were run individually for each pairwise comparison of oxygen condition in both sponge species. These analyses were not completed for the two symbiont types because GO terms were not available for these microbes.

Phylogenetic positions of nitroreductase proteins

To determine if nitroreductases in sponge symbionts are related to genes that may be involved in oxygen production in N. maritimus, nitroreductase proteins from the Thaumarchaeota symbionts of H. stellifera and E. sp. 2 were extracted using annotations assigned by egg-NOG. Phylogenetic position of nitroreductase proteins from H. stellifera and E. sp. 2 were investigated by aligning them to other known nitroreductase proteins from Genbank using Mafft v7 [START_REF] Katoh | MAFFT Multiple Sequence Alignment Software Version 7: Improvements in Performance and Usability[END_REF]. The ModelTest-NG v.0.1.7 function [START_REF] Edler | raxmlGUI 2.0: A Graphical Interface and Toolkit for Phylogenetic Analyses Using RAxML[END_REF]) implemented within RaxmlGUI 2.0 was used to calculate the evolutionary model. According to AIC, LG+I+G was chosen as the best-fitting model. Maximum Likelihood (ML) phylogeny was then reconstructed using raxmlGUI v.2.0.6 [START_REF] Edler | raxmlGUI 2.0: A Graphical Interface and Toolkit for Phylogenetic Analyses Using RAxML[END_REF] with 100 (rapid) bootstraps. The resulting protein tree was rendered by FigTree v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).

Results

Sponge differential expression

There were three pairwise comparisons for significantly differentially expressed genes -hypoxia vs. normoxia, anoxia vs. normoxia, and anoxia vs. hypoxia. The number of significantly differentially expressed genes in these comparisons for sponges, their mitochondria and their symbionts are summarized in Table 1. In the sponge Eurypon sp. 2, 737, 1,376 and 1,090 genes were significantly differentially expressed under hypoxia (vs. normoxia), anoxia (vs. normoxia), and anoxia vs. hypoxia, respectively (p < 0.1, Table 1). The differentially expressed genes under each oxygen condition ranged between 1.5-2.8% of the genes of the Eurypon sp. 2 reference transcriptome, and 41-48% of these genes are annotated, which reflects the overall annotation percentage of 39.7% for the whole transcriptome (Supplemental Table 2). In hypoxia, the majority of these differentially expressed genes were upregulated (68%). Under anoxia, on the other hand, the majority of genes were downregulated (60%). In anoxia vs. hypoxia, 57% of genes were upregulated in anoxia compared to hypoxia in Eurypon sp. 2 (Table 1). Table 1. Number of significantly differentially expressed genes (DEGs) under deoxygenation in the sponge holobionts summarized per reference genome or transcriptome. The direction (up or down) within each paired oxygen condition comparison is relative to normoxia or in the case of anoxia vs. hypoxia, hypoxia. The number of annotated genes is shown in parentheses. In the sponge H. stellifera, 523, 284, and 2,154 genes were significantly differentially expressed in hypoxia, anoxia, and in anoxia vs. hypoxia, respectively (p < 0.1, Table 1). Differentially expressed genes in the different oxygen conditions corresponded to 0.42-3.4% of the total number of genes in the H. stellifera transcriptome. Approximately 25-37% of the differentially expressed genes were annotated, which, as in the other sponge species, roughly corresponded to the global percentage (37.0%) of genes annotated within the H. stellifera reference transcriptome. Under hypoxia in this species, the vast majority of differentially expressed genes (89%) were upregulated. As in Eurypon sp. 2, the majority of differentially expressed genes (84%) were downregulated in anoxia in H. stellifera. The majority of differentially expressed genes between anoxia and hypoxia were upregulated in anoxia (71%, Table 1).

Hypoxia

0 (0) 0 (0) 0 (0) 1 (1) 0 (0) H. stellifera mitochondria 14 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) Eurypon sp. 2 Thaumarchaeota 1,556 19 (14) 1 (1) 3 (3) 0 (0) 0 (0) 0 (0)

H. stellifera Thaumarchaeota

KOG enrichments in sponges

The global gene expression patterns of Lough Hyne sponges in response to deoxygenation were compared across the two species by assessing up and down regulation of KOG classes (broad functional annotations assigned to genes). Within these sponges, KOG expression in response to hypoxia and anoxia (versus normoxia) were largely distinct in both species (Figure 1A). However, KOG expression in H. stellifera in hypoxia vs. anoxia was significantly positively correlated, i.e. similar, to that of Eurypon sp. 2 under hypoxia (r = 0.43, p < 0.05, Figure 1E), but global expression patterns were only similar between the two species in this single case. ). Unhighlighted columns represent external datasets. Up or down regulation with respect to normoxia is indicated in 'Hypox' (hypoxia) and 'Anox' (anoxia) columns, and differential regulation with respect to hypoxia is shown in 'Hypox vs. Anox' columns, i.e. upregulated KOGs (red delta ranks) are upregulated in anoxia with respect to hypoxia. Within each column on the heatmap, KOGs that are significantly enriched (FDRadjusted p < .05) are marked with asterisks. KOG enrichments from the sponge Tethya wilhelma after four days exposure to hypoxia (Tw Hypox Long) and one hour exposure to anoxia (Tw Anox Short), the coral Acropora millepora after three days of exposure to heat stress (Am Coral Heat), Caenorhabditis elegans in a dauer state (Ce Dauer), and the midge Sitodiplosis mosellana in diapause (Sm Diapause) from external studies were included for comparative purposes. B-E. Significant (p < 0.05) correlations between the KOG enrichments for different species and oxygen levels from this study. F-J. Significant correlations between KOG enrichments in this study and those of external datasets. All other possible comparisons between the various datasets were not significant. Within Eurypon sp. 2, KOG class expression patterns in hypoxia vs. anoxia were significantly positively correlated to those of hypoxia (r = 0.71, p < 0.01, Figure 1B). Similarly, in H. stellifera, KOG enrichment in hypoxia vs. anoxia was positively correlated to that of hypoxia (r = 0.72, p < 0.001, Figure 1C) and anoxia (r = 0.46, p < 0.05, Figure 1D). All significant correlations between datasets are shown in (Figure 1B-J), and all relationships that were not significantly correlated are not shown. At the level of individual KOGs, the KOG class 'replication recombination and repair' (RRR) was significantly enriched in all pairwise comparisons of oxygen conditions for both sponge species (Padj < 0.05), though in different directions depending on oxygen. It was significantly upregulated in anoxia and anoxia vs. hypoxia but downregulated in hypoxia in both species (Figure 1A). In Eurypon sp. 2, five genes with RRR KOG annotations were significantly upregulated in hypoxia, e.g. HELZ and HELZ2 transcription associated helicases, and 11 genes were significantly downregulated (p < 0.1, Supplemental Figure 4A), including five genes from the minichromosome maintenance (MCM) protein complex. In anoxia, two genes were significantly upregulated within RRR in this species -'siz regulation of ARF protein signal transduction' and DNA polymerase epsilon catalytic subunit (p < 0.1, Supplemental Figure 4B). The majority of significantly differentially expressed RRR genes (five) were downregulated in anoxia (p < 0.1, Supplemental Figure 4B). The majority of RRR genes (nine) were also downregulated in anoxia vs. hypoxia, with only two significantly upregulated in anoxia (Supplemental Figure 4B). The genes NBN and RFC1 were downregulated in anoxia compared to both normoxia and hypoxia (Supplemental Figure 4B,C). Similarly, the gene 'siz regulation of ARF protein signal transduction' was upregulated in anoxia compared to both normoxia and hypoxia (Supplemental Figure 4B,C).

In H. stellifera, three genes were significantly differentially expressed in the RRR KOG class in hypoxia, and all were upregulated (p < 0.1, Supplemental Figure 4D). No genes within this class were significantly differentially expressed in anoxia (vs. normoxia). Eleven genes were upregulated and eight genes were downregulated in anoxia vs. hypoxia (Supplemental Figure 4E). Similar to Eurypon sp. 2, HELZ2 was downregulated in anoxia vs. hypoxia in H. stellifera. In contrast to Eurypon sp. 2 a distinct MCM gene (MCM6) was downregulated in anoxia vs. hypoxia. Additionally, four genes from the structural maintenance of chromosomes protein 6 (SMC6) complex were downregulated in anoxia vs. hypoxia. Like RRR, the KOG class 'energy production and conversion' (EPC) was significantly enriched (Padj < 0.05) between all conditions except for anoxia vs. hypoxia in H. stellifera (Padj = 0.07). EPC was significantly upregulated in hypoxia in H. stellifera but downregulated in all remaining conditions in both sponge species (Figure 1A). Within Eurypon sp. 2 specifically, four genes were significantly downregulated and twelve genes were upregulated within the EPC class in hypoxia (p < 0.1, Supplemental Figure 5A). Various nuclear genes for processes in the mitochondria were differentially expressed in hypoxia. SLC25A27 (a mitochondrial membrane uncoupling protein), ndufs8 (a NADH dehydrogenase core subunit in the mitochondrial membrane -complex I), and Surf1 (part of the complex that regulates cytochrome c oxidase assembly) were downregulated in hypoxia (Supplemental Figure 5A). However, other nuclear mitochondrial genes were upregulated, including four subunits of ATP synthase and two additional genes with annotations corresponding to ATP synthesis (Supplemental Figure 5A). In anoxia, fourteen genes were upregulated and five genes were downregulated in the EPC class (Supplemental Figure 5B). As in hypoxia, ndufs8 was downregulated, but no other significant DEGs were shared across the two oxygen levels. The following genes were also downregulated in anoxia: mitochondrial sulfite oxidase (SUOX), malate dehydrogenase (MDH1), and two genes belonging to the ATP gaunido phosphatase family (Supplemental Figure 5B). Conversely, four genes involved in phosphoenolpyruvate carboxykinase activity were significantly upregulated in anoxia (p < 0.1, Supplemental Figure 5B). Additionally, three mitochondrial membrane transporter proteins from the SLC family, two CYB5R4 cytochrome reductase genes, a mitochondrial carrier protein, and fdx1 cytochrome reductase were upregulated in anoxia. In anoxia vs. hypoxia in Eurypon sp. 2, ten genes were upregulated in anoxia (Supplemental Figure 5C). Upregulated genes included all of the same genes that were significantly upregulated in anoxia, except for the phosphoenolpyruvate carboxykinase associated genes. Six genes were significantly downregulated in anoxia vs. hypoxia. Of these, three were also downregulated in anoxia (vs. normoxia) and two genes were significantly upregulated in anoxia (Supplemental Figure 5A-C). It is noteworthy, however, that all of the upregulated genes in anoxia vs. normoxia and hypoxia were only strongly upregulated in one of the two anoxia replicates (DC56, Supplemental Figure 5 B,C). In H. stellifera, eight genes were upregulated and two genes were downregulated in hypoxia. As in Eurypon sp. 2, the ATP5A1 ATP synthase subunit was upregulated in hypoxia. An additional ATP synthase subunit (ATP6V0C) that was not significantly upregulated in Eurypon sp. 2 was upregulated in H. stellifera in hypoxia (p < 0.1, Supplemental Figure 5D). Mitochondrial sulfide:quinone oxidoreductase (SQRDL), three genes associated with phosphoenolpyruvate carboxykinase, and a histone demethylase (KDM6A) were also significantly upregulated in hypoxia in H. stellifera. Conversely, PCK1 (also associated with phosphoenolpyruvate carboxykinase) and MDH1 were significantly downregulated in hypoxia (p < 0.1, Supplemental Figure 5D). In anoxia in H. stellifera, ATP5A1 was significantly upregulated, and ATP6V0C, PCK1, and MDH1 were significantly downregulated (Supplemental Figure 5E). Comparing anoxia vs. hypoxia, eight genes were upregulated and five genes were downregulated (Supplemental Figure 5F). As in hypoxia, SQRDL, four genes associated with phosphoenolpyruvate carboxykinase activity, and KDM6A were upregulated in anoxia vs. hypoxia. Additionally, two unshared genes, including the mitochondrial Succinyl-CoA:3-ketoacid coenzyme A transferase 1 (OXCT1) and a mitochondrial carrier protein, were upregulated in anoxia vs. hypoxia (Supplemental Figure 5F). Five genes were downregulated in anoxia vs. hypoxia in H. stellifera. Two of these four only exhibited high expression levels in one hypoxic sample (DC105) that drove the observed difference (Supplemental Figure 5F). The other three downregulated genes in anoxia vs. hypoxia were mitochondrial creatine kinase (CKMT2), NAD P transhydrogenase (NNT), and proto-oncogene tyrosine-protein kinase (Src). The 'cytoskeleton' KOG class was also significantly enriched (Padj < 0.01) in all but one oxygen conditionanoxia in Eurypon sp. 2 (Padj = 0.1). Interestingly, enrichment patterns in this KOG were species specific, being significantly upregulated in hypoxia and anoxia vs. hypoxia in Eurypon sp. 2 and downregulated in all conditions for H. stellifera. The other twenty KOG classes were only significantly enriched in Lough Hyne sponges under a few (or no) different oxygen conditions and patterns of significant enrichment were generally condition and species-specific (Figure 1A). With the exception of EPC, KOGs pertaining to metabolic processes were rarely significantly enriched. Expression of the KOG classes 'carbohydrate metabolism' and 'amino acid transport and metabolism' were not significantly changed under any oxygen condition (Padj > 0.05). There was some significant enrichment of metabolic KOGs, but only in Eurypon sp. 2. 'Coenzyme transport and metabolism' and 'inorganic ion transport and metabolism' were both downregulated in anoxia, and 'nucleic acid transport metabolism' was downregulated in hypoxia vs. anoxia. For both these conditions (anoxia and hypoxia vs. anoxia), 'Lipid transport and metabolism' and 'Secondary metabolites biosynthesis, transport and catabolism' were significantly downregulated. However, the transcription of the majority of genes labeled with metabolic KOGs remained unchanged under anoxia and hypoxia. Comparing the Lough Hyne dataset with external datasets of invertebrate KOG expression in response to stress, only two positive correlations were found. In Eurypon sp. 2, expression patterns under anoxia vs. hypoxia were comparable to those of T. wilhelma under short-term anoxia (r = 0.44, p < 0.05, Figure 1H). The hypoxic response of Eurypon sp. 2 resembled that of diapause in S. mosellana (r = 0.68, p < 0.001, Figure 1F), but this was the only correlation found between Lough Hyne sponge data and diapause or dauer states. Importantly, KOG enrichments under anoxia (versus normoxia) were not positively correlated to any other external datasets for either sponge species. Three significant negative correlations, which corresponded to opposing KOG expression patterns, were observed between Lough Hyne sponges and external datasets. Hypoxia responses in H. stellifera were negatively correlated with those of A. millepora under heat stress (r = -0.56, p < 0.01, Figure 1G). The long exposure of T. wilhelma to hypoxia was negatively correlated to both anoxia (r = -0.64, p < 0.01, Figure 1J) and hypoxia vs. anoxia (r = -0.61, p < 0.01, Figure 1I) in Eurypon sp. 2. GO term enrichments in sponges Generally, significantly enriched GO terms between the different treatments were not shared between the two sponge species (Supplemental Figures 123). The strongest upregulated GO enrichments (p < 0.01) in hypoxia in Eurypon sp. 2 were cellular signaling and 'lipoprotein particle receptor binding' (Supplemental Figure 1A). For H. stellifera, on the other hand, the strongest upregulated enrichments (p < 0.01) in hypoxia included the terms 'protein tyrosine kinase binding,' 'molecular adaptor,' 'structural molecule,' and 'structural constituent of ribosome'; whereas the downregulation of the term 'motor' was strong (p < 0.01, Supplemental Figure 1B). Under anoxia, strong upregulation of the 'proteasome binding' term was observed in Eurypon sp. 2 (p < 0.01), and the terms 'phosphoenolpyruvate carboxykinase', 'carbohydrate kinase,' monovalent inorganic cation transmembrane transporter,' and 'sodium ion transmembrane transporter' were highly downregulated (p < 0.01, Supplemental Figure 2A). In H. stellifera under anoxia, the terms 'damaged DNA binding' and 'calcium-dependent phospholipid binding' were highly upregulated, and four terms associated with various aspects of the cytoskeleton were downregulated (p < 0.01, Supplemental Figure 2B). In Eurypon sp. 2 under anoxia vs. hypoxia, the terms 'lipoprotein receptor binding' and 'proteasome binding' were strongly upregulated in anoxia, but no terms were strongly downregulated in anoxia versus hypoxia (p < 0.01, Supplemental Figure 3A). The former term was also upregulated under anoxia vs. normoxia. Conversely, the strongest enrichments in H. stellifera in anoxia vs. hypoxia were downregulated in anoxia and were composed of three of the four cytoskeletal binding GO terms that were also downregulated in anoxia vs. normoxia (see above, p < 0.01, Supplemental Figure 3B), but no strong upregulation of GO terms was observed.

Considering GO enrichments at lower significance levels (0.01 < p < 0.1), within Eurypon sp. 2, the GO terms 'proteasome binding', 'lipoprotein particle receptor binding,' and 'Hsp70 binding' (HSP: Heat shock protein) were significantly upregulated in all cases of decreased oxygen albeit at varied significance levels (p < 0.05, Supplemental Figures 1-3A). Between only two levels, the terms 'transmembrane signaling receptor,' 'secondary active transmembrane transporter,' and 'molecular transducer' were upregulated in hypoxia but downregulated in anoxia. However, the term lipoprotein particle receptor binding was downregulated in anoxia and upregulated in hypoxia (p < 0.1, Supplemental Figures 1,2A). Additionally, the term 'structural constituent of ribosome' was upregulated in anoxia but downregulated in anoxia vs. hypoxia (p < 0.1, Supplemental Figures 1,3A). No additional GO terms other than the three present in all oxygen levels for Eurypon sp. 2 were enriched in both the pairwise comparison of anoxia and anoxia vs. hypoxia. No GO terms were significantly up or downregulated in all three levels of oxygen conditions for H. stellifera. However, some GO terms were significantly enriched in two levels. The term 'ATPase, coupled to transmembrane movement of ions, rotational mechanism' was downregulated in hypoxia but upregulated in anoxia (p < 0.1, Supplemental Figures 1B,2B respectively). The GO terms 'I-SMAD binding', thioesterase binding,' and 'ubiquitin-like protein ligase binding' were all upregulated in hypoxia (vs. normoxia) and anoxia vs. hypoxia, while the 'motor' term was downregulated (p < 0.1, Supplemental Figures 1B,3B). There were only two shared GO enrichments between both sponge species for any given oxygen condition. The terms 'protein tyrosine kinase binding' and 'structural constituent of the ribosome' were upregulated in both species in hypoxia (p < 0.1, Supplemental Figure 1), although they were enriched to different degrees of significance. No significantly enriched GO terms were shared between species in normoxia vs. anoxia or hypoxia vs. anoxia. At the gene level, significantly differentially expressed heat shock proteins (Hsps) in both sponges are shown in Supplemental Figure 8. In Eurypon sp. 2, one Hsp gene (DNAJC8/Hsp40) was downregulated in hypoxia (p < 0.01), but four Hsp genes (HSPA8, two genes for HSP90B1, and Hsp70) were significantly upregulated in hypoxia (p < 0.1, Supplemental Figure 8A). Additionally, HSPA8 and C11orf73 (Hsp70 binding protein) were upregulated in anoxia vs normoxia (p < 0.05, Supplemental Figure 8A). Another Hsp70 gene was also significantly downregulated in anoxia vs. hypoxia (p < 0.05, Supplemental Figure 8A). In H. stellifera, HSP90A1 was upregulated in hypoxia and anoxia (p < 0.1, Supplemental Figure 8B). Two additional Hsp genes (both for HSP90B1) were also upregulated in anoxia. Finally, the two HSP90B1 genes were upregulated in anoxia vs. hypoxia, but HSPA5 was downregulated in anoxia vs. hypoxia (p < 0.05, Supplemental Figure 8B).

Mitochondrial gene expression

All 14 coding genes from the sponge mitochondria were transcribed in both species and in all oxygen conditions. There was only one mitochondrial gene (nad2) that was significantly differentially expressed but only in Eurypon sp. 2. This gene was significantly upregulated (p < 0.1) in anoxia vs hypoxia (Supplemental Figure 6A, Table 1). Similar counts for each individual mitochondrial gene across both sponge species and all oxygen conditions were found. For instance, counts for atp6 were consistently twice as high as counts for atp8 (Supplemental Figure 6).

Prokaryotic differential expression Thaumarchaeota gene expression

For the Thaumarchaeota associated with Eurypon sp. 2, 20, 3 and 0 genes were differentially expressed in hypoxia, anoxia, and hypoxia vs. anoxia, respectively (p < 0.1, Table 1). These genes accounted for 0-1.3% of this Thaumarchaeota genome, and 75-100% of the differentially expressed genes were annotated (Table 1), compared to 69.2% annotated in the entire reference genome. All but one of the differentially expressed genes in hypoxia were upregulated (Table 1, Figure 2A). The downregulated gene belonged to the sodium solute symporter family (SSF, TC 2.A.21, Figure 2A). The upregulated genes in hypoxia included three stress response genes (COG0589 Universal stress protein [UspA], small heat shock protein [HSP20], and DnaJ molecular chaperone homology domain), two CBS domains, a member of the lactate/malate dehydrogenase (LDH MDH) superfamily, a tetratricopeptide (TPR) repeat, two transcription regulators, and a PII signal transduction protein. Figure 2. A-D. All significantly differentially expressed genes for the Thaumarchaeota in both sponge species (left: E. sp 2, right: H. stellifera. For Thaumarchaeota associated with E. sp. 2, only two pairwise comparisons of oxygen conditions had significant differences: A. hypoxia vs. normoxia and C. anoxia vs. normoxia. In symbionts associated with H. stellifera, there were significantly differentially expressed genes in hypoxia vs. normoxia (B) and anoxia vs. hypoxia (D). There were no significant differences in the three genes involved in ammonia oxidation (AmoA, B and C) in either symbiont, but their heatmaps are shown in E and F for E. sp 2 and H. stellifera, respectively.

Only three genes were significantly differentially expressed in anoxia in Eurypon sp. 2's Thaumarchaeota, and all were upregulated (Table 1, Figure 2C). The upregulated genes in anoxia were a UbiD family decarboxylase, CAAX protease self-immunity, and the same TPR repeat that was also upregulated in hypoxia. There were no significantly differentially expressed genes between anoxia and hypoxia in this Thaumarchaeota (Table 1). The genes necessary for ammonia oxidation, i.e. subunits A, B and C of ammonia monooxygenase (Amo), were present in this Thaumarchaeota, but they were not differentially expressed under deoxygenation (Figure 2E). In the Thaumarchaeota symbiont of H. stellifera, 17, 0 and 28 genes were significantly differentially expressed in hypoxia, anoxia and hypoxia vs. anoxia, respectively. Similar to the symbiont of Eurypon sp. 2, all 17 differentially expressed genes in the H. stellifera symbiont were upregulated in hypoxia. These genes included the stress response genes 'universal stress protein, an HSP20 family protein, and two copies of UspA. HSP20 and UspA were also upregulated in the Eurypon sp. 2 holobiont. Similarly, a PII signal transduction protein, CBS domain, and two transcription regulators (one identical: TBP stabilizer, and one distinct to H. stellifera symbionts) were upregulated in both Thaumarchaeota species in hypoxia. Unlike in Eurypon sp. 2, the Thaumarchaeota in H. stellifera upregulated a helix turn helix ASNC type domain and a Lrp/AsnC ligand binding domain (Figure 2B). No genes were significantly differentially expressed in anoxia in this Thaumarchaeon. As in the other holobiont species, Thaumarchaeota genes for ammonia oxidation were not significantly differentially expressed with H. stellifera symbionts (Figure 2F). Twenty-seven genes were upregulated and one gene (K141_339108_19, unannotated) was downregulated in anoxia vs. hypoxia in the Thaumarchaeota within H. stellifera (Figure 2D). Most of the upregulated genes in anoxia vs. hypoxia were also upregulated in hypoxia. Upregulated genes in anoxia vs. hypoxia unique to the H. stellifera holobiont included additional transcriptional regulators, an integrase domain, and a putative nitroreductase.The phylogenies of this upregulated nitroreductase and all nitroreductases from the studied holobionts as well as publicly available nitroreductase sequences from other Thaumarchaeota are shown in Supplemental Figure 9. For both Thaumarchaeota, only two COGs ('signal transduction mechanisms' and 'transcription') were significantly enriched (Padj < 0.05, Supplemental Figure 7A). These two COGs were both downregulated in hypoxia and anoxia vs. hypoxia within H. stellifera symbionts. All other COG classes were comparable in their expression levels regardless of holobiont species or oxygen concentration (Supplemental Figure 7A).

Gammaproteobacteria gene expression

The Gammaproteobacteria associated with Eurypon sp. 2 only exhibited significant differential expression under hypoxia but not in any other oxygen condition (Table 1). In hypoxia, three genes were upregulated: type 1 restriction-modification system, flavin adenine dinucleotide binding, and signal sequence binding (Figure 3A). The 29 downregulated genes in hypoxia included two genes involved in sulfur processing (k212078_2, k237243_12), cytochrome c (class I), a gene associated with the assembly of the periplasmic c-type cytochromes (k561840_2), a NQR complex gene, various genes involved in pore formation, a Rhodaneselike gene, glutamine synthetase (type I), and YbiS YcfS YnhG family protein. In the Gammaproteobacteria associated with H. stellifera, genes were only significantly differentially expressed between anoxia and hypoxia (436 genes, p <0.1, Table 1, Figure 3B,C). Twenty genes were upregulated in anoxia vs. hypoxia and all were annotated (Figure 3B). Upregulated genes in anoxia vs. hypoxia included five genes associated with the ribosome (k141_198365_15, k141_198365_8, k141_198365_17, and k141_245017_3), various genes for molecular transport, a transmembrane proton pump (k141_226857_1), succinyl-CoA synthetase, and formate dehydrogenase (alpha subunit). There were also 416 genes that were significantly downregulated in anoxia vs. hypoxia (Table 1). The top ten upregulated genes (lowest p values) are exhibited in Figure 3C. Gammaproteobacterial COG expression is shown in Supplemental Figure 7B. Only one COG class, 'translation, ribosomal structure and biogenesis,' was significantly enriched in Gammaproteobacteria under the different oxygen conditions (Padj < 0.05). This COG class was downregulated in Gammaproteobacteria associated with Eurypon sp. 2 under hypoxia as well as those of H. stellifera in anoxia vs. hypoxia (Figure 7B).

Discussion

Sponge holobiont gene expression was significantly changed by exposure to deoxygenation in situ, though responses were generally species-specific, dependent on oxygen level, and varied across holobiont members. The potential functions of this differential expression are discussed below for each member of the two holobionts: sponge, sponge mitochondria, Thaumarchaeota and Gammaproteobacteria. The implications of this study for past and future oceans are considered.

Differential expression in sponges during deoxygenation

Cellular stress responses were observed in the expression of HSPs, a group of stress response genes that have been widely studied in marine invertebrates. These proteins are molecular chaperones involved in protein folding, unfolding, transport and the assembly of protein complexes, and they are present in nearly all organisms [START_REF] Li | Heat-Shock Proteins[END_REF]. Both Lough Hyne sponges upregulated HSP genes in response to hypoxia and anoxia (vs normoxia), though there was downregulation of some HSP genes in deoxygenation as well. Eurypon sp. 2 and H. stellifera both upregulated Hsp90 genes in hypoxia, but they upregulated different HSPs in anoxia: HSPA8 and Hsp90 for Eurypon sp. 2 and H. stellifera, respectively. It should be noted that Hsp90 upregulation in H. stellifera in anoxia was driven by only one replicate. However, in the pairwise comparison with the highest replication in this study, hypoxia vs. normoxia in Eurypon sp. 2, Hsp90 genes were significantly upregulated in hypoxia. Hsp90 expression was also higher in anoxia than normoxia for Eurypon sp. 2, but this was not significant. Based on these trends, it is possible that increased replication would result in a clearer upregulated signal from Hsp90 genes, since false negatives in gene expression analyses decrease with more replication [START_REF] Schurch | How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use?[END_REF]. Eurypon sp. 2 also upregulated HSPA8 and Hsp70 in hypoxia but downregulated Hsp40. In long term hypoxia, T. wilhelma exhibited a similar HSP expression pattern, upregulating Hsp90 but downregulating Hsp70 [START_REF] Mills | The Last Common Ancestor of Animals Lacked the HIF Pathway and Respired in Low-Oxygen Environments[END_REF]. Sponges under heat stress also upregulate Hsp70 [START_REF] López-Legentil | Bleaching and Stress in Coral Reef Ecosystems: hsp70 Expression by the Giant Barrel Sponge Xestospongia Muta[END_REF][START_REF] Guzman | Gene Expression Dynamics Accompanying the Sponge Thermal Stress Response[END_REF][START_REF] Webster | A Complex Life Cycle in a Warming Planet: Gene Expression in Thermally Stressed Sponges[END_REF]) and HSP90 [START_REF] Guzman | Gene Expression Dynamics Accompanying the Sponge Thermal Stress Response[END_REF]. Pollutants such as lead and zinc can also stimulate the upregulation of HSPs in sponges [START_REF] Efremova | Heat Shock Protein Hsp70 Expression and DNA Damage in Baikalian Sponges Exposed to Model Pollutants and Wastewater from Baikalsk Pulp and Paper Plant[END_REF][START_REF] Schröder | Stress Response in Baikalian Sponges Exposed to Pollutants[END_REF]. HSPs are therefore crucial to sponge stress response. HSP expression patterns and other key transcriptomic responses of the sponge holobionts are outlined in Figure 4. In metazoans with a complete HIF pathway, Hsp90 also serves to stabilize HIFa and therefore plays an important role in hypoxia tolerance (Semenza 2007). In the absence of HIFa in sponges, Hsp90 may play a greater role in oxygen homeostasis given that it is upregulated under hypoxia in all sponges that have been examined transcriptomically to date, i.e. Eurypon sp. 2, H. stellifera and T. wilhelma, and that Hsp90 is upregulated in H. stellifera under anoxia. Hsp90 could stabilize currently unknown proteins that are involved in sponge oxygen regulation or modulate transcription factors other than HIFa. Further research is needed to determine the exact role of Hsp90 in sponge oxygen homeostasis, but it presents a good candidate for a gene underlying sponge tolerance to deoxygenation.

Patterns of species and oxygen-condition-dependent expression of GO terms were observed in sponges in response to deoxygenation. GO terms such as 'proteasome binding' and 'Hsp70 binding' were upregulated in all cases of decreased oxygen in Eurypon sp. 2, suggesting that the species-specific response to deoxygenation in Eurypon sp. 2 may center around repairing or degrading damaged proteins. In H. stellifera, on the other hand, upregulated GO terms related to DNA repair in anoxia, including 'damaged DNA binding', 'DNA binding,' and 'catalysis of DNA and RNA,' so DNA repair may be more central to the deoxygenation tolerance of this species. Both protein and DNA repair upregulation might be a response to reactive oxygen species (ROS), which can become more abundant under hypoxia due to imbalances in the mitochondrial electron transport chain [START_REF] Görlach | Reactive oxygen species, nutrition, hypoxia and diseases: Problems solved?[END_REF]. If there is any oxygen present in the sponges under environmental anoxia, e.g. oxygen produced by Thaumarchaeota (see [START_REF] Kraft | Oxygen and Nitrogen Production by an Ammonia-Oxidizing Archaeon[END_REF], Thaumarchaeota gene expression below), then ROSs could cause cellular stress in anoxia. Otherwise, the presence of sulfide in the water under anoxia [START_REF] Schuster | Effects of Seasonal Anoxia on the Microbial Community Structure in Demosponges in a Marine Lake in Lough Hyne, Ireland[END_REF] or ROSs other than superoxide (O2 . ) could damage proteins and DNA in anoxic conditions.

The EPC KOG class was significantly downregulated in all oxygen conditions in both sponges except for hypoxia and hypoxia vs. anoxia in H. stellifera, when it was significantly upregulated and unchanged, respectively (Figure 1,4). Within the EPC class, however, there were speciesspecific responses to anoxia and hypoxia at the level of individual differentially expressed genes. Most nuclear ATP synthase genes were upregulated in hypoxia in Eurypon sp. 2 but were unchanged under anoxia. Similarly, two ATP synthase subunits (ATP5A1 and ATP6V0C) were upregulated in hypoxia in H. stellifera. In anoxia, ATP5A1 was also upregulated, but ATP6V0C was downregulated. Additionally, downregulation of SLC25A27 in Eurypon sp. 2 in hypoxia may increase the coupling across the mitochondrial membrane to improve ATP synthesis. We note that the electron transport chain genes Ndufs8 and Surf1 were also downregulated in hypoxia Eurypon sp. 2. The NNT gene, also involved in electron transport, was downregulated in anoxia vs. hypoxia in H. stellifera. Nevertheless, most of the genes associated with the respiration electron transport chain were not differentially expressed, and electron transport likely continued under deoxygenation. Sponges may therefore have increased ATP production under hypoxia but not under anoxia. In a similar way, mollusks may upregulate mitochondrial electron transport system activity and ATP synthesis capacity to improve respiration efficiency under hypoxia [START_REF] Hochachka | Mechanism, Origin, and Evolution of Anoxia Tolerance in Animals[END_REF].

Glycolysis, which can produce ATP and NADH without oxygen (Stryer 1995), was expected to increase under anoxia. However, relatively few genes directly or indirectly involved in glycolysis were differentially expressed in our sponges under deoxygenation. In Eurypon sp. 2, pantothenate kinase 1 (PANK1), ATP-dependent 6-phosphofructokinase (PFKL) and Phosphoenolpyruvate carboxykinase 2 (PCK2) were significantly upregulated in anoxia, and only aldehyde dehydrogenase 3 family, member a2 (ALDH3A2) was significantly upregulated in hypoxia (see Supplemental Reports, Supplemental Figure 5). In H. stellifera, acetyl-coenzyme A synthetase 1 (ACSS1) was significantly downregulated in hypoxia, but PCK1 and pyruvate kinase (PKM) were upregulated in hypoxia. In anoxia for H. stellifera, beta-enolase (ENO3) and PCK1 were upregulated. Based on these results, it is unlikely that glycolysis rates increased under anoxia and hypoxia in either species. It is more probable that glycolysis continued under anoxia and hypoxia at normal rates in both species, since gene expression in the other ~70 genes associated with glycolysis was not significantly changed under deoxygenation (Figure 4). Either way, the sponges would have at least some access to anaerobically generated ATP through glycolysis during environmental deoxygenation.

In the KOG class 'replication, recombination and repair' (RRR), there was significant upregulation in anoxia and anoxia vs. hypoxia, but downregulation in hypoxia (vs. normoxia) for both sponge species. In Eurypon sp. 2, many MCM genes involved in the initiation of DNA replication were downregulated in hypoxia but remained unchanged in anoxia. Certain aspects of DNA replication, e.g. the DNA polymerase epsilon catalytic subunit, may even be upregulated in anoxia. The pattern was less clear in H. stellifera, but MCM genes were downregulated in anoxia compared to hypoxia (Figure 4). These data indicated that DNA replication may decrease under hypoxia but remain stable during anoxia, possibly due to increased oxidative stress under hypoxia compared to normoxia (see above). This pattern of a 'stronger' response to hypoxia than anoxia mirrored that of nuclear ATP synthesis genes.

Between the two study species, there was only one significant, positive correlation of KOG class expression between hypoxia in Eurypon sp. 2 and hypoxia vs. anoxia in H. stellifera, demonstrating the general species-specific responses of most KOGs to deoxygenation. Within a single sponge species, no correlations were seen between KOG enrichments in anoxia or hypoxia (vs. normoxia), highlighting the transcriptomic differences between responses to hypoxia and anoxia. KOG enrichments in hypoxia versus anoxia, on the other hand, were positively correlated with those of anoxia in Eurypon sp. 2 and both anoxia and hypoxia in H. stellifera.

When quantitatively compared to transcriptomic datasets from the sponge T. wilhelma using KOGMWU, transcriptomic responses in Lough Hyne sponges were only positively correlated between T. wilhelma under short-term anoxia and Eurypon sp. 2 in hypoxia vs. anoxia (Figure 1H). This positive correlation could indicate similar transcription patterns demonstrated by T. wilhelma and Eurypon sp. 2 transitioning from hypoxia to anoxia. However, no further similarities between the responses of T. wilhelma and Lough Hyne sponges to deoxygenation were observed. In fact, the long-term exposure of T. wilhelma to hypoxia was negatively correlated to KOG enrichments for both anoxia vs. hypoxia and normoxia vs. anoxia in Eurypon sp. 2 (Figure 1I, J), meaning that KOG expression was generally opposite between these two species for these conditions in terms of up and down regulation. The physiological responses of T. wilhelma and Eurypon sp. 2 to deoxygenation are presumably very different. These results demonstrate that large-scale transcriptomic responses to deoxygenation depend on both the studied species and oxygen concentration.

Based on KOG correlation analyses, the hypothesis that sponges are dormant, i.e. temporarily ceasing growth, development, and physical activity, under anoxia was not supported (Figure 4). Sponge gene expression under anoxia did not correlate to either diapause or a dauer state in C. elegans [START_REF] Sinha | Divergent Gene Expression in the Conserved Dauer Stage of the Nematodes Pristionchus Pacificus and Caenorhabditis Elegans[END_REF] or S. mosellana [START_REF] Gong | Global Transcriptome Analysis of Orange Wheat Blossom Midge, Sitodiplosis Mosellana (Gehin) (Diptera: Cecidomyiidae) to Identify Candidate Transcripts Regulating Diapause[END_REF], respectively. Although these dormancy datasets come from distantly related invertebrates, [START_REF] Strader | Red Fluorescence in Coral Larvae Is Associated with a Diapause-like State[END_REF] demonstrated a correlation between the KOG expression patterns of a specific phenotype of coral larvae and diapause in S. mosellana, hypothesizing that this potential functional analog of dormancy could serve to increase larval dispersal. Therefore, the sponges of Lough Hyne could only be considered dormant under anoxia if their transcriptional method of dormancy was distinct from that of an insect, nematode and coral, which is possible but not parsimonious given that animals are monophyletic [START_REF] Feuda | Improved Modeling of Compositional Heterogeneity Supports Sponges as Sister to All Other Animals[END_REF]). Nevertheless, it was inconclusive whether or not sponges are pumping water or feeding during anoxia, but sponges were definitively pumping under hypoxia [START_REF] Schuster | Effects of Seasonal Anoxia on the Microbial Community Structure in Demosponges in a Marine Lake in Lough Hyne, Ireland[END_REF]. Sponge pumping and filtering activity remain to be quantified under anoxia and hypoxia in situ.

Paradoxically, Eurypon sp. 2 responses to hypoxia were positively correlated to diapause in S. mosellana (Figure 1F,5). Therefore, the hypoxic response in Eurypon sp. 2 species, rather than the anoxic response, was similar to diapause in terms of global KOG expression patterns. The RRR kog was also downregulated in Eurypon sp. 2 under hypoxia as well as in S. mosellana under diapause, which may partially explain the correlation between these datasets. Despite the correlation, the EPC expression patterns between the two species were opposing. EPC was significantly downregulated in hypoxia in Eurypon sp. 2 and upregulated in diapause in S. Mosselana. These KOG enrichments are calculated from thousands of genes, so any individual gene within a KOG class might still be upregulated even if the KOG class is significantly downregulated. Therefore, it is possible that ATP production increased in hypoxia based on the upregulation of ATP synthase genes, in spite of the downregulation of the EPC KOG. Diapause can be associated with the reduction of cell division (MacRae 2010), which may occur in sponges under hypoxia but not anoxia based on MCM gene expression. In comparison, KOG enrichments under anoxia were not positively correlated to any external datasets for either sponge species (Figure 1), indicating that sponge transcriptomic responses to anoxia are relatively novel compared to those of other taxa in response to various factors in the existing literature [START_REF] Dixon | Coral Reefs. Genomic Determinants of Coral Heat Tolerance across Latitudes[END_REF][START_REF] Strader | Red Fluorescence in Coral Larvae Is Associated with a Diapause-like State[END_REF].

Sponge mitochondrial gene expression during deoxygenation

A high degree of transcriptional stability under deoxygenation was observed in sponge mitochondria. Only one of the 14 mitochondrial protein-coding genes, nad2, was significantly differentially expressed in the different oxygen conditions, where it was upregulated in anoxia vs. hypoxia in Eurypon sp 2 (Figure 4). Both of the sponge species had higher expression of nad2 in anoxia than in normoxia or hypoxia, but this difference was only significant in Eurypon sp. 2. Surprisingly little is known about differential gene expression in mitochondria, even within humans. In the limited studies available, mitochondrial gene expression was significantly different between human cancer cells and normal cells [START_REF] Reznik | Mitochondrial respiratory gene expression is suppressed in many cancers[END_REF] as well as between different populations of healthy pancreas tissue [START_REF] Medini | Mitochondrial gene expression in single cells shape pancreatic beta cells´subpopulations and explain variation in insulin pathway[END_REF]. In these studies, nad2 was significantly downregulated in most studied cancers, though it was still significantly upregulated in a few cancer types, and nad2 was also significantly differentially expressed between pancreas cell populations. Based on the general upregulation of nuclear genes associated with respiration in hypoxia, mitochondrial gene upregulation was expected in hypoxia when compared to normoxia and anoxia, but the opposite was true for nad2 in Eurypon sp. 2 in the current study. The other six nad genes, which code for parts of the same membrane respiratory chain complex (complex 1, NADH dehydrogenase), however, were not significantly differentially expressed in either species. In both human studies, when one mitochondrial gene was significantly differentially expressed, most other mitochondrial genes were significantly differentially expressed in the same direction [START_REF] Reznik | Mitochondrial respiratory gene expression is suppressed in many cancers[END_REF][START_REF] Medini | Mitochondrial gene expression in single cells shape pancreatic beta cells´subpopulations and explain variation in insulin pathway[END_REF]. The physiological meaning (if any) of this single, significantly differentially expressed mitochondrial gene in sponges remains to be determined. Compared to that of humans, even less is known about mitochondrial gene expression in marine invertebrates. Mitochondrial gene expression is usually excluded from differential expression analyses due to the differences in library preparation potentially needed for sequencing mitochondrial RNA and/or the exclusion of mitochondrial RNA from reference transcriptomes (sensu Meyer 2016). This dearth is partly due to the nature of mitochondrial transcription, which generates a single continuous transcript. Following transcription, however, mitochondrial transcripts are broken up into individual coding mRNA, tRNA or rRNA (D'Souza and Minczuk 2018). After the breakup, the levels of mRNAs from specific genes are controlled post-transcriptionally (D'Souza and Minczuk 2018). It is therefore possible to have different counts between the different mitochondrial genes despite their single transcript origin (Supplemental Figure 6). The consistent 'expression' of transcripts across oxygen conditions for each mitochondrial protein-coding gene following continuous transcription are likely to be necessary for the maintenance of mitochondrial function. It remains unclear how sponge mitochondria maintain this transcription from an energetic standpoint under anoxia. Within the mitochondrial electron transport chain, oxygen is used as the terminal electron acceptor and is crucial for the generation of ATP through oxidative phosphorylation. Theoretically, anoxia should cause electron transport and oxidative phosphorylation in the mitochondria to cease [START_REF] Chandel | Cellular Oxygen Sensing by Mitochondria: Old Questions, New Insight[END_REF]. If the mitochondria are transcribing solely to maintain oxidative phosphorylation machinery until oxygen returns, then they would need to use at least enough energy (~2 ATP per base) to transcribe for months at a time without oxygen [START_REF] Lynch | The Bioenergetic Costs of a Gene[END_REF], and this energetic requirement would apply to the transcription of sponge nuclear genes as well. Glycolysis may provide some anaerobic ATP production, even if it is not upregulated in anoxia or hypoxia. Alternatively, to maintain ATP production through aerobic respiration, the sponges of Lough Hyne might have access to oxygen produced by their Thaumarchaeota symbionts, a topic explored in more detail below.

Microbiome gene expression during deoxygenation

Thaumarchaeota gene expression

Thaumarchaeota symbionts from both sponge species exhibited significant differential gene expression under hypoxia. Only the Thaumarchaeota from Eurypon sp. 2 (OTU1) exhibited differential gene expression under anoxia, and only those of H. stellifera (OTU3) differentially expressed genes under anoxia vs. hypoxia. Though there were some species-specific responses, overall gene expression in both Thaumarchaeota was similar, supporting the hypothesis from [START_REF] Schuster | Effects of Seasonal Anoxia on the Microbial Community Structure in Demosponges in a Marine Lake in Lough Hyne, Ireland[END_REF] that these taxonomically similar symbionts between two sponge species might have similar functional roles.

Relatively few genes involved in anaerobic metabolism were expressed in Thaumarchaeota symbionts under any oxygen condition. Lactate/malate dehydrogenase, which could be involved in fermentation [START_REF] Bayer | Microbial Strategies for Survival in the Glass Sponge Vazella Pourtalesii[END_REF], was upregulated in Eurypon sp. 2 Thaumarchaeota, but only in hypoxia. Genes for the other anaerobic pathways suggested by Bayer et al. 2020, i.e. succinate dehydrogenase/ fumarate reductase, nitrite reductase and sulfite reductase, were present in both Thaumarchaeota genomes in the present study but were not significantly differentially expressed under different oxygen conditions. The genes for nitrite reductase exhibited low expression across all samples (0-20 counts) in both Thaumarcheaota, and so did the dehydrogenase/ fumarate reductase gene in the Thaumarcheaota of H. stellifera. Therefore, these genes with low expression may not be used under the conditions of this study. Only three genes in the Thaumarchaeota of Eurypon sp. 2 under anoxia were significantly differentially expressed, and none of these were involved in anaerobic metabolism. No genes were differentially expressed in H. stellifera Thaumarchaeota under anoxia. Therefore, there was no transcriptomic evidence for the upregulation of anaerobic metabolism within Thaumarchaeota in anoxia.

Transcription of ammonium monooxygenase (AMO) genes did not significantly change under anoxia or hypoxia, even though the ammonia oxidation pathway is oxygen dependent. Like the sponge mitochondria, the Thaumarchaeota are therefore either transcribing AMO genes to be ready when oxygen returns, or they are obtaining trace oxygen for this pathway by synthesizing oxygen themselves [START_REF] Kraft | Oxygen and Nitrogen Production by an Ammonia-Oxidizing Archaeon[END_REF]. If sponge Thaumarchaeota symbionts are also able to produce oxygen under anoxic conditions like their relative N. maritimus, then they should express amo genes under anoxia as well as normoxia. The proposed pathway for oxygen production in N. maritimus involves the reduction of nitrite to nitric oxide that is dismutated into oxygen and nitrous oxide [START_REF] Kraft | Oxygen and Nitrogen Production by an Ammonia-Oxidizing Archaeon[END_REF]. Water sampled for nitrite concentrations in Lough Hyne using established methods [START_REF] Grasshoff | Methods of Seawater Analysis[END_REF]) revealed that nitrite concentrations under anoxia ranged from 0.2-1.2 μM (Supplemental Table 3), increasing with depth below the thermocline. Therefore, the conditions in Lough Hyne would theoretically allow for Thaumarchaeota to produce oxygen under anoxia.

AsnC family proteins and nitroreductase expression was upregulated during oxygen production in N. maritimus [START_REF] Kraft | Oxygen Production by an Ammonia-Oxidizing Archaeon[END_REF]. In Archaea, AsnC proteins can regulate metabolic processes, including the switch between aerobic and anaerobic metabolisms [START_REF] Lemmens | Transcription Regulators in Archaea: Homologies and Differences with Bacterial Regulators[END_REF]). Similarly, nitroreductases can have a wide range of metabolic functions [START_REF] Akiva | Evolutionary and Molecular Foundations of Multiple Contemporary Functions of the Nitroreductase Superfamily[END_REF].

29

Thaumarchaeota within H. stellifera upregulated a Lrp/AsnC ligand binding domain under hypoxia and in anoxia vs. hypoxia. Another AsnC type protein was also upregulated under hypoxia, and a nitroreductase was upregulated under anoxia vs. hypoxia. However, these proteins did not directly correspond to those upregulated in N. maritimus during oxygen production (Supplemental Figure 9), and the latter proteins are not definitively known to be involved in oxygen production, since the pathway remains unknown.

In hypoxia, both Thaumarchaeota upregulated the stress response genes USpA and a protein from the HSP20 family (Figure 4). It was surprising that Thaumarchaeota may be stressed by hypoxia since free-living thaumarchaeal Nitrosopumilales can be found in high abundances in oxygen-minimum zones (OMZs) where oxygen is either at low concentrations or below detection [START_REF] Löscher | Production of Oceanic Nitrous Oxide by Ammonia-Oxidizing Archaea[END_REF], though some Thaumarchaeota OTUs decrease in abundance in suboxic and anoxic zones [START_REF] Sollai | A combined lipidomic and 16S rRNA gene amplicon sequencing approach reveals archaeal sources of intact polar lipids in the stratified Black Sea water column[END_REF]. Members of the genus Nitrosopumilus, however, increase their absolute abundance in suboxic regions, but then decrease in abundance with depth in the anoxic zone [START_REF] Sollai | A combined lipidomic and 16S rRNA gene amplicon sequencing approach reveals archaeal sources of intact polar lipids in the stratified Black Sea water column[END_REF], indicating a preference for suboxic environments. If suboxia is also the preferred niche of sponge Thaumarchaeota symbionts, the upregulation of stress response genes in response to hypoxia may correspond to a rise in intracellular oxidative stress following an increase in metabolic activity in hypoxia, like their sponge hosts.

Gammaproteobacteria gene expression

Gene expression in the Gammaproteobacteria was species and oxygen-condition dependent. In Eurypon sp. 2, Gammaproteobacteria downregulated genes in hypoxia that were associated with respiration and energy production, including cytochrome c and associated genes [START_REF] Welchen | Cytochrome c, a hub linking energy, redox, stress and signaling pathways in mitochondria and other cell compartments[END_REF]. The Gammaproteobacteria also downregulated key metabolic genes in hypoxia such as glutamine synthetase and structural genes like YbiS YcfS YnhG family protein, which is involved in cell wall production [START_REF] Magnet | Identification of the L,D-transpeptidases for peptidoglycan cross-linking in Escherichia coli[END_REF]. So unlike sponges or Thaumarchaeota, Gammaproteobacteria in Eurypon sp. 2 (OTU2) may decrease metabolic activity under hypoxia, though this is likely not the case in H. stellifera symbionts (OTU7).

Gammaproteobacteria in H. stellifera downregulated hundreds of genes in anoxia vs. hypoxia, including USpA, so, like other members of the holobiont, hypoxia requires a stronger transcriptomic response and, therefore, may be more stressful than anoxia for these Gammaproteobacterial symbionts. The few upregulated genes in anoxia vs. hypoxia included genes for energy production, e.g. a transmembrane proton pump, succinyl-CoA synthetase (part of the citric acid cycle), and formate dehydrogenase (alpha subunit), although other genes associated with energy production were downregulated. Formate dehydrogenase can function in both aerobic respiration and in anaerobic respiration [START_REF] Mika | Formate Dehydrogenasea Versatile Enzyme in Changing Environments[END_REF], but this enzyme was not significantly upregulated in anoxia compared to normoxia or in the Gammaproteobacteria of Eurypon sp. 2 under any condition. So, it is unlikely that this anaerobic pathway is active in sponge-associated Gammaproteobacteria under anoxia.

Gammaproteobacteria symbionts were previously hypothesized to oxidize sulfide under anoxia [START_REF] Schuster | Effects of Seasonal Anoxia on the Microbial Community Structure in Demosponges in a Marine Lake in Lough Hyne, Ireland[END_REF]. Although they are not present in OTU2, OTU7 from H. stellifera has the sulfur oxidation genes soxA, soxY, and soxZ, which could oxidize sulfide, as well as the nitrate reductase gene NAP. However, the sox genes were not expressed in any samples in the current study, and NAP was expressed at very low levels across all samples (2-28 counts). Moreover, neither symbiont exhibited significant differential gene expression in anoxia versus normoxia. Consequently, there was very little transcriptomic evidence that the anaerobic pathways in these Gammaproteobacteria were activated during environmental anoxia.

Summary of potential adaptations of sponge holobionts to low oxygen

Metabolically, sponge transcription under anoxia was surprisingly similar to normoxia throughout the holobiont members, and there was no evidence of dormancy within the holobionts in anoxia (Figure 4). No clear anaerobic respiration responses were observed in the transcriptomes of sponges, their mitochondria or their two microbial symbionts under environmental anoxia. The limited significant differential expression of glycolysis genes could indicate that this metabolic process continues to comparable degrees under normoxia, hypoxia, and anoxia. Therefore, in the sponges and their microbes, glycolysis could potentially provide ATP without oxygen (Figure 4). Hypoxia, rather than anoxia, elicited upregulation in energy production machinery such as ATP synthase (Figure 4). These potential metabolic increases suggest that within the hypoxic limits studied here (1.30-3.56 mg L -1 ), sponges perform optimally. Indeed, these limits are above those at which mortality or substantial physiological changes have been observed for sponges (Micaroni et al. 2021;[START_REF] Kumala | Spatial and temporal anoxia in single-osculum Halichondria panicea demosponge explants studied with planar optodes[END_REF][START_REF] Mills | Oxygen Requirements of the Earliest Animals[END_REF]. This increased energy production could be stressful, however, and the diversion of resources could limit cell division in hypoxia (Figure 4). Somewhere between the lower limit of hypoxia (1.37 mg L -1 ) and the full anoxia studied here, sponges may no longer attempt to increase energy efficiency, and therefore elevated transcription of ATP synthase genes was not observed in anoxia.

For both the sponges and their Thaumarcheaota symbionts, cellular stress responses were evident under deoxygenation. Both sponge species upregulated HSP genes under anoxia and hypoxia. These genes may be crucial for sponge deoxygenation tolerance in these species (Figure 4). We hypothesize that Hsp90 plays a role in sponge oxygen homeostasis, but this needs to be tested. The Thaumarchaeota symbionts upregulated HSP20 and other stress response genes in anoxia and hypoxia for H. stellifera and hypoxia within E. sp. 2 (Figure 4). In hypoxia, archaeal and host upregulation of stress response genes may have been a response to oxidative stress following increased metabolic activities. In anoxia, cellular stress response might result from the presence of sulfide. Additionally, post-transcriptional, metabolic changes may play a role in deoxygenation tolerance, and proteomic and metabolomic changes should be examined in future studies.

Given the largely unchanged metabolic gene expression patterns under hypoxia or anoxia in these sponges, the production of oxygen by Thaumarchaeota could also be crucial to holobiont survival (noted with a question mark in Figure 4). Although any oxygen produced would be at extremely low concentrations, it would have a short path of diffusion to the sponge cells, as Thaumarchaeota symbionts maintain high cell densities within sponges [START_REF] Moeller | Characterization of a Thaumarchaeal Symbiont That Drives Incomplete Nitrification in the Tropical Sponge Ianthella basta[END_REF]). In addition, mitochondria in general can utilize very low concentrations of oxygen (500 nM) without changing their oxygen consumption rate, and cells can still respire below that threshold albeit at a lower rate [START_REF] Chance | Reaction of Oxygen with the Respiratory Chain in Cells and Tissues[END_REF][START_REF] Longmuir | Respiration Rate of Rat-Liver Cells at Low Oxygen Concentrations[END_REF][START_REF] Wilson | The Oxygen Dependence of Cellular Energy Metabolism[END_REF]. In culture, oxygen generated by Thaumarchaeota reached 100-200 nM [START_REF] Kraft | Oxygen and Nitrogen Production by an Ammonia-Oxidizing Archaeon[END_REF]. If sponge Thaumarchaeota produce oxygen at similar levels under anoxia, cellular respiration could continue in sponges under anoxia, though likely at lower respiration rates when compared to normoxia. However, given that the genes responsible are still unknown, the production of oxygen by sponge Thaumarchaeota needs to be verified experimentally.

Relevance to future ocean deoxygenation

Sponges are predicted to be the 'winners' over other benthic organisms in the future due to their tolerance of ocean warming and acidification associated with climate change [START_REF] Bell | Sponges to Be Winners under Near-Future Climate Scenarios[END_REF]. Indeed, sponges have already increased in abundance in some locations while corals have decreased due to altered carbon and nutrient cycling [START_REF] Pawlik | A Vicious Circle? Altered Carbon and Nutrient Cycling May Explain the Low Resilience of Caribbean Coral Reefs[END_REF]. The additional stressor of deoxygenation that comes with climate change and coastal eutrophication has been relatively overlooked compared to ocean warming and acidification. In other marine taxa, a meta-analysis demonstrated that the incorporation of deoxygenation as a factor along with warming could substantially decrease survival times of marine organisms compared to warming alone (Vaquer-Sunyer and Duarte 2011). Since sponges may be more tolerant than other taxa to changes in all three stressors of temperature, pH and oxygen, they could indeed become more abundant in a changing ocean. However, even sponges can be lethally affected by warming, acidification [START_REF] Bennett | Interactive Effects of Temperature and pCO2 on Sponges: From the Cradle to the Grave[END_REF]) and deoxygenation (Micaroni et al. 2021) in scenarios that might occur with climate change, since tolerance to these stressors is speciesspecific. So, surviving sponge species in a warmed, acidified, and deoxygenated ocean might be abundant, but their abundance, and diversity would likely be lower than at present. The effects of all three of these stressors together need to be understood on a much wider range of organisms to accurately predict the composition of future oceans. More importantly, the causes and effects of anthropogenic climate change need to be ameliorated as soon as possible to avoid such diversity loss.

Relevance for animal evolution

Assuming that the first metazoans were more like an adult sponge than a ctenophore, based on morphology [START_REF] Nielsen | Six Major Steps in Animal Evolution: Are We Derived Sponge Larvae?[END_REF]) and phylogeny [START_REF] Feuda | Improved Modeling of Compositional Heterogeneity Supports Sponges as Sister to All Other Animals[END_REF], then some species of early sponges or proto-animals could have tolerated low environmental oxygen levels and at least some periods of benthic anoxia and euxinia in the Neoproterozoic Era (~538.8-1,000 million years ago) when early animals, including sponges, first evolved [START_REF] Dohrmann | Dating Early Animal Evolution Using Phylogenomic Data[END_REF]. Lough Hyne sponges are a valuable model system for understanding early animal evolution. Early sponges could have survived under low and variable oxygen levels using mechanisms like those of Lough Hyne sponges. Therefore, contrary to longstanding hypotheses, the evolution of animals probably did not require high and stable levels of environmental oxygen. also thank the Nordcee lab technicians of the Department of Biology at University of Southern Denmark for constant support in the laboratory.

Figures Dessinées, d'Après Nature". Paris: Chez Deterville. 2: 1-296, pls. 6-18., available online at

Figure 1 .

 1 Figure1. KOG expression profiles of sponges under varied oxygen concentrations compared to each other and external datasets. A. Clustered heatmap (based on correlation) of KOG class (rows) enrichments for upregulated or downregulated genes. Highlighted columns correspond to data generated in this study (Purple: Hymeraphia stellifera [Hs], Cyan: Eurypon sp. 2 [Es]). Unhighlighted columns represent external datasets. Up or down regulation with respect to normoxia is indicated in 'Hypox' (hypoxia) and 'Anox' (anoxia) columns, and differential regulation with respect to hypoxia is shown in 'Hypox vs. Anox' columns, i.e. upregulated KOGs (red delta ranks) are upregulated in anoxia with respect to hypoxia. Within each column on the heatmap, KOGs that are significantly enriched (FDRadjusted p < .05) are marked with asterisks. KOG enrichments from the sponge Tethya wilhelma after four days exposure to hypoxia (Tw Hypox Long) and one hour exposure to anoxia (Tw Anox Short), the coral Acropora millepora after three days of exposure to heat stress (Am Coral Heat), Caenorhabditis elegans in a dauer state (Ce Dauer), and the midge Sitodiplosis mosellana in diapause (Sm Diapause) from external studies were included for comparative purposes. B-E. Significant (p < 0.05) correlations between the KOG enrichments for different species and oxygen levels from this study. F-J. Significant correlations between KOG enrichments in this study and those of external datasets. All other possible comparisons between the various datasets were not significant. Within Eurypon sp. 2, KOG class expression patterns in hypoxia vs. anoxia were significantly positively correlated to those of hypoxia (r = 0.71, p < 0.01, Figure1B). Similarly, in H. stellifera, KOG enrichment in hypoxia vs. anoxia was positively correlated to that of hypoxia (r = 0.72, p < 0.001, Figure1C) and anoxia (r = 0.46, p < 0.05, Figure1D). All significant correlations between datasets are shown in (Figure1B-J), and all relationships that were not significantly correlated are not shown. At the level of individual KOGs, the KOG class 'replication recombination and repair' (RRR) was significantly enriched in all pairwise comparisons of oxygen conditions for both sponge species (Padj < 0.05), though in different directions depending on oxygen. It was significantly upregulated in anoxia and anoxia vs. hypoxia but downregulated in hypoxia in both species (Figure1A). In Eurypon sp. 2, five genes with RRR KOG annotations were significantly upregulated in hypoxia, e.g. HELZ and HELZ2 transcription associated helicases, and 11 genes were significantly downregulated (p < 0.1, Supplemental Figure4A), including five genes from the minichromosome maintenance (MCM) protein complex. In anoxia, two genes were significantly upregulated within RRR in this species -'siz regulation of ARF protein signal transduction' and DNA polymerase epsilon catalytic subunit (p < 0.1, Supplemental Figure4B). The majority of significantly differentially expressed RRR genes (five) were downregulated in anoxia (p < 0.1, Supplemental Figure4B). The majority of RRR genes (nine) were also downregulated in anoxia vs. hypoxia, with only two significantly upregulated in anoxia (Supplemental Figure4B). The genes NBN and RFC1 were downregulated in anoxia compared to both normoxia and hypoxia (Supplemental Figure4B,C). Similarly, the gene 'siz regulation of ARF protein signal transduction' was upregulated in anoxia compared to both normoxia and hypoxia (Supplemental Figure4B,C).

Figure 3 .

 3 Figure 3. All significantly differentially expressed genes for the Gammaproteobacteria in both sponge species (A: Eurypon sp 2, B: H. stellifera. For Gammaproteobacteria associated with Eurypon sp. 2, there were only significant differences between hypoxia vs. normoxia (A). In symbionts associated with H. stellifera, there were only significantly differentially expressed genes in anoxia vs. hypoxia (B,C). All downregulated genes in anoxia vs. hypoxia are shown in B. There were 402 genes that were downregulated in anoxia vs. hypoxia in Gammaproteobacteria from H. stellifera, so only the top 10 based on p value are shown (C).

Figure 4 .

 4 Figure 4. Summary of sponge holobiont responses to deoxygenation presented in a tabular format with species on the left (Es: Eurypon sp. 2, Hs: H. stellifera) and oxygen condition on the top (hypoxia, anoxia, and hypoxia vs. normoxia). The major components of the holobiont are illustrated as follows: sponge cell (tan), mitochondria (red), Gammaproteobacteria (green), and Thaumarchaeota (purple). Upregulated and downregulated genes and KOGs are indicated with red and blue arrows, respectively. Unchanged gene expression (GE) is noted by the = sign. Correlations (corr.) with external dormancy datasets are indicated along with the direction of the correlation (+ or -).

Acknowledgements

Field work at Lough Hyne was carried out with the permission of the National Parks and Wildlife Service of Ireland. We thank Luke Harman and Allen Whittaker of University College Cork. We

Funding

This project was funded by Villum Fonden grant no. 16518. B.K. was funded by the Villum Young Investigator Grant 25491.

Data, scripts, code, and supplementary information availability

Raw sequence data as well as genome and transcriptome assemblies are available at NCBI under BioProject number PRJNA893197 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA893197). Scripts, code, and additional data are available here: https://zenodo.org/badge/latestdoi/573040118. Supplemental figures and tables are available on BioRxiv under 'Supplementary material.'

Conflict of interest disclosure

The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.