Elucidating gene expression adaptation of phylogenetically divergent coral holobionts under heat stress

As coral reefs struggle to survive under climate change, it is crucial to know whether they have the capacity to withstand changing conditions, particularly increasing seawater temperatures. Thermal tolerance requires the integrative response of the different components of the coral holobiont (coral host, algal photosymbiont, and associated microbiome). Here, using a controlled thermal stress experiment across three divergent Caribbean coral species, we attempt to dissect holobiont member metatranscriptome responses from coral taxa with different sensitivities to heat stress and use phylogenetic ANOVA to study the evolution of gene expression adaptation. We show that coral response to heat stress is a complex trait derived from multiple interactions among holobiont members. We identify host and photosymbiont genes that exhibit lineage-specific expression level adaptation and uncover potential roles for bacterial associates in supplementing the metabolic needs of the coral-photosymbiont duo during heat stress. Our results stress the importance of integrative and comparative approaches across a wide range of species to better understand coral survival under the predicted rise in sea surface temperatures.


Symbiodiniaceae physiology reveals a distinctive response to heat stress
After 9 days of heat exposure O. faveolata-Symbiodinium A3 and P. clivosa-B. faviinorum showed a drastic reduction in symbiont density (t-test; p=0.0005866; p=0.01263 respectively), in comparison with control organisms maintained at 28ºC. In contrast, S. radians-Breviolum B5 exhibited a non-significant reduction in symbiont density after the heat stress exposure (Fig.  1D, Supplementary Table 3,4).
A decrease of PSII maximum quantum yield efficiency (Fv/Fm) during heat stress (34ºC) indicates significant increases in light stress for Symbiodiniaceae in hospite during the elevated temperature treatment, which may or may not result in photosynthesis impairment 1,2 , but certainly increases the costs of repair and maintenance of photosynthetic activity. Fv/Fm revealed a statistically significant decrease in Symbiodinium A3 after 2 days of heat stress, whereas Fv/Fm in B. faviinorum and Breviolum B5 decreased significantly after 3 days ( Fig  1C). After 9 days, a significant decline in Fv/Fm on the heat-treated samples showed major differences between coral species (ANOVA, p=4.59e-06, Supplementary Table 9, Fig. 1C). By the end of the experiment, maximum quantum yield in B. faviinorum dropped abruptly, in Symbiodinium A3 it decreased to a lesser degree, and in Breviolum B5 it only diminished slightly (Fig. 1C).
The large values of a*sym determined for Symbiodinium A3 in O. faveolata, respond to the low number of Symbiodiniaceae cells that this species had at the time this experiment was performed (November 2008). This low symbiont content does not correspond with the large pigmentation and symbiont species that O. faveolata tends to associate with in this area of the Caribbean 2 and elsewhere 3 . The reef lagoon of Puerto Morelos experienced two consecutive events of elevated temperatures during the previous summer, and O. faveolata was found particularly paled during the summer of 2007 4 . Therefore, the process of recovery after both bleaching events could explain the low number of photosymbionts that O. faveolata harbored in our experimental study in comparison with other analysis performed in the same area 2 .
Based on ITS2 rRNA Symbiodiniaceae sequencing, we found that after 9 days the photosymbiont composition did not change between control and heat-treated colonies, P. clivosa hosted Breviolum faviinorum and Breviolum B5 were associated with S. radians whereas O. faveolata associated with Symbiodinium A3.
A set of genes are upregulated in P. clivosa when compared to the ortholog profile expression from the other species. For example, during the temperature challenge, the over expression of the TNFR 8 and 10 is less prominent in O. faveolata and higher in P. clivosa (Supplementary Data 1). When considering the Faviina expression profiles as reference to assess the ortholog expression in S. radians, we found that a vast majority of genes is not even differentially expressed, followed by a set of DEGs that aren't statistically significant and in some cases with an opposite transcriptional profile ( Supplementary Fig.2 A-B). A gene encoding for an allene oxide synthase-lipoxygenase is significantly overexpressed in P. clivosa when compared to S. radians, where it is under expressed (Supplementary Fig.2C). This protein plays a role in jasmonic acid biosynthesis in plants and has been studied in octocorals due to its catalase activity 5 . S. radians significant DEG profiles in general are not conserved with the Faviina orthologues ( Supplementary Fig.2C).
As not all the DEGs in each species were part of the KEGG enrichment analyses, we identified a subset of those DEGs previously described in coral transcriptomic studies of heat stress and bleaching. These DEGs involved in heat shock response, oxidative stress, immunity, ion transport, cytoskeleton and extracellular matrix, cell adhesion, biomineralization, protein folding and processing, apoptosis, and with miscellaneous activity 6 were found in all three coral hosts. The expression profiles of these DEGs are not consistent with a conserved heat stress response across scleractinian corals. In addition to previously described bleaching gene markers (reviewed in Louis et al., 2017 66 ; Maor-Landaw and Levy, 2016 6 ), we identified to the best of our knowledge a novel subset of DEGs present in all coral hosts that to our knowledge has not been implicated before as a part of heat stress or bleaching ( Supplementary Fig.3D). These shared DEGs did not elicit the same transcriptional profiles across taxa. For example, NFX1 type zinc finger contains protein 1 which is overexpressed in O. faveolata and displays an opposite trend in S. radians. This protein has been recently implicated in C. elegans as a factor required for epigenetic inheritance to guarantee an equilibrated amplification of small RNA signals and it has been hypothesized that it may establish epigenetic programs that allow a rapid adaptation to environmental perturbations in the wild 7 . Further studies on different coral populations and species to understand the molecular mechanisms of epigenetic inheritance and their role of this protein during bleaching are warranted.  11 identified this transporter homolog sequences in Symbiodinium. Our orthology analyses also found this transporter in Breviolum B5 and in the Symbiodinium A3 ( Fig. 3B-C and Supplementary Data 1). Symbiodinaceae in hospite relies heavily on host ammonium uptake and translocated as the main nitrogen source instead of nitrate 12 . The transport of nitrate via NRT2 seems unlikely since a structural analysis predicted its ligand as monoacylglycerol, and since the 3D structure resembles a bacterial glycerol-3-phosphate transporter 11 . Glycerol amounts and its release during osmotic stress in Symbiodinium in culture and during heat stress in hospite (A. aspera 6 and 8 d) have been documented and suggested as a metabolite conferring thermotolerance 13,14 . If monoacylglycerol (MAG) is transported by NRT2, during elevated temperature this transporter is under-expressed to increase MAG accumulation. Evidence from our evolutionary analyses indicates a divergent expression between lineages of the monoacylglycerol lipase (MAGL), an enzyme that cleaves this molecule in order to produce glycerol. Increased expression of MAGL under heat stress in the two Breviolum species when compared with the Symbiodinium A3 (Fig. 5) indicates an increase of glycerol production which may be also involved in lipid signaling during increasing temperature.

