Corals were collected under the Israel Nature and National Parks Protection Authority permit number 2018/41994. Coral fragments (2 cm long) from three Stylophora pistillata colonies (28 fragments per colony) were collected at a depth of 4–5 meters from the coral nursery of the Inter-University Institute (IUI) in the Gulf of Eilat Red Sea (28.6929°N, 34.7299°E) in February 2019 and immediately transported to a controlled environment aquarium system at the Leon H. Charny School of Marine Sciences in the University of Haifa, Israel. The fragments were glued onto pin-like bases (see Figure 1 in Zaquin et al., 2019) and allowed to recover for two weeks under simulated ambient conditions (light parameters, temperature, salinity, pH, calcium, magnesium, phosphate and nitrate), replicating the IUI coral nursery conditions that were tested biweekly and adjusted to standards as needed. Dissolved oxygen (DO) levels of ~5.5 mg L^-1, simulating the average normoxic conditions at the site of collection (Drake et al., 2021), were maintained via continuous aeration. Water changes of 10% were performed on a weekly basis using artificial seawater (Red Sea Salt, Red Sea Ltd.) and all corals were fed identical food twice a week (Microvore, Brightwell Aquatics, Fort Payne, AL).

After an acclimation period of two weeks in the same aquarium system, four identical 20 L aquaria were established, two aquaria for each treatment (deoxygenation and control). Fragments from the three colonies were placed in each aquarium, so that in each aquarium all three colonies (i.e., genotypes) were represented, and acclimated again for 5 days ( Figure S2A ). The aquaria were designed with airtight lids containing holes fitted for nitrogen bubbling, a water pump, a multi-parameter measuring device, and a water inlet/outlet. A water pump was used to ensure constant water circulation and to maintain homogenous DO levels. DO concentrations were measured every 15 minutes with a multi-parameter measuring device (WTW Multi 3420, Xylem analytics, Germany). Reduced DO concentrations (deoxygenation) in the experimental aquaria were generated by bubbling the water with nitrogen (Haas et al., 2014), (Johnson et al., 2021b) and kept at a level of 0.34 mg L^-1 (± 0.16 mg L^-1) that is close to the minimum level of 0.5-0.2 mg L^-1 observed in water-quality surveys (Altieri et al., 2017) (Johnson et al., 2021a) and in laboratory experiments ( Figure S3 ) (Altieri et al., 2017). By manipulating O₂ levels alone, pathways of acclimation to this specific stressor might be resolved (Johnson et al., 2021b). However, as deoxygenation is also accompanied with alterations in CO₂ and pH (Klein et al., 2017), our laboratory conditions did not address acclimatization of S. pistillata to natural deoxygenation.

For RNA extractions, three fragments (one from each colony) from each treatment (deoxygenation and control) were sampled at five time points: three time points at night, after 12 hours, one week, and two weeks of incubation in deoxygenation; and two time points at day, after one week and two weeks of incubation in deoxygenation ( Figure S2B ). In order to represent both the colony (i.e., genotype) and the aquarium at each sampling point, the sampling of colonies within treatment and within control was randomized but the following criteria were considered: the same colony combination (e.g., colonies 1 and 2) will not be sampled from the same aquarium at consecutive times, and at all time-point both aquaria per treatment were sampled. Therefore, all three colonies were represented at each sampling time point. Hence, controlling for colony effect was enabled ( Figure S2B, C ). The fragments were snap-frozen with liquid nitrogen, and stored at -80°C. For DNA extractions, the same fragments sampled for RNA after one week and two weeks were used (n=6 for control and deoxygenation, night and day sampling combined). Additionally, four fragments (three from different colonies and one randomly chosen) from each treatment (control and deoxygenation) were snap-frozen and used for DNA extractions (n=10 in total for DNA extractions). In order to be consistent with other studies, the time points for the physiological measurements were 1 week (Altieri et al., 2017) and 11 days (Johnson et al., 2021b), reporting prolonged deoxygenation, in addition to our 14 day deoxygenation experiment. For calcification rates fragments from the three different colonies were sampled at the same time of day from each of the control and deoxygenation treated aquaria (n=3 for 1 week, n=4 for 11 days and for 2 weeks; see Figure S2 ). The same fragments were used for protein extraction, symbiont density counts, and chlorophyll-a quantification (see below). The extensive time series design restricted the number of fragments used for each measurement, resulting in replicate numbers of n=3 for transcriptom analysis and n=3 or 4 for calcification rates, protein extraction, algal symbiont density, and chlorophyll-a quantification. However, this replicate number is in accordance with other similar works [e.g., (Johnson et al., 2021b), (Alderdice et al., 2020)] and is considered sufficient to acquire information.

