Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

RNAi based transcriptome suggests genes potentially regulated by HSF1 in the Pacific oyster Crassostrea gigas under thermal stress

BMC Genomics volume 20, Article number: 639 (2019) Cite this article

Abstract

Background

The Pacific oyster Crassostrea gigas is an important fishery resource that is sensitive to temperature fluctuations. Thus, it has evolved a protection mechanism against heat stress by increasing the expression of the gene coding for heat shock protein (HSP) 70 under elevated temperatures. In other animals, heat shock response is a transcriptional response driven by the heat shock transcription factor 1 (HSF1) and thermal stress can trigger HSP70 expression to protect the organism via HSF1. However, the regulatory relationship between HSF1 and HSP remains unclear in Pacific oyster. Therefore, in the present study, we examined the transcriptomic response of several to thermal stress following HSF1 interference.

Results

We identified 150 genes responsive to heat shock including seven HSP genes, six of which belonging to the group of 17 HSP genes enriched in response to heat shock, according to weighted gene co-expression network analysis (WGCNA). The other gene was enriched in the module correlated with HSF1 interference. In addition, we found 48 and 47 genes that were upregulated and downregulated by HSF1 in response to heat shock, respectively. In the upregulated genes, we identified one HSP70 potentially regulated by HSF1 in response to heat shock. Furthermore, based on differentially expressed genes and WGCNA analyses, we found that the hypoxia signaling pathway was enriched under heat shock conditions. Five genes were then selected to detect dynamic changes through time. The results suggested that gene expression was correlated with HSF1 expression. The regulation of HSP70 by HSF1 was preliminarily confirmed by binding site predictions and by a dual luciferase assay.

Conclusions

Our results revealed that the expression of HSP70 and HSP20 was initially triggered after 2 h of heat shock, and one of the HSP70 genes was potentially regulated by HSF1. From these results, it is evident that not all heat-inducible genes were triggered simultaneously in response to heat shock stress. Overall, the results revealed a possible HSF1–HSP regulatory relationship in Pacific oyster, providing valuable information on the mechanisms of thermal tolerance in this commercially important oyster.

Background

Temperature is an important environmental factor affecting the physiological metabolism [1], distribution [2], and growth [3] of organisms. Naturally distributed in the intertidal zone, Pacific oyster (Crassostrea gigas) are prone to stress and seasonal temperature fluctuations [4]. The species has become a valuable model for studying the adaptability of marine mollusks because its reference genome (GCA_000297895.1, GenBank) [4] is available. Due to its ability to adapt to a wide range of environmental conditions and high market value, C. gigas is an aquaculture resource in many countries. Temperature, an important abiotic factor, is related to C. gigas mass mortality during summer as individuals become more susceptible to infection by bacteria and viruses [5, 6], and high temperature hinders their development [7]. Tolerance to thermal stress requires compensation for the metabolic energy used through molecular regulatory processes.

The physiological effects of temperature fluctuations on the growth and metabolism of Pacific oyster have been well studied [8, 9]. Previous studies have focused on heat shock tolerance accumulation, detection of physiological indicators, and analysis of gene expression under different temperature treatments in terms of transcriptome [10,11,12,13,14,15] and proteome [16] changes. Many heat responsive genes were reported in these studies. Notably, HSP genes, which function as molecular chaperones, play an important role in heat shock response. Bivalves such as the Pacific oyster [4], the pearl oyster, Pinctada fucata [17], the golden mussel, Limnoperna fortunei [18], and the Yesso scallop, Patinopecten yessoensis [19] underwent massive HSP70 gene family expansion compared with other animals [11], according to recent genomic and transcriptomic surveys. In addition, HSP70 genes are particularly expanded and highly inducible during heat shock in oysters, suggesting that these genes play important roles in heat and other stress factors adaptation in oysters [20] and other bivalves [17,18,19, 21, 22].

Although heat responsive genes have been widely studied, their regulatory mechanisms are not fully understood. Heat shock transcription factor 1 (HSF1) can form trimers in response to thermal stress and it regulates the expression of heat shock response genes (especially HSP genes) by binding to heat shock elements (HSEs) in mammals [23, 24]. The isoforms of HSF1 in Pacific oyster were described by Kawabe and Yokoyama [25], and studies on the functions of HSF1 in other species provided insights into the molecular interactions that occur during heat shock. For example, the HSF1–HSP pathway has been described in the roundworm Caenorhabditis elegans [26, 27], in the fruit fly Drosophila melanogaster [28,29,30], and in mammals [31, 32]. Thus, we hypothesized that a similar pathway may exist in Pacific oyster.

RNA interference (RNAi) and sequencing are widely used to study the mechanisms of various life processes in C. elegans [27], D. melanogaster [33], and mammal cells [34]. Weighted gene co-expression network analysis (WGCNA), which can help to identify associations between genes and treatments, has also been widely used to identify hub genes of specific biological processes. In the present study, untreated and HSF1-interfered Pacific oyster individuals were subjected to a heat shock treatment at 35 °C for 2 h. In addition, oyster RNAi was used to identify genes potentially regulated by HSF1. Heat shock- responsive transcripts were revealed through transcriptome analysis the pathway in which HSP70 is regulated by HSF1 under heat shock conditions was illustrated. Using a dual-luciferase reporter assay and quantitative real-time PCR (qRT-PCR), we preliminarily validated the regulatory relationship and the expression patterns of enriched genes.

Results

RNA interference experiment

An initial experiment on heat shock and HSF1 interference treatment times showed that HSF1 was significantly highly expressed compared with the control after 2-h heat shock treatment. Therefore, 2 h was considered the adequate time for heat shock treatment in further experiments (Additional file 1: Figure S1) and 48 h was the adequate time to detect RNAi effects of HSF1 (for gene expression levels see Additional file 1: Figure S2).

Results of the sequencing analysis for HSF1 gene expression and protein content of the RNAi oysters with and without heat shock are shown in Fig. 1. qRT-PCR results showed no significant difference in HSF1 mRNA expression between the RNAi group and the RNAi&HS group (p > 0.05), while both groups showed significantly lower expression of HSF1 compared with that of the control group (both p < 0.001) (Fig. 1a). Western blotting was used to detect the expression of HSF1 protein in each group, with β-tubulin as the control (Fig. 1b), and protein content was measured using a gray scale (Fig. 1c). Protein levels were significantly different between the RNAi and the control groups (p < 0.05). Overall, these data indicated that our RNAi treatment was effective.

Fig. 1
figure 1

Effects of HSF1 gene silencing on HSF1 expression in C. gigas gills. The control represents group 1, RNAi&HS represents group 4 (heat shock following RNA interference), RNAi represents group 3 (RNA interference) a HSF1 mRNA expression detected by qRT-PCR (n = 15). b HSF1 protein content detected by western blotting (n = 2). Each band represents the protein content of an individual, with the upper band representing the content of the HSF1 protein and the lower band representing the content of β-tubulin. c Statistical analysis of HSF1 protein content detected by western blotting (n = 2). The error bar denotes the standard error of the mean; **, p < 0.05; *** p < 0.001

Full size image

Transcriptome data analysis