Symbiodiniaceae expression profiles reveal differential regulation in key metabolic and cellular pathways
Symbiodinium A3 exhibits a general stress response of downregulation of metabolic and cellular processes, such as ion transporters (zupT), exocytosis/vesicular transport (VPS35) (Fig.  3B). VPS (vacuolar protein sorting proteins) and syntaxins have been hypothesized to be involved in symbiont shuffling and symbiont exit from the host during coral bleaching 15 . Interestingly, a slightly different transcript coding for VPS35 displays an opposite trend being over expressed (logFC =8.4), in a similar fashion as VPS53A, highlighting the complexity of the transcriptional regulation of the SNARE processes during host-symbiont cellular interactions. Signal transduction regulates key cellular and metabolic pathways during heat stress. This cascade, however, is not well characterized in Symbiodiniaceae. We found a GABA receptor overexpressed solely in Symbiodinium A3-O. faveolata (absent in the other Symbiodiniaceae transcriptomes) (logFC= 2.44). This receptor is a component of a heterodimeric G-protein coupled receptor for GABA. This metabolite has been involved in pH and osmoregulation, used as an antioxidant and cell signaling 16 . Recently the receptor and GABA have emerged as part of a larger module in the response to biotic and abiotic stress in plants and diatoms 17 . We hypothesize that this receptor may be key to regulate and orchestrate Symbiodiniaceae stress response.
Another gene overexpressed and present uniquely in Symbiodinium A3 is a glutamyl tRNA (Gln) amidotransferase (GATA, logFC= 2.12). During heat stress translation machinery is prone to errors, and misacylated tRNAs may occur 18 . GATA upregulation likely ensures the proper tRNA loading ultimately to maintain the protein synthesis.
In the profile expression of thermotolerant symbionts, Breviolum B5 mounted a reduced specific response that seems effective to contend with the heat stress. In general, this strategy consists in the downregulation of genes involved in ion and amino acid transport. Decreased expression of ZIP transporters (Zinc (Zn 2+ )-Iron (Fe 2+ ) Permease (ZIP) Family) in this symbiont hosted by S. radians as well as in O. faveolata symbionts highlight the importance of metal homeostasis during heat stress. DeSalvo, et al. 19 reported a homolog zinc transporter gene downregulated in the Acropora palmata associated Symbiodinium during bleaching induced by darkness. Zinc is an essential trace metal since it is involved in many biological processes such as redox biochemistry, DNA stability, gene expression, autophagy and apoptosis [20][21][22] . Carbonic anhydrase has a zinc cofactor that catalyzes the reversible conversion of HCO3 and CO2, this CO2 is then concentrated at the active site of Rubisco in order to be fixed during the Calvin cycle. Furthermore, zinc accumulation and uptake are greater in the host and Symbiodinium in symbiosis than aposymbiotic anemones 23 concordantly with the differential expression of different zinc transporters upregulated in symbiotic anemones 24 . An impairment between host and symbiont zinc absorption may occur during heat stress and bleaching, photosymbionts may accumulate metals and nutrients to sustain the energetic budget imposed by the symbiosis.
Over expressed flavodoxin has been documented in diatoms during iron limitation, to overcome the iron expense imposed by the synthesis of chloroplast ferredoxin 25 . We reported a unique flavodoxin gene present and downregulated in Breviolum B5, which indicates that during elevated temperature the thermotolerant symbiont is not iron limited. We corroborate this notion by observing the constitutively expression of ferredoxin genes irrespectively of the temperature increase. Furthermore, in silico prediction of this flavodoxin cellular location indicates likely to reside in the cytoplasm, suggesting a specific role other than electron transference during photosynthesis.
B-hexosaminidase is upregulated as well in Breviolum B5. The specific role of this enzyme during heat stress or symbiosis is unknown, but it has been shown to have diverse functions. In the protist Entamoeaba hystolitica, it acts as a pathogenicity factor by cleaving the terminal N-acetyl-galactosamine (GalNAc) residues from the intestinal epithelial surface after lectinmediated binding 26 . In mycorrhizal symbioses, the GalNAc residues residing in the cell wall have been described as a symbiont associated molecular patterns recognized by the plant. These residues may function as a mediator of the host cell-photosymbiont cell recognition and attachment. Lundgren, et al. 27 identified this gene as a pre-selection candidate where Single Nucleotide Polymorphisms (SNPs) correlated against a temperature gradient on different Pocillopora damicornis populations across the Great Barrier Reef, suggesting that this gene likely plays a role in heat stress tolerance within coral populations.
During heat stress, alternative metabolic pathways may be activated to sustain the symbiosis nutritional demand. 6-phosphogluconolactonase, an enzyme of the pentose phosphate pathway is over expressed in the thermotolerant Breviolum B5. There is some evidence that this enzyme is involved in commensal and pathogenic Enterococcus faecalis interactions with the tobacco hornworm Manduca sexta 28 . The pentose phosphate pathway is a source for the synthesis of aromatic amino acids and nucleotides. This pathway is key for the robust response mounted by Breviolum B5 since it is a metabolic redox sensor that generates the reducing equivalent NADPH, preventing oxidative stress by maintaining the pool of reduced glutathione (GSH), via glutathione reductase (GR). A previous study with Symbiodinaceae species in culture showed that a thermosensitive Breviolum B1 on a sub-lethal temperature increases the activity of GR, contributing to the availability of free glutathione notwithstanding even with increased ROS 29 . While this stress tolerance mechanism malfunctions under excessive heat in Breviolum B1 species, a thermotolerant Fugacium F1 species had an increase on GR activity irrespectively of the heat treatment, highlighting the importance of a redox balancing mechanism conferring thermotolerance 29 . Taken together, these results highlight the unique trajectories followed by each photosymbiont species and the key strategies to sustain symbiosis under adverse temperature increases.

