Crassostrea hongkongensis specimens were obtained from an oyster farm in Zhanjiang, Guangdong Province, China. The healthy adult oysters which will quickly close and spray seawater upon stimulation with stirring the seawater were then aerated in sand-filtered natural seawater in a 500-L tank at 23° ± 2°C for 2 weeks prior to the experiment. The oysters were fed twice daily with microalgae, Tetraselmis suecica and Isochrysis galbana. For thermal stress treatment, 100 oysters were subjected to heat shock for 1 h in a 100-L tank with pre-warmed (37°C) seawater followed by transfer to ambient seawater (23° ± 2°C). The seawater is static during the treatment with salinity at 20 ppt and pH at 7.2. Previous studies indicate that 1 h thermal exposure is stressful to this species (Lang et al., 2009). The control group was maintained in ambient seawater (23° ± 2°C) throughout the study. After 24 h, 10 oysters were randomly sampled from thermal-stressed and control groups.

Gill tissue of 10 oysters from each of the two groups was sampled in this study. We chose to study gill tissue as it is a large organ with high surface area exposed to the environment and has been shown to be perturbed by thermal stress (Lang et al., 2009). Total RNA was extracted from gill suspensions using 1 mL TRIzol Reagent (Invitrogen, Carlsbad, CA) per 100 mg of tissue following the manufacturer's instructions. RNA concentrations were measured using a Nanodrop 2000c Spectrophotometer (Thermo Scientific, USA) and the quality of the RNA was analyzed by Bioanalyzer 2100.

DGE Libraries were prepared using the Illumina Digital Gene Expression Tag Profiling Kit, according to the manufacturer's protocol. Briefly, the total RNA samples were treated with deoxyribonuclease I (DNase I) to remove any possible DNA contamination. Then 8 μg of total RNA was treated with Oligo (dT) magnetic beads to purify mRNA. Purified mRNA was mixed with fragmentation buffer, to yield short fragments. The first strand cDNA was synthesized from the cleaved RNA fragments using random primers. RNase H and DNA polymerase I were subsequently added to remove the RNA fragments and synthesize the second strand of cDNA, respectively. The cDNA fragments were then end repaired, followed by 3′ single nucleotide A (adenine) addition and ligation of adapters. Finally, products were purified and amplified by PCR to create the sequencing libraries. The libraries were quantified using Agilent 2100 Bioanalyzer and ABI Step One Plus Real-Time PCR System and then sequenced by Illumina HiSeq 2000 sequencer at the Beijing Genomics Institute (BGI, Shenzhen, China). The sequencing data obtained were deposited in NCBI Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE85707, accession number: GSE85707).

The raw image outputs from the sequencer were transformed into base calls using the Illunina pipeline. Prior to mapping the DGE tags, empty tags (no tag sequence between the adaptors), adaptors, low quality tags (tags with unknown sequences “N”), tags with a copy number of 1 (likely sequencing error), and tags that were too long or too short were excluded from the analysis, leaving behind clean tags of 21 nt in length. The clean tags were mapped to the transcriptome of C. hongkongensis sequence (Tong et al., 2015) using Bowtie (Langmead et al., 2009). Only tags with a perfect match or 1 nt mismatch were further considered and annotated. To monitor mapping events on both strands, tags mapping to both sense and antisense sequences were included in the data analysis. All tags mapped to the reference sequences were filtered, and the remaining clean tags were designated as unambiguous clean tags. Gene expression was normalized using the Reads Per kb per Million (RPKM) method (Mortazavi et al., 2008).