The mRNA of four groups subjected to different treatments [group 1: control; group 2: heat shock treatment (HS); group 3: RNAi; and group 4: RNAi&HS], and each consisting of five randomly assigned individuals with three biological replicates was sequenced in the Illumina HiSeq X Ten platform (Illumina, CA, USA). The high-quality 150 bp paired-end reads of 12 RNA-sequencing (RNA-seq) samples were mapped to the C. gigas genome (quality information for the 12 samples is shown in Additional file 2: Table S1). The high-quality reads of each sample were sequence-aligned with the C. gigas genome, and the efficiency of the alignment varied from 75.18 to 76.51%. Based on the alignment results, 7873 new genes were identified (new genes are listed in Additional file 3: Table S2), of which 6212 were functionally annotated in five databases (Additional file 4: Table S3). We normalized each condition (groups 1–3) to the control in order to determine fold changes in the relative abundance of mRNA. Differentially expressed genes (DEGs) [false discovery rate (FDR) ≤0.05, fold change (FC) ≥1.5] for groups 2, 3, and 4 compared with the control group are designated gene sets 1, 2, and 3, respectively, throughout this manuscript (Additional file 1: Figure S3). A complete list of the significant genes altered in response to each condition, after normalization to the control (150, 173, and 215 genes were identified for gene sets 1, 2, and 3, respectively), is provided in Additional file 5: Table S4 (FDR ≤0.05, FC ≥1.5). The DEGs were annotated using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (gene annotation is presented in Additional file 5: Table S4; annotation statistics are presented in Additional file 6: Table S5). To test whether high temperature and HSF1 levels could affect the expressions of C. gigas housekeeping genes, we examined the changes in the expressions of elongation factor 1α (EF-1α), glyceraldehyde 3-phosphate dehydrogenase (GADPH), β-actin, β-tubulin, 18 s rRNA, and 28 s rRNA among treatments, but detected no significant differences (Additional file 1: Figure S4).

Heat shock-responsive transcripts

According to GO classification, the 150 genes responsive to heat shock in gene set 1, i.e., genes that were differentially expressed between groups 1 and 2 (FDR ≤0.05, FC ≥1.5), were annotated in the following terms: “single-organism process” (GO:0044699, 37 genes) and “cellular process” (GO:0009987, 34 genes) for biological processes; “cell part” (GO:0044464, 25 genes) for cellular components; and “binding” (GO:0005488, 28 genes) and “catalytic activity” (GO:0003824, 21 genes) for molecular functions. Five HSP70 and two HSP20 genes were enriched compared with the control group (Table 1). Pathways that were significantly enriched (p < 0.05) under heat shock conditions are shown in Table 2, and they include apoptosis, protein processing in endoplasmic reticulum, ubiquitin-mediated proteolysis, the hypoxia-inducible factor (HIF)-1 signaling pathway, and the vascular endothelial growth factor (VEGF) signaling pathway.

Table 1 HSP in differentially expressed genes of heat shocked and control groups
Full size table
Table 2 Enriched KEGG pathways in response to heat shock response compared with the control
Full size table

DEGs potentially regulated by HSF1 response to thermal stress

The genes potentially regulated by HSF1 in the early stage of heat shock were unique to gene set 1, after excluding those that overlapped with gene sets 2 and 3 (group details are presented in Additional file 1: Figure S3, S5). The 95 DEGs potentially regulated by HSF1 included 48 upregulated genes (pink area in Fig. 2a) and the most abundant GO terms were the same as those for the heat shock treatment (Fig. 2b). The KEGG annotation of these DEGs (Fig. 2c, d) revealed that protein processing in endoplasmic reticulum, endocytosis, spliceosome, and apoptosis pathways were enriched. One HSP70 gene (CGI_10002387), which was annotated in “posttranslational modification, protein turnover chaperones”, belonged to the DEGs unique to gene set 1.

Fig. 2
figure 2

The DEGs upregulated by HSF1 in response to heat shock (FDR < 0.05). a Venn diagram showing the overlap of three groups (groups 2–4) compared with group 1 (control). The pink area indicates genes that were upregulated by HSF1 in response to heat shock. b Gene Ontology classification of the DEGs upregulated by HSF1 in response to heat shock. The red bars represent biological processes; the green bars represent the cellular component; and the blue bars represent the molecular functions. c Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched for the HSF1-regulated genes, which were dependent on heat shock. d KEGG classification of the HSF1-regulated genes dependent on heat shock. Each row on the left side represents a pathway, and the collection to which the pathway belongs is displayed on the right side

Full size image

Among the 95 DEGs unique to gene set 1, 47 were downregulated (pink area in Fig. 3a), and the most abundant GO terms were the same as those for the heat shock treatment (Fig. 3b). Five genes were annotated in KEGG, of which two were related to environmental information processing (Fig. 3c).

Fig. 3
figure 3

DEGs downregulated by HSF1 in response to heat shock (FDR < 0.05). a Venn diagram showing the overlapping groups (groups 2–4) compared with group 1 (control). The pink area indicates genes that were downregulated by HSF1 in response to heat shock. b Gene Ontology classification of the DEGs downregulated by HSF1 in response to heat shock. c Kyoto Encyclopedia of Genes and Genomes classification of the HSF1 downregulated genes that were dependent on heat shock

Full size image

Co-expression network

To find the genes responsive to HS and RNAi&HS treatments, we identified the co-expressed genes by applying WGCNA [35, 36] across all groups. The 2758 genes that fulfilled the stringent criteria (FPKM ≥1, module similarity threshold ≥0.25) belonged to 14 co-expression modules (14 different colors in Fig. 4) containing from 40 (cyan) to 499 (turquoise) genes (see details in Additional file 7: Table S6). The 14 distinct modules identified (Fig. 4a) are displayed as a single sample type to show the relationships between modules and sample (Fig. 4b). Although the three biological replicates of each condition belonged to different modules, these modules were correlated with each other (Additional file 1: Figure S6; the Eigengene adjacency heatmap shows that the black, purple, and salmon modules were highly correlated with each other, and that the blue module was positively correlated with the salmon module).

Fig. 4
figure 4

Weighted gene coexpression network (WGCNA) analysis of genes in all samples. a Hierarchical cluster tree representing the coexpression modules identified by WGCNA. Each leaf in the tree represents one gene. The major tree branches constitute 14 modules [fragments per kilo base of exon per million reads mapped (FPKM) ≥1, module similarity threshold ≥0.25], labeled with different colors. b Module–sample association. Each row corresponds to a module and it is labeled with a different color as in (a). Each column corresponds to a sample. The color of each cell at the row–column intersection indicates the correlation coefficient between the module and the sample. Each cell represents the correlation value between the module and each sample; the upper number is the correlation value, and the number underneath that is the p-value (significant at p < 0.05)

Full size image

Under the heat shock treatment, 14 HSP70 and three HSP20 were enriched according to WGCNA (Additional file 8: Table S7). Three modules (tan, yellow, and blue) were strongly correlated with heat shock conditions. The tan module was highly correlated to the heat shock treatment sample (HC-1) and the genes were highly correlated with those of the tan module (r = 0.98, p < 0.05). The genes classified by GO analysis as “binding” in the molecular function and as “metabolic process” in the biological process categories were the most abundant. Pathway enrichment was observed for signaling transport, amino acid metabolism, carbohydrate metabolism, xenobiotics biodegradation and metabolism, and other catabolism processes. Notably, the HIF-1 signaling pathway (Ko04066) was included in the environmental information process (p < 0.05). Both the yellow and blue modules were enriched for metabolism processing co-expression and were related to the heat shock treatment samples (HC-2 and HC-3, respectively). For the yellow module, lipid metabolism, ubiquitin-mediated proteolysis, and phagosome formation were the most abundant processes. Furthermore, the tricarboxylic acid (TCA) cycle and amino acid metabolism were the most abundant metabolic processes. The genes in the blue module were annotated in six kinds of processes. The most abundant were involved in genetic information processing such as RNA transport and protein processing.

Genes within group 4 were identified in the black, purple, green, and salmon modules, in which some heat shock protein related genes were included (Table 3). The salmon module was also associated with heat shock (Fig. 4b), and HSP genes were included in this module.

Table 3 HSP genes in salmon, black, purple, and green modules
Full size table

Hub genes were selected by their high eigengene connectivity values (KME > 0.7) and showed the most connections in the network (i.e., the six genes with the highest KME among the four modules; modules related to HS and RNAi&HS are shown in Additional file 9: Table S8). The correlation network of the tan, yellow, blue, and salmon modules and their hub genes are shown in Fig. 5.

Fig. 5
figure 5