To understand the effect of prolonged deoxygenation on coral growth, calcification rates were determined using the total alkalinity (TA) anomaly method with minor modifications, as described (Zaquin et al., 2019). In brief, each fragment was sealed in a specimen cup filled with 0.22 mm of filtered artificial seawater collected from the aquarium system holding the fragments. For deoxygenation, the filtered artificial seawater was bubbled with nitrogen (Haas et al., 2014; Johnson et al., 2021b) to achieve hypoxic O₂ levels of 0.34 mg L^-1 before sealing the cups. The cups were kept in their respective aquarium systems for 1.5 h. Water samples were collected for TA at the beginning of the incubation and at the end of the incubation after 1.5 h from each specific cup. Water samples were filtered to remove debris and measured with an automatic alkalinity titrator (855 Robotic Titrosampler, Metrohm, Switzerland). The calcification rates [µmol CaCO₃ cm^-2 h^-1] were then calculated using the equation proposed (Schneider and Erez, 2006):

where ΔTA [µmol g^-1] is the difference in alkalinity between the beginning and the end of the incubation period; V is the recorded volume of the artificial seawater in the respective

incubated cup [ml]; 1.028 is the density of seawater in the northern Gulf of Eilat [g ml^-1]; and SA is the surface area of the coral (cm²) estimated as described (Veal et al., 2010); T is the duration of the incubation [hours].

To quantify coral tissue loss, protein concentrations were calculated and normalized to skeletal surface area (cm²). To assess the symbiont-host-system, symbiont density was calculated and normalized to skeletal surface area (cm²). To assess symbiont physiology, chlorophyll-a concentration per symbiont cell was calculated.

After measuring calcification rates, the same fragments were frozen in liquid nitrogen and kept at -80°C. Coral tissue was removed by an airbrush connected to a reservoir of phosphate buffer saline (PBS) solution filtered through a 0.22 mm filter, and the skeletons were kept for further analysis. The extracted tissue was mechanically disrupted using an electrical homogenizer (HOG-160-1/2, MRC Labs, Israel) for 3 × 10 s at the maximum power setting of 10. The density of zooxanthellae in the homogenate was determined from three replicate counts of the homogenized seawater slurry using a haemocytometer (Improved Neubauer, BOECO Germany) under a microscope and normalized to coral surface area (Veal et al., 2010). The homogenate was centrifuged at 5000 g for 5 min at 4°C to separate the debris and the zooxanthellae cells from the coral host tissue. A protease inhibitor cocktail (cat. G652A, Promega, United States) was added to the supernatant with the host tissue and sonicated for 3 × 30 s (Ultrasonic cell crusher, MRC Labs, Israel). The protein concentration of the coral host was determined using a fluorometric BCA protein kit (Pierce BCA, United States) following the manufacturer protocol. A PerkinElmer (2300 EnSpireR, United States) plate reader was used to determine the total protein concentration with a 540 nm wavelength emission. To measure the chlorophyll-a (chl-a) concentrations of the zooxanthellae, 1 ml of 90% cold acetone was added to the pellet. The pellet was homogenized and incubated overnight at 4°C in the dark. The solution was filtered into a glass cuvette through a 0.22 mm syringe filter (Zaquin et al., 2019). A NanoDrop (Thermo-Fisher, United States) was used for spectrophotometric measurements at wavelengths of 630, 647, 664, and 691 nm, and the light absorbance results were used to calculate the chl-a concentration based on the following equation (Ritchie, 2008):

