Transcriptomic Analysis of Thermally Stressed Symbiodinium Reveals Differential Expression of Stress and Metabolism Genes

Endosymbioses between dinoflagellate algae (Symbiodinium sp.) and scleractinian coral species form the foundation of coral reef ecosystems. The coral symbiosis is highly susceptible to elevated temperatures, resulting in coral bleaching, where the algal symbiont is released from host cells. This experiment aimed to determine the transcriptional changes in cultured Symbiodinium, to better understand the response of cellular mechanisms under future temperature conditions. Cultures were exposed to elevated temperatures (average 31°C) or control conditions (24.5°C) for a period of 28 days. Whole transcriptome sequencing of Symbiodinium cells on days 4, 19, and 28 were used to identify differentially expressed genes under thermal stress. A large number of genes representing 37.01% of the transcriptome (∼23,654 unique genes, FDR < 0.05) with differential expression were detected at no less than one of the time points. Consistent with previous studies of Symbiodinium gene expression, fold changes across the transcriptome were low, with 92.49% differentially expressed genes at ≤2-fold change. The transcriptional response included differential expression of genes encoding stress response components such as the antioxidant network and molecular chaperones, cellular components such as core photosynthesis machinery, integral light-harvesting protein complexes and enzymes such as fatty acid desaturases. Differential expression of genes encoding glyoxylate cycle enzymes were also found, representing the first report of this in Symbiodinium. As photosynthate transfer from Symbiodinium to coral hosts provides up to 90% of a coral’s daily energy requirements, the implications of altered metabolic processes from exposure to thermal stress found in this study on coral-Symbiodinium associations are unknown and should be considered when assessing the stability of the symbiotic relationship under future climate conditions.