Correlation networks and genes with high KME values in the modules. a, c, e, and g represent tan, yellow, blue, and salmon modules, respectively. The larger the circle, the wider the lines between the genes and the higher the degree of connectivity between the gene and the module. b, d, f, and h represent the hub genes of the tan, yellow, blue, and salmon modules, respectively

Full size image

Validation of the key responsive genes by qRT-PCR and dual luciferase assay

Five genes, including three HSP70 (HSP70–2387, HSP70–2823, and HSP70–2594), one HSP20 (HSP20–4164), and one HIF signaling pathway-related (CGI_10022835) genes associated with heat shock, were selected for qRT-PCR analysis (materials treated at 35 °C for 0, 0.25, 0.5, 1, 2, 6, 12, and 24 h from the start of the experiment) and then compared with the RNA-seq data. The expression of the HIF signaling pathway-related gene changed under HSF1 interference after 1 h of heat stress (Fig. 6a). The response of HSP genes was delayed but expression levels were higher than that of HSF1 when continuously subjected to heat shock (Fig. 6a). The RNA-seq data showed that 2 h of heat shock significantly upregulated HSP genes and the expression level of each responsive gene under HSF1 interference within 2 h of heat shock (Fig. 6b). The fold changes observed in these genes revealed that the mRNA expression of HSP70–2387 changed the most when the organism was exposed to heat shock and had the highest interference level compared with that of other genes (Fig. 6b). This gene was one of the enriched genes unique to gene set 1 (after excluding gene sets 2 and 3), indicating that it might be regulated by HSF1. Thus, it was chosen to study the regulatory relationship between HSF1 and HSP70. The predicted HSEs of HSP70–2387 are listed in Fig. 6c, and the details of binding sequencing are listed in Fig. 6d. Dual luciferase assay results are shown in Fig. 6e. The transcription of the different isoforms of the HSF1 gene was significantly different, but both isoforms (HSF1a and HSF1d) were positively activated.

Fig. 6
figure 6

Quantification of enriched genes and validation of the regulatory relationship between HSF1 and HSP70–2387 genes. a Quantification of the expression of six genes at different sampling times under 35 °C heat stress treatment (n = 5). The X axis represents the sampling time from the start of the experiment (ck represents 0 h from the start, and sampling times were 0, 1/4, 1/2, 1. 2, 6, 12, and 24 h). The Y axis on the left represents fold change of the genes presented in the histogram while the right axis represents fold change of the genes presented in the line chart. The different colors represent the fold changes of each gene. b FPKM values and the fold change of enriched genes. Control represents group 1, HS represents group 2, and RNAi&HS represents group 4. The X axis represents different genes. The Y axis on the left represents FPKM of HS (heat shock) and control groups, while the right axis represents HS/control (FPKM of heat shock divided by FPKM of control) and HS/RNAi&HS (FPKM of heat shock divided by FPKM of RNAi&HS). c Sequencing of the HSP70–2387 gene and the binding site of HSF1 predicted by AnimalTFDB 3.0. The red, blue, and black frames represent the predicted heat shock elements. d Relative luciferase activity of the HSF1 to the HSP70–2387 genes. The control had no HSF1 plasmid added to the system, and HSF1a and HSF1d represent the plasmids that can produce proteins containing the HSF1a and HSF1d isoforms of HSF1, respectively. The error bar denotes the standard error of the mean; **, p < 0.05; *** p < 0.001

Full size image

Discussion

Oysters are exposed to fluctuating temperatures that can lead to physiological stress owing to their intertidal habitat. As they can trigger the expression of stress-related genes for protection from heat stress, oysters might have regulatory capabilities [4, 8, 9, 16, 37]. Heat shock response genes were expressed under heat shock treatment at 35 °C for 2 h. Two HSP20 and five HSP70 genes were enriched among the 150 DEGs. The results of the WGCNA analysis suggested that the HSPs were highly associated with the heat shock treatment. Three HSP20 and 14 HSP70 genes were enriched in WGCNA analysis. Although some HSP70 genes can be induced by a specific stress, there is little change in the expression level of heat shock cognates under stress or control treatments. The expression of four HSP70 genes (HSP70–10647, HSP70–2594, HSP70–10646, and HSP70–2823) was 20 times higher under 35 °C than under the control treatment (Additional file 8: Table S7). Four of the five heat-induced genes reported by Zhang [4] were consistent with our results, but the expression of another heat-induced gene reported in Zhang’s study (CGI_10003417) was about five times higher under heat shock treatment than under the control treatment. The four heat-induced genes (Fig. 6a) exhibited increasing mRNA levels with heat shock treatment. Zhu [11] revealed the dynamic expression pattern of thermal stress-responsive candidate genes under heat shock (35 °C), and that the recovery stages trigger heat shock responses. Two of the four genes (HSP70–2594 and HSP70–2823) were also reported by Zhu [11], who considered that these genes were highly induced in response to 35 °C heat shock. All these results suggest that the highly induced genes revealed in our study are credible.

The expression of HSP70–2387 was over 45 times higher under heat shock conditions than under control conditions (Fig. 6b). Heat shock proteins are molecular chaperones that play important roles in protein folding and refolding, stabilization of denatured proteins, and protein transport under heat or other stresses [38,39,40,41]. In Pacific oyster, HSPs protect the organism from heat stress, and the high expression levels of HSP70 genes indicate that they may contribute more to this protective function. These results suggest that the heat-induced genes (mostly HSP70s) may play important roles in response to heat shock conditions. Several reports have shown that HSPs can be triggered by thermal stress. In Pacific oyster, most individuals exposed to 37 °C for 1 h could survive a subsequent heat treatment (43–44 °C, 1 h), and the expression of a 69 kDa HSP was induced by heat shock, suggesting that these proteins play an important role in oyster survival under heat shock [10]. Heat shock proteins and genes involved in cellular homeostasis were the most highly expressed when C. gigas were kept at 25 °C for 24 days and then analyzed by RNA-seq [42]. The study also showed that gills were more sensitive than the mantle, and therefore we choose the gills for transcriptome analysis. Furthermore, Zhang [16] examined the proteome in gills under high temperature (38 °C, 1 h) and found that HSPs, including eight HSP20 genes, three HSP70 genes, and one small HSP, HSP9/12, were enriched. Our findings are consistent with these results, as HSP20 and HSP70 were initially triggered under thermal stress.

Notably, the HIF-signaling pathway and VEGF signaling pathway contained enriched DEGs in the KEGG annotation. One of the DEGs (CGI_10022835) was annotated as “angiopoietin 4”. Angiopoietins are a family of vascular growth factors that play a role in angiogenesis. It was reported that angiogenesis can be triggered during thermal stress in zebrafish [43]. Previous studies have revealed that angiopoietin 4 and VEGF were expressed at higher levels during hypoxia in U87 glioblastoma cells [44]. Angiogenic genes, such as VEGF, are regulated by HIF during hypoxic conditions [45]. Baird [30] showed that HIF1 could control heat shock pathways by regulating HSF1 in adult flies (D. melanogaster). In addition, Kawabe and Yokoyama [25] reported the HIF-HSF pathway in Pacific oyster. According to a previous study, intermittent hypoxia occurs when Pacific oyster is exposed to heat shock [46], and the genes involved in aerobic energy metabolism are downregulated in response to the thermal stress [47]. A possible explanation for this might be that the HIF1 signaling pathway is involved in the heat shock response process. In addition, some processes were enriched in the heat shock module, such as the TCA cycle, lipid metabolism, amino acid metabolism, and signal transduction. These results are supported by those of previous studies where metabolic processes were involved in heat stress responses in Mytilus [48], C. gigas [14, 15], and zebrafish [43]. Under heat shock conditions, these metabolic processes were regulated to compensate for the energy consumed to produce HSPs [48, 49].

