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² 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^® 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^® 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.

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 (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 dark-adapted 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^® Inc., Lincoln, NE, USA).

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.

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^®-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 on-column 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.

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¹ 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 sequence-length and library-size (Li et al., 2015; 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² (Charif and Lobry, 2007), venn diagrams were drawn with the ‘VennDiagram’ package³ (Chen and Boutros, 2011), using R version 3.3.2⁴ (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), molecular function (Supplementary Figure S2), and cellular component (Supplementary 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⁵ (Kolde, 2015).

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.

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 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).

Antioxidant genes important in stress responses were detected with differential expression in thermally stressed cells. Superoxide dismutase (SOD) transcripts displayed increased expression throughout the experiment. Manganese SOD (MnSOD) expression were significantly increased at day 19 (comp78421_c0, 1.166 up, p < 0.03), whereas, copper/zinc SOD (CuZnSOD) were significantly upregulated at day 4 (comp47575_c0, 1.243 up, p < 0.03) and for transcript comp28011_c0 at day 19 (1.187 up, p < 0.02) and day 28 (1.136 up, p < 0.03; Figure 3A; Supplementary Table S3). Additionally, a transcript encoding both ubiquitin and nickel superoxide dismutase (NiSOD) domains (comp24219_c0) were downregulated on days 19 (1.304 down, p < 0.04) and 28 (1.197 down, p < 0.005; Figure 3A; Supplementary Table S3). A further three transcripts encoding NiSODs and two transcripts encoding MnSODs were annotated (Supplementary Table S5), though no significant expression changes were detected under the conditions used. Two transcripts encoding catalase peroxidase (KatG) exhibited significantly increased expressed in thermally stressed cells (comp80428_c0, day 4, 1.215 up, p < 0.001 and comp61565_c0, day 28, 1.110 up, p < 0.03; Figure 3A; Supplementary Table S3). Nine transcripts encoding ascorbate peroxidases (APX), four peroxiredoxin (Prx) genes, 27 thioredoxin (Trx) genes, and 10 glutathione S-transferase (GST) transcripts were detected with significantly different expression throughout the thermal stress (Figure 3A; Supplementary Table S3).

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 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).

Transcripts encoding anti-apoptosis proteins were also differentially expressed in thermally stressed cells. Three transcripts encoding inhibitors of apoptosis (IAP; Baculoviral IAP repeat-containing protein isoforms) were detected with increased expression on day 4 (BIRC2; comp13582_c0, 1.596 up, p < 0.02) and day 19 (BIRC3; comp12814_c0, 1.862 up, p < 0.04 and comp40288_c0, 1.482 up, p < 0.02; Figure 3B; Supplementary Table S3). One further apoptosis inhibitor 1 (IAP1, comp115062_c0) displayed decreased expression on day 19 (1.183 down, p < 0.02) and day 28 (1.197 down, p < 0.03; Figure 3B; Supplementary Table S3). Transcripts encoding suppressors of apoptosis, Bax inhibitor 1 (BI1) family genes, were detected in the Symbiodinium transcriptome including 13 protein lifeguard isoforms (LFG1–LFG4; Figure 3B; Supplementary Table S3). In thermally stressed cells nine BI1 family genes were detected with differential expression including three LFG1, three LFG2, two LFG3, and a BI1-like protein (Figure 3B; Supplementary Table S3). Interesting three transcripts encoding the LFG4 isoform did not exhibit altered expression under the conditions used (Supplementary Table S5).

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,5-bisphosphate 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 peridinin-chlorophyll 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). Fifty-four 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).

Genes encoding enzymes for fatty acid desaturation and fatty acid β-oxidation were detected with differential expression in Symbiodinium cells exposed to thermal stress. Seven delta-fatty acid desaturase transcripts were detected with differential expression in thermally stressed cells (Figure 4B). Five contigs were annotated as delta-5 desaturases (isoforms fad5A, fad5B, and fad5C) and two were annotated as palmitoyl-monogalactosyldiacylglycerol delta-7 desaturase. Throughout the experiment two palmitoyl-monogalactosyldiacylglycerol delta-7 desaturase (fad7) contigs were downregulated (comp72598_c0, day 4, 1.399 down, p < 0.04; comp21243_c0, day 19, 1.565 down, p < 0.01, day 28, 1.210 down, p < 0.01) whereas, three of the delta-5 desaturase contigs were significantly upregulated at days 19 and 28 (Figure 4B; Supplementary Table S3). Representatives of all genes involved in the β-oxidation of fatty acids were differentially expressed over the course of the experiment. Six transcripts encoding acyl-CoA dehydrogenases (ACADs) were detected with differential expression over the course of the experiment, with five displaying significantly increased fold changes on day 28 (Figure 4B). Transcripts encoding enoyl-CoA hydratases (three ECH), fatty acid oxidation complex subunit alpha (one FADJ), and peroxisomal bifunctional enzymes (1 MFEA and 2 ECHP) were all detected with increased expression in thermally stressed cells (Figure 4B; Supplementary Table S3). One transcript, encoding a probable 3-hydroxyacyl-CoA dehydrogenase (HCDH), was detected with decreased expression on day 28 (comp107722_c0, 1.102 down, p < 0.02; Figure 4B). Three transcripts encoding β-ketothiolases (FADA), responsible for the final cleavage step of the β-oxidation pathway 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 sub-family 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). Twenty-two 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).

Differential regulation of cellular component biosynthesis were detected in thermally stressed Symbiodinium cells. Transcripts encoding lipid biosynthetic acyl carrier protein (acp; comp62787_c0), CDP-diacylglycerol-serine O-phosphatidyltransferase (pss; comp44746_c0), and two phosphatidylserine decarboxylase proenzymes (psd) exhibited decreased expression in thermally stressed cells (Figure 4B; Supplementary Table S3). Four magnesium-chelatase subunit H (chlH) transcripts were detected with significantly increased expression over the experiment, whereas, two magnesium-chelatase subunit I (chlI) transcripts exhibited decreased expression on day 28 (comp46484_c0, 1.173 down, p < 0.001 and comp39765_c0, 1.181 down, p < 0.05; Figure 4B; Supplementary Table S3). Additionally, one transcript encoding a magnesium-chelatase subunit D (chlD; comp84001_c0; Supplementary Table S5), displayed no changes in expression under the conditions used here (FDR < 0.05).

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.

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., 2014, 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 3-phosphate 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₂ fixation via the Calvin cycle under thermal stress requires further investigation.

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 (Anderson et al., 2004; 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 (Anderson et al., 2004; 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 MEI2-like 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 (Anderson et al., 2004; 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 post-transcriptional 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.

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 manganese-stabilizing protein (or PsbO protein) encoded by psbO, and PS I core protein psaA (P₇₀₀ 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₆f complex, photosystem I (PS I), ATP synthase, cytochrome c₆, 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₆f complex or cyclic electron flow (CEF) through PS I and the cytochrome b₆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₇₀₀ 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 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.

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.