To identify DEGs, the number of raw clean tags in each library was normalized to RPKM. To examine expression differences between the thermal treatment group and the control group, the RPKM of the two groups were compared. Identification of the DEGs were performed following a rigorous algorithm described previously (Benjamini and Yekutieli, 2001). False discovery rate (FDR) ≤ 0.001 and log2 fold change ≥ 1 were used as cut-off to determine significant differentially expressed genes (DEGs). All DEGs were further subjected to Gene Ontology (GO) and KEGG Orthology (KO) enrichment analysis by mapping the DEGs to the GO (http://www.geneontology.org/) and KEGG (http://www.genome.jp/kegg/) databases followed by hypergeometric tests. A 5% significance level following the Bonferroni correction was used as the threshold for considering enriched terms for DEGs.

The oyster samples from another independent treatment groups were used for quantitative real-time PCR (qRT-PCR) assays. 5 oysters were randomly sampled for each group. The qRT-PCR analysis was performed on selected DEGs using Lightcycler480 (Roche, Basel, Switzerland). All reactions were performed using a SYBR Green MasterMix (SYBR Premix EX Taq™, TaKaRa) according to the manufacturer's instructions. EF1α gene was used as the internal control. Primers used for qRT-PCR are listed in Table 1. PCR conditions were as follows: initial denaturation at 95°C for 3 min; followed by 45 amplification cycles of 95°C for 15 s, 55°C for 15 s, 72°C for 10 s; 85°C for 20 s for signal collection in each cycle. Melting curve analysis was performed at the end of each PCR reaction to confirm specificity of the amplicons. Reactions were performed in triplicate, and the relative expression levels were calculated using the 2^−ΔΔCT method (Livak and Schmittgen, 2001). Significance of the difference in the expression levels of selected genes between the control and treatment were estimated using Student's t-test.

DGE analysis is a powerful tool to identify and quantify gene expression at the whole genome level and could provide comprehensive transcript data for studying differentially expressed genes (Mardis, 2008). Recently, DGE techniques have been successfully utilized to study the molecular mechanism of thermal adaptation in Chlorostoma funebralis, Cryptolaemus montrouzieri, and Delia antiqua (Gleason and Burton, 2015; Zhang Y. H. et al., 2015; Hao et al., 2016). To identify genes involved in thermal stress tolerance in the Hong Kong oyster, a cDNA library was constructed from the gill tissues of 10 pooled oysters subjected to a 1 h heat shock in pre-warmed (37°C) seawater. The control library was also prepared using gills of 10 pooled oysters from a group maintained in ambient seawater (23°C). A total of 5,908,582 and 6,178,057 reads were generated from the thermal stress and the control libraries, respectively (Table 2). Following adapter trimming and clean-up, 4,387,274 (74.25% of total reads) and 4,635,084 (75.02% of total reads) reads could be successfully mapped to the reference transcriptome for the treatment and control, respectively (Figure S1). Among these, 3,400,605 (thermal stress) and 3,691,555 (control) reads matched perfectly to the reference.

Sequencing saturation data can reflect the sequencing quality to some extent. Saturation analysis of the clean reads showed that the genes that the clean tags were mapped to were saturated when the tag count approached 3 million (Figure S2). Therefore, the sequencing depth used in our study is sufficient to obtain appropriate transcriptome coverage. Additionally, the unigene coverage, defined by the number of clean reads aligned to a particular reference transcript, was analyzed for the transcriptome reference for each sample. As showed in Figure S3, the distributions of the unigene coverage obtained in the two libraries are more or less uniform, with a >50% coverage obtained for 72 and 76% of the total number of genes for thermal stress and control groups, respectively. This implies that the DGE technology is sensitive enough to detect transcriptomic changes in our model species and can be exploited for gene expression studies at the genome wide scale.

To analyze the gene expression changes associated with thermal stress in oysters, a rigorous algorithm was utilized to identify DEGs from the normalized DGE datasets using pairwise comparison between the two groups (thermal stress and control). To obtain the authentic DEGs, 10 pooled oysters were used to reduce the variation among individuals. DGE analysis of the two libraries detected transcripts of 20,824 and 20,751 genes under thermal stress and control conditions, respectively (Table S1). The transcripts showing at least 2-fold change (FDR≤0.001 and |log2-Ratio (thermal stress/control)|≥1) in the thermal stress library compared to the control library were considered to be differentially regulated. Using this cut-off criterion, a total of 605 genes were found to be differentially regulated (DEGs) in oysters exposed to thermal stress in comparison to the controls (Table S2). Among these DEGs, 378 genes were up-regulated and 227 were down-regulated in response to thermal stress (Figure 1, Table 3).

To gain a better insight into the possible mechanisms of thermal stress response all the 605 differentially expressed genes were mapped to the terms in GO database. The main GO terms included “cellular component,” “molecular function,” and “biological process.” In this study, 110 differentially expressed genes were mapped to a GO ID and 9 functional groups were found to be significantly enriched (Bonferroni-corrected P ≤ 0.05; Table S4) (Figure 2). The genes in the enriched functional groups participated in various metabolic and cellular functions including endoplasmic reticulum membrane (GO:0005789), nuclear outer membrane-endoplasmic reticulum membrane network (GO:0042175), endoplasmic reticulum (GO:0005783), endoplasmic reticulum part (GO:0044432), purine nucleoside binding (GO:0001883), ribonucleoside binding (GO:0032549), purine ribonucleoside binding (GO:0032550), nucleoside binding (GO:0001882) and small molecule binding (GO:0036094). These functional groups were thus inferred to participate in thermal stress response of the oyster.

To further investigate the function of differentially expressed genes during thermal stress, KEGG based pathway analysis of the DEGs were performed. The DEGs were mapped to annotated genes in the KEGG database to identify significantly enriched genes involved in metabolic and/or signal transduction pathways. A BLAST search against the KEGG database indicated that 605 DEGs were assigned to 150 KEGG pathways, and 18 of these pathways were significantly enriched with a cut-off Q ≤ 0.05 (Table S4). A variety of signaling pathways including protein processing in endoplasmic reticulum, apoptosis, ubiquitin mediated proteolysis, endocytosis, spliceosome, RIG-I-like receptor signaling pathway and MAPK signaling pathway were found to be significantly affected (Figure 3). Out of the 18 significantly enriched pathways, protein processing in endoplasmic reticulum was the most affected metabolic pathway, with 40 unigenes mapping to the genes of this pathway (Figure 4). Amongst these 40 unigenes, 25 (mapped to 16 genes) are the main players of unfolded protein response (UPR) and ER-associated degradation (ERAD) processes. After thermal stress, several major components of the UPR and ERAD process were up-regulated, including BiP, PERK, Derlin, HSP40, HSP70, NEF, sHSF, and Prg1. Misfolded proteins resulting from thermal stress are retained within the ER lumen in complex with the above mentioned molecular chaperones, and then directed for degradation through the proteasome in the ERAD process (Fewell and Brodsky, 2009). Results from our pathway analysis are consistent with this mechanism of misfolded protein removal under thermal stress. Corresponding with our results, previous reports in Oncorhynchus tshawytscha and Stylophora pistillata suggest that the genes involved in protein processing in endoplasmic reticulum are significantly up-regulated under heat stress (Maor-Landaw et al., 2014; Tomalty et al., 2015). Thus, our transcriptomic study indicated that protein processing in endoplasmic reticulum may be involved in thermal stress response in the Hong Kong oysters which are vital for survival of the organism in such unfavorable conditions.

To validate the digital gene expression (DGE) results, 8 up-regulated (Heat shock protein 70, HSPbeta-1, DNA polymerase, C-type lectin, apoptosis inhibitor protein, cAMP-responsive element-binding protein-like 2, Rad17, metallothionein) and 4 down-regulated (PATS1, Alkaline phosphatase, Hemicentin-1, C1q1) genes representing different biological processes were selected for qRT-PCR analysis (Figure 5). The results indicated a strong correlation (R2 = 0.908) between the expression fold changes detected by DGE and qRT-PCR (Figure 6). This finding is comparable with other studies analyzing correlations between RNA-seq and qRT-PCR results (Camarena et al., 2010; Zhang et al., 2013; Zhai et al., 2016). The qRT-PCR validation results confirmed the accuracy and reliability of the expression changes detected by DGE analysis. Thus, the DEGs functional enrichment analysis could be used to make reasonable deductions. We hypothesize that these DGEs could play an important role in the response to thermal stress in the Hong Kong oyster.

Thermal stress is one of the most prominent forms of abiotic stress which increases the susceptibility of oysters to diseases, thereby affecting oyster production in the aquaculture industry (Hofmann and Todgham, 2010; Shiel et al., 2015). Therefore, functions related to thermal stress response form the major focus of this study. Gene expression analysis is a useful method to discover genes involved in the regulation of thermal stress. Techniques for gene expression analysis include northern blots, qRT-PCR, DNA microarray and RNA-seq DGE). Northern blots and qRT-PCR are good for detecting whether a single gene is being expressed. With DNA microarrays, it could only use the available network knowledge for investigation of possible pathways (Wheelan et al., 2008). And the advantages of DGE technique in profiling lower expressed genes and measuring greater dynamic range have made it an attractive alternative to other next-generation sequencing techniques. To generate an information platform for thermal stress response, representative transcripts showing differential expression of 2 fold or more were analyzed under control and thermal stress condition in C. hongkongensis. Functional analysis for these DEGs showed that genes related to protein folding (molecular chaperones or HSP family), antioxidants, immune response, and cytoskeleton was differentially expressed between the thermal treatment and control group respectively.

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.