Several reports have also shown that RNAi can be an efficient method for studying gene functions [50,51,52]. We chose three combinations of small interfering RNAs for RNAi to avoid activating the interferon system [53,54,55,56]. Overall, knockdown of the HSF1 gene allowed examining gene functions when the oyster was under heat stress. By performing western blotting and qRT-PCR, we found that the levels of the HSF1 protein and mRNA were significantly different compared with the control, suggesting that these materials could be used for HSF1 function analysis. By analyzing the transcriptome data, we found which candidate genes could be regulated by HSF1 under heat stress conditions. The results showed that the most abundant terms enriched in DEGs were also enriched under the heat shock treatment. Thus, HSF1 might regulate the genes related to these functions. For the KEGG pathway analysis, most of the upregulated genes were annotated in: metabolism, cellular processes, and environmental information processing; for example, DNA repair and metabolism. In a previous study, HSF1 was found to transactivate genes involved in heat shock as cytoprotective proteins [32]. The results of the present study indicate that the “antioxidant activity” process can be triggered as part of the HSF1 regulation in response to heat shock. This finding is supported by the results of Madeira et al. [57] regarding two common fish species living in an environment commonly exposed to high temperatures.

Although HSP70 genes are expanded in Pacific oyster, it remains unclear which genes are heat-inducible and how these genes are regulated. The expression level of HSP70–2387 changed the most under heat shock conditions (Fig. 6b), and that of the other four genes (HSP70–2823, HSP70–2594, HSP20–4164, and CGI_10022835) were almost the same with or without HSF1 interference under heat shock conditions, suggesting that HSP70–2387 is more likely to be regulated by HSF1 than the other four genes. Two isoforms of HSF1 can be expressed under heat shock and be transferred to the nucleus where they form trimers and bind the HSE to regulate the expression of target genes. The dual luciferase assay has been used for examining regulatory relationships [58] and, in our study, it showed that HSF1 (two isoforms) could activate HSP70–2387 (CGI_10002387) in HEK293T cells. Increases in the expressions of the five HSP genes through time were detected by qRT-PCR, and previous studies have shown that HSP genes are involved in the heat shock process [4, 59]. Gene CGI_10002387 was considered to be regulated by HSF1 and it exhibited the same trend as HSF1, albeit a delay in the mRNA level, suggesting that the expression levels determined by the different methods are reliable. The HSF–HSP pathway has been reported in C. elegans [27] and in some mammals [31].

Although many genes could potentially be regulated by HSF1, only one HSP70 gene was enriched in the unique gene set, therefore not providing the most informative result for understanding the HSF1-HSP regulatory relationship. According to our analyses, several reasons might account for this result. Firstly, the heat shock condition applied was not so extreme, and therefore, gene responses might not have been determined as significant by our statistical analyses. A study on HSP70 overexpression in D. melanogaster larvae showed that although extra HSP70 provides additional protection against the immediate damage from heat stress, abnormally high HSP70 concentrations can decrease growth, development, and survival to adulthood [60]. Two hours of heat shock at 35 °C might be a common condition for Pacific oyster, and thus, not all HSP70s would participate in protection owing to the high energy demands of the process. These results suggest that a trade-off between energy and cell protection might exist for response to heat shock stress. Although few genes were enriched, they provided insights into the initial response of the organism to heat shock conditions. Thus, the results of the present study are relevant for studying the molecular-level reaction and the regulatory mechanism underlying response to heat shock.

Secondly, the levels of HSF1 protein and mRNA might have been further reduced, but not abolished, by the RNAi treatment. The abundance of HSF1 gene may be related to incomplete interference (~ 40%). Although 40% interference suppressed the target gene, the protein level also decreased; however, it should be adequate for the confirmation of gene knockdown effects. Therefore, although there are limitations to the present study, our results provide insights into gene function. With the development of technology, more advanced methods of RNAi can be used with model organisms, such as nematodes, to further study the gene function of mollusks.

Conclusions

The present study was designed to determine the effects of heat shock stress and the regulatory mechanisms of HSF1 when the Pacific oyster is exposed to elevated temperatures. The DEG and WGCNA analyses revealed that seven HSP genes and 17 HSP genes were enriched during heat shock, respectively. These included HSP70 and HSP20 genes only, suggesting that these genes were initially triggered under thermal stress. In addition, 48 genes were upregulated and 47 genes were downregulated by HSF1, and the regulation of one HSP70 gene that was enriched in the upregulated gene set by HSF1 was preliminarily validated. The dual luciferase assay allowed examining the regulatory relationship between HSF1 and the enriched HSP70 gene. In addition, the DEG and WGCNA analyses revealed that the HIF signaling pathway was associated with heat shock. These findings will be of interest for future studies aiming to understand the mechanisms underlying thermal tolerance in Pacific oyster.

Methods

Experimental materials

Pacific oyster (Crassostrea gigas) was cultured at Shentanggou (36°21′N, 120°41′E), Qingdao, Shandong province, China. After separating the oysters into monomers, 7-month-old growing oysters of intermediate size were selected for the experiment. These individuals were 58.6 ± 1.3 mm in height and 23.4 ± 1.6 g in weight (values are mean ± standard errors of the mean). During the acclimation period, individuals were fed spirulina powder, and the water temperature was maintained at an average of 12 ± 1 °C for 1 week (water was changed daily).

After the acclimation period, the group of individuals subject to this treatment only was used as the control group (group 1). Individuals that were then subjected to heat shock treatment at 35 °C for 2 h (group 2) were used for detecting heat shock responsive genes. The RNAi-treated group (group 3) and the group treated with RNAi and heat shock (group 4) were used for detecting genes regulated by HSF1 in response to heat shock by comparing the genes expressed in these groups with genes expressed in the other two groups. A diagram of the experimental design is presented in Additional file 1: Figure S5. For the experimental design, we normalized each condition to the control (group 1). The difference between groups 2 and 4 is RNAi of HSF1, which means that without the HSF1 RNAi, gene set 1 should be similar to gene set 3. Therefore, the overlap of gene sets 1 and 3 means that the genes were regulated independently of HSF1 gene in response to heat shock. For gene sets 1 and 2, the genes that overlap indicate that gene expression changed under either heat shock or HSF1 interference. The overlap of gene sets 2 and 3 indicates that the genes regulated by HSF1 were independent from heat shock. Overall, the genes unique to gene set 1, after excluding those that overlapped with gene sets 2 and 3, represented the one regulated by HSF1 in response to heat shock.

RNAi and heat shock conditions

Small interference RNA (siRNA) was synthesized by GenePharma (Shanghai, China) and used for RNAi assays (details of the sequence are shown in Additional file 10: Table S9). After acclimation, individuals were anesthetized according to Suquet [61] (500 g MgCl2 + 5 L seawater + 5 L freshwater) and then randomly divided into four groups, the siRNA group (n = 42), the diethyl pyrocarbonate (DEPC) water group (n = 42), the negative control (NC) group (n = 15), and the control group (n = 42). The muscles of the individuals in the former three groups were injected with 100 μL of 10 μg/100 μL siRNA, DEPC water, or 10 μg/100 μL NC strands, and then the individuals in the four groups were equally divided into two sub-groups for the following treatments. One sub-group was placed into sea water at 35 ± 0.7 °C for 2 h, 46 h after the start of the experiment. At 48 h, their gills were sampled and frozen in liquid nitrogen until used for the qRT-PCR. The time and duration of the heat shock and selected siRNA were optimized in the experiment (see Additional file 1: Figure S1, S2). The control group was used to detect whether the RNAi decreased HSF1 expression.

The qRT-PCR was performed in the ABI 7500 Fast Real-Time PCR System (Applied Biosystems, CA, USA) using the SYBR Green real-time PCR mix (Takara, Japan). The primers for gene expression detection were from Kawabe [25] (HSF1 total) and Li [62] (EF-1a).

Immunoblotting and protein quantification