To assess the effect of deoxygenation on the photophysiology of symbionts, the photochemical efficiency of the symbionts in the coral host was measured. Photochemical efficiency of the algal symbiont was assessed with quantum yield of chlorophyll-a fluorescence measured with an underwater diving-PAM fluorometer (Walz, Effeltrich, Germany) as previously described (Johnson et al., 2021b) with slight modifications: saturation intensity=7 and measuring light intensity=6. Measurements of the maximum (Fv/Fm) quantum yield of photosystem II were taken at one week, 11 days and two weeks of incubation under deoxygenation. Coral fragments from each treatment were randomly chosen (representing all three colonies) and were measured after two hours of the dark period. All measurements were conducted in triplicates per each coral fragment in order to minimize variations.

Total DNA from frozen coral fragments was isolated using the Wizard^® SV genomic DNA purification system tissue protocol (Promega, WI, United States) according to the manufacturer’s instructions. DNA concentration and quality were evaluated using the NanoDrop 2000c (Thermo Scientific, United States).

Changes in the amount of mitochondrial DNA were determined relative to the DNA of a nuclear single-copy gene (that has two copies/diploid genome) (Nicklas et al., 2004). qRT–PCR primers were designed using the software Primer Express 3.0 (Applied Biosystems, CA, USA). Primers were planned for the S. Pistillata mitochondrial gene Cox1 (AB_441231) (F-CCTCTTGCTAGTATTCAAGCACAC;R-TTTATTAAAAGACACACCCGGGGCTCG) and the nuclear Timeless gene (XM_022932799) (F-CCAGAAGATGAAGAGATGGCAGCAG; R-TAATACTTCAAAAGCCAGCAGTATGG). Transcript levels were determined by qRT–PCR using the 7500 Fast Real-Time PCR system (Applied Biosystems). Thermal cycling was 95°C for 10 min, followed by 40 cycles of 95°C for 15 s, and 60°C for 100 s, then a final extension at 60°C for 1 min. Triplicate reactions were carried out in 20 μL volumes consisting of: 1× FastStart SYBR Green Master Mix (Applied Biosystems), 0.25 μM of each forward/reverse primer set, 20ng DNA, and PCR grade water. The difference between CT values of Cox1 and Timeless in each sample was used as a measure of relative abundance of the mitochondrial genome (Nicklas et al., 2004).

Statistical analysis of calcification rates, zooxanthellae cell numbers, chlorophyll-a and Cox1 transcription levels between the control and deoxygenation group at each sampling point was done in R. For each measurement, the Shapiro-Wilk test for normality and the Levent test of homogeneity of variance were performed. Where normality and homogeneity of variance were confirmed, Two-way ANOVA test were performed (calcification rates, zooxanthellae cell numbers and chlorophyll -a). For the Cox1 transcription measurements normal distribution could not be assumed (Shpiro-Wilk W=0.825, P<0.05). Therefore, the Aligned Rank Transform Analysis of Variance (ART ANOVA) was done with R package ARTool (Mansouri et al., 2004). This test allows for non-parametric, rank based factorial analysis of variance.