INTRODUCTION
Unicellular dinoflagellates (genus Symbiodinium) form symbiotic relationships with reef-building corals and other marine invertebrates. The success of coral reef ecosystems is due to mutualistic nutrient exchange between host and endosymbionts (Yellowlees et al., 2008). Exposure to stressors (e.g., elevated temperature) has been attributed as causing coral bleaching, the dysfunction of the symbiotic relationship, resulting in the expulsion of Symbiodinium from the coral host. Coral bleaching is the release of either the Symbiodinium cells from host tissue or the loss of their photosynthetic pigments (Iglesias-Prieto et al., 1992). Depending on the degree of the bleaching event the result may vary, with the host being recolonised by Symbiodinium, disease outbreak or widespread coral mortality (Hoegh-Guldberg, 1999).
Experimentation on Symbiodinium and the coral holobiont has focused on many environmental factors implicated in the onset of coral bleaching including elevated seawater temperatures, acidification, eutrophication (nutrient stress), and disease. The effect of high sea-surface temperatures have been a key focus due to mass coral bleaching events [∼42% GBR reefs bleached in 1998 and ∼54% reefs bleached in 2002 (Berkelmans et al., 2004)], attributed to global climate change (Hoegh-Guldberg, 1999) with the 1998 bleaching event coinciding with an El Niño Southern Oscillation event (Bruno et al., 2001;Fujise et al., 2014). Modeling of bleaching patterns have shown that short periods of high temperature are highly stressful to corals (Berkelmans et al., 2004) with studies emulating these acute conditions in an attempt to understand mechanisms of coral bleaching (Iglesias-Prieto et al., 1992;Ralph et al., 2001Ralph et al., , 2005Takahashi et al., 2008). These studies though extremely valuable in providing an insight into bleaching processes, employ experimental conditions that are not reflective of future longterm predicted sea-surface temperatures or coral bleaching induced by moderate thermal stress over long periods of time (Berkelmans et al., 2004;Fujise et al., 2014;Ainsworth et al., 2016). Additionally mechanisms of thermal acclimation within Symbiodinium are also unknown, with only two recent studies investigating the effect of moderate thermal stress on photobleaching (Takahashi et al., 2013;Fujise et al., 2014).
Dinoflagellates have a number of cellular traits and features that make them unique. Dinoflagellates have large nuclear genomes, the chromosomes remain permanently condensed through the cell cycle and contain highly expressed genes with elevated copy numbers and tandem repeats (Hackett et al., 2004). Approximately half of the dinoflagellates are photosynthetic, having acquired a variety of plastids via endosymbiotic events (Delwiche, 1999;Hackett et al., 2004). In general, the plastids are surrounded by three envelope membranes and have unique chloroplast genome structure, having been reduced to single gene minicircles, with the majority of genes transferred to the nucleus (Hackett et al., 2004;Barbrook et al., 2014;Mungpakdee et al., 2014).
The genus Symbiodinium is a taxonomically diverse species complex, divided into nine phylogenetically distinct clades (A-I; Pochon and Gates, 2010). Additionally intra-cladal diversity exists subdividing the genus further. Associations between Symbiodinium exist with many taxa including ciliates, platyhelminthes, and a variety of marine invertebrates such as cnidarians, molluscs, poriferans, and foraminiferans (Baker, 2003;Coffroth and Santos, 2005;Stat et al., 2008;Venn et al., 2008;. Symbiodinium associations may either exist as mixed populations, at intra-clade or inter-clade level, or as host-symbiont specific interactions (Baker, 2003;Coffroth and Santos, 2005;Pochon and Gates, 2010). Mixed Symbiodinium populations may occur in different proportions, with a single dominant species with multiple species detected at low abundances (Baker, 2003). Therefore, understanding the different Symbiodinium strains may provide an insight into the complexity of symbiotic interactions that are observed. Development of high throughput sequencing technologies has seen many advances in Symbiodinium genomics and transcriptomics [summarized in Shinzato et al. (2014)]. These include publication of three Symbiodinium nuclear genomes, Symbiodinium minutum (clade B1; Shoguchi et al., 2013), the Symbiodinium kawagutii (clade F1) nuclear genome (Lin et al., 2015), and the Symbiodinium microadriaticum (clade A1) nuclear genome (Aranda et al., 2016) and a further 13 sequencing projects representing organelle genomes and transcriptomes of various Symbiodinium sp. (clades A-D and F; Supplementary Table S1). The Symbiodinium minutum (clade B) genome has identified a large number of protein-coding genes (41,925), attributed to duplication events due to the presence of highly repetitive gene clusters (Shoguchi et al., 2013;Shinzato et al., 2014). Further, expansions of regulator of chromosome condensation family protein (RCC1) provide a possible molecular basis for permanently condensed chromatin observed in dinoflagellates (Shoguchi et al., 2013;Shinzato et al., 2014). Other expanded multi-copy gene families identified include ion channel proteins and the chlorophyll a/b-binding proteins (lhcb, PF00504; Shoguchi et al., 2013). Study of the S. minutum mitochondrial genome  has shown that it is greatly expanded and fragmented, whereas, the plastid genome is greatly reduced with most (all but 14 located on DNA minicircles) plastid-associated genes being transferred to the nucleus . Mechanisms for RNA editing in Symbiodinium have also been revealed from studying the plastid genome . Symbiodinium transcriptomes and EST datasets have been published using different sequencing technologies representing various clades that associate with a variety of hosts (summarized in Leggat et al., 2011b;Rosic et al., 2014;Shinzato et al., 2014;Xiang et al., 2015;Levin et al., 2016). Due to the interest in studying establishment of symbiosis, gene expression studies have focused on genes associated with these molecular mechanisms (Voolstra et al., 2009a;Mohamed et al., 2016).
Symbiodinium genome and transcriptome data have allowed for studies of specific gene families of interest . Integral light-harvesting complex (LHC) gene families [chlorophyll a-chlorophyll c 2 -peridinin protein complex (acpPC)] in Symbiodinium have been investigated due to their unique gene arrangement (annotated as chlorophyll a/b-binding proteins in some instances; Takahashi et al., 2008;Boldt et al., 2012;Maruyama et al., 2015;Xiang et al., 2015). Gene-mining has revealed high diversification of the integral light-harvesting gene family (acpPC) in Symbiodinium species (Boldt et al., 2012;Maruyama et al., 2015), supporting theories of intra-and intergenic duplication from ancestral LHC gene(s) possibly of red algal origin (Engelken et al., 2010;Boldt et al., 2012;Maruyama et al., 2015). The functional purpose for the highly redundant multicopy gene family has been hypothesized as photoprotective (Boldt et al., 2012;Maruyama et al., 2015;Xiang et al., 2015), and differential expression of acpPC genes in Symbiodinium has been detected under thermal stress in targeted gene expression studies (Takahashi et al., 2008;Gierz et al., 2016). This study investigates the expression of Symbiodinium sp. (clade F) lightharvesting acpPC genes under thermal stress at the transcriptome level.
Whilst these next-generation sequencing projects have generated large quantities of data and provide basal reference transcriptomes, some of these studies lack comparison between stress and non-stress conditions, which is one of the aims of this study. Recently, the transcriptional response of cultured Symbiodinium (strain SSB01) maintained under different light intensities with different growth conditions (Xiang et al., 2015), and of two type C1 Symbiodinium cultures exposed to a 13-days thermal stress (32 • C) (Levin et al., 2016) have been published, and are providing insights into Symbiodinium transcriptomes under stress conditions. Previous studies of Symbiodinium using targeted gene expression analysis (quantitative-PCR), demonstrate that alteration of gene expression generally occurs at low fold changes, with regulation hypothesized to occur instead through translational or post-translational mechanisms (Leggat et al., 2011a;Rosic et al., 2011;Krueger et al., 2015;Gierz et al., 2016). This study employed Illumina RNA-Seq and used four biological replicates at each time point per condition to ensure a robust analysis. Samples were selected for sequencing to encapsulate stages of exposure to elevated temperatures at beginning (day 4), middle (day 19), and end (day 28) of the experimental period. This experimental design has allowed us to investigate the potential effects of exposure to elevated temperatures on Symbiodinium molecular processes.

Culture Conditions and Experimental Design
Cultures of Symbiodinium sp. [clade F (ITS2 rDNA region; Supplementary Table S2)] were obtained from the Australian Institute of Marine Science (AIMS). Cells were grown in ASP-8A media (Blank, 1987), in 75 cm 2 vented tissue culture flasks at 25 • C. Cultures were maintained in light cabinets (Thermoline Scientific refrigerated incubator, Sanyo) at an irradiance of 80-100 µmol quanta m −2 s −1 measured using a LI-193SA Underwater Spherical Quantum Sensor with a LI-250A Light Meter (LI-COR R Inc., Lincoln, NE, USA).
Over a period of 28 days cultures of Symbiodinium were exposed to control temperatures (24.5-25 • C) or thermal treatment temperatures of approximately (30-31.5 • C; Figure 1A). In both treatment and controls daily fluctuations in temperature reflect 12:12 h light:dark photoperiod. Temperatures in each incubator were recorded every 10 min with HOBO R temperature/alarm pendant data loggers (Onset, MA, USA). To maintain an adequate number of cells throughout the experiment 16 flasks were designated to each treatment with five replicate culture flasks sampled at each time point per treatment (Supplementary Figure S1). Samples were taken at 10 h into the light photoperiod.

Symbiodinium Density and Chlorophyll Pigment Analysis
On days 1, 4, 7, 19, and 28, five replicate culture flasks were taken from each treatment for cell number approximation and pigment quantification. Cell numbers were determined using a Neubauer haemocytometer, with replicate cell counts performed (n = 4). Symbiodinium cells were pelleted by centrifugation at 4500 × g for 3 min for chlorophyll a and c quantification. Chlorophylls were extracted in 90% acetone for 20 h in the dark at 4 • C and pigment content quantified using the equations of Jeffrey and Humphrey (1975).

Imaging-Pulse Amplitude-Modulated Fluorometry
Imaging-Pulse amplitude-modulated (PAM) fluorometry (MAXI Imaging-PAM and ImagingWin software, Walz, Effeltrich, Germany) were used to measure photosynthetic efficiency of Symbiodinium cultures. The preprogrammed Induction Curve + Recovery kinetic recording type were used to examine the ability of Symbiodinium to dissipate excess light energy and recover from light stress after exposure to elevated temperature. Symbiodinium cells were aseptically transferred to 15 mL falcon tubes and pelleted by centrifugation at 4,500 × g for 3 min. Cell pellets were resuspended in approximately 250 microliters ASP-8A media and transferred to a 96-well plate and darkadapted for 20 min prior to analysis. Following dark-adaption, minimum chlorophyll fluorescence (F o ) were determined using blue measuring light (intensity 2), and maximum chlorophyll fluorescence (F m ) were determined by applying a pulse (0.72 s) of saturating light (intensity 5, ∼2800 µmol quanta m −2 s −1 ) allowing calculation of the dark-adapted maximal quantum yield of PS II (F v /F m ). For the Induction Curve, actinic illumination (254 µmol quanta m −2 s −1 , intensity 6) was switched on and 15 saturating pulses of photosynthetically active radiation [∼2800 µmol quanta m −2 s −1 (intensity 5, 0.72 s)] were applied at 20 s intervals for 5 min. During the Recovery phase, a further 16 saturation pulses were applied within a 14 min period without actinic illumination, time between each pulse exponentially increased. Imaging-PAM fluorometry were used to determine photo-kinetic parameters, such as the maximal quantum yield of PS II (F v /F m ), effective quantum yield at the end of the induction curve, and non-photochemical quenching (NPQ) at the beginning of the recovery phase. Light levels were measured using a LI-190SA Quantum Sensor with a LI-250A Light Meter (LI-COR R Inc., Lincoln, NE, USA).

Data Analysis
Statistics software package SPSS statistics (v19.0, IBM, USA) were used for all statistical analyses of physiological parameters. A generalized linear model with 'day' and 'treatment' as main effects and 'day × treatment' as an interaction were used for pairwise comparison of cell density, chlorophyll a and c, and imaging-PAM with sequential Bonferroni post hoc test to determine significant differences between control and treatment. The generalized linear model approach was chosen as samples taken at each time point were considered independent.

RNA Isolation and Sequencing
For total RNA isolation, 15 mL of cells were pelleted by centrifugation at 4,500 × g for 3 min. Cells were then transferred to a screw cape tube and centrifuged at 8,050 × g for 3 min. Pellets were snap frozen in liquid nitrogen and stored at −80 • C. Total RNA were extracted using the RNeasy Plant Mini Kit (Qiagen, USA). Symbiodinium cells were first lysed using the FastPrep R -24 sample preparation system (MP Biomedicals, Australia). Four hundred and fifty microliters of Buffer RLT containing 1% β-mercaptoethanol were used to resuspend cells and were transferred to a lysing matrix D tube (MP Biomedicals, Australia), and shaken three times for 40 s at 4.5 M s −1 to lyse the cells. Total RNAs were isolated from cells using the Purification of Total RNA from plant cells and tissues protocol. The optional oncolumn DNase Digestion was performed using the RNase-Free DNase set (Qiagen, USA) as per the manufacturer's protocol.
Concentrations of isolated total RNAs were checked using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and quality were assessed using a Bioanalyzer (Agilent) prior to library generation. RNA-Seq libraries were prepared by the Australian Genome Research Facility (AGRF, Melbourne, Australia), using the Illumina TruSeq RNA sample preparation kit V2 (Illumina) and the standard Illumina protocols. Multiplexed sequencing were performed on the 24 libraries by AGRF on an Illumina HiSeq 2000 platform on two lanes, generating over 370 million 100 bp single-end reads.

RNA-Seq Analysis
Image analyses were performed in real time by the HiSeq Control Software (HCS) v1.4.8 and Real Time Analysis (RTA) v1.18.61, running on the instrument computer. The Illumina CASAVA (Consensus Assessment of Sequence and Variation) 1.8.2 pipeline was used to generate the sequence data. Sequencing reads were filtered according to their multiplexing tags and multiplexing tags were removed. The sequenced library were mapped against a Symbiodinium reference transcriptome generated from the same culture (Bobeszko et al., unpublished). Briefly, for the Symbiodinium reference transcriptome (BioProject number PRJNA371519), paired-end reads were trimmed, removing sequence adaptors and low-quality regions using libngs 1 with a minimum quality of 20 bp and a minimum size of 75 bp. Trimmed reads were then assembled with Trinity (Grabherr et al., 2011) and the resulting assembly clustered using CD-HIT-EST (Fu et al., 2012) using 90% similarity and a word size of 8. TransDecoder (Haas et al., 2013) and Blast2GO (Conesa et al., 2005) were used to predict protein-coding sequences.
The single-end sequenced library were mapped using the ArrayStar application and QSeq module of the DNAStar Lasergene Genomics Suite (Version 11; DNASTAR, Inc., Madison, WI, USA) with all parameters set to defaults. QSeq parameters were set to default, read counts were normalized via RPKM (reads per kilo base of exon model per million mapped reads; Mortazavi et al., 2008) and processed genes defined as 'use sequences as genes.' The RPKM method standardizes molar concentration of transcripts by determining transcript length (in kilobases) and the read abundances by dividing each read count by the library-size (in millions) to normalize (Mortazavi et al., 2008). The RPKM normalization method is accepted in RNA-Seq analysis as it removes technical biases introduced by sequencelength and library-size Conesa et al., 2016), and is suitable for the single-end reads generated from sequencing. Statistically significant expression changes between the control and treatment data sets were determined using the student's t-test with multiple test correction by Benjamini-Hochberg false discovery rate (FDR < 0.05). The results for each time point comparison [day 4,9,471 DEG;day 19,12,701 DEG;and day 28,13,269 DEG (Figures 2A,B)], uniquely differentially expressed at any time point (23,654 DEG) and differentially expressed at every time point (2,798 DEG) were exported and saved in Microsoft Excel (Microsoft, Redmond, WA, USA) and used for subsequent analysis.
Nucleotide fasta files for candidate transcripts were generated using the 'SeqinR' package 2 (Charif and Lobry, 2007), venn diagrams were drawn with the 'VennDiagram' package 3 (Chen and Boutros, 2011), using R version 3.3.2 4 (R Core Team, 2014). Distribution of gene ontology (GO) terms within the sequences that were differentially expressed at all time point (2,798 DEG) were determined using Blast2GO (v3.3.5;Conesa et al., 2005). Briefly, nucleotide sequences were imported into Blast2GO and GO terms generated using default parameters for the blast (blastx), mapping and annotation steps (Conesa and Götz, 2008). To visualize GO term distributions, combined GO annotation graphs were produced (default settings) and used to generate graph level pie charts of biological processes ( Figure 2C Figure S3), sequences were filtered by GO terms and a cut-off ontology level of 3 were applied, slices smaller than 2% were grouped in the "other" slice. To further analyze the distribution of GO terms within the biological process category data were split into three sets, significantly upregulated at all time points (1,428 DEG), significantly downregulated at all time points (1,331 DEG), or significantly expressed (up and down) at all time points (39 DEG; Supplementary Figure S4). Sequences with GO terms contributing to biological process categories were selected/filtered and used to further identify genes and pathways differentially expressed over all time points. Heatmaps were drawn with the 'pheatmap' package 5 (Kolde, 2015).