Individuals in groups 1 (control), 3 (RNAi), and 4 (RNAi and heat shock) were used for HSF1 protein detection (n = 2 per group). Their gills were firstly ground in liquid nitrogen, and proteins were extracted using a cell lysis buffer (Beyotime, Jiangsu, China). Proteins were quantified using the Bradford assay and incubated with 2× protein sodium dodecyl sulfate polyacrylamide gel electrophoresis buffer (Takara) at 100 °C for 3 min. Proteins were analyzed by western blotting using the anti-HSF1 monoclonal antibody (unpublished data) and β-tubulin antibody (CWbiotech, China) as controls. Band intensity was quantified using ImageJ software (http://imagej.nih.gov/ij/).

RNA preparation and sequencing

The gills were also used for total RNA extraction with an RNAprep Kit (Tiangen, China) according to the manufacturer’s instructions. The quality and quantity of RNA were assessed using 1.5% agarose gel electrophoresis and the NanoDrop 2000c UV-Vis Spectrophotometer (Thermo Fisher Scientific, MA, USA), respectively.

According to the analysis of HSF1 expression, 15 individuals of each group (three replicates of five individuals) were selected for RNA-seq. RNA integrity was assessed using the RNA Nano 6000 Assay Kit and the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). One microgram of RNA per sample was used for generating sequencing libraries with the NEBNext Ultra™ RNA Library Prep Kit for Illumina (NEB, MA, USA) following the manufacturer’s recommendations. Index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5×). First strand complementary DNA (cDNA) was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase. Second strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of the 3′ ends of DNA fragments, NEBNext Adaptors with a hairpin loop structure were ligated to prepare for hybridization. To preferentially select cDNA fragments of 240 bp, the library fragments were purified with the AMPure XP system (Beckman Coulter, CA, USA). Then, 3 μl USER enzyme (NEB) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C. The PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers, and the index (X) primer. Finally, PCR products were purified (AMPure XP system) and library quality was assessed using the Agilent Bioanalyzer 2100 system (Agilent Technologies). Clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v4-cBot-HS kit (Illumina, CA, USA) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq X Ten platform and paired-end reads were generated.

RNA-seq data analysis

The quality control of raw data in FASTQ file format was firstly processed through in-house scripts. The scripts were used to filter out adapter sequences, duplicated sequences, ambiguous reads (‘N’), and low-quality reads. At the same time, Q20, Q30, GC-content, and sequence duplication were calculated for the high-quality data. All downstream analyses were based on the high-quality data. Software HISAT2 [63] was used to map the reads on C. gigas genome. Based on the C. gigas genome sequence, mapped reads were spliced using StringTie [64] software and compared with the original genomic annotation information to find the transcripts which have never been annotated before in the reference genome (C. gigas genome) in order to identify new transcripts and new genes, thus complementing the original genome annotation information. New genes were screened by filtering out sequences encoding short peptide chains (< 50 amino acid residues) or containing only a single exon (in order to ensure the result of new gene discovering was reliable, we filtered out the sequence which contained only one exon according to the observation that most eukaryotic genes have introns [65]). Gene function was annotated based on the National Center for Biotechnology Information non-redundant, Protein family, Swiss-Prot, KEGG Orthology (KO), and GO databases. The fragments per kilo base of exon per million reads mapped (FPKM) [66] were used to estimate gene expression levels.