Total RNA was extracted from all 30 fragments using TRIzol reagent, Pure Link RNA kit and Pure Link DNase (Invitrogen, CA, United States) according to manufacturer’s instructions. Purified RNA samples were analyzed using a Nano Drop 1000c spectrophotometer (Thermo Scientific, United States) to assess RNA quantity and the 2200 Tape Station (Agilent, United States) to assess RNA quality (RIN > 8.5). RNA-Seq libraries were prepared using an in-house protocol at the Weizmann Institute of Science. RNA-Seq single end (SE) reads were adapter trimmed using cutadapt 1.5 (https://cutadapt.readthedocs.io), and quality filtered with Trimmomatic 0.30 (Bolger et al., 2014). In order to detect the species origin of symbiont sequences, high-quality reads were mapped to the selected proteomes databases of Symbiodinium/SAR and cnidarian/metazoa, using Diamond (Buchfink et al., 2015). Top hits almost exclusively belonged to robust corals (mainly S. pistillata), and Symbiodinium “clade A” species (mainly Symbiodinium microadriaticum). We further aligned the reads to the host genome assembly (NCBI GCA_002571385.1) and the symbiont genome assembly (NCBI GCA_001939145.1) using STAR 2.5 (Voolstra et al., 2017). For different samples, 82-86% of the reads were concordantly aligned to S. pistillata genome, and 0.7-3.8% to S. microadriaticum. For the host genes, DE analysis was conducted using the Generalized Linear Model (GLM) framework of DESeq2 (Love et al., 2014). Three factors were considered: Treatment (H, C), Diurnal (Day: zt6; Night: zt18), Temporal (12 hour, 1 week, 2 weeks), in the following two formulas in DESeq2: (a) “~ Treatment + Diurnal + Temporal + Treatment : Diurnal + Treatment : Temporal + colony”; (b) “~ x + colony”, where the x factor is a string combination of the factors: Treatment, Diurnal, Temporal (e.g., “Control_Night_day1”, “Deoxygenation_Night_day1”, etc…). Colony was added as a factor to all formula in order to control for the colony effect. The first formula models the additive effect of three factors plus interaction effects of treatment with the diurnal factor, and with the Temporal factor ( Data S1 and S2 ). Heatmap results are based on FPKM (Fragments per kb million), and multidimensional scaling ordination (NMDS) were conducted using Vegan tool in R. All significant results shown here have Benjamini and Hochberg (BH) adjusted p-values<0.05.

Proteome sequences of the tested metazoan species ( Data S3 ) were clustered using Orthofinder 2.4.0 (Emms and Kelly, 2015), to give ~16,000 orthogroups ( Data S3 ). Each orthogroup represents a cluster of protein sequences that share sequence similarity. That similarity indicates that these genes are descended from a single gene in the last common ancestor of a clade of species. For each orthogroup, sequences were aligned using MAFFT (Rozewicki et al., 2019), and poorly aligned sequences/sites were removed. Gene trees were constructed with iqtree (Nguyen et al., 2014), using the model option “LG+G4”, namely: an improved general amino acid replacement matrix (10.1093/molbev/msn067), plus discrete Gamma distribution of 4 rates across sites categories. Iqtree provides binary unrooted trees, and we further rooted the tree and estimated orthology relationships between proteins (one:one, one:many, many:many) using Orthofinder pipeline (Emms and Kelly, 2019). In order to estimate orthology in Orthofinder, the well-established Opisthokonta species tree was considered as reference. We also used reciprocal Blast between S. pistillata and homo sapiens, in order to find closely-related homologous genes that were not identified as orthologous genes. We also considered closely related cnidarians vs. Deuterostomia genes belonging to relatively small Orthogroups. Coral gene symbols were determined based on human 1:1 orthology, or otherwise we used the prefix of the gene family symbol. For selected gene trees (specifically PHD4 and HIFs family), we further tested that aligned S. pistillata sequences are not highly segmented. For that, the S. pistillata members of these gene families were aligned by multiple alignment using MAFFT (Rozewicki et al., 2019). The alignment was manually inspected. It is noted here that the proteome data used for S. pistillata for orthology analysis are directly related to the transcriptome data used for gene expression quantification.

Genes’ biological terms were assigned based on Uniprot S. pistillata (www.uniprot.org), KEGG S. pistillata (www.kegg.jp), S. pistillata Trinotate annotations (Bryant et al., 2017), and Uniprot S. microadriaticum databases. Enrichment analysis was conducted in Bioconductor GOSeq (Subramanian et al., 2005), as described previously (Malik et al., 2021).

We wanted to study the response of S. pistillata and its symbionts to severe prolonged deoxygenation conditions, not hitherto exerted under laboratory conditions yet relevant to acute deoxygenation events recorded on coral reefs. Therefore, in this experiment, coral fragments were subjected to over two weeks of deoxygenation where the oxygen level was 0.34 mg L^-1 (± 0.16 mg L^-1) ( Figure S3 ).

The morphology and physiology of the coral and its symbionts showed high tolerance to the prolonged deoxygenation. All the coral fragments survived two weeks of severe deoxygenation conditions, without apparent morphological changes ( Figure 1A ). After one week of deoxygenation, coral protein amount was unaffected by treatment (n=3 for 1 week; n=4 for 11 days and for 2 weeks p=0.063) ( Figure 1B ). Coral calcification rates showed a trend of reduction that was not significant (n=3 for 1 week; n=4 for 11 days and for 2 weeks, p =0.295), resuming control levels after two weeks ( Figure 1C ). Symbiont numbers showed a trend of increase after two weeks of deoxygenation that was not significant (n=3 for 1 week; n=4 for 11 days and for 2 weeks, p =0.277) ( Figure 1D ). The chlorophyll-a content per symbiont cell showed a trend of decrease after two weeks of deoxygenation that was not significant (n=3 for 1 week; n=4 for 11 days and for 2 weeks p=0.174) ( Figure 1E ). The photochemical efficiency (Fv/Fm) of the symbionts decreased by 4% after two weeks of deoxygenation. Although this decrease was statistically significant (n=12 for 1 week and 11 days; n=9 for 2 weeks p= 0.0019) ( Figure 1F ), we cannot assume it has biological significance as decreases of 40-50% after three to five days are reported for deoxygenation stress (Johnson et al., 2021b). These results indicate that S. pistillata, the symbiont and the host-symbiont system are tolerant to a prolonged reduction in oxygen levels.

To point to molecular mechanisms that may enable the tolerance of S. pistillata to the deoxygenation conditions applied here, we measured the corals’ diurnal transcriptional profiles (midnight, zeitgeber time 18; and midday, zeitgeber time 6) at different time points (12 h, 1 week, and 2 weeks) during the prolonged deoxygenation ( Figure S2B ). Gene expression was quantified based on mapping sequence reads to the S. pistillata genome assembly (see Methods) and identified 24,577 expressed putative genes (hereafter termed “genes”).

To identify the general trends of gene expression in control and deoxygenation conditions through time, we performed a non-metric multidimensional scaling (NMDS, Figure 2A ). The NMDS clearly separates the samples of night deoxygenation from the rest of the samples ( Figure 2A ). The most distinct samples are the night-time samples from 12h and one week of deoxygenation. After two weeks of deoxygenation, two of night-time deoxygenation samples are close to the day-time deoxygenation samples, which are relatively close to the controls. These results point towards a dynamic acclimation to the deoxygenation conditions over the two weeks of the experiment. Genes that were differentially expressed (DE) under the deoxygenation treatment were identified by comparing deoxygenation to control (from here on this comparison is referred to as HvsC) at each sampling time point. Differences were considered significant for Benjamini-Hochberg adjusted P<0.05 (see Experimental procedures for details on annotations and statistics).

The number of genes showing significant differential expression (DE, p<0.05) in HvsC at each time point was higher at night compared to the number of DE genes during the day ( Figure 2B ). The number of DE genes at night decreases after two weeks ( Figure 2B ) in agreement with the trend shown by the NMDS results ( Figure 2A ). Thus, based on both NMDS and HvsC DE genes, the transcriptional response to deoxygenation was more extensive at night than during the day, and was gradually reduced with time.

Then, we used functional enrichment analysis to illuminate pathways significantly affected by deoxygenation. To put our results in an evolutionary context, we investigated the expression of conserved hypoxia-related orthologous genes. We also looked for coral-specific genes that include acquisition of components from the symbionts and stinging cell formation.

We sought to better understand the effect of deoxygenation at night-time, when both coral and algae are respiring and oxygen levels are extremely low, as compared to day-time, when photosynthesis occurs and oxygen is produced (Kühl et al., 1995), (Linsmayer et al., 2020). To that end, we identified DE genes that their response to deoxygenation was significantly different between the night-time points and the two day-time points (HvsC x Day vs Night interaction model, Figure 3A ). We divided these genes into two sets, genes that their HvsC was up-regulated at night compared to the day (color coded in orange in Figure 3A ), and genes that their HvsC was down-regulated at night (color coded in purple in Figure 3A ). The set of genes that showed significant up-regulation of HvsC at night compared to the day, was larger (1101 genes) than the set that showed significant HvsC down-regulation at night compared to the day (442 genes) ( Figure 3A , Data S1 ). We then looked for statistically significant overrepresentation of specific functional groups within these two sets, and identified more functional groups within night-up-regulated genes (15 groups) rather than night-down-regulated ones (two groups) ( Table 1 , Data S2 ).

Within the night-up-regulated genes, enrichment was found in genes related to mechanisms of energy production by the mitochondria, that are part of the highly conserved hypoxia response in metazoans (e.g., GO terms of glycolysis/gluconeogenesis, fatty acid beta-oxidation and oxidative phosphorylation), suggesting enhanced mitochondria activity. The night-time up-regulation of glycolysis rather than gluconeogenesis is supported by the expression of major genes in this pathway ( Figure 3C (1)-(3)). Specifically, phosphofructokinase (PFK), the rate-limiting enzyme for glycolysis and a direct HIF target gene is significantly up-regulated in HvsC at all night-time points and at day-time after 1 week, and its expression is higher at night ( Figure 3C (1)). Up-regulation of oxidative phosphorylation genes, including subunits of the electron transport chain (ETC, represented in Table 1 by the GO term: mitochondrial respiratory chain complex III) is exemplified by the expression of Succinate dehydrogenase [ubiquinone] iron-sulfur subunit, mitochondrial (SDHB) gene. SDHB was significantly up-regulated at all night-time points ( Figure 3C (4)). Enhanced antioxidant activity is indicated by the enrichment in glutathione transferase activity GO term ( Table 1 ) and the expression of anti-oxidant catalase genes that was up- regulated at night-time (( Figure 3C , (5) and (6)). Also, DNA replication and repair were enriched with up-regulated genes, suggesting enhanced maintenance of cellular integrity ( Table 1 ).

The observed enrichment in genes related to mitochondria dynamics (i.e., mitochondrion organization), suggests possible changes in the number and composition of the mitochondria. These changes could include increasing mitochondria numbers (e.g., by replication of the mitochondrial genome as well as synthesis of new mitochondrial components) or conversely, reduction in mitochondria numbers (e.g., by disassembly of mitochondria or by apoptosis/mitophagy). The expression of central genes of apoptosis/mitophagy show that the anti-apoptotic gene styBcl-2 was up-regulated at night time ( Figure 3C , (7)) and the pro-apoptotic, executioner of apoptosis gene styCaspase-3 was down-regulated ( Figure 3C , (9)). Also, BNIP, the major mitophagy-inducing and HIF target gene was not DE ( Figure 3B , (8)). Accordingly, these pathways of reducing cells/tissues/mitochondria numbers were not induced. To better understand the effect of deoxygenation on mitochondria dynamics we searched for mitochondria-related functional enrichment in the HvsC comparisons for each time point. An increase in the number of mitochondria is suggested by the enrichment of up-regulated genes in the “mitochondria biogenesis” functional group at night-time at one and two weeks and at day-time at two weeks ( Figure S4 ). Additionally, after two weeks at day-time the functional group “mitophagy” was enriched with down-regulated genes ( Figure S4 ). Therefore, we studied the mitochondrial DNA copies to establish how mitochondria numbers were affected by deoxygenation. By using quantitative real-time PCR we confirmed a significant increase in mitochondria DNA copy numbers after one to two weeks of deoxygenation ( Figure 3D ), suggesting mitochondria numbers are increased.

Taken together these results suggest that at night-time, deoxygenation imposes great demand for energy production, possibly provided by increasing mitochondria numbers.

We wanted to investigate the molecular acclimation to prolonged deoxygenation as compared to the immediate acute response. To do that we identified DE genes that their response to deoxygenation (HvsC) was significantly different between the 1^st night and the 1^st week night of deoxygenation, and between the 1^st night and the 2^nd week-night of deoxygenation (temporal interactions, Figure 3B ). We divided these genes into two sets, genes that their HvsC was up-regulated over time in the 2^nd week or 1^st week compared to the 1^st night (color coded in orange in Figure 3B ), and genes that their HvsC was down-regulated over time (purple coded in Figure 3B ). We then looked for statistically significant overrepresentation of specific functional groups within these sets ( Table 2 , Data S2 ).

Genes that their HvsC was down-regulated over time were related to metabolism, e.g., glycolysis, valine, leucine and isoleucine degradation, arabinose metabolic process and phagosome. This suggests that the acute metabolic-related response to deoxygenation is followed by metabolic acclimation to prolonged deoxygenation. For example, the expression of PFK, PDK and LDHA subsided over time ( Figure 3B (1)-(3)). At night-time up-regulation of transcription at 1^st week and of the mismatch repair pathway at 2^nd week of deoxygenation as compared to the first 12 hours was observed. Overall, the molecular acclimation to deoxygenation is reflected in the decrease of the metabolic response and maintenance of DNA correction. Also, after two weeks gene expression patterns at Day and Night are more similar than at one week ( Figure 3C , 2W), as observed for the whole transcriptome analysis ( Figure 2A ).

We next wanted to test whether the conserved components of the hypoxia pathway of S. pistillata respond to deoxygenation in a similar temporal sequence as in other metazoans. This sequence includes down-regulation of PHDs, stabilization and up-regulation of HIF-α, HIF-α association with HIF-β and HIF complex formation, HIF activation of its target genes (Shaw, 2008; Lendahl et al., 2009; Fuhrmann and Brüne, 2017) ( Figure S1 ). Specifically, we focused on PHDs orthologs, since up regulation of HIF-α and activation of its target genes were identified in corals, but only one PHD2 was identified and its expression was up-regulated rather than down-regulated (Alderdice et al., 2021), (Alderdice et al., 2020). To do that, we identified the relevant orthologous genes in S. pistillata genome, established their evolutionary position and then followed their transcription in HvsC (see methods Data S5 ).

The PHD1-3 family gene tree indicates that Cnidaria have only two PHD1-3 orthologs ( Figure S5A ). One of them (PHD3 homolog) was up-regulated at all time points, and the change in its expression level was remarkably high ( Figure S5B ). In contrast, the PHD4 tree included a single human PHD4 gene, and many S. pistillata orthologs, comprising 16 groups in total ( Figure 4A ). This expansion is evident only in Cnidaria containing photosynthetic symbionts (i.e., S. pistillata, Acropora sp. and Nematostella sp.) but not in Medusozoa and Myxozoa. Deoxygenation induced an immediate significant down-regulation of most PHD4 orthologs ( Figure 4B ): out of 20 orthologs, 14 were down-regulated in HvsC, while only two were up-regulated. These results suggest that corals and Nematostealla had expanded their PHD4-like gene repertoire that may function as the main oxygen sensors, being down-regulated under deoxygenation conditions.

The HIF gene family tree is composed of six highly conserved groups, including HIF-α, NPAS, AHR, Clock, HIF-β and ARNTL- like genes (in accordance with (Mills et al., 2018)). S. pistillata genome includes orthologes of all of these key genes ( Figure 5A ). Orthologous genes of HIF-α, AHR, NPSA, HIF-β and ARNT were significantly up-regulated in deoxygenation, especially in the night time points ( Figure 5B ). Importantly, the up-regulation of HIF-α at 12 hours subsided at two weeks ( Figure S6 ). The up-regulation of some of the circadian rhythm genes suggests that these genes participate in the response to hypoxia in corals.

Next, we identified 30 genes that are orthologous to the conserved hypoxia target genes activated by HIF ( Figure S7A ). Most of these genes show significant up-regulation under deoxygenation conditions, throughout the two weeks ( Figure S7B ). Taken together these results suggest that the conserved sequential core hypoxia signature, from PHD down-regulation, through HIF up-regulation and the up-regulation of HIF target genes, is conserved and functional in corals, possibly enabling survival under deoxygenation.

Increasing mitochondria numbers and biogenesis requires cell membrane components such as sterols while increasing glycolysis requires glycogen and starch sources (Kopp et al., 2015; Hambleton et al., 2019). Our studies reveal an immediate up-regulation of most genes that enhance acquisition of sterols (NPCs) and glycogen (Ppp1r3b and Gbs-76A) from the symbiont at all time points ( Figure 6A “symbiosis-derived components”, and Figure 6B (1)-(6)). The enhanced acquisition of carbohydrates is also indicated by the enrichment of up-regulated genes of “starch binding” GO term at all night-time points and at one-week day-time ( Figure 6A ). Conversely, after two weeks at night, “nematocyst” GO term was enriched with down-regulated genes ( Figure 6A ). Indeed, a significant down-regulation of the expression of nematogalactin and minicollagen-1, essential for nematocyst formation (Beckmann and Özbek, 2012) is observed at one week, and becomes more prominent over time ( Figure 6B (7), (8)), supporting down-regulation of nematocyst formation. Nematocysts are essential for predation, defense and for territorial aggression (Yosef et al., 2020), therefore these results could indicate that prolonged deoxygenation enhance the corals’ dependence on symbiosis-derived components and reduce its predation efforts.