RNA-Seq Reveals Extensive Transcriptional Response to Heat Stress in the Stony Coral Galaxea fascicularis

Galaxea fascicularis, a stony coral belonging to family Oculinidae, is widely distributed in Red Sea, the Gulf of Aden and large areas of the Indo-Pacific oceans. So far there is a lack of gene expression knowledge concerning this massive coral. In the present study, G. fascicularis was subjected to heat stress at 32.0 ± 0.5°C in the lab, we found that the density of symbiotic zooxanthellae decreased significantly; meanwhile apparent bleaching and tissue lysing were observed at 10 h and 18 h after heat stress. The transcriptome responses were investigated in the stony coral G. fascicularis during heat bleaching using RNA-seq. A total of 42,028 coral genes were assembled from over 439 million reads. Gene expressions were compared at 10 and 18 h after heat stress. The significantly upregulated genes found in the Control_10h vs. Heat_10h comparison, presented mainly in GO terms related with DNA integration and unfolded protein response; and for the Control_18h vs. Heat_18h comparison, the GO terms include DNA integration. In addition, comparison between groups of Control_10h vs. Heat_10h and Control_18h vs. Heat_18h revealed that 125 genes were significantly upregulated in common between the two groups, whereas 21 genes were significantly downregulated in common, all these differentially expressed genes were found to be involved in stress response, DNA integration and unfolded protein response. Taken together, our results suggest that high temperature could activate the stress response at the early stage, and subsequently induce the bleaching and lysing through DNA integration and unfolded protein response, which are able to disrupt the balance of coral-zooxanthella symbiosis in the stony coral G. fascicularis.


INTRODUCTION
Stony corals play a crucial role in forming the structure of coral reefs, one of the most important marine ecosystems with immense biodiversity in tropical and subtropical oceans (Zhang et al., 2015). They form complex mutualistic symbiosis with unicellular photosynthetic zooxanthellae, allowing them to proliferate even in the oligotrophic conditions. In the coral-zooxanthella symbiosis, stony coral provides carbon dioxide and inorganic nutrients for engaged photosynthetic ooxanthellae, while zooxanthellae supply stony coral with oxygen and organic nutrients (Shinzato et al., 2011). However, in recent years the increased surface seawater temperature due to the global warming and climate change has led to severe coral bleaching events, and the stony corals are being eliminated at unsustainable rates (Swain et al., 2016). During the bleaching process, stony corals either lose their symbiotic zooxanthellae or the zooxanthellae' photosynthetic pigments are degraded (Coles and Brown, 2003). The resulting disruption of the coralzooxanthella symbiosis is accompanied with increased mortality and decreased adaptation of coral (Goreau and Macfarlane, 1990).
The mechanism of coral bleaching has been investigated in earlier studies through the coral-zooxanthella symbiosis. For the coral host, its physiological responses to heat stress were observed, including photoprotective mechanisms, fluorescent pigments changes, molecular chaperone mechanism, apoptosis regulation and redox equilibrium (Hayes and King, 1995;Salih et al., 2000;Lesser, 2004). These responses were thought to be triggered by excess oxygen free radical and resulting oxidative damage in the symbiont during exposure to elevated temperatures, which further induce the coral bleaching (Lesser, 1997;Downs et al., 2002). For symbiotic zooxanthellae, their genotype and density are thought to be related with the susceptibility of stony coral to bleaching. Certain genotypic zooxanthellae with thermotolerant ability can be endowed with the strong resistance of coral host to heat stress (Rowan et al., 1997;Hawkins et al., 2014). For example, Symbiodinium thermophilum sp. nov. has been considered as a thermotolerant symbiotic zooxanthella, which allows the survival of its coral host in the world's warmest sea, the Persian/Arabian Gulf (Hume et al., 2015). In addition, high-density symbiotic zooxanthellae were also considered to be able to increase the susceptibility of stony corals to bleaching (Cunning and Baker, 2012). However, coral bleaching is a systemic reaction of coral-zooxanthella symbiont to heat stress, its complex mechanism is still not clear.
Galaxea fascicularis, a massive stony coral belonging to family Oculinidae, is widely distributed in Red Sea, the Gulf of Aden and large areas of the Indo-Pacific oceans (Veron, 2000). In common with other corals containing zooxanthellae, G. fascicularis is likely stressed by warmer seawater due to global climate change. However, most previous studies have focused on the heat stress of branching corals, subsequently there is a lack of gene expression knowledge concerning massive corals that are known to be less susceptible to bleaching comparing to branching corals (Maor-Landaw and Levy, 2016). In the present study, the transcriptome response in the massive coral G. fascicularis under heat stress was investigated using the next generation sequencing technology (NGS). As a system biology approach, the NGS has been successfully applied to greatly enrich the transcriptome resource of corals (Polato et al., 2011;Barshis et al., 2013;Sun et al., 2013;Shinzato et al., 2014;Vidal-Dupiol et al., 2014;Anderson et al., 2016), and hopefully it will pave a new way to further reveal the molecular mechanism of coral bleaching and environmental adaption. The purposes of this study were (1) to identify the transcripts from coral G. fascicularis, (2) to survey the differentially expressed genes of coral under heat stress, and (3) to investigate the potential physiological changes resulted from differentially expressed genes, so as to better understand the molecular mechanism triggering the thermal bleaching in the stony coral G. fascicularis.