Differential gene expression analysis of selected groups was performed using DESeq2 [67], which provided statistical routines for determining differential expression in digital gene expression data using a model based on a negative binomial distribution. The resulting p-values were adjusted using the Benjamini and Hochberg approach for controlling FDR. Genes with adjusted p-values < 0.05 found by DESeq2 were considered differentially expressed. Venny 2.1 (http://bioinfogp.cnb.csic.es/tools/venny/index.html) was used to construct Venn diagrams with the significantly altered mRNAs for each condition (FDR < 0.05) compared to the control group. Gene ontology is an internationally standardized gene-function classification system consisting of terms that provide a global representation of gene functions using a controlled vocabulary. The GO enrichment analysis was performed using the GOseq R package [68]. Significantly enriched metabolic and signal transduction pathways containing DEGs were determined using KEGG pathway enrichment analysis and compared with the whole genome background. Software KOBAS [69] was used to test the statistical enrichment of DEGs in the KEGG pathways.

Identification of coexpression modules

The R package WGCNA [36] was used to identify modules of highly correlated genes based on the FPKM data. Genes with low FPKM [FPKM < 1 for 12 samples (group 1: NC-1, NC-2, NC-3; group 2: HC-1, HC-2, HC-3; group 3: NR-1, NR-2, NR-3; group 4: HR-1, HR-2, HR-3)] or a module similarity threshold < 0.25 were filtered out. The resulting adjacency matrix was then converted to a topological overlap (TO) matrix with the TO similarity algorithm and genes were hierarchically clustered based on TO similarity. The Dynamic Tree Cut package (https://cran.r-project.org/web/packages/dynamicTreeCut/) was used to cut the hierarchal clustering tree, and modules were defined as the branches from the tree cutting. Modules with fewer than 30 genes were merged into their closest larger neighbor module. Each module was summarized by the first principal component of their scaled module expression profiles (referred to as module eigengene [ME]). Module eigengene-based connectivity (KME) of a gene to a given module was calculated as Pearson correlation coefficient (PCC) between the expression levels (FPKMs) of the gene and the ME of the module using the KME algorithm. Finally, genes were reassigned using the module Merge Using KME Algorithm to ensure that each gene had the highest KME in its own assigned module.

Visualization of hub genes

Genes with the highest degree of connectivity within a module are referred to as hub genes [36]. The top 150 connections (based on topological overlap) among the genes in each module ranked by KME were visualized by VisANT [70]. The GO term enrichment analysis of the WGCNA-identified coexpression modules was performed using a modified Fisher’s exact test in Blast2GO (https://www.blast2go.com/) software (FDR < 0.05). The GO annotations for coexpressed modules were obtained from the C. gigas genome (GenBank accession No., GCA_000297895.1). For each module, only protein-coding genes were subject to GO enrichment analysis, and the most specific biological processes were reported.

Validation of RNA-seq data by qRT-PCR and dual luciferase reporter assay

The qRT-PCR was performed for individuals subjected to heat shock at 35 °C and sampled at 15 and 30 min, and at 1, 2, 6, 12, and 24 h. Five individuals were sampled each time. After sampling, total RNA was prepared using an RNAprep Pure Tissue Kit (Tiangen) following the manufacturer’s protocol. The quality and quantity of prepared RNA were assessed using 1.5% agarose gel electrophoresis and the NanoDrop 2000c UV-Vis Spectrophotometer (Thermo Fisher Scientific). Each RNA sample was reverse transcribed using a cDNA synthesis kit (Takara) for gene expression detection by qRT-PCR. The primer used for HSF-total was from a previous study [25]. We used the housekeeping gene EF-1a as a control gene. The primers used are listed in Additional file 11: Table S10. The qRT-PCR was performed as described in a previous section.

The dual-luciferase assay has been widely used in cell lines to rapidly and accurately determine the activity of a given promoter. The Dual-Glo® Luciferase Assay System (Promega, WI, USA) was used to estimate the effect of HSF1 on the mRNA level of the target gene in HEK293T cells following the manufacturer’s instructions. The potential binding site, the area of the HSP70 promoter where HSF1 can bind, was predicted by AnimalTFDB 3.0, following Guo [71]. The region ~ 2 kb upstream of the HSP70 gene (CGI_10002387) was cloned into the pGL3-basic luciferase reporter plasmid, and the HSF1a and HSF1d fragments were constructed to pCMV-N-Myc plasmids. Cells were co-transfected with varying amounts of pRL-TK, pGL3-basic-2387, pCMV-N-Myc-HSF1a/pCMV-N-Myc-HSF1d, and pCMV-N-Myc plasmids and reporter gene activity was measured using the Dual-Luciferase Reporter Assay System 24 h after transfection. Fluorescence was detected using the Varioskan Flash Multimode Reader (Thermo Fisher Scientific). According to the fluorescence value of each group, we estimated the regulatory activity of HSF1.

Statistical analyses

Statistical analyses were performed using the software SPSS 16.0 (IBM, NY, USA). One-way analysis of variance (ANOVA) followed by Duncan’s test was used to evaluate the expression levels of HSF1 and control genes. The results were significant at p < 0.05.

Availability of data and materials

The transcriptome datasets supporting the results of this study are available in the NCBI BioProject repository (Accession Number PRJNA516762, https://www.ncbi.nlm.nih.gov/bioproject/PRJNA516762) and in the SRA database (https://www.ncbi.nlm.nih.gov/sra) (Accession Number SRP181984).

Abbreviations

DEGs:

Differentially expressed genes

FDR:

False discovery rate

FPKM:

Fragments per kilo base of exon per million reads mapped

GO:

Gene ontology

HIF :

Hypoxia induced factor

HS:

Heat shock

HSC:

Heat shock cognate

HSEs:

Heat shock element

HSF1 :

Heat shock transcription factor 1

HSP :

Heat shock protein

KEGG:

Kyoto Encyclopedia of Genes and Genomes

KO:

KEGG orthology

ME:

Module eigengene

NC:

Negative control

qRT-PCR:

Quantitative real-time polymerase chain reaction

RNAi:

RNA interference

RNA-seq:

RNA-sequencing

SEM:

Standard error of means

siRNA:

Small interference RNA

WGCNA:

Weighted gene co-expression network analysis

References

  1. Portner HO, Farrell AP. ECOLOGY physiology and climate change. Science. 2008;322(5902):690–2.

    Article  PubMed  Google Scholar 

  2. Viña J. Biochemical adaptation: mechanism and process in physiological evolution. Biochem Mol Biol Educ. 2002;30(3):215–6.

    Article  Google Scholar 

  3. Pörtner HO, Berdal B, Blust R, Brix O, Colosimo A, De Wachter B, Giuliani A, Johansen T, Fischer T, Knust R, et al. Climate induced temperature effects on growth performance, fecundity and recruitment in marine fish: developing a hypothesis for cause and effect relationships in Atlantic cod (Gadus morhua) and common eelpout (Zoarces viviparus). Cont Shelf Res. 2001;21(18):1975–97.

    Article  Google Scholar 

  4. Zhang GF, Fang XD, Guo XM, Li L, Luo RB, Xu F, Yang PC, Zhang LL, Wang XT, Qi HG, et al. The oyster genome reveals stress adaptation and complexity of shell formation. Nature. 2012;490(7418):49–54.

    Article  CAS  PubMed  Google Scholar 

  5. Ugalde SC, Preston J, Ogier E, Crawford C. Analysis of farm management strategies following herpesvirus (OsHV-1) disease outbreaks in Pacific oysters in Tasmania, Australia. Aquaculture. 2018;495:179–86.

    Article  Google Scholar 

  6. Zhang H, Wang H, Chen H, Wang M, Zhou Z, Qiu L, Wang L, Song L. The transcriptional response of the Pacific oyster Crassostrea gigas under simultaneous bacterial and heat stresses. Dev Comp Immunol. 2019;94:1–10.

    Article  CAS  PubMed  Google Scholar 

  7. Lang RP, Bayne CJ, Camara MD, Cunningham C, Jenny MJ, Langdon CJ. Transcriptome profiling of selectively bred Pacific oyster Crassostrea gigas families that differ in tolerance of heat shock. Mar Biotechnol. 2009;11(5):650–68.

    Article  CAS  Google Scholar 

  8. Ximing G, Yan H, Linlin Z, Christophe L, Aude J. Immune and stress responses in oysters with insights on adaptation. Fish Shellfish Immunol. 2015;46(1):107–19.

    Article  CAS  Google Scholar 

  9. Zhang G, Li L, Meng J, Qi H, Qu T, Xu F, Zhang L. Molecular basis for adaptation of oysters to stressful marine intertidal environments. Ann Rev Anim Biosci. 2015;4(1):357.

    Article  CAS  Google Scholar 

  10. Clegg JS, Uhlinger KR, Jackson SA, Cherr GN, Rifkin E, Friedman CS. Induced thermotolerance and the heat shock protein-70 family in the Pacific oyster Crassostrea gigas. Mol Mar Biol Biotechnol. 1998;7(1):21–30.

    CAS  Google Scholar 

  11. Zhu Q, Zhang L, Li L, Que H, Zhang G. Expression characterization of stress genes under high and low temperature stresses in the Pacific oyster, Crassostrea gigas. Mar Biotechnol (NY). 2016;18(2):176–88.

    Article  CAS  Google Scholar 

  12. Hamdoun AM, Cheney DP, Cherr GN. Phenotypic plasticity of HSP70 and HSP70 gene expression in the Pacific oyster (Crassostrea gigas): implications for thermal limits and induction of thermal tolerance. Biol Bull. 2003;205(2):160–9.

    Article  CAS  PubMed  Google Scholar 

  13. Kim J-H, Jeong S-Y, Kim P-J, Dahms H-U, Han K-N. Bio-effect-monitoring of long-term thermal wastes on the oyster, Crassostrea gigas, using heat shock proteins. Mar Pollut Bull. 2017;119(1):359–64.

    Article  CAS  PubMed  Google Scholar 

  14. Lim HJ, Kim BM, Hwang IJ, Lee JS, Choi IY, Kim YJ, Rhee JS. Thermal stress induces a distinct transcriptome profile in the Pacific oyster Crassostrea gigas. Comp Biochem Physiol Part D Genomics Proteomics. 2016;19:62–70.

    Article  CAS  PubMed  Google Scholar 

  15. Kim B-M, Kim K, Choi I-Y, Rhee J-S. Transcriptome response of the Pacific oyster, Crassostrea gigas susceptible to thermal stress: a comparison with the response of tolerant oyster. Mol Cell Toxicol. 2017;13(1):105–13.

    Article  CAS  Google Scholar 

  16. Zhang Y, Sun J, Mu H, Li J, Zhang Y, Xu F, Xiang Z, Qian PY, Qiu JW, Yu Z. Proteomic basis of stress responses in the gills of the Pacific oyster Crassostrea gigas. J Proteome Res. 2015;14(1):304–17.

    Article  CAS  PubMed  Google Scholar 

  17. Takeuchi T, Koyanagi R, Gyoja F, Kanda M, Hisata K, Fujie M, Goto H, Yamasaki S, Nagai K, Morino Y, et al. Bivalve-specific gene expansion in the pearl oyster genome: implications of adaptation to a sessile lifestyle. Zool Lett. 2016;2.

  18. Uliano-Silva M, Americo JA, Brindeiro R, Dondero F, Prosdocimi F, Rebelo MD. Gene Discovery through Transcriptome Sequencing for the Invasive Mussel Limnoperna fortunei. PLoS One. 2014;9(7):e102973.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  19. Cheng J, Xun XG, Kong YF, Wang SY, Yang ZH, Li YJ, Kong DX, Wang S, Zhang LL, Hu XL, et al. Hsp70 gene expansions in the scallop Patinopecten yessoensis and their expression regulation after exposure to the toxic dinoflagellate Alexandrium catenella. Fish Shellfish Immunol. 2016;58:266–73.

    Article  CAS  PubMed  Google Scholar 

  20. Zhang G, Li L, Meng J, Qi H, Qu T, Xu F, Zhang L. Molecular basis for adaptation of oysters to stressful marine intertidal environments. Ann Rev Anim Biosci. 2016;4(1):357–81.

    Article  CAS  Google Scholar 

  21. Hofmann G, Somero G. Evidence for protein damage at environmental temperatures: seasonal changes in levels of ubiquitin conjugates and hsp70 in the intertidal mussel Mytilus trossulus. J Exp Biol. 1995;198(7):1509.

    CAS  PubMed  Google Scholar 

  22. Sun J, Zhang Y, Xu T, Zhang Y, Mu H, Zhang Y, Lan Y, Fields CJ, Hui JHL, Zhang W, et al. Adaptation to deep-sea chemosynthetic environments as revealed by mussel genomes. Nat Ecol Evol. 2017;1:0121.

    Article  Google Scholar 

  23. Gomezpastor R, Burchfiel ET, Thiele DJ. Regulation of heat shock transcription factors and their roles in physiology and disease. Nat Rev Mol Cell Biol. 2018;19(1):4–19.

    Article  CAS  Google Scholar 

  24. Ortner V, Ludwig A, Riegel E, Dunzinger S, Czerny T. An artificial HSE promoter for efficient and selective detection of heat shock pathway activity. Cell Stress Chaperones. 2015;20(2):277–88.

    Article  CAS  PubMed  Google Scholar 

  25. Kawabe S, Yokoyama Y. Novel isoforms of heat shock transcription factor 1 are induced by hypoxia in the Pacific oyster Crassostrea gigas. J Exp Zool A Ecol Genet Physiol. 2011;315(7):394–407.

    Article  CAS  PubMed  Google Scholar 

  26. Walker GA, Thompson FJ, Brawley A, Scanlon T, Devaney E. Heat shock factor functions at the convergence of the stress response and developmental pathways in Caenorhabditis elegans. Faseb J. 2003;17(11):1960.

    Article  CAS  PubMed  Google Scholar 

  27. Jessica B, Stephanie M, Lu Y, Cheng F, Westerheide SD. The genome-wide role of HSF-1 in the regulation of gene expression in Caenorhabditis elegans. BMC Genomics. 2016;17(1):559.

    Article  CAS  Google Scholar 

  28. Clos J, Westwood JT, Becker PB, Wilson S, Lambert K, Wu C. Molecular cloning and expression of a hexameric Drosophila heat shock factor subject to negative regulation. Cell. 1990;63(5):1085–97.

    Article  CAS  PubMed  Google Scholar 

  29. Tian SB, Haney RA, Feder ME. Phylogeny Disambiguates the Evolution of Heat-Shock cis-Regulatory Elements in Drosophila. PLoS One. 2010;5(5):266–73.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  30. Baird NA, Turnbull DW, Johnson EA. Induction of the heat shock pathway during hypoxia requires regulation of heat shock factor by hypoxia-inducible factor-1. J Biol Chem. 2006;281(50):38675–81.

    Article  CAS  PubMed  Google Scholar 

  31. Mahat DB, Salamanca HH, Duarte FM, Danko CG, Lis JT. Mammalian heat shock response and mechanisms underlying its genome-wide transcriptional regulation. Mol Cell. 2016;62(1):63–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Dayalan Naidu S, Dinkova-Kostova AT. Regulation of the mammalian heat shock factor 1. FEBS J. 2017;284(11):1606–27.

    Article  CAS  PubMed  Google Scholar 

  33. Boutros M, Kiger AA, Armknecht S, Kerr K, Hild M, Koch B, Haas SA, Consortium HFA, Paro R, Perrimon N. Genome-wide RNAi analysis of growth and viability in Drosophila cells. Science. 2004;303(5659):832.

    Article  CAS  PubMed  Google Scholar 

  34. Cullen LM, Arndt GM. Genome-wide screening for gene function using RNAi in mammalian cells. Immunol Cell Biol. 2005;83(3):217–23.

    Article  CAS  PubMed  Google Scholar 

  35. Horvath S, Zhang B, Carlson M, Lu KV, Zhu S, Felciano RM, Laurance MF, Zhao W, Qi S, Chen Z, et al. Analysis of oncogenic signaling networks in glioblastoma identifies ASPM as a molecular target. Proc Natl Acad Sci. 2006;103(46):17402.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9(1):559.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Zhang LL, Hou R, Su HL, Hu XL, Wang S, Bao ZM. Network Analysis of Oyster Transcriptome Revealed a Cascade of Cellular Responses during Recovery after Heat Shock. PLoS One. 2012;7(4):e35484.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Frydman J. Folding of newly translated proteins in vivo: the role of molecular chaperones. Annu Rev Biochem. 2001;70(1):603–47.

    Article  CAS  PubMed  Google Scholar 

  39. Hartl FU, Hayer-Hartl M. Molecular chaperones in the cytosol: from nascent chain to folded protein. Science. 2002;295(5561):1852.

    Article  CAS  PubMed  Google Scholar 

  40. Young JC, Agashe VR, Siegers K, Hartl FU. Pathways of chaperone-mediated protein folding in the cytosol. Nat Rev Mol Cell Biol. 2004;5:781.

    Article  CAS  PubMed  Google Scholar 

  41. Sørensen JG, Kristensen TN, Loeschcke V. The evolutionary and ecological role of heat shock proteins. Ecol Lett. 2003;6(11):1025–37.

    Article  Google Scholar 

  42. Meistertzheim A-L, Tanguy A, Moraga D, Thébault M-T. Identification of differentially expressed genes of the Pacific oyster Crassostrea gigas exposed to prolonged thermal stress. FEBS J. 2007;274(24):6392–402.

    Article  CAS  PubMed  Google Scholar 

  43. Scott GR, Johnston IA. Temperature during embryonic development has persistent effects on thermal acclimation capacity in zebrafish. Proc Natl Acad Sci U S A. 2012;109(35):14247–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Lund EL, Høg A, Olsen MWB, Hansen LT, Engelholm SA, Kristjansen PEG. Differential regulation of VEGF, HIF1α and angiopoietin-1, −2 and −4 by hypoxia and ionizing radiation in human glioblastoma. Int J Cancer. 2004;108(6):833–8.

    Article  CAS  PubMed  Google Scholar 

  45. Bracken CP, Whitelaw ML, Peet DJ. The hypoxia-inducible factors: key transcriptional regulators of hypoxic responses. Cell Mol Life Sci. 2003;60(7):1376–93.

    Article  CAS  PubMed  Google Scholar 

  46. Pörtner HO. Oxygen- and capacity-limitation of thermal tolerance: a matrix for integrating climate-related stressor effects in marine ecosystems. J Exp Biol. 2010;213(6):881–93.

    Article  PubMed  Google Scholar 

  47. Tomanek L. Proteomics to study adaptations in marine organisms to environmental stress. J Proteome. 2014;105:92–106.

    Article  CAS  Google Scholar 

  48. Lars T, Zuzow MJ. The proteomic response of the mussel congeners Mytilus galloprovincialis and M. trossulus to acute heat stress: implications for thermal tolerance limits and metabolic costs of thermal stress. J Exp Biol. 2010;213(20):3559–74.

    Article  CAS  Google Scholar 

  49. Sokolova IM, Frederich M, Bagwe R, Lannig G, Sukhotin AA. Energy homeostasis as an integrative tool for assessing limits of environmental stress tolerance in aquatic invertebrates. Mar Environ Res. 2012;79(4):1–15.

    Article  CAS  PubMed  Google Scholar 

  50. Fabioux C, Corporeau C, Quillien V, Favrel P, Huvet A. In vivo RNA interference in oyster--vasa silencing inhibits germ cell development. FEBS J. 2009;276(9):2566–73.

    Article  CAS  PubMed  Google Scholar 

  51. Huvet A, Fleury E, Corporeau C, Quillien V, Daniel JY, Riviere G, Boudry P, Fabioux C. In vivo RNA interference of a gonad-specific transforming growth factor-β in the Pacific oyster Crassostrea gigas. Mar Biotechnol. 2012;14(4):402–10.

    Article  CAS  Google Scholar 

  52. Bo YJ, Min SK, Mi YC, Lee SJ, Park MA, Jin WK, Choi SH, Jeong HD, Kim KH. Simultaneous and systemic Knock-down of big Defensin 1 and 2 gene expression in the Pacific oyster Crassostrea gigas using long double-stranded RNA-mediated RNA interference. Fishe Aquat Sci. 2014;17(3):377–80.

    Google Scholar 

  53. Elbashir SM, Harborth J, Lendeckel W, Yalcin A, Weber K, Tuschl T. Duplexes of 21-nucleotide RNAs mediate RNA interference in cultured mammalian cells. Nature. 2001;411:494.

    Article  CAS  PubMed  Google Scholar 

  54. Sledz CA, Holko M, de Veer MJ, Silverman RH, Williams BRG. Activation of the interferon system by short-interfering RNAs. Nat Cell Biol. 2003;5:834.

    Article  CAS  PubMed  Google Scholar 

  55. Persengiev SP, Zhu XC, Green MR. Nonspecific, concentration-dependent stimulation and repression of mammalian gene expression by small interfering RNAs (siRNAs). RNA-Publ RNA Soc. 2004;10(1):12–8.

    Article  CAS  Google Scholar 

  56. Bridge AJ, Pebernard S, Ducraux A, Nicoulaz A-L, Iggo R. Induction of an interferon response by RNAi vectors in mammalian cells. Nat Genet. 2003;34:263.

    Article  CAS  PubMed  Google Scholar 

  57. Madeira C, Mendonça V, Leal MC, Flores AAV, Cabral HN, Diniz MS, Vinagre C. Thermal stress, thermal safety margins and acclimation capacity in tropical shallow waters—an experimental approach testing multiple end-points in two common fish. Ecol Indic. 2017;81:146–58.

    Article  Google Scholar 

  58. Wang T, Meng J, Li L, Zhang G. Characterization of CgHIFalpha-like, a novel bHLH-PAS transcription factor family member, and its role under hypoxia stress in the Pacific oyster Crassostrea gigas. PLoS One. 2016;11(11):e0166057.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  59. Weibezahn J, Schlieker C, Tessarz P, Mogk A, Bukau B. Novel insights into the mechanism of chaperone-assisted protein disaggregation. In: Biological Chemistry, vol. 386; 2005. p. 739.

    Google Scholar 

  60. Krebs RA, Feder ME. Deleterious consequences of Hsp70 overexpression in Drosophila melanogaster larvae. Cell Stress Chaperones. 1997;2(1):60–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Suquet M, de Kermoysan G, Araya RG, Queau I, Lebrun L, Le Souchu P, Mingant C. Anesthesia in Pacific oyster, Crassostrea gigas. Aquat Living Resour. 2009;22(1):29–34.

    Article  CAS  Google Scholar 

  62. Li A, Li L, Song K, Wang W, Zhang G. Temperature, energy metabolism, and adaptive divergence in two oyster subspecies. Ecol Evol. 2017;7(16):6151–62.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Xing D, Wang Y, Hamilton M, Ben-Hur A, Reddy ASN. Transcriptome-wide identification of RNA targets of Arabidopsis SERINE/ARGININE-RICH45 uncovers the unexpected roles of this RNA binding protein in RNA processing. Plant Cell. 2015;27(12):3294.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Florea L, Song L, Salzberg SL. Thousands of exon skipping events differentiate among splicing patterns in sixteen human tissues. F1000Res. 2013;2:188.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  67. Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26(1):136–8.

    Article  PubMed  CAS  Google Scholar 

  68. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  69. Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, Kong L, Gao G, Li C-Y, Wei L. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011;39(suppl_2):W316–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Hu Z, Mellor J, DeLisi C. Analyzing networks with VisANT. Curr Protoc Bioinformatics. 2004;8(1):8.8.1–8.8.24.

    Article  Google Scholar 

  71. Hu H, Miao Y-R, Jia L-H, Yu Q-Y, Zhang Q, Guo A-Y. AnimalTFDB 3.0: a comprehensive resource for annotation and prediction of animal transcription factors. Nucleic Acids Res. 2019;47(D1):D33–8.

    Article  PubMed  CAS  Google Scholar 

Download references

Acknowledgements

The authors thank all members of the laboratory for helping to sample the oyster and valuable discussion.

Funding

This research was supported National Key R&D Program of China (2018YFD0900304, in design of the study and collection, analysis and interpretation of data and in writing the manuscript), the National Natural Science Fundation of China (31572620, in polish the manuscript) and the Earmarked Fund for Modern Agro-industry Technology Research System (CARS-49, in sample collection).

Author information

Authors and Affiliations

  1. Key Laboratory of Experimental Marine Biology, Institute of Oceanology, Chinese Academy of Sciences, Qingdao, 266071, China

    Youli Liu, Li Li, Baoyu Huang, Wei Wang & Guofan Zhang

  2. Laboratory for Marine Biology and Biotechnology, Qingdao National Laboratory for Marine Science and Technology, Qingdao, 266071, China

    Youli Liu, Wei Wang & Guofan Zhang

  3. University of Chinese Academy of Sciences, Beijing, 100039, China

    Youli Liu

  4. Laboratory for Marine Fisheries Science and Food Production Processes, Qingdao National Laboratory for Marine Science and Technology, Qingdao, 266071, China

    Li Li & Baoyu Huang

  5. Center for Ocean Mega-Science, Chinese Academy of Sciences, Qingdao, 266071, China

    Li Li, Baoyu Huang, Wei Wang & Guofan Zhang

  6. National and Local Joint Engineering Laboratory of Ecological Mariculture, Qingdao, 266071, China

    Li Li, Baoyu Huang, Wei Wang & Guofan Zhang

Authors
  1. Youli Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Li Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Baoyu Huang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Wei Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Guofan Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

LL conceived the project. YLL, LL and GFZ participated in the design of the study and discussed the results. RNA extraction, content measurement of HSF1 and validation the results were conducted by YLL. YLL carried out bioinformatics analyses. YLL, LL and BYH wrote and polished the manuscript. WW collected the experimental materials. All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Li Li.

Ethics declarations

Ethics approval and consent to participate

The oysters used in this study were marine-cultured animal and were collected from a local farm in Qingdao, China. There are no specific permissions required for the sampling of these oysters. All experiments were conducted with approval from the Experimental Animal Ethics Committee, Institute of Oceanology, Chinese Academy of Sciences, China.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Figure S1 qRT-PCR of total HSF-1 gene during heat shock (35 °C for 24 h). Figure S2. Small RNA interference strands were selected and tested separately by measuring the expression level of HSF-1. Figure S3. Scheme for RNA-seq data normalization. Figure S4. Housekeeping genes share a similar expression profile between each RNA-seq treatment condition. Figure S5. Experimental scheme for RNA-seq. Figure S6. The eigengene adjacency heat map of 14 modules. (DOCX 413 kb)

Additional file 2:

Table S1. Quality information of 12 sample data output statistics. (XLSX 10 kb)

Additional file 3:

Table S2. New gene of Crassostrea gigas discovered after filtering (the parameters we applied with StringTie are: -p 6 --merge -F 0.1 -T 0.1 -i -eB -p 6). (XLSX 3207 kb)

Additional file 4:

Table S3. The number of new genes that are annotated in each database. (XLSX 9 kb)

Additional file 5:

Table S4. Significant genes altered in response to each condition after normalization to the control. (XLSX 115 kb)

Additional file 6:

Table S5 Number of differentially expressed genes (DEGs) in each database. (XLSX 9 kb)

Additional file 7:

Table S6. Number of genes of the 14 modules. (XLSX 9 kb)

Additional file 8:

Table S7. Heat shock proteins (HSPs) in the tan, yellow, blue, and salmon modules. (XLSX 12 kb)

Additional file 9:

Table S8. Top six genes with the highest KME values in the tan, yellow, blue, green, black, purple, and salmon modules. (XLSX 10 kb)

Additional file 10:

Table S9. SiRNA of the RNA interference assay. (XLSX 9 kb)

Additional file 11:

Table S10. Primer for validation of the genes enriched by RNA sequencing. (XLSX 9 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Liu, Y., Li, L., Huang, B. et al. RNAi based transcriptome suggests genes potentially regulated by HSF1 in the Pacific oyster Crassostrea gigas under thermal stress. BMC Genomics 20, 639 (2019). https://doi.org/10.1186/s12864-019-6003-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-019-6003-8

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us