Baseline expression adaptation among coral and symbiont species changing during heat stress (extended)
We identified orthologues common across all three host and symbiont species (1,579 and 1,419 respectively). There are candidates for expression level adaptation in the O. faveolata-P. clivosa host lineage (Fig. 4), such as COPD and COPE (coatomer subunit delta and epsilon) including genes putatively related to cycling between the endoplasmic reticulum (ER) and Golgi apparatus during heat stress, increasing the influx of damaged or misfolded proteins into the ER. Similar results have been found in the symbiotic anemone Exaiptasia pallida, where coatomer protein was greatly elevated in heat-shocked conditions (33.5ºC~24h) 30 . This gene was already expressed in control O. faveolata and P. clivosa samples, which may reflect that these species live at a very fragile threshold of stress. In S. radians, COP (delta/epsilon) is not expressed either before or after stress, possibly indicating life at a closer to optimum temperature.
The transcriptional trajectory of a catecholamine receptor for octopamine (OCTR) is conserved in the Faviina lineage. The susceptible and intermediate species P. clivosa and O. faveolata expressed the octopamine receptor, which is not expressed in S. radians during heat stress. Under this environmental challenge, corals may experience starvation due to photosymbiont cell loss 31 . In other invertebrates such as worms, flies and snails; octopamine, an analogue of norepinephrine, is related to starvation 32,33 . Upon starvation, an increase in octopamine also triggers lipolysis to balance energy homeostasis 32 . Upon octopamine receptor stimulation, cAMP Response Element-Binding Protein (CREB) is activated. Our analysis detected CREB3 as a candidate for adaptive response. This gene is upregulated in P. clivosa and O. faveolata. During ER stress response, CREB3 activates transcription of unfolded protein response (UPR) target genes 34 and it may also be involved in regulating vesicle trafficking from the ER to the Golgi 35,36 . We hypothesize that CREB3 may be activated via OCTR and could be part of putative conserved starvation and UPR modules in corals.
From the Symbiodiniaceae EVE-R analysis, we identified several divergent gene expression candidates for the lineage Breviolum when compared to the Symbiodinium A3 outgroup symbionts. Breviolum and Symbiodinium sensu stricto are thought to have diverged 147.3 MYA, whereas the Breviolum radiation is less than 3.4 MYA 37 (Fig. 1A). The elicited transcriptional profiles in response to heat stress was unique in Breviolum spp. A life history trait that may contribute to a better coupling of the symbiosis is photosymbiont mode of transmission. While O. faveolata and P. clivosa acquire Symbiodinaceae from the environment through horizontal transmission, S. radians acquires its photosymbiont through vertical transmission 38 , leading to a stronger selective pressure imposed by the host environment.
HSP20 is an example of a gene candidate for Symbiodiniaceae expression-level adaptation (Fig. 5). HSP20 is involved in protein folding and heat stress and its expression is conserved in the two Breviolum species in hospite under control conditions. The expression in B. faviinorum is increased significantly in response to a temperature increase, indicating a maximum limit to the damage threshold. However, in Breviolum B5, expression of HSP20 remains at the same basal level of expression even after heat treatment.