MATERIALS AND METHODS
Coral Sampling, Identification, and Ethical Statement Four G. fascicularis colonies were collected from the fringing reefs in Luhuitou (18 • 12 ′ 45 ′′ N, 109 • 28 ′ 29 ′′ E), Sanya City, Hainan Province, China. Sample collection for this study was approved by the Management Office of Sanya National Coral Reef Nature Reserve (China). Three colonies of G. fascicularis were sectioned into 12 fragments (each was sectioned into four fragments, ∼25 cm 2 /fragment), preparing for the following heat stress treatment. The colonies were maintained in an aquarium containing 25∼26 • C filtered seawater (same as the in situ surface water temperature, SST). Four fluorescent bulbs (Philips T5HO Activiva Active 54 W) were used as light sources and the corals were subjected to a 12 h light/12 h dark cycle for 2 weeks to acclimatize to the aquarium conditions. All corals were allowed to recover from the fragmentation until the onset of the heat treatment experiment.
G. fascicularis is classified into soft (S) and hard (H) types based on nematocyst morphology (Hidaka, 1992), and this morphological characteristics is correlated with the length of the non-coding region between the mitochondrial genes cytb and nad2 (Watanabe et al., 2005). The soft type has a 290 bp indel (Smt-long, SL), while the hard type does not (H-mt-short, HS). The nematocysts from G. fascicularis colonies were examined under microscope, and the non-coding region was amplified using PCR with the primers (188-F2: 5 ′ -TCCTGTAGAATAGGGTATAC-3 ′ ) and (188-R2: 5 ′ -TTTGCCTTTCCGTATCCACCAT-3 ′ ). All colonies used in the present study were confirmed to be SL type.

Heat Stress Treatment
An aquarium containing heated seawater (32.0 ± 0.5 • C) was prepared and maintained using heating rods (Fish Baby TM ). As described above, three colonies of G. fascicularis were sectioned into 12 fragments (each was sectioned into four fragments, ∼25 cm 2 /fragment), now six fragments (two from each colony) were transferred into the heated aquarium, while the other six were kept in the 25.0 ± 0.5 • C seawater aquarium as the control. At 10 and 18 h post heat stress, those fragments were sampled randomly with three biological repetitions (one fragment from each colony). Four groups of samples were collected, including Control_10h, Heat_10h, Control_18h and Heat_18h groups (n = 3). Meanwhile, the fourth coral colony was sampled in the control aquarium at 0 h, which was employed as Control_0h group. Each sample was stored immediately in liquid nitrogen for RNA extraction later.

Symbiotic Zooxanthellae Dynamics
According to the methods described in Drew (1972), the densities of symbiotic zooxanthellae were measured in G. fascicularis at 0, 10, and 18 h after heat stress, respectively. Briefly, three neighboring polyps were cut from the colony and their dimension was calculated using a square calculation paper. The polyps were decalcified with 5% HCl for about 10 min and then the soft tissue was homogenized with lysis solution (150 mM NaCl, 1% Triton X-100, 0.1% SDS, 50 mM Tris-HCl, pH 7.4). The number of zooxanthella in the homogenate was counted using a hemocytometer (Neubauer improved).