Data Deposition
The Illumina sequenced read data reported in this article have been deposited into the National Center for Biotechnology Information (NCBI) Sequence Read Archive under the accession number SRA467551, which is associated with BioProject number PRJNA342240.
FIGURE 2 | Thermal stress-induced differential gene expression. (A) Venn diagram illustrates the differential expression of 35,441 genes (FDR < 0.05) of Symbiodinium sp., after exposure to thermal stress for 4, 19, and 28 days. Venn diagram generated using the VennDiagram package in R (B) Illustration of the distribution of upregulated or downregulated genes (FDR < 0.05) at days sampled following exposure to thermal stress. (C) Visualization of the distribution of biological process gene ontology (GO) classifications for the 2,798 genes differentially expressed at all time points in Symbiodinium exposed to thermal stress (FDR < 0.05). GO annotation graph produced using Blast2GO, GO categories displayed at ontology level 3 and slices smaller than 2% grouped into the 'other' term, numbers displayed represent the number of sequences assigned to each ontology category.

Physiological Responses of Symbiodinium to Thermal Stress
Symbiodinium cell densities were significantly decreased in thermally stressed cultures from day 4 onward of the experiment ( Figure 1B). Flasks sampled on day 1 were sampled again on day 14, day 4 again on day 19, and day 7 again on day 28.
Maximum quantum yield of photosynthesis, F v /F m , were measured on sampling days throughout the experiment. For cultures maintained at control temperatures dark-adapted yield ranged between 0.536 and 0.616 (average 0.577, standard error ± 0.003; Figure 1C). Analysis of the dark-adapted yield found that there were a significant interaction between treatment and day (p < 0.001, df = 5; Figure 1C) and for days (p < 0.001, df = 5) and treatments (p < 0.001, df = 1). Sequential Bonferroni post hoc analysis found that F v /F m decreased in the heated treatment on day 14 (p < 0.05), day 19 (p < 0.01), and day 28 (p < 0.01) of the experiment ( Figure 1C). Symbiodinium effective quantum yield at the end of the induction phase were also determined ( Figure 1D). This measurement is taken after the cells were exposed to a phase of actinic light with periodic saturating pulses. Cells maintained at ∼31 • C only exhibited decreased effective quantum yield on day 4 (p < 0.01; Figure 1D). NPQ a proxy for cell protective mechanisms was measured over the course of the experiment. Analysis of NPQ found that there were significant effects in the interaction between day and treatment (p < 0.001, df = 5) and between days (p < 0.001, df = 5; Figure 1E). NPQ values did vary slightly from controls in thermally stressed cells, though it only significantly decreased on day 19 (p < 0.001; Figure 1E).
Chlorophyll a content per Symbiodinium cell were increased in treatment cells compared to control cells over the experimental period ( Figure 1F). Analysis of chlorophyll a content found that there were significant interaction day × treatment (p < 0.001, df = 5; Figure 1F) and differences between day (p < 0.001, df = 5) and treatment (p < 0.001, df = 1). Though chlorophyll c content were increased in thermally stressed cells no significant differences were found between control and treatment cells ( Figure 1G). No significant differences were found in the ratio of chlorophyll c to chlorophyll a in Symbiodinium cells throughout the experiment ( Figure 1H).