Intrinsic expression adaptation in response to heat stress among corals and symbiont species (extended)
Apoptosis genes have been reported in the response to heat stress in corals 9,39,40 . In our study, BAX (p=0.024), caspase 3 (p=0.000056), and caspase 8 (p=0.003) were found to have a divergent expression difference across conditions, which can lead to different heat susceptibilities trajectories due to the differential triggering of apoptotic transduction. P. clivosa expressed these transcripts more abundantly relative to O. faveolata and S. radians.
EVEReSt Symbiodinaceae divergent genes are involved in sugar transport, DNA repair, antioxidant response, lipid metabolism, photorespiration and photosynthesis. During heat stress GLO5, a glycolate oxidase increases slightly in expression. This enzyme, involved in photorespiration, has been proposed to act as a photoprotective mechanism during thermal stress where it is key in diverting excessive excitation energy preventing downstream ROS production 41,42 .
Along with photorespiration, all of these molecules release NH4 which in excess is toxic to the cell. NH4 can be assimilated via the GS-GOGAT pathway, though in plants an increase in ammonia inhibits glutamine synthase (GS) 43 promoting photorespiration during thermal stress 44 . In agreement to these observations, we found evidence for the organic nitrogen accumulation via allantoin. It has been proposed that during stress in plants, allantoins are part of a conserved mechanism to retain NH4 in a non-toxic form 45 . An allantoin permease transcript is upregulated in B. faviinorum and slightly in Breviolum B5. Allantoins may be used as a global adaptive strategy shared among photosynthetic eukaryotes during the inhibition of GS-GOGAT cycle.
Our analysis detected a key gene to regulating the symbiont response to thermal and oxidative stress. SIR2 is a dual function protein that belongs to a family of NAD dependent deacetylases and ADP-ribosyltransferases. It acts as a central metabolic energy sensor of NAD+. Increments in this metabolite cause SIR2 increase [46][47][48] , and in yeast 49 this increases deacetylase histones and key transcription factors involved in the antioxidant and heat shock and UPR response regulation such as catalase, SOD, and HSF1 50-52 . There is a phylogenetically conserved signal of SIR2 expression in Breviolum. Both Breviolum B5 and B. faviinorum express this gene, while there is a divergent lack of expression in Symbiodinium A3. The Breviolum species have opposite photochemical efficiencies in response to a temperature increase, yet the expression of this essential key metabolic sensor/regulator is conserved.
O. faveolata can harbor multiple photosymbiont species at once from different Symbiodiniaceae genera enabling switching under elevated temperatures 53,54 . A holobiont strategy during stressful periods may rely on symbiosis disruption with a subsequent association with another photosymbiont genus. Rather than a pigment content adjustment, the strategy observed in A3 relies on the reduction of Symbiodinium cells. We found an expression divergent strategy mounted by the Symbiodinium A3 in hospite remaining cells. Thioredoxin is upregulated in those Symbiodinium A3 cells and not in Breviolum spp. Cells, highlighting an exceeding ROS production during heat stress. This protein is used to maintain an antioxidant response by reducing-oxidizing its cysteine residues redirecting the electron flow scavenging ROS.
A divergent gene involved in symbiosis maintenance expressed higher in Breviolum B5 but not in the other symbiont species is the α1,2 mannosyltransferase MNN23. In Candida albicans, this protein has been involved in fibrin length and in the immune recognition on the host 55 . In this yeast, the mannans resulting from cleaving are displayed on the cell walls, where they are recognized by the host innate immune system 56,57 including Toll and C-type lectins macrophage mannose receptors 58 .