The Construction and Deep Sequencing of Transcriptome Libraries
Total RNA was isolated from each coral sample using Trizol reagent (Invitrogen) according to the manufacture's protocol. The total RNA was quantified by Nanodrop 2000 (Thermo Scientific) at 260/280 nm (ratio > 2.0), and its integrity was checked with Agilent 2100 Bioanalyzer (Agilent Technologies). One paired-end fragment library (2 × 150 bp, Control_0h group) and twelve single-end fragment libraries (50 bp, Control_10h, Heat_10h, Control_18h, and Heat_18h groups) were constructed and sequenced on the Illumina Hiseq4000 platform according to the manufacturer's instructions (BGI, Shenzhen, China). The raw sequencing reads had been submitted to NCBI Short Read Archive under the accession number SRP083089.

Transcript Assembly and Differentiation
After the evaluation of sequence quality and the removal of low quality reads, the clean reads in the all thirteen fragment library were used to assembly the transcripts and genes using the Trinity software (http://trinityrnaseq.github.io/; Haas et al., 2013). The coral and zooxanthellae transcripts were further differentiated as described previously with minor modifications (Yuan et al., 2017). Briefly, the assembled transcripts were aligned to the protein database through BLASTX algorithm (E < 0.00001), which was constructed using the protein sequences from coral Acropora digitifera and zooxanthellae Symbiodinium minutum and Symbiodinium kawagutii. The transcripts encoding proteins which share higher sequence similarity with S. minutum or S. kawagutii proteins were regarded as zooxanthellae-derived ones, while the transcripts sharing higher similarity with A. digitifera proteins were referred as the ones from coral G. fascicularis. The coral transcripts and corresponding genes were used in the subsequent reads mapping and function annotation.

The Identification of Differentially Expressed Genes (DEGs)
The assembled coral transcripts served as the reference sequence in the reads mapping. The alignment of all reads in the twelve single-end fragment libraries was performed using TopHat software (http://ccb.jhu.edu/software/tophat/index.shtml), while Cufflinks software (http://cole-trapnell-lab.github.io/cufflinks/) was used to estimate gene expression abundance and identify the differentially expressed genes between two groups (Trapnell et al., 2012). The differential expression analysis was presented graphically using the CummeRbund package included in the Bioconductor project (http://www.bioconductor.org/packages/ release/bioc/html/cummeRbund.html).

Function Annotation of Assembled Coral Transcripts
The possible coding regions and corresponding encoded proteins of assembled coral genes were identified and retrieved using TransDecoder software (http://transdecoder.github. io/). The sequences of these genes and proteins were aligned by local BLASTX and BLASTP programs to the SwissProtdatabases (max_target_seqs = 1). The alignment results were further parsed by Trinotate software (http:// trinotate.github.io/) for function annotation and GO term assignment.
The GO enrichment analysis of differentially expressed gene was implemented by the hypergeometric test with FDR < 0.01. The differentially expressed genes were selected as the test set while all assembled coral genes were taken as the reference set. The significantly overrepresented GO terms were calculated from the test set, and displayed as a network using BiNGO plug-in to Cytoscape software (http://cytoscape.org/; Maere et al., 2005).

Coral Bleaching/Lysing and Density Dynamics of Symbiotic Zooxanthellae during Heat Stress Treatment
Tissue lysing began to appear at 10 h after heat stress treatment, meanwhile cup edges of a few polyps and coenosarcs turned white (Figure 1). At 18 h after heat stress treatment, the coenosarcs became lysed and degraded, and the tentacles were dissolved into dark brown mucus, while no apparent bleaching and lysing were observed in the control group during the whole treatment process (Figure 1). Meanwhile, the density of symbiotic zooxanthellae decreased from 4.20 ± 0.95 × 10 6 cell cm −2 at 0 h to 3.18 ± 0.85 × 10 6 cell cm −2 at 10 h after heat stress in heat stress group (that of control group is 4.14 ± 0.95 × 10 6 cell cm −2 ), and continued declining to 1.41 ± 0.19 × 10 6 cell cm −2 at 18 h (3.80 ± 1.16 × 10 6 cell cm −2 in the control group). The repeated ANOVAs showed that zooxanthellae densities changed significantly with temperature [F (25, 32) = 8.75, p < 0.05] but not over time [F (10, 18) = 3.45, p > 0.05, Supplementary Figure 1].

Construction and Sequencing of Transcriptome Libraries
A total of 50,400,754 paired-end reads with length of 2 × 150 bp was obtained from the Control_0h transcriptome library. Meanwhile, twelve single-end transcriptome libraries were constructed for the four groups, namely Control_10h, Heat_10h, Control_18h, and Heat_18h, with three biological repetitions for each group. To compare the expression of coral genes after heat stress, 12 libraries were sequenced to the saturated level. After the filtering out of low-quality and adaptor sequences, totally 439,876,050 single-end reads were obtained with length of 50 bp. The reads numbers of each single-end transcriptome library were shown in Table 1.

Assembly and Differentiation of the Transcripts
A total of 333,107 assembled transcripts (N50 = 1,557 bp) was obtained. These transcripts were aligned to the BLAST protein database derived from A. digitifera, S. minutum, and S. kawagutii. The encoding proteins of 77,986 transcripts (42,028 genes, transcript N50 = 2,503 bp, gene N50 = 2,313 bp) shared higher homology with those of A. digitifera, were referred as the transcripts of coral G. fascicularis, while 66,509 transcripts (50,028 genes, transcript N50 = 1,787 bp, gene N50 = 1,736 bp) whose encoded protein sharing higher homology with that of S. minutum or S. kawagutii were considered as the transcripts of symbiotic zooxanthellae ( Table 2).

Abundance Estimation of Gene Expression
The number of mapped reads in 12 libraries ranged from 6,672,958 to 16,436,429, and the mapping rates ranged from 18.2 to 44.8% ( Table 1). The constant coefficients of variation of gene expression levels in four groups of samples were relatively small (Supplementary Figure 2), and the PCA plot defined the integral difference of gene expression abundance among four sample groups (Figure 2)

Identification of Differentially Expressed Genes (DEGs)
After library calibration, the expression levels of 42,028 coral genes were compared between the control and heat stress groups. Total of 4,425 DEGs were obtained in the comparisons, the number and information of differently expressed genes are shown in Figure 3 and Supplementary

GO Overrepresentation Analysis of Significantly Upregulated Genes
There were 383 and 394 significantly upregulated genes in the comparison of Control_10h vs. Heat_10h and Control_18h vs. Heat_18h groups, respectively. Moreover, 125 common significantly upregulated genes were observed in these two comparisons (Supplementary Figure 3). Base on the annotated GO terms of all coral genes (Supplementary Table 3), the GO overrepresentation analysis of those genes was completed at multiple GO levels in the Biological Process category.
For the 383 significantly upregulated genes in the comparison of Control_10h vs. Heat_10h groups, there were eight major overrepresented GO terms, including DNA integration  Table 5). FIGURE 1 | Bleaching and lysing in G. fascicularis during heat stress treatment. No apparent bleaching and lysing were observed in control groups (C_10h and C_18h), while in group H_10h (10 h after heat stress treatment), some phenomena of bleaching and lysing initiation began to appear, the cup edges of a few polyps and the coenosarcs turned white (yellow ellipse). As the heat stress continued, in group H_18h (18 h after heat stress treatment), the coenosarcs became lysed and degraded and the tentacles were dissolved into dark brown mucus (red ellipse).
Frontiers in Genetics | www.frontiersin.org Seven overrepresented GO terms in the Biological Process category were observed for the common significantly upregulated genes, including DNA integration (GO:0015074), DNA metabolic process (GO:0006259), response to heat (GO:0009408, Table 3), response to unfolded protein (GO:0006986), protein folding (GO:0006457), response to protein stimulus (GO:0051789) and response to temperature stimulus (GO:0009266) (Figure 5 and Supplementary Table 6).

GO Overrepresentation Analysis of Significantly Downregulated Genes
From the comparisons of Control_10h vs. Heat_10h and Control_18h vs. Heat_18h groups, 186 and 309 significantly downregulated genes were observed, respectively. No GO term was overrepresented significantly for the two significantly downregulated gene lists.
Only one GO term (pantothenate metabolic process, GO:0015939) was overrepresented for significantly downregulated genes in the comparison of Control_10h vs. Heat_10h groups (Supplementary Table 7). For the 309 significantly downregulated genes in the comparison of Control_18h vs. Heat_18h, 6 overrepresented GO terms were observed including positive regulation of epithelial cell migration (GO:0010634), positive regulation of epithelial to There were 21 common significantly downregulated genes in the two comparisons (Supplementary Figure 4), and only one overrepresented GO term (pantothenate metabolic process, GO:0015939) for these common significantly downregulated genes was observed (Supplementary Table 9).

DISCUSSION
Coral bleaching can cause corals to lose their vital symbiotic zooxanthellae, therefore threatening coral reefs worldwide. To assess the potential response of the stony coral G. fascicularis to high temperature, its phenotypes were observed and the densities of symbiotic zooxanthellae were determined at 10 and 18 h after heat stress. In the present study, the density of symbiotic zooxanthellae decreased significantly during 10∼18 h period after heat stress, and reached the lowest level after 18 h. This demonstrates that an acute heat stress (32.0 • C) could induce a significant decrease of symbiotic zooxanthellae in the stony coral G. fascicularis. The similar decrease of symbiotic zooxanthellae in the stony coral G. fascicularis under heat stress was also reported in a previous study (Bhagooli and Hidaka, 2004). Furthermore, the result was consistent with the present observation that the cup edges of a few polyps and the coenosarcs began bleaching at 10 h after heat stress. Therefore, our results demonstrate that heat stress could induce the collapse of coral-zooxanthella symbiosis  and the expulsion of symbiotic zooxanthellae, resulting in the bleaching of the stony coral G. fascicularis. The reason for the symbiosis collapsing after heat stress could also attribute to the changes in the internal environment of the stony coral, such as the increase of reactive oxygen species (ROS) level (Lesser, 1997), the activation of cell apoptosis and necrosis, etc. (Tchernov et al., 2011). The continued decrease in symbiotic zooxanthellae after acute heat stress could cause coral bleaching, and even cell death in the stony coral G. fascicularis.
To systematically explore the molecular mechanism of coral bleaching, the transcriptome libraries of G. fascicularis were sequenced and analyzed. In the present study, the transcriptome of the stony coral G. fascicularis was first assembled and 333,107 transcripts were obtained. To distinguish these coral transcripts, homologous search was employed using BLASTx method, which was also used in other studies on stony coral symbionts (Vidal-Dupiol et al., 2013;Shinzato et al., 2014). There were 77,986 transcripts from 42,028 genes whose encoding protein sharing higher homology with A. digitifera protein, indicating that the majority of these transcripts were encoded by the G. fascicularis genome and could be used in the subsequent analysis of gene expressions. Then, 569 and 703 DEGs were observed in the comparisons of Heat_10h vs. Control_10h and Heat_18h vs. Control_18h groups, respectively. Our data showed that the transcriptome of the stony coral G. fascicularis could be altered under heat stress. RNA-seq technology has been widely employed to reveal the mechanism underlying heat stress and bleaching in other coral species, such as Acropora hyacinthus (Seneca and Palumbi, 2015), Acropora millepora (Kaniewska et al., 2015), Orbicella faveolata (Pinzón et al., 2015), and Pocillopora damicornis (Vidal-Dupiol et al., 2014), and identified a number of DEGs.
In this study, three biological replicates in each group were analyzed and the DEG analysis showed the consistent results. In addition, those heat stress-induced GO terms found in other scleractinian corals [e.g., protein folding (GO:0006457) (Meyer et al., 2011), response to heat (GO:0009408) Louis et al., 2017], were also overrepresented in our study. Therefore, the present RNA-seq result should be reliable, yet need to be further validated by qRT-PCR in future studies.
To understand the early systemic response of the stony coral G. fascicularis to heat stress, the GO overrepresentation analysis of these differentially expressed genes was conducted in the Control_10h vs. Heat_10h comparison. Eight GO terms were found overrepresented, which were related with unfolded protein response and DNA integration. The unfolded protein response was also observed in the response to darkness in Acropora palmata and Montastraea faveolata (DeSalvo et al., 2011). The unfolded protein response, a cellular stress response related to the endoplasmic reticulum, was found to be conserved amongst most eukaryotes (Walter and Ron, 2011;Hetz, 2012). The unfolded protein response might be activated by misfolded proteins in the stony coral G. fascicularis under heat stress, which could  result from excessive oxygen free radicals observed in other corals during exposure to elevated temperatures (Lesser, 1997;Downs et al., 2002). This physiological process could attenuate the negative effect of heat stress on the stony coral G. fascicularis. However, it might also lead to cell apoptosis and death if the heat disruption was prolonged (Skalka and Katz, 2005;Ron and Walter, 2007). Intriguingly, only one overrepresented GO term was observed for significantly downregulated genes at 10 h after heat stress, which was related with pantothenate metabolic process. Pantothenate, as a kind of vitamin, is used mainly in the synthesis of coenzyme A, which is involved in signal transduction, enzyme activation and deactivation through acylation and acetylation. Deficiency of pantothenate could result in the disorders of the nervous, gastrointestinal and immune systems (Smith and Song, 1996). We speculated that the stony coral G. fascicularis might be susceptible to pathogenic bacteria under heat stress owing to the downregulation of immune response through pantothenate deficiency, this was also supported by the report that the infection of pathogen Vibrio coralliilyticus would lead to the bleaching and lysis of the coral P. damicornis under high temperature (Ben-Haim et al., 2003). This result suggests that the heat stress could induce DNA integration and unfolded protein response during the early stage to attenuate the negative effects of heat stress and maintain the homeostasis.
The GO overrepresentation analysis of the differentially expressed genes at prolonged 18 h after heat stress was conducted to better understand the developing mechanism of heat bleaching in the stony coral G. fascicularis. The overrepresented GO terms response to heat demonstrated that the heat stress response might function throughout the heat treatment. The unfolded or misfolded proteins owing to heat stress response were refolded and renatured generally via molecular chaperones (such as heat shock proteins), and the unfolded protein response was triggered immediately to degrade denatured proteins (Feder and Hofmann, 1999). The prolonged unfolded protein response could trigger cell apoptosis and death, and has been suggested as a universal molecular mechanism underlying the coral bleaching (Dunn et al., 2007;Tchernov et al., 2011). However, the TNF (tumor necrosis factor) was also thought to activate the cell apoptosis and induce the bleaching of coral A. digitifera (Quistad et al., 2014), indicating that the activator of apoptosis and cell death in coral bleaching is different according to species and stressors. Due to the sustained pantothenate deficiency, the negative regulation of immune response would increase the susceptibility to pathogenic microorganisms in stony coral under heat stress and accelerate the bleaching (Ben-Haim et al., 2003;Pinzón et al., 2015). Therefore, more severe bleaching phenomena appeared in the stony coral G. fascicularis at prolonged 18 h after heat stress (Figure 1). Furthermore, the lysis of coral tissues after heat stress could attribute to the suppression of the differentiation, development and migration epithelial cell in the stony coral G. fascicularis, which was demonstrated by the overrepresented GO terms for the significantly downregulated genes at 18 h after heat stress. Taken together, our results reveal the induction of unfolded protein response, DNA integration and immune response regulation during the heat bleaching in the stony coral G. fascicularis. In future studies, we will target the main DEGs in those overrepresented GO terms and verify their expression and function during coral bleaching.

ACCESSION CODE
The raw transcriptome read data of G. fascicularis have been deposited into NCBI Short Read Archive (SRA) under accession number SRP083089.

AUTHOR CONTRIBUTIONS
YWa: Conceived and designed the research, and wrote the paper; ZZ and TX: Performed the transcriptomic analysis, interpreted the data, and wrote the paper together; YWa, TX, and DS: Conducted the sampling; JH, DS, TX, YWu, LC, and JW: Completed the experiments.

ACKNOWLEDGMENTS
The authors were grateful to the staffs, Chuanlian Wu, Jingjing Tian and Yaxing Liu from Management Office of Sanya National Coral Reef Nature Reserve for coral sampling, and to all the laboratory members for continuous technical advice and helpful discussion.