Differential Gene Expression at a Pre-Bleaching Temperature Threshold
Differential gene expression in Symbiodinium in response to elevated temperature were determined at days 4, 19, and 28. Overall, the number of unique DEGs detected throughout the thermal stress accounts for 37.01% of the transcriptome (FDR < 0.05). Though a large number of differentially regulated genes were detected, only 2.78% (1,776 contigs) were unique genes with ≥2-fold change in expression. Comparison of gene expression between all time points identified 2,798 common transcripts differentially expressed (both up and down) under elevated temperature at all time points (Figure 2A) and these were further analyzed. Gene OntologyGO analysis of these common transcripts identified biological processes including 409 genes encoding proteins for cellular metabolic processes, 204 genes encoding proteins involved in cellular component organization, and 133 genes encoding proteins associated with response to stress ( Figure 2C).
Transcripts encoding 41 heat shock proteins (HSPs; HSP90, HSP70, HSP20), heat shock transcription factors (HSTF), and molecular chaperones (DNAJ) were detected with differential expression through all time points. In the downregulated data set, nine transcripts encoding HSPs, DNAJs, and heat shock-related proteins (HRPs) were detected with significantly decreased expression at all three time points (Figure 3A; Supplementary Table S3). In the upregulated data set, two HSP70 transcripts (comp71407_c0 and comp76430_c0) were detected with significantly increased expression at all time points (Figure 3A; Supplementary Table S3). Further, chloroplast targeted HSPs exhibited no difference in expression patterns compared with cytosolic associated HSPs.
DNA damage repair and proteasomal degradation pathways were differentially regulated in thermally stressed Symbiodinium cells. Nine DNA repair RAD proteins (RAD5, RAD16, and RAD23-1) were annotated with differential expression over the course of the experiment ( Figure 3B; Supplementary Table S3). A single RAD50 DNA repair transcript (comp72013_c0) displayed significantly increased expression at all time points (Figure 3B; Supplementary Table S3). DNA photolyases (PHR) and cryptochrome transcripts [cryptochrome DASH (CRYD)] and cryptochrome 2 (CRY2) showed varied expression, six transcripts were detected with significantly increased expression, and six transcripts displayed significantly decreased expression over the experiment (Figure 3B; Supplementary  Table S3). Ubiquitin proteasome pathway (UPP) components were detected with significant fold changes, including ubiquitination enzymes for conjugation (E1, E2, and E3 enzymes) and deubiquitination (DUBs) and ubiquitin-like modifiers (SUMO, NEDD8, ISG15, APG8, and APG12; Supplementary Table S4). Though a large number (>260) of UPP enzymes and modifiers displayed significant fold changes (Supplementary Table S4), only five transcripts were detected with significantly decreased expression at all time points and 15 transcripts displayed significantly increased FIGURE 3 | Heatmap illustration of differentially expressed stress response genes (FDR < 0.05) in Symbiodinium exposed to thermal stress at days 4, 19, and 28. Data are expressed as fold-changes relative to control; only significant data are shown (p < 0.05), non-significant data denoted as white boxes. (A) Differential expression of antioxidant defenses (Continued)
(B) Differential expression of stress related transcripts including genes encoding DNA damage repair proteins, selected ubiquitin proteasome pathway components, metacaspases, and anti-apoptosis proteins. fold change at all time points ( Figure 3B; Supplementary  Table S3). Apoptosis-like transcripts were detected in thermally stressed Symbiodinium cells. Seventeen transcripts encoding three metacaspase 1 isoforms (MCA1, MCA1A, and MCA1B) were detected with differential expression in thermally stressed cells (Figure 3B; Supplementary Table S3). Three detected metacaspase contigs (two MCA1 isoforms and one MCA1B isoform) were significantly downregulated at all time points (Figure 3B; Supplementary Table S3). One transcript (comp71271_c0) encoding a MCA1A isoform displayed significantly increased expression at all time points (Figure 3B; Supplementary Table S3). A further 13 metacaspase transcripts were detected with differential expression over the experiment, 10 transcripts were significantly upregulated at at least one time point, and three were significantly downregulated at at least one time point (FDR < 0.05; Figure 3B; Supplementary Table S3). Six transcripts encoding apoptosis-inducing factor homologs (AIFB and AIFM3) were detected, one AIFM3 transcript and three AIFB transcripts displayed increased gene expression and two displayed decreased expression across the experiment though not all time points were significantly different to controls ( Figure 3B; Supplementary Table S3).

Photosynthesis Related Genes
Nine transcripts encoding polypeptide subunits of PS II (psbC, psbF, psbH, psbK, psbO, psbP, and psbY) were detected with significant changes in expression ( Figure 4A). Genes encoding PS II D1 protein, PS II D2 protein, and CP47 reaction center protein (encoded on plastid mini-circles by psbA, psbD, and psbB, respectively) were annotated in the transcriptome (Supplementary Table S5), though no significant fold changes were detected under the experimental conditions used (FDR < 0.05). One transcript encoding psbF (comp42202_c0), which encodes cytochrome b559 β forming part of the reaction center core of PS II, displayed significant downregulation at all time points (Figure 4A). Seven transcripts encoding polypeptide subunits of PS I (psaA, psaC, psaD, psaF, psaJ, and psaL) were detected with significant changes in expression ( Figure 4A). Plastid minicircle genes encoding integral membrane peptide subunits of PS I were annotated in the transcriptome, expression of the PS I P700 chlorophyll a apoprotein A1 gene (psaA) were decreased at day 28 (3.306 down fold change, p < 0.02; Figure 4A), no significant differences in expression were detected for the PS I P700 chlorophyll a apoprotein A2 (psaB; Supplementary Table S5). Differential expression of transcripts encoding ferredoxin-NADP+ reductase [petH (comp42084_c0)] and ferredoxin [petF (comp68408_c0 and comp35855_c0)] were detected with significant differences in expression under thermal stress ( Figure 4A; Supplementary  Table S3).
Plastid-targeted genes were differentially expressed in thermally stressed Symbiodinium cells. For example nine transcripts encoding the unique Form II ribulose 1,5bisphosphate carboxylase/oxygenase (RuBisCo) enzyme (rbcL, large subunit) were differentially expressed under thermal stress ( Figure 4A; Supplementary Table S3). Expression of RuBisCo transcripts were temporally varied, six isoforms displayed significantly increased expression at day 28 (largest 1.465 up fold change). Additionally, a single RuBisCo isoform (comp68158_c0) with significant differential expression at all time points were identified ( Figure 4A; Supplementary Table S3).
Genes involved in photoprotective mechanisms in Symbiodinium including NPQ were differentially expressed under thermal stress. In Symbiodinium, the xanthophyll cycle involves the de-epoxidation and epoxidation reactions of diadinoxanthin/diatoxanthin for energy dissipation and to regulate the amount of energy reaching the photosystem reaction centers. Three violaxanthin de-epoxidase (vde) transcripts were detected with significant changes in expression at time points throughout the thermal stress ( Figure 4A; Supplementary Table S3). Two zeaxanthin epoxidase (zep) transcripts (comp79868_c0 and comp88413_c0), displayed significantly increased fold changes through all time points (Figure 4A; Supplementary Table S3).
Expression of the nuclear-encoded LHC proteins responsible for enhancing light transfer to core photosystems and photoprotection were determined in thermally stressed cells. A transcript encoding the extrinsic water-soluble LHC peridininchlorophyll a-binding protein (PCP) was detected in the Symbiodinium transcriptome (comp80938_c0; Supplementary  Table S5), though no significant fold changes were detected under the experimental conditions used (FDR < 0.05). Fiftyfour transcripts encoding the Symbiodinium specific intrinsic membrane-bound LHC isoforms (acpPC) were detected with significant fold changes during the thermal stress (PFAM ID PF00504.16). Blast annotations shown, (cb, chlorophyll binding protein; ccac, caroteno-chlorophyll a-c binding protein; fcp, fucoxanthin-chlorophyll a-c binding protein; lh18, LHC I protein; li818, chlorophyll a-b binding protein l1818; Figure 4A; Supplementary Table S3). Of these 54 intrinsic LHC transcripts, 14 were detected with significant fold changes at all time points (13 were upregulated and one transcript displayed mixed expression; Figure 4B). The remaining 40 light-harvesting protein complex transcripts detected displayed predominately increased expression at all time points (85% up fold change), though not all were significantly different to controls (Figure 4A; Supplementary Table S3). Additionally, a further 29 transcripts annotated in the transcriptome as light-harvesting protein complexes displayed no change in expression under the experimental conditions (Supplementary Table S5).
were detected with significantly increased expression on day 28 ( Figure 4B). One of the β-ketothiolases (comp79906_c0) was annotated as a peroxisomal associated transcript, and displayed significantly increased expression on day 19 (1.276 up, p < 0.01) and day 28 (1.169 up, p < 0.03; Figure 4B). Further, three transcripts encoding peroxisome membrane proteins (PMP34), two peroxisome adenine nucleotide carrier isoforms (PNC1 and PNC2), and two peroxisome ATP-binding cassette subfamily D members were all detected with significant increases in expression under these conditions (Supplementary Table S5). Glyoxylate cycle and gluconeogenic pathway enzymes were detected with differential expression in thermally stressed cells. Transcripts encoding enzymes of the glyoxylate cycle including citrate synthase (CS), aconitase (acnB), isocitrate lyase (aceA), malate synthase (aceB), and malate dehydrogenase (MDH2) were detected with differential expression in thermally stressed Symbiodinium cells ( Figure 4B). Specifically, four transcripts encoding isocitrate lyases (aceA) were detected with significant changes in expression under thermal stress ( Figure 4B; Supplementary Table S3). Two of the isocitrate lyase contigs displayed decreased expression on day 4, whereas on day 28 three of the isocitrate lyase contigs displayed significantly increased fold change compared to controls ( Figure 4B). Two malate synthase (aceB) transcripts were detected with decreased expression, comp85614_c0 on day 4 (1.171 down, p < 0.04) and comp96309_c0 on day 19 (1.259 down, p < 0.002), and day 28 (1.263 down, p < 0.008; Figure 4B). Three transcripts encoding succinate dehydrogenase flavoprotein subunits were also detected with significantly increased expression (comp25945_c0 on day 28, 1.294 up, p < 0.009), and comp72308_c0 on day 19, 1.433 up, p < 0.02 and day 28, 1.287 up, p < 0.004, and comp337315_c0 on day 19, 31.456 up, p < 0.02 and day 28 25.270 up, p < 0.03 ( Figure 4B). In addition, three contigs encoding phosphoenolpyruvate carboxykinase (PEPCK) were detected with differential expression in thermally stressed cells (Figure 4B; Supplementary Table S3).
Serine/threonine-protein kinases are crucial components of diverse signaling pathways and for regulation of cell proliferation, meiosis, and apoptosis. One hundred and seventy-seven transcripts representing more than 20 serine/threonine-protein kinase families were detected with significant changes in expression in thermally stressed cells (Supplementary Table S3). Nine transcripts encoding three classes of aurora kinases [Aurora-A, Aurora-B, and Aurora-C (AurK)] and two aurora-related kinases (ARKs) were detected with significant fold changes throughout the experiment (Figure 4B; Supplementary Table S3). Twentytwo never-in-mitosis A serine/threonine kinase (Nek) transcripts representing seven Nek families were detected with differential expression (Figure 4B; Supplementary Table S3). One transcript encoding a cyclin-dependent kinase (CDK5, comp29198_c0) was detected with significantly decreased expression at all time points (Figure 4B; Supplementary  Table S3).
Meiosis-specific and meiosis-related genes previously annotated in Symbiodinium (Chi et al., 2014;Levin et al., 2016) were detected in the transcriptome and in the repertoire of differentially expressed transcripts. Eight meiosis-specific transcripts representing four genes (three Dmc1, two Hop2, one Mnd1, and two Msh4) were annotated with differential expression in thermally stressed cells (Figure 4C; Supplementary Table S3). Twenty meiosis-related transcripts representing 17 genes were also annotated with differential expression in thermally stressed Symbiodinium cells (Figure 4C; Supplementary Table S3). Additionally, 52 transcripts encoding protein MEI2 and MEI2-like isoforms were detected with differential expression (Figure 4C; Supplementary Table S3). MEI2 genes have been annotated in Schizosaccharomyces pombe and MEI2-like proteins have been annotated in Arabidopsis thaliana and Oryza sativa subsp. japonica. Two radial spoke head homologs (Rsph1) were detected with differential expression (Figure 4C; Supplementary  Table S3), previously detected around the chromosomes during metaphase in male gametes undergoing meiotic division. However, Rsph1 has also been implicated in axenomal central pair regulating dynein activity for flagella and cilia movement.