Glutathione and glyoxylate QMFs in coral holobionts
Genes involved in the glutathione biosynthesis have been found upregulated as a part of an antioxidant response during heat stress in corals 39,59 and Symbiodiniaceae 60 species. An increased glutathione metabolism QMF expression after heat stress is observed in different holobiont members in each holobiont. In the O. faveolata holobiont the increase occurs in the host, whereas in S. radians glutathione module QMF increases in its photosymbiont Breviolum B5 counterpart. In contrast, the associated bacterial community of P. clivosa increased this activity during the heat treatment (Fig. 5A). The fact that both P. clivosa host and its B. faviinorum photosymbiont did not change glutathione expression, highlights the possibility of the thermosensitive coral species limited capacity to activate the biosynthesis of this key antioxidant during heat stress. Glyoxylate metabolism follows a consistent pattern where all the QMF holobiont members from the thermotolerant species are higher (S. radians) when compared to the intermediate and susceptible holobionts (O. faveolata and P. clivosa). This pathway has been described as part of corals strategy to contend with the starvation imposed by heat stress, where fatty acids are converted to carbohydrates 61 and may be key in the resistance to heat stress documented in S. radians holobiont thermotolerance.

Transcriptome sequencing, assembly statistics and annotation
Eighteen cDNA libraries, nine from coral species at the control tank and nine from corals subjected to the treatment tank were constructed and sequenced on the Illumina HiSeq 2500 platform. The number of high-quality, cleaned reads per coral species considering the three control and three treatment replicates yielded 925,234,952 total reads for O. faveolata, 914,039,660 total reads for P. clivosa, and 812,003,378 total reads for S. radians 150-bp, paired end cDNAs (Supplementary Fig.1 Three distinct coral holobiont metatranscriptomes were obtained by combining these highquality, cleaned reads of the three control and three treatment replicates per species. Each coral species was assembled separately. The total number of transcripts assembled for the O. faveolata holobiont was 1,297,972 (N50=1,845; largest transcript=51,329 nucleotides), whereas the P. clivosa holobiont result in 1,255,827 (N50=1,965; largest transcript=41,481), and although the total number of reads for S. radians was the lowest, the transcriptome assembly yielded the highest number of transcripts among species with a total of 2,058,347 (N50=1,405; largest transcript=41,407). Although the total number of reads was less in S. radians, the assembled metatranscriptome was bigger than for the other two species. The higher proportion of bacterial transcripts in S. radians when compared to the other species, revealed an intrinsic diverse bacterial community irrespective of the sequencing depth ( Supplementary  Fig.1A).
To quantify gene abundance through pseudoaligments employing kallisto 63 , transcripts were prior separated per coral host, photosymbiont and bacterial associated communities by each coral species. We obtained 640,070 genes for the O. faveolata host, 175,772 genes from the associated Symbiodinium A3, and 4,928 bacterial genes. For the P. clivosa holobiont 428,458, 166,020, and 2,976 genes were retrieved from the host, B. faviinorum and bacteria respectively. S. radians host genes assessed were 475,058, whereas 155,089 were obtained from Breviolum B%, and 55,109 genes were retrieved for the associated bacterial communities inhabiting this holobiont species.
The quality of each transcriptome for the host and symbiont were evaluated by using the Benchmarking Universal Single-Copy Ortholog (BUSCO) assessment tool 64 . BUSCO analysis of the three coral hosts using a set of 954 conserved metazoan single copy orthologues as a reference, revealed P. clivosa host transcriptome 233 complete and single copy metazoan orthologues and 507 duplicated orthologues, resulting in 77.9% of complete orthologues. The closely related host species, O. faveolata resulted in 322 complete and single copy metazoan orthologues and 395 duplicated orthologues which constitutes 75.2% complete orthologues. S. radians host single copy orthologues were 301 and 404 duplicated orthologues, 73.89% of complete orthologues (Supplementary Fig.1B).
Symbiodiniaceae BUSCO analysis was performed by using a set of 255 eukaryote single copy orthologues as a reference. In contrast to coral hosts completeness, symbiont transcriptomes were not as complete, which has been previously reported for other Symbiodiniaceae transcriptomes 65 . This may be reflecting a limited repertoire of genes expressed in hospite or a scarce representation of Symbiodiniaceae sequences in the BUSCO dataset. Breviolum faviinorum resulted in 108 single copy eukaryote orthologues, and 28 duplicated orthologues, 53.4% of complete orthologues. Breviolum B5 transcriptome yielded 62.3% of complete orthologues, from which 123 were single copy, and 36 were duplicated copy eukaryote orthologues. Symbiodinium A3 is the most complete transcriptome with 70.2% of complete orthologues (77 single copy, 102 duplicated copy eukaryote orthologues) ( Supplementary  Fig.1B).

Supplementary Figures
Supplementary Figure 1. Transcriptome statistics. Total number of reads vs. total number of assembled transcripts. The total number of reads from all the samples including control and experimental samples contributed to the metatranscriptome assemblies per coral and photosymbiont species is represented. S. radians (blue), O. faveolata (yellow), P. clivosa (red) (A). Cumulative percentage of orthologues inferred from the BUSCO search from the Metazoan database against the host transcriptomes from three coral species (left) and using Eukarya database against the photosymbiont transcriptomes from the three respective Symbiodiniaceae species (right) (B). Complete orthologues can be either single-copy (S) or duplicated (D); incomplete orthologues are considered fragmented (F), if orthologues from databases, they are marked as missing (M). Source data are provided as a Source Data file.     Table 7), the rarefaction curves for each coral species at both control (28ºC) and treatment (34ºC), reached their asymptotes consistently. This shows that the read depth was likely adequate. Source data are provided as a Source Data file.

Supplementary Tables
Supplementary Table 1