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²/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 (S-mt-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.

An aquarium containing heated seawater (32.0 ± 0.5°C) was prepared and maintained using heating rods (Fish Baby™). As described above, three colonies of G. fascicularis were sectioned into 12 fragments (each was sectioned into four fragments, ~25 cm²/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.

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

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.

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

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

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⁶ cell cm⁻² at 0 h to 3.18 ± 0.85 × 10⁶ cell cm⁻² at 10 h after heat stress in heat stress group (that of control group is 4.14 ± 0.95 × 10⁶ cell cm⁻²), and continued declining to 1.41 ± 0.19 × 10⁶ cell cm⁻² at 18 h (3.80 ± 1.16 × 10⁶ cell cm⁻² 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].

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.

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

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). The expression levels of all genes are shown in the Supplementary Table 1. A similar number of expressed genes were detected in the four groups: Control_10h (33,799), Heat_10h (33,440), Control_18h (33,248) and Heat_18h (32,829).

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 Table 2. Among six DEG lists, there were 569 and 703 DEGs observed in the comparisons of Control_10h vs. Heat_10h and Control_18h vs. Heat_18h groups, respectively.

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 (GO:0015074), DNA metabolic process (GO:0006259), protein folding (GO:0006457), DNA recombination (GO:0006310), macromolecule metabolic process (GO:0043170), cellular macromolecule metabolic process (GO:0044260), nucleic acid metabolic process (GO:0090304) and transposition (GO:0032196) (Figure 4 and Supplementary Table 4). Only one overrepresented GO term, DNA integration (GO:0015074), was observed for the 394 significantly upregulated genes in the comparison of Control_18h vs. Heat_18h groups (Supplementary Table 5).

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

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 mesenchymal transition (GO:0010718), positive regulation of cell morphogenesis involved in differentiation (GO:0010770), regulation of epithelial to mesenchymal transition (GO:0010717), regulation of epithelial cell migration (GO:0010632) and positive regulation of cell development (GO:0010720) (Figure 6 and Supplementary Table 8).

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

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.