DISCUSSION
Predictions for future climate conditions estimate the seasurface temperature to rise between 1.1 and 6.4 • C depending on emissions scenarios (Solomon et al., 2007). The implications of this change for coral reefs and the marine ecosystems they support are unknown. However, this prolonged exposure of coral reefs to elevated temperatures may have catastrophic consequences and result in the loss of many species. Through examining the physiological response of Symbiodinium cultures to elevated temperatures, and linking this to molecular processes that are altered under thermal stress we may begin to understand how reefs may be affected in the future. In this study, we exposed cultured Symbiodinium sp. (clade F) to thermal stress (28 days, ∼30-31.5 • C; Figure 1A) and generated a library of contigs that represent the transcriptome under future temperature conditions.
Analysis of the Symbiodinium transcriptome revealed more than 37% were differentially expressed under thermal stress (FDR < 0.05). A large number of DEGs (23,654 unique contigs) were detected, of which 2,798 were differentially expressed at all time points (Figure 2A). The majority of DEGs [21,878 genes (92.49%)], displayed a ≤2-fold change in expression. This is reflective of previous targeted expression studies in Symbiodinium where relatively small changes in gene expression were determined (Leggat et al., 2011a;Rosic et al., 2011;McGinley et al., 2012;Ogawa et al., 2013;Gierz et al., 2016), due to this it has been hypothesized that translational or post-translational regulation may be critical in Symbiodinium cellular responses. Biological process GO visualization revealed that distribution of combined time point DEGs (2,798 transcripts) were relatively even in assigned terms between those up and downregulated transcripts (Supplementary Figure S4). Under these conditions regulation of molecular processes in this Symbiodinium strain is variable, with components of many pathways exhibiting dissimilar expression patterns.
Over the course of the experiment temperature significantly effected Symbiodinium density and photosynthetic efficiency. Previously, growth rates of six cultured Symbiodinium strains (representing clades A, B, C, D, and F) with various thermal tolerances where determined at different temperatures (25, 30, and 33 • C; Karim et al., 2015). Showing that increased temperatures can alter growth rates and photosynthetic efficiency differently, resulting in the classification of three categories depending on the thermal tolerance of the strain (Karim et al., 2015). In this study, Symbiodinium cell density were decreased from day 4 compared to controls (Figure 1B), and the photosynthetic ability of thermally stressed cells were maintained until day 14 ( Figure 1C). The slight but significantly depressed dark-adapted yield (day 14, 19, and 28; Figure 1C) indicates that cells exposed to elevated temperatures were exhibiting altered photosynthetic efficiencies in response to thermal stress. However, as complete loss of photosynthetic efficiency did not occur (Figure 1C), this strain is categorized as photosynthetically tolerant though growth response were highly sensitive to elevated temperature (Karim et al., 2015).
Analysis of chlorophyll pigments in Symbiodinium found significantly increased chlorophyll content in a manner seen previously (McBride et al., 2009;Ogawa et al., 2013;Gierz et al., 2016). Chlorophyll a concentration were significantly higher in thermally stressed cells on days 19 and 28 ( Figure 1F). Comparison of growth rates and chlorophyll content in Symbiodinium californium found cultures exhibiting low growth rates (incubated at 5, 10, and 30 • C) also showed an increased chlorophyll a content, whereas those actively growing had reduced chlorophyll a content (McBride et al., 2009). In phytoplankton, variation of chlorophyll a-specific absorption has been attributed to packaging of chlorophylls within different pigment-protein complexes (chlorophyll a-chlorophyll c-peridinin (ACP) versus PSI; Bissett et al., 1997). It is possible that the increased chlorophyll a content observed in thermally stressed Symbiodinium cells may be due to changes in the specific pigment-protein complexes within the chloroplasts. Analysis of chlorophyll c content ( Figure 1G) and the ratio of chlorophyll c to chlorophyll a ( Figure 1H) found that there were no significant differences between controls and thermally stressed cells. Biosynthesis of chlorophyll involves the ATPdependent insertion of a magnesium ion into protoporphyrin IX, catalyzed by three magnesium-chelatase subunits (Von Wettstein et al., 1995). Six transcripts encoding two of the three magnesium-chelatase subunits (two chlI and four chlH) were differentially expressed in thermally stressed cells ( Figure 4B; Supplementary Table S3). However, expression of these subunits were not consistent with chlI subunits detected with significant decreased fold changes, no changes detected in a chlD subunit (Supplementary Table S5), and four chlH subunit with significantly increased expression during exposure to increased temperatures ( Figure 4B; Supplementary  Table S3). The implications for the variable expression of these subunits critical for chlorophyll biosynthesis requires further investigation.
Measurements of photosynthetic ability are often employed as indicators of Symbiodinium cell health (Buxton et al., 2012;Hill and Takahashi, 2014). In this study, Induction Curve + Recovery kinetic recording type were used to determine the ability of cells to respond to light stress. Significant decreases in effective quantum yield at the end of the induction phase were detected on day 4 ( Figure 1D) and may be indicative of an early photosynthetic response to elevated temperatures. Throughout the remainder of the experiment no significant difference in effective quantum yield between controls and thermally stressed cells were detected after the recovery phase, though values were slightly depressed in treatment cells (Figure 1D). NPQ were found to be significantly decreased in treated cells on day 19, and may indicate that the photoprotective mechanisms of cells were impacted by exposure to thermal stress ( Figure 1E). Though little changes were observed in NPQ, we detected xanthophyll cycle enzyme genes (VDE and ZEP) with significantly different expression ( Figure 4A; Supplementary Table S3). VDE and ZEP are responsible for the epoxidation and de-epoxidation of dinoxanthin/diadinoxanthin as a photoprotective mechanism by dissipating excess energy.

Differential Expression of the Symbiodinium Antioxidant Network
Photobleaching in Symbiodinium induced by thermal stress and high solar irradiance has been linked to oxidative damage resulting from the production of reactive oxygen species (ROS; Murata et al., 2007;Takahashi et al., 2008;Krueger et al., 2014). ROS can be deleterious to cells resulting in oxidative damage to lipids, proteins, and DNA but may also function as second messengers for signal transduction (Lesser, 2006). Defenses against ROS in Symbiodinium and other photoautotrophs include an antioxidant network of enzymes such as SODs, catalases, and peroxidases as well as non-enzymatic antioxidants such as glutathione, Trx, Prx, and carotenoids (via the xanthophyll cycle; Lesser, 2006;Bayer et al., 2012;Krueger et al., 2014Krueger et al., , 2015. Comparison of Symbiodinium types has shown that the antioxidant network and the antioxidant response to thermal stress can vary between strains of different thermal tolerance with implications for photosynthesis and cell viability (Krueger et al., 2014). Within Symbiodinium cells exposed to thermal stress transcripts encoding CuZnSOD, MnSOD, KatG, and ZEP were all significantly upregulated whereas, expression of APX and VDE transcripts varied (Figures 3A and 4A; Supplementary  Table S3). Within the Symbiodinium transcriptome assembly, three transcripts encoding prokaryotic-like NiSOD and one transcript encoding both ubiquitin and NiSOD domains were also identified. Genes containing ubiquitin/NiSOD and NiSOD domains have also been identified in the antioxidant gene repertoire of Symbiodinium sp. CassKB8 and Symbiodinium sp. Mf1.05b (Bayer et al., 2012) and two C1 type Symbiodinium (Levin et al., 2016) and may represent genes acquired from prokaryotes by lateral gene transfer. Within Symbiodinium sp. CassKB8 and Symbiodinium sp. Mf1.05b a high number of Trx domain containing genes [106 and 73 Trx genes, respectively (PF00085.14)] were identified (Bayer et al., 2012), within this transcriptome assembly 85 Trx domain containing genes were identified with 27 Trx genes differentially expressed under thermal stress (Figure 3A; Supplementary Table S3). Trx superfamily proteins (Trx and Trx-like proteins) have roles in the oxidative stress response by regulating the redox state, aid in the repair of damaged photosystems (Nishiyama et al., 2011) and regulate many photosynthetic enzymes in plants (including Calvin cycle enzymes such as glyceraldehyde 3phosphate dehydrogenase, phosphoribulokinase, and RuBisCo activase; Hisabori et al., 2007). Differential expression of components of the antioxidant network implicated in protecting cells from ROS and regulating photosynthetic processes were detected here in thermally stressed Symbiodinium. Given that we did not observe a loss in photosynthetic ability of cells, but did observe a shift in metabolism, and don't know the mechanism of regulation of RuBsiCo form II in Symbiodinium the effect of cell redox state on photosynthesis and CO 2 fixation via the Calvin cycle under thermal stress requires further investigation.

Cell Cycle in Thermally Stressed Symbiodinium
The life cycle of Symbiodinium has predominately been considered asexual, deduced by studying morphological transitions and direct observations of mitotic growth. During the vegetative growth phase, haploid cells undergo a diel cycle of mitosis (Fitt and Trench, 1983;Santos and Coffroth, 2003). Progression of the cell cycle in Symbiodinium has been shown to halt when fatty acid syntheses were inhibited by addition of cerulenin [interpreted from decreased free fatty acid and phosphatidylethanolamine (PE) content] (Wang et al., 2013). In prokaryotic and eukaryotic cells, PEs are structural components of membranes and de novo synthesis of PE occurs via the CDP-ethanolamine pathway a branch of the Kennedy pathway (Gibellini and Smith, 2010). Mutant Escherichia coli cells deficient in PE due to defective CDP-ethanolamine pathway genes (pss and psd), showed reduced transcription of flagella genes, resulting in decreased motility and chemotaxis compared to wild type cells (Shi et al., 1993). In this study, components of the CDP-ethanolamine pathway were differentially expressed in thermally stressed cells, one CDP-diacylglycerol-serine O-phosphatidyltransferase transcript (encoded by pss) displayed significantly decreased expression on day 19 (comp44746_c0) and two phosphatidylserine decarboxylase transcripts (encoded by psd) were significantly decreased (comp151799_c0 on day 4 and comp44050_c0 on days 19 and 28; Figure 4B; Supplementary Table S3). Therefore, in thermally stressed Symbiodinium, decreased expression of CDP-ethanolamine pathway genes for PE synthesis could impact the cell cycle due to reduced glycerophospholipid content available for cellular processes.
Mitotic kinases implicated in cell cycle regulation and DNA damage responses were identified with differential expression in thermally stressed Symbiodinium. For example, aurora kinases regulate cell proliferation by controlling M-phase events, such as mitotic spindle attachment (Aurora-A; Nigg, 2001) and NimA-related kinases (Neks) have roles in cell cycle control regulating establishment of the mitotic spindle, chromosome condensation, response to DNA damage, and flagella/cilia development (Nigg, 2001;Fry et al., 2012). Additionally cell division control protein homologs and cyclin-dependent kinases were also detected with differential expression in thermally stressed cells (Figure 4B; Supplementary Table S3). Expression of these cell cycle regulators were not consistent across the gene families identified. Of the nine differentially expressed transcripts encoding nek1 genes in thermally stressed cells, four contigs were significantly increased and three were significantly decreased on day 28 (Figure 4B; Supplementary Table S3). Therefore, although a large number of mitotic kinases exhibited differential expression, linking the expression of these cell cycle regulatory proteins in thermally stressed cells to the observed physiological response is not feasible with the current data.
Documentation of Symbiodinium sexual recombination has been difficult, as cytological evidence of karyogamy has not been found (Chi et al., 2014). However, investigation of population genetic patterns and advances in genomic data has potentially identified a collection of meiotic genes in Symbiodinium (Chi et al., 2014;Wilkinson et al., 2015;Levin et al., 2016). In this study, we identified a number of these meiosis-specific and meiosis-related transcripts (Figure 4C; Supplementary  Table S3), providing further support for the occurrence of sexual recombination in Symbiodinium. Additionally the occurrence of sexual recombination may be restricted to specific conditions, i.e., to free-living symbionts, or be inactivated by symbiotic conditions (Chi et al., 2014), reducing the opportunities for cytological observation. Recently, analysis of the transcriptional response of two type C1 Symbiodinium to heat stress (32 • C), identified upregulation of two mutS homolog genes (Msh4 and Msh5) and no significant changes in other Msh genes (Msh1, 2, 3, and 6), suggesting these genes are essential to meiosis lending support to adaption (Levin et al., 2016). In this study, two transcripts annotated as Msh4 were significantly downregulated on day 19 (comp151677_c0 and comp191937_c0), one transcript annotated as Msh2 exhibited significant upregulation on days 4 and 19 (comp79008_c0), whereas, three transcripts annotated as Msh5 (comp169914_c0, comp201605_c0, and comp32251_c0) and an Msh6 transcript (comp43641_c0) displayed no change in expression (Supplementary Tables S1 and S2). In this study exposure of Symbiodinium to thermal stress resulted in reduced growth rate potentially inducing a cell cycle phase conducive to meiotic-like division. However, differential expression of these genes may also be for repair of DNA damage or the stabilization of DNA in thermally stressed Symbiodinium.
Genes encoding RNA binding proteins containing the RNA recognition motif (RRM) were detected with differential expression. In thermally stressed Symbiodinium cells, RRM containing genes were detected with similarity to the MEI2 gene identified in S. pombe and the MEI2-like gene families identified in A. thaliana and O. sativa. In S. pombe, accumulation and localization of the MEI2 protein to meiRNA results in pre-meiotic DNA synthesis and entry into meiosis I Jeffares et al., 2004). MEI2-like genes have been identified in various eukaryotes including Chlamydomonas reinhardtii, with single copies identified in ascomycete fungi, alveolates, and entamoebidae and gene families identified in plants, however, they are not found in metazoans Jeffares et al., 2004). Analysis of conserved orthologs between clade C3 and A1 Symbiodinium identified an RNA binding protein MEI2 homolog with hits to the alveolate Plasmodium falciparum genome (Voolstra et al., 2009b). Recently, two MEI2like proteins (MEI2-like 2 and MEI2-like 4) were identified with differential expression between species in the comparison of clade B Symbiodinium species (Parkinson et al., 2016). The MEI2-like genes share conserved RRM domains with MEI2 genes Jeffares et al., 2004), though functionally are believed to be involved in cell development during vegetative growth in A. thaliana as well as regulating meiosis (Kaur et al., 2006). RNA binding proteins have also been implicated in chromatin organization and remodeling (Kaur et al., 2006) and MEI2-like gene knockout in Plasmodium yoelii prevents a posttranscriptional regulatory mechanism inhibiting liver schizont stage maturation (Dankwa et al., 2016). The function of these RNA binding proteins in thermally stressed Symbiodinium is unclear, but with 52 transcripts displaying differential expression further studies are needed.
The molecular response to stress can vary with many pathways in place to minimize damage and re-establish cellular homeostasis. The type and degree of stress can ultimately alter the fate of the cell with various cellular functions targeted during stress responses such as cell cycle control, molecular chaperoning, protein repair, protein degradation, and DNA repair, however, if cellular function cannot be regained cell death (apoptosis) may occur (Kültz, 2005). Here within Symbiodinium exposed to thermal stress DEGs encoding RAD DNA repair proteins, DNA photolyases, UPP components, molecular chaperones (HSPs and DNAJ), pro-apoptosis (metacaspases and apoptosis-inducing factors) and anti-apoptosis [IAPs and Bax1 (protein lifeguard) were differentially expressed; Figure 3B; Supplementary Table S3]. Previous studies of Symbiodinium have identified HSPs (HSP70 and HSP90) and HRPs (Leggat et al., 2007;Rosic et al., 2014) and quantified their expression in response to stress (Leggat et al., 2011a;Rosic et al., 2011;Barshis et al., 2014), however, under the various conditions used expression patterns of these molecular chaperones varied. Exposure of Symbiodinium to thermal stress, elicited components of stress response pathways detected here using gene expression analysis (Figures 3A,B; Supplementary  Table S3), which were reflected by decreased cell growth within the physiological parameters measured. By understanding the roles of molecular chaperones in maintaining cellular functions including protein folding and the effect of the proteasomal repair and degradation pathways we may improve our understanding of the stress response of Symbiodinium cells.

Photosynthesis in Thermally Stressed Symbiodinium
Understanding of the photosynthetic machinery gene homologs and their organization has recently advanced with sequencing of the Symbiodinium chloroplast and nuclear genomes (Shoguchi et al., 2013;Barbrook et al., 2014;Mungpakdee et al., 2014). Previously, characterization of Symbiodinium photosystem subunit proteins and genes has examined PS II core proteins psbA (D1 protein) and psbD (D2 protein), PS II manganesestabilizing protein (or PsbO protein) encoded by psbO, and PS I core protein psaA (P 700 protein; Iglesias-Prieto and Trench, 1997;Warner et al., 1999;Iida et al., 2008;Takahashi et al., 2008;McGinley et al., 2012;Castillo-Medina et al., 2013;Gierz et al., 2016). Photoinhibition of PS II and decreased expression of PS II D1 protein (D1 protein content and psbA gene expression) have been observed in thermally stressed Symbiodinium cells (Takahashi et al., 2008;McGinley et al., 2012;Gierz et al., 2016). In this study, exposure of Symbiodinium to thermal stress resulted in a slight decrease in dark-adapted yield ( Figure 1C), a measurement of the PSII activity, though PS II core proteins encoded by psbA and psbD showed no significant changes under these conditions. However, transcripts encoding various subunits of photosystem II (PS II), the cytochrome b 6 f complex, photosystem I (PS I), ATP synthase, cytochrome c 6 , phycocyanin beta, ferredoxin (FRX), and ferredoxin-NADP (+) reductase (FNR) were detected with differential expression under thermal stress conditions (Figure 4A; Supplementary Table S3). Notably four transcripts encoding two extrinsic proteins of the PS II oxygen-evolving complex (psbO and psbP) were significantly upregulated in thermally stressed cells (Figure 4A; Supplementary Table S3). PsbO homologs are found in higher plants, green algae, red algae, diatoms and cyanobacteria and stabilize the Mn cluster and may be critical for the recruitment and assembly of PS II (Ifuku and Noguchi, 2016). PsbP homologs have been annotated in plants, green algae, and cyanobacteria, have high calcium ion binding affinity and may aid in stabilizing the PS II-LHC II supercomplexes in higher plants (Ifuku et al., 2011;Ifuku and Noguchi, 2016).
Biosynthesis of photosynthetic complexes is a controlled process relying on synthesis, insertion, and coordination of each subunit for successful assembly (Rochaix, 2011). Assembly of PS I in C. reinhardtii relies on the insertion of the scaffold anchor protein PsaB, followed by PsaA forming the chlorophyll a-protein complex I and finally by PsaC, after which the other subunits are incorporated (Rochaix, 2011). Photosynthetic ATP synthesis relies on the generation of a proton gradient, either by linear electron flow from PS II to PS I (producing reduced NADP) through the cytochrome b 6 f complex or cyclic electron flow (CEF) through PS I and the cytochrome b 6 f complex (Rochaix, 2011). Additionally, CEF around PSI in Symbiodinium has been demonstrated to increase under short moderate heat stress, proposed to alleviate photoinhibition by dissipating excess energy (qE; Aihara et al., 2016). In this study, expression of the PS I P 700 gene (psaA) varied over the experiment but significantly decreased on day 28 ( Figure 4A; Supplementary  Table S3). Previous targeted studies of psaA across multiple strains of Symbiodinium also found thermal stress resulted in decreased gene expression (McGinley et al., 2012). This change in psaA expression could therefore impair assembly of new PS I complexes in Symbiodinium exposed to thermal stress, potentially disrupting photoprotection via CEF and the synthesis of ATP and reduction of NADP, which are required for cellular metabolic processes including carbon fixation via the Calvin cycle.
In Symbiodinium, the integral LHC family has been studied due to their proposed functional roles of enhancing light capture and photo-protection by dissipating excess energy (Iglesias-Prieto et al., 1993;Takahashi et al., 2008). Intrinsic LHCs have therefore been implicated in the stress response of Symbiodinium cells (Takahashi et al., 2008;Maruyama et al., 2015;Xiang et al., 2015;Gierz et al., 2016). High diversification of the integral light-harvesting gene family (acpPC) in Symbiodinium has been shown following analysis of the S. minutum genome (Maruyama et al., 2015) and analysis of the Symbiodinium sp. C3 acpPC gene repertoire (Boldt et al., 2012). The expression of acpPCs have been characterized in Symbiodinium between two strains of varied thermal tolerance (Takahashi et al., 2008), within a coral host under thermal stress (Gierz et al., 2016) and under light stress (Xiang et al., 2015). In this study, we identified 54 differentially expressed transcripts encoding acpPCs in thermally stressed cells (Supplementary Table S3), of these LHC genes 14 were significantly expressed at every time point throughout the experiment (Figure 4B). The functional associations of Symbiodinium LHCs cannot be determined from this study, though the varied expression of acpPC genes following exposure to thermal stress may be indicative of specific gene function.

Fatty Acid Desaturases
Fatty acids have been quantified in Symbiodinium in attempts to link thermal tolerance to thylakoid membrane lipid composition (Tchernov et al., 2004) and to develop lipid biomarkers for stress (Kneeland et al., 2013). Analysis of lipid content has shown that thermally tolerant species possess different ratios of polyunsaturated fatty acid (PUFA) compared to thermally sensitive strains (Tchernov et al., 2004). However, in multiple clades the lipid composition of whole cells versus enriched chloroplast fractions (Díaz-Almeyda et al., 2011) and the lipid profile (Kneeland et al., 2013) has shown that PUFA desaturation cannot be used to estimate thermal sensitivity of Symbiodinium. Fractionation of chloroplasts has shown that between clades lipid composition can differ, and desaturation and isomerization of these fatty acids can alter the melting point of thylakoid membranes, with increased membrane fluidity observed in thermally stressed S. microadriaticum A1 (Díaz-Almeyda et al., 2011). Further studies of Symbiodinium sp. type C1 and subtype D1 identified decreases in the desaturation ratio and in the fatty acid-to-sterol ratio in cells incubated above 30 • C (Kneeland et al., 2013). However, as total fatty acids were saponified, the sources of change in the lipid profile (storage lipids versus membrane lipids) could not be discerned (Kneeland et al., 2013). In this study Symbiodinium lipid content were not quantified, so the lipid profiles cannot be estimated. However, DEGs encoding fatty acid desaturase were detected in thermally stressed Symbiodinium cells (Figure 4B; Supplementary Table S3). Recently, transcriptome analysis of multiple Symbiodinium clade B strains revealed differences in expression of fatty acid metabolism and biosynthesis pathway genes potentially related to membrane composition, energy storage, and varied growth rates between species (Parkinson et al., 2016).
Previously, in clades C and D type Symbiodinium, orthologs of the palmitoyl-monogalactosyldiacylglycerol delta-7 desaturase were annotated in each assembly, with significantly elevated d N /d s along the clade D lineage (Ladner et al., 2012). As described in the clades C and D analysis, these orthologs are very similar to the palmitoyl-monogalactosyldiacylglycerol delta-7 desaturase (ADS3) in A. thaliana which are involved in the desaturation of hexadecatrienoic acid (16:3 7,10,13 ) (Ladner et al., 2012). As mentioned previously, the lipid composition of thylakoids determined that the ratios of fatty acids (C16, C18, and C22) can change under thermal stress influencing thylakoid membrane fluidity (Díaz-Almeyda et al., 2011). In this study, two classes of fatty acid delta desaturases involved in PUFA biosynthesis were detected with differential expression in Symbiodinium exposed to thermal stress ( Figure 4B). Without having quantified the lipid content of fractionated thylakoids in this study we cannot identify if the increased delta-5 desaturase activity and decreased delta-7 desaturase activity were specific for membrane lipids or storage lipids. Though the differential expression of fatty acid desaturase enzymes detected has not previously been reported, and future studies in Symbiodinium may benefit from pairing assays of both fatty acid desaturase expression and lipid content quantification.

Lipid Catabolism in Thermally Stressed Symbiodinium
Analysis of DEGs detected multiple transcripts encoding the four enzymes of the β-oxidation pathway and two enzymes of the glyoxylate cycle ( Figure 4B). DEGs encoding the four enzymes (acyl-CoA dehydrogenase, enoyl-CoA hydratase, HCDH, and β-ketothiolase) of the fatty acid β-oxidation pathway ( Figure 4B; Supplementary Table S3) were detected in thermally stressed Symbiodinium. Additionally enzymes targeted for both mitochondrial and peroxisomal β-oxidation were annotated, in mammalian cells β-oxidation occurs in both the mitochondria and peroxisomes, whereas, in plant cells β-oxidation is restricted to peroxisomes (Poirier et al., 2006). The peroxisomal β-oxidation pathway in mammalian cells is used for very-long-chain fatty acids and their derivatives that are otherwise slowly degraded in the mitochondria (Poirier et al., 2006), and in plants has been linked to the metabolism of fatty acids from membrane lipids supplying acetyl-CoA to the glyoxylate cycle (Cornah and Smith, 2002). As mentioned previously, fatty acid content of Symbiodinium sp. C1 and S. microadriaticum A1 thylakoid membranes were significantly modified under thermal stress (Díaz-Almeyda et al., 2011), it is possible that in this study we are detecting DEGs of the β-oxidation pathway for the modification of membrane lipid content by removing free fatty acids in thermally stressed cells or the use of storage lipids.
In plants, bacteria, protists, and fungi degradation of fatty acids by β-oxidation generates acetyl-CoA, which may then enter the glyoxylate cycle to produce substrates for gluconeogenesis (Cornah and Smith, 2002). Glyoxylate cycle enzymes have been identified in yeast and plants under starvation (Cornah and Smith, 2002), in a number of corals from genome and transcriptome analyses (Meyer et al., 2009;DeSalvo et al., 2010;Kenkel et al., 2013;Shinzato et al., 2014) in dinoflagellates (Karenia brevis and Amphidinium carterae; Bachvaroff et al., 2004;Butterfield et al., 2013), and recently in the S. kawagutii genome (Lin et al., 2015). The glyoxylate cycle in plants and filamentous fungi occurs in glyoxysomes, which have also been identified in Symbiodinium by transmission electron microscopy (Sammarco and Strychar, 2013). It is also possible that the peroxisomal β-oxidation transcripts identified may instead be glyoxysomal-associated enzymes in Symbiodinium. Two enzymes of the glyoxylate cycle, isocitrate lyase (aceA) and malate synthase (aceB; Kornberg and Madsen, 1958) were detected with differential expression in thermally stressed Symbiodinium (Figure 4B; Supplementary Table S3). Three of the four differentially expressed transcripts encoding isocitrate lyase orthologs were significantly upregulated on day 28 (Figure 4B), this enzyme cleaves isocitrate into glyoxylate and succinate (Kornberg and Madsen, 1958). However, the contigs encoding malate synthase were downregulated in thermally stressed cells (Figure 4B; Supplementary Table S3). In gluconeogenesis, the conversion of malate/oxaloacetate to phosphoenolpyruvate is catalyzed by PEPCK (Pilkis and Granner, 1992). In thermally stressed Symbiodinium cells we detected contigs encoding orthologs of Dictyostelium discoideum and Cucumis sativus PEPCKs (Figure 4B; Supplementary Table S3). Of the three PEPCKs annotated all displayed significantly decreased expression on day 4 and the two D. discoideum orthologs were later significantly upregulated on day 28 (Figure 4B; Supplementary Table S3), indicating that the gluconeogenic pathway were in use in thermally stressed Symbiodinium on day 28. This shift to gluconeogenic metabolism could be a mechanism to reduce free fatty acids in cells from membrane remodeling or potentially be indicative of mobilization of fatty acid stores due to inhibition of photosynthesis induced by thermal stress.
Success of the coral-dinoflagellate symbiosis is largely reliant on the exchange of metabolic products (Davy et al., 2012). Though still unclear due to the difficulties of studying the symbiotic relationship, the translocation of photosynthates from symbionts to host cells (including glycerol, glucose, alanine, and organic acids such as citrate, succinate, fumarate, malate, and glycolate) is estimated to account for up to 60% of photosynthetically fixed carbon (Davy et al., 2012). Although the behavior of isolated Symbiodinium can vary from that of cells in symbiosis (Sutton and Hoegh-Guldberg, 1990;Davy et al., 2012), the implications of altered metabolism from exposure to thermal stress may ultimately influence the stability of the symbiotic relationship.
This study provides an overview of the transcriptome response of Symbiodinium exposed to thermal stress and highlights differential expression of key genes. This study is the first of its kind to employ a moderate thermal stress regime, for a period of 28 days reflective of future temperature conditions, and provides an assessment of physiological parameters that are paired with RNA-Seq analysis. Our results provide a basis for further studies as the transcriptome analysis provides documentation of differentially expressed genes in Symbiodinium exposed to thermal stress for an extended time period. Surprisingly, despite small fold changes a large proportion (23,654 genes) of the transcriptome exhibited altered expression. The longitudinal approach used here has also allowed us to identify genes that display consistently altered expression and those that are only transiently differentially expressed. Differential expression of key stress, photosynthesis, metabolism, and cell cycle genes were detected. Differential expression of glyoxylate cycle enzymes reported here, represents the first instance of this in Symbiodinium. The implications for the change in Symbiodinium metabolism under extended thermal stress and the effect this may have on Symbiodinium-host interactions is unknown, though future studies investigating impacts of extended thermal stress should aim to incorporating metabolomics.