PacBio Isoform Sequencing and Illumina RNA Sequencing Provide Novel Insights on Responses to Acute Heat Stress in Apostichopus japonicus Coelomocytes

Significant increases in global sea surface temperatures are expected with climate change and may cause a serious challenge for marine organisms cultured in aquatic environments that are characterized by short and long-term fluctuations in water temperatures. Apostichopus japonicus, a sea cucumber with high nutritional value and pharmacological properties, is an important economic species that is widely raised in aquaculture in China. In recent years, continuous extreme high temperatures (up to 30°C) have occurred frequently in summer leading to mass mortality of sea cucumbers cultured in semi-open shallow regions seriously restricting the sustainable development of sea cucumber aquaculture. In the present study, we combined RNA-seq and PacBio single-molecular real-time (SMRT) sequencing technology to unveil the potential mechanisms of response to acute heat stress in A. japonicus coelomocytes. A total of 1,375 differentially expressed genes (DEGs) were identified in a comparison of control and 48 h heat stress (HS) groups. Pathway enrichment analysis indicated that nine important pathways induced by HS were significantly enriched (q-value < 0.05) and mostly fell into four classes: folding, sorting, and degradation, immune and infectious diseases, signal transduction, and post-transcriptional regulation. Among them, all 41 genes connected with protein processing in endoplasmic reticulum were significantly up-regulated, and 12 of these were selected and validated via qPCR. Furthermore, changes in alternative splicing (AS) were also identified in sea cucumbers following HS. A total of 1,224 and 1,251 differential alternative splicing (DAS) events were identified using splice junction counts (JC only) and reads on target and junction counts (JCEC) as the input for rMATS in CO-HS comparison. We further found that the RNA splicing-related genes were enriched in the spliceosome pathway and showed DAS in control versus heat-stressed animals. In particular, we compared and confirmed that the hsfs1 gene, the master regulator of the heat shock response, showed differentially spliced exons in response to HS. This is the first comprehensive study showing that transcriptional and post-transcriptional (AS) controls are involved in the acute heat stress response of sea cucumber coelomocytes and provides novel insights into the molecular mechanisms of echinoderm adaptation to environmental stress.

Significant increases in global sea surface temperatures are expected with climate change and may cause a serious challenge for marine organisms cultured in aquatic environments that are characterized by short and long-term fluctuations in water temperatures. Apostichopus japonicus, a sea cucumber with high nutritional value and pharmacological properties, is an important economic species that is widely raised in aquaculture in China. In recent years, continuous extreme high temperatures (up to 30 • C) have occurred frequently in summer leading to mass mortality of sea cucumbers cultured in semi-open shallow regions seriously restricting the sustainable development of sea cucumber aquaculture. In the present study, we combined RNA-seq and PacBio single-molecular real-time (SMRT) sequencing technology to unveil the potential mechanisms of response to acute heat stress in A. japonicus coelomocytes. A total of 1,375 differentially expressed genes (DEGs) were identified in a comparison of control and 48 h heat stress (HS) groups. Pathway enrichment analysis indicated that nine important pathways induced by HS were significantly enriched (q-value < 0.05) and mostly fell into four classes: folding, sorting, and degradation, immune and infectious diseases, signal transduction, and post-transcriptional regulation. Among them, all 41 genes connected with protein processing in endoplasmic reticulum were significantly up-regulated, and 12 of these were selected and validated via qPCR. Furthermore, changes in alternative splicing (AS) were also identified in sea cucumbers following HS. A total of 1,224 and 1,251 differential alternative splicing (DAS) events were identified using splice junction counts (JC only) and reads on target and junction counts (JCEC) as the input for rMATS in CO-HS comparison. We further found that the RNA splicing-related genes were enriched in the spliceosome pathway and showed DAS in control versus heat-stressed animals. In particular, we compared and confirmed that the hsfs1 gene, the master regulator of the heat shock response, showed differentially INTRODUCTION Climate change has been predicted to cause a significant increase in global sea surface temperatures (Collins et al., 2013) and alter amplitudes of diurnal temperature fluctuations (Easterling et al., 2000;Meehl et al., 2000). Over the past several decades, rising sea water temperature has impinged upon all biochemical and physiological processes of marine organisms (Somero, 2010). Heat waves and seasonal variation linked with global warming are causing frequent fluctuations in water temperature in aquatic environments, causing further challenges for marine organisms, particularly those in shallow water locations (Frölicher et al., 2018;Oliver et al., 2018).
Sea cucumbers play an important role in marine ecosystem functionality through nutrient recycling, bioturbation, benthic productivity, and meiofauna population dynamics (Schneider et al., 2011;MacTavish et al., 2012;Oh et al., 2017). To date, more than 1,250 species of sea cucumbers have been identified in benthic areas of both shallow and deep seas world-wide and approximately 20 of them are known to be edible (Jo and Park, 2016). The sea cucumber Apostichopus japonicus, has a high nutritional value and important pharmacological properties. The species has long been popular as a food and folk medicine in Asian markets and is widely cultured in China as one of the most important economic species (Chang et al., 2009;Bordbar et al., 2011). As a poikilothermic marine animal and a temperate species, A. japonicus is very sensitive to temperature variation and shows a limited ability to adapt to environmental temperature changes (Xu et al., 2016). Over a long evolutionary process, sea cucumbers have acquired a high temperature-induced special physiological behavior, estivation, to aid survival and reproduction (Yang et al., 2005). Since the late summer of 2013, continuous extreme high temperatures (up to 30 • C) occurred frequently during the summer in Shandong and Hebei provinces of China and have resulted in mass mortality of sea cucumbers cultured in semi-open regions, seriously restricting the sustainable development of sea cucumber aquaculture .
In previous studies, heat stress (HS) responses of A. japonicus were investigated and it was demonstrated that acute and/or long-term exposure to HS results in extensive changes in gene transcription (Meng et al., 2009;Zhao et al., 2011;Xu et al., 2014Xu et al., , 2016Xu et al., , 2018Li and Xu, 2018;Huo et al., 2019;Xu and Wang, 2019). Although transcriptome analyses have greatly broadened our knowledge of the HS response in sea cucumbers, the available data are still not enough because of the limited annotations and different sampling times and tissues reported in different studies. Most previous transcriptome analyses of the HS response in sea cucumbers were performed using intestine, respiratory tree, or muscle tissues Huo et al., 2019;Li et al., 2019). However, information on the HS responses of coelomocytes at the transcriptional level are still limited. Coelomocytes are suspended in coelomic fluid and function similarly to the hemocytes of other invertebrates. Previous studies have reported that immune function in sea cucumbers is mainly executed by coelomocytes, and that humoral immunity is based on a variety of macromolecules in the coelomic cavity that are secreted by coelomocytes Ramírez-Gómez et al., 2010;Xue et al., 2015).
In addition to transcriptional changes, post-transcriptional modifications generated by alternative splicing (AS) also play important roles in the HS response (Fujikake et al., 2005;Liu et al., 2013;Tian et al., 2019). HS-induced AS has been suggested to be a common feature among heat shock factor (hsf ) transcripts, being reported from Drosophila melanogaster (Fujikake et al., 2005), zebrafish (Danio rerio) (Råbergh et al., 2000), and catfish (Ictalurus punctatus) (Tan et al., 2019). However, to date, there have been no studies of the potential role of AS in the HS response of sea cucumbers. Therefore, the strategy of applying "Next-generation" sequencing (NGS) technology and single-molecular real-time (SMRT) sequencing combined with genome reference to analyze AS under HS was worth exploring.
Currently, SMRT sequencing has been confirmed as a valuable tool for the improvement and identification of alternative isoforms produced by AS with higher accuracy by optimizing gene structures of genome annotation. In the present study, we aimed to investigate transcriptional and post-transcriptional changes generated by AS in response to HS in sea cucumber coelomocytes. Both Pacific Biosciences (PacBio) SMRT and Illumina-based sequencing were performed to construct transcriptomes from control and 48 h heat-shocked sea cucumbers. Our results help to elucidate the underlying mechanisms of the acute HS response in sea cucumber coelomocytes and provide important clues for further study of the molecular physiology of A. japonicus.

Sample Collection for PacBio Sequel Sequencing
One hundred and twenty A. japonicus adults (100-120 g) were collected from an aquaculture farm in Weihai (Shandong, China) in April, 2019. The animals were acclimated in a seawater aquarium system at 15 • C and 30 salinity for 2 weeks and fed with an artificial food compound (30% gulfweed powder, 20% murina powder, and 50% sea mud) once a day. All animals were fasted for 2 days prior to experiments.
After acclimation, animals were randomly divided into control (CO, 15 • C) or experimental groups (48 h HS, 30 • C) with three parallel tanks for each group. Dissolved oxygen (DO) was maintained at 8 mg/L in both groups. Following exposure, coelomocytes were collected immediately and frozen rapidly in liquid nitrogen from all groups. All samples were stored at −80 • C for subsequent analysis. All animal protocols for care, experimentation, and sacrifice had the approval of the Experimental Animal Ethics Committee of the Ocean University of China.

RNA Extraction
Total RNA was isolated from coelomocytes using Trizol RNA isolation reagent (TaKaRa, Japan) following the manufacturer's instructions and digested with RNase-free DNase I (TaKaRa, Japan) to remove genomic DNA. RNA degradation and contamination were monitored using 1% agarose gels. The purity and concentration of RNA was measured using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Rockland, DE, United States) with OD 260/280 reading. Accurate quantification and integrity of RNA was further determined using a Qubit RNA Detection Kit with a Qubit 2.0 Fluorometer (Invitrogen, CA, United States) and Agilent 2100 Bioanalyzer (Agilent Technologies, CA, United States).

PacBio Library Construction and Sequencing
The total RNA samples from control (n = 3) and 48 h heat stressed (n = 3) A. japonicus coelomocytes were pooled and sequenced using the Pacific Biosciences Isoform sequencing (Iso-Seq) platform. The mRNA with polyA tails was enriched by use of Oligo (dT) magnetic beads and reverse transcribed into cDNA using a Clontech SMARTer PCR cDNA Synthesis Kit (TaKaRa, Japan). Then optimized PCR reactions were used to generate double-stranded cDNA, followed by size selection (fragments larger than 4 kb) using a Blue Pippin TM Size-Selection System (Sage Science, MA, United States) as described by PacBio (PN 100-092-800-03). The generated cDNA was then re-amplified using PCR to obtain sufficient cDNA. Full-length cDNA damage/terminal repair, ligation of the SMRT dumbbelltype adapters and SMRTbell template preparation were then carried out. The SMRTbell template was bound with sequencing primer and polymerase and sequenced on the PacBio sequel platform according to the manufacturer's protocol. The Iso-seq datasets generated are available from the NCBI Sequence Read Archive database (SRA) with accession number PRJNA785124.

SMRT Sequencing Data Processing
The raw SMRT sequencing reads were processed using PacBio SMRTlink v5.0 software. After removing the adapter and lowquality data, the circular consensus sequence (CCS) reads were generated from subreads.bam files and further classified as fulllength (FL) and non-full-length (NFL) reads based on following parameters: whether the sequence contains both 5 primer, 3 primer, and poly A tail signals. Subsequently, the full-length nonchimeric (FLNC) reads were clustered by Iterative Clustering for Error Correction (ICE) software to generate the cluster consensus isoforms. To obtain FL polished consensus sequences, Quiver software was used to finally obtain the FL polished high quality consensus sequences (accuracy 99%), which were applied for gene structure optimization of the A. japonicus reference genome (NCBI BioProject ID PRJNA757529).

Open Reading Frame Prediction and Gene Structure Optimization
The open reading frames (ORFs) were analyzed using ANGEL (Shimizu et al., 2006) software for the isoform sequences to obtain the coding sequences (CDS), protein sequences, and UTR sequences. To optimize gene structure of the A. japonicus reference genome, each isoform from FL polished high-quality consensus sequences was compared with the existing gene models of A. japonicus genome annotation. Known isoforms were 5 UTR or 3 UTR extended or shortened. Other isoforms were further classified into span isoforms, exonic overlap isoforms and intron isoforms based on their exon structure (splicing junctions).

RNA-seq Library Construction and Sequencing
Total RNA from nine individuals per group (CO and HS) was pooled into three biological replicates (three individuals each) and a total of six libraries were constructed. The Illumina libraries were generated using NEBNext R Ultra TM RNA Library Prep Kit (E7530L) (NEB, United States) following manufacturer's recommendations, and index codes were added to attribute sequences to each sample. The library preparations were sequenced on an Illumina Novaseq6000 platform and pairedend reads were generated. Then, the raw reads were processed using the NGS QC Toolkitv2.3.3 (Patel and Jain, 2012). The first five bases from the 5 end of the read were trimmed, and low quality bases > 20% or ambiguous bases > 1% were removed. Cleaned reads were obtained for the following analyses. The RNA-seq datasets generated are available from the NCBI SRA with accession number PRJNA687597.

Differentially Expressed Genes, Enrichment Analysis, and Gene Set Enrichment Analysis
The clean reads from the above RNA-Seq data were mapped to the optimized reference genome using HISAT2.2.4 with "-rnastrandness RF" and other parameters as default (Kim et al., 2015). The mapped reads of each sample were assembled using StringTie v1.3.1 in a reference-based approach and then counted and subsequently normalized to fragments per kilobase of transcript per million (FPKM) values to quantify expression abundance. The DESeq R package (1.10.1) (Love et al., 2014) was used to compare the differential expression between CO and HS groups. The P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR). An FDR < 0.05 and absolute fold change 1 were set as the threshold for significantly differentially expressed genes (DEGs).
Gene ontology (GO) analysis of the DEGs was implemented by the GOseq R packages (V1.10.0) (Young et al., 2010) based on the Wallenius non-central hypergeometric distribution. KOBAS software version 3.0.0  was used to test the statistical enrichment of DEGs in kyoto encyclopedia of genes and genomes (KEGG) pathways. A q-value < 0.05 was considered as a significant change.
Gene set enrichment analysis (GSEA) was performed using software GSEA and MSigDB (Subramanian et al., 2005) to identify whether a set of genes in specific pathways showed significant differences in two groups. Briefly, we input gene expression matrix and rank genes by the SinaltoNoise normalization method. | NES| > 1, NOM p-val < 0.05 and FDR q-val < 0.25 were considered as significant.

Alternative Splicing Analysis
Alternative Splicing analysis was based on the optimized reference genome. Identification and quantification of AS events including skipped exon (SE), alternative 5 splice site (A5SS), alternative 3 splice site (A3SS), mutually exclusive exon (MXE), and retained intron (RI) were carried out using rMATS (replicate multivariate analysis of transcript splicing) (Shen et al., 2014) (version 4.0.1) 1 with default parameters. Choosing an FDR of <0.05, significant splicing events were determined.

Validation Experiments
Quantitative real-time PCR (qPCR) analysis was used to verify the differentially expressed transcripts. Total RNA was isolated from coelomocytes of A. japonicus from the CO and HS groups as described as above. cDNA was synthesized using the PrimeScript TM RT reagent kit with gDNA Eraser (TaKaRa, Japan, Cat. No. RR047A) following the manufacturer's instructions. All transcript-specific primers for qPCR are listed in Supplementary  Table 1. SYBR R Premix Ex Taq (Tli RNaseH Plus) was used for qPCR (TaKaRa, Japan, Cat. No. RR420A). Each PCR reaction system consisted of 10.0 µl SYBR R Premix Ex Taq (Tli RNaseH Plus), 0.4 µl of each forward and reverse primer (10 µM), 0.4 µl ROX Reference Dye, 2 µl cDNA, and 6.8 µl dH 2 O in a total volume of 20 µl. Then, qPCR was performed on a StepOnePlus TM Real-Time PCR machine (Applied Biosystems, CA, United States). The standard procedure for amplification was as follows: 95 • C for 30 s and 40 cycles of 95 • C for 5 s, 57 • C for 30 s. β-actin and β-tubulin were used as control genes for internal standardization as previous reported Xu et al., 2018;Huo et al., 2019). The 2 − CT method was applied to analyze transcript expression levels. One-way ANOVA was conducted followed by Duncan's multiple tests with SPSS 19.0 software. Significant difference was set as P 0.05.
To verify the accuracy of AS events, total RNA was reversetranscribed into cDNA using a PrimeScript RT reagent kit (TaKaRa, Japan, Cat. No. RR047A) following the manufacturer's instructions. Transcript-specific primers were designed to 1 http://rnaseq-mats.sourceforge.net/index.html span the predicted splicing events using Primer 6 software (Supplementary Table 1). PCR conditions were 95 • C for 2 min and 35 cycles of 95 • C for 15 s, Tm for 15 s, 72 • C for 3 min, and 72 • C for 5-10 min that depended on the predicted product size. PCR products were monitored on 1% agarose gels stained by GelStain and verified by Sanger sequencing.

Apostichopus japonicus Transcriptome Sequencing Using PacBio Iso-Seq and Reference Genome Optimization
A total of 8,949,487 subreads with an average length of 2,443 were generated. After clustering, 479,677 CCSs were obtained, of which 263,473 full length non-chimeric (FLNC) sequences were further classified containing 5 primer, 3 primer, and poly (A) tails (Table 1). To further improve the sequence quality, the NFL reads were used to polish the above cluster consensus isoforms using Quiver software to obtain FL polished high quality consensus sequences (accuracy 99%). The length distribution of the subreads, CCSs, FLNC reads and polished consensus reads are shown in Supplementary Figure 1. A total of 20,162 polished high quality consensus sequences were applied for gene structure optimization of the A. japonicus reference genome, and 380 genes in the reference genome were optimized with detailed information shown in Supplementary Table 2.

Differential Alternative Splicing Events Analysis in Response to Heat Stress
In the present study, DAS events were identified by comparing control and HS groups. A total of 1,224 and 1,251 DAS events (Figure 1) were identified using the splice junction counts (JC only) and reads on target and junction counts (JCEC) as the input for rMATS. Skipping exon (SE) dominated in both count methods (61.68 and 61.23%, respectively). Detailed information of DAS events is shown in Supplementary Table 3.

Differentially Expressed Genes Analysis in Response to Heat Stress
A total of 1375 DEGs were identified between the control and 48 h HS groups. As shown in Figure 2, 671 genes were  up-regulated and 704 genes were down-regulated; FDR was <0.05 and | log2fold-change| 1). The expression of DEGs with FDR < 0.01 and | log2fold-change| 2.5 from the six samples were clustered into two distinct groups by hierarchical clustering (Supplementary Figure 2). Additional details of the DEGs are presented in Supplementary Table 4.

GO, Kyoto Encyclopedia of Genes and Genomes Enrichment Analysis of DEGs and GSEA
Gene ontology enrichment analysis showed that 543 DEGs were significantly enriched in 83 GO categories (Supplementary Table 5) (q-value < 0.05) including 127 DEGs assigned to "response to stress" (BP, GO:0006950), 77 DEGs assigned to "immune response and immune system process" (BP, GO:0006955 and GO: 0002376, respectively), and 23 DEGS assigned to "protein folding" (BP, GO:0006457). The top 20 GO terms are shown in Figure 3. KEGG enrichment analysis showed that 9 pathways were significantly enriched (q-value < 0.05) after 48 h of HS including "Protein processing in endoplasmic reticulum, " "Prostate cancer, " "PI3K-Akt signaling pathway, " "Antigen processing and presentation, " "Parathyroid hormone synthesis, secretion, and action, " "Spliceosome, " "Systemic lupus erythematosus, " "Thyroid hormone synthesis, " and "Herpes simplex infection" (Figure 4 and Supplementary Table 6  The results of GSEA further showed that two significantly enriched pathways were confirmed to be activated in response to HS: "Protein processing in endoplasmic reticulum" and "Antigen processing and presentation" (Figure 5).

Validation of RNA-seq Results and Representative DAS Gene
To validate the accuracy of the transcriptome quantification, 12 up-regulated and 7 down-regulated genes (all from the enriched KEGG pathways) were selected and quantified using qRT-PCR. The results of quantitative real-time PCR showed that all upand down-regulated genes showed the same trends as the highthroughput sequencing data (Supplementary Figure 3). Hsf1 was selected as the representative DAS gene, being the only heat shock transcription factor gene that was significantly alternatively spliced under HS. The differentially spliced exon was located at exon 10 ( Figure 6) resulting in two hsf1 transcript isoforms: Ajhsf1-1 and Ajhsf1-2. In the control group of sea cucumbers, the inclusion level of exon 10 was 86.3% but, by contrast, the exon inclusion level was only 7.3% in sea cucumbers under HS (Figure 6 and Supplementary Table 3).

DISCUSSION
Seasonal fluctuations in water temperature are unavoidable in semi-open aquaculture areas and, with global warming, it is predicted that suboptimal high temperatures may become more frequent and severe in coastal areas (Frölicher et al., 2018;Oliver et al., 2018). In recent years, extreme weather caused by high temperatures (up to 30 • C) has been considered to be an added factor for sea cucumber mortality syndrome in the summer . The HS response of sea cucumbers is a rapid and complicated process, potentially involving numerous genes and with differential responses to HS by different tissues Huo et al., 2019;Li et al., 2019). Coelomocytes, the major component of the non-specific defense system of sea cucumbers, are involved in a series of immune responses (Xue et al., 2015). It is believed that the molecular basis of the HS response in coelomocytes is more complicated than those in intestine, muscle, or respiratory tree. The current study presents the first analysis of transcriptional and post-transcriptional changes generated by AS of genes in response to HS (48 h at 30 • C) in sea cucumber coelomocytes.
Overall, a total of 1,375 DEGs identified in response to HS were significantly enriched in nine pathways including protein processing in endoplasmic reticulum, PI3K-AKT signaling pathway, antigen processing and presentation, and the spliceosome. Two of these were confirmed by GSEA to be activated by HS: (a) protein processing in endoplasmic reticulum and (b) antigen processing and presentation.

Heat Shock Response and Endoplasmic Reticulum Stress
The HSP family has been widely confirmed to play an important role as "assistants" that support proper folding, assembly, or disassembly of proteins and protein complexes and promote the proper refolding of abnormal or damaged proteins (Parsell and Lindquist, 1993;Velichko et al., 2013). Synthesis of heat shock proteins under HS conditions is one of the most important molecular responses to stress and is a highly conserved mechanism across phylogeny that helps to maintain cell homeostasis and protect cells and tissues from structural damage in most organisms, including aquatic ectotherms (Roberts et al., 2010). Expression of HSPs in animals exposed to HS has been well characterized in a wide variety of aquatic animals (Meng et al., 2009;Zhao et al., 2011;Xu et al., 2014Xu et al., , 2016. As observed in the present study, significant overexpression of several HSPs was detected in response to HS in sea cucumber coelomocytes. These included the inducible form and its constitutive counterpart of the HSP90 family members, hsp90ab1, hsp90b1, and hsp90aa1 (alias hsp90-alpha). The HSP90 protein family is essential to the folding and stabilization of cellular proteins and is implicated in a wide variety of cell processes as a response to stress (Jackson, 2013). These findings are in line with previous HS studies in aquatic animals including sea cucumbers (Liu et al., 2013;Xu et al., 2018;Huo et al., 2019;Li et al., 2019;Swirplies et al., 2019). Interestingly, sea cucumbers exposed to a moderate high temperature (26 • C) showed no changes in hsp90 expression after 24-72 h HS (Huo et al., 2019), which differs from our results, but suggests that the expression of hsp90 may vary in different tissues or has a trigger temperature somewhere in between 26 and 30 • C.
We also identified two cofactors, CDC37 and AHSA1, that were significantly up-regulated after 48 h HS. CDC37 is an accessory factor for HSP90, directing HSP90 to protein kinases as substrates and can also function as a molecular chaperone FIGURE 3 | Circular map of the top 20 gene ontology (GO) classifications of differentially expressed genes (DEGs) between 48 h_HS and CO groups. The first (outer) circle shows different Ontology: Biological Process (yellow) and Cellular Component (green). The second circle (orange) shows the background gene numbers enriched in each GO term, the higher the number, the longer the bar; the depth of color represents the q-value, the smaller of the q-value, the deeper the color. The third circle represents the number of up-regulated (dark purple) and down-regulated (light purple) DEGs enriched in each GO term. The fourth (inner, white) circle shows the rich factor value of each GO term (the number of DEGs enriched in this GO term divided by all the numbers of background genes enriched in this GO term). In the fourth circle, each grid of the background grid line represents 0.1. FIGURE 4 | Circular map of top 20 KEGG pathways enriched in the differentially expressed genes (DEGs) between 48 h_HS and CO groups. The first (outer) circle: different colors represent different KEGG A-classes: genetic information processing, environmental information processing, cellular processes, organismal systems, and human diseases. The second circle: the number shows the background gene numbers enriched in this pathway, the higher the number, the longer of the bar; depth of color represents the q-value, the smaller the q-value, the deeper of the color. The third circle: dark purple bar represents the number of up-regulated DEGs enriched in this pathway, and the light purple bar represents the number of down-regulated DEGs enriched in this pathway. The fourth (inner) circle shows the rich factor value of each pathway (the number of DEGs enriched in each pathway divided by all the numbers of background genes enriched in this pathway). In the fourth circle, each grid of the background grid line represents 0.1.  itself, interacting with other HSP90 co-chaperones to support protein folding (MacLean and Picard, 2003). This current finding is similar to the results of a HS study in rainbow trout (Huang et al., 2018). AHSA1 is also a HSP90 chaperone and stimulates the ATPase activity of HSP90 (Shao et al., 2016). It has been reported that HSP90 and AHSA1 are up-regulated in parallel under HS in sea cucumber intestine , which is consistent with our present study and suggests the potential cooperation of these two genes in response to HS.
In addition, seven members of the HSP40 family were identified and all were up-regulated in our present study, including DNAJA1, DNAJA2, DNAJB4, DNAJB6, DNAJB9, DNAJC3, and DNAJC7. These are all involved in protein translation, folding, unfolding, translocation, and degradation, primarily by stimulating the ATPase activity of HSP70 chaperones (Qiu et al., 2006). HSP70 is usually considered to be the central actor of the HS response, and elevated hsp70 mRNA expression upon HS has been reported in numerous aquatic animals including sea cucumber (Roberts et al., 2010;Li et al., 2019). Thus, in the present study, inducible forms of hsp70 were likely involved in cellular protein folding and remodeling processes within the coelomocytes of stressexposed sea cucumbers. Strong expression of these HSP genes induced by HS suggest that these family members contribute to "thermotolerance" for sea cucumber cells by providing temporary resistance to stress until potential longer-term adaptive mechanisms can be acquired. HSP70 also cooperates with other cellular chaperones including HSP90 and HSP40 and small heat shock proteins (e.g., HSP26) together constituting a dynamic and functionally versatile chaperone network for protein folding, unfolding, regulation, targeting, aggregation, and disaggregation, in response to HS (Rosenzweig et al., 2019).
In addition, we also identified the most significantly enriched KEGG pathway in the coelomocytes of sea cucumbers in response to HS. This was "protein processing in the endoplasmic reticulum (ER)" and all 41 genes belonging to this pathway were upregulated in comparison with controls, which is consistent with observations in other aquatic animals and other tissues in sea cucumber in response to HS (Guo et al., 2016;Qian and Xue, 2016;Li et al., 2017;Huang et al., 2018;Xu et al., 2018). It is well known that the ER has strict quality-control systems that ensure the correct folding of newly synthesized proteins and contributes importantly to maintaining cell homeostasis. An accumulation of unfolded or misfolded proteins in the ER induces the unfolded protein response (UPR), which is an evolutionarily conserved stress-detecting system used to adapt to changing environmental or stress conditions and, thereby, re-establish normal ER function (Raghunath et al., 2018). Our previous study reported the key ER chaperone response induced by HS in intestinal cells of sea cucumbers (Wang et al., 2021). A marked signature of the UPR stress response is the coordinated transcriptional up-regulation of resident ER chaperones such as the 78 kDa glucose-regulated protein (Grp78), heat shock protein 90 kDa beta (Grp94), and PDI family members (Lee, 2001;Ma and Hendershot, 2002;Dunlop et al., 2018). The overexpression of ER chaperones in our present study suggested that they may help to correct the folding and assembly of nascent proteins and prevent the transport of proteins or protein subunits that are folded incorrectly in sea cucumbers in response to HS. A second characteristic of the UPR is inhibition of protein synthesis under stress conditions (Prostko et al., 1992;Brostrom et al., 1996;Ma and Hendershot, 2002). In the present study, two eukaryotic translation initiation factor 2alpha kinase (EIF2AK) proteins (EIF2AK1 and EIF2AK3) were significantly overexpressed in our HS group. EIF2AKs share a common physiological role of phosphorylating eIF2α to inhibit protein synthesis at the level of translation initiation in response to various stress conditions, including oxidative stress, heme deficiency, osmotic shock, and heat shock (Krishna and Kumar, 2018). This indicates another conserved regulatory mechanism present in sea cucumbers and activated in response to HS.
Through a process called ER quality control, proteins that do not mature properly are retained in the ER and eventually targeted for ER-associated degradation (ERAD), which is the third major pathway of the UPR and disposes of misfolded proteins from the ER lumen in response to stress (Brodsky et al., 1999;Nishitoh, 2012;Dunlop et al., 2018). In A. japonicus, we also observed a significant up-regulation of endoplasmin in response to HS, an important protein of the ERAD process (Christianson et al., 2008). ERAD also involves chaperones of the HSP family that act in concert with ubiquitin ligases (Goldberg, 2003).
Hence, in response to HS, the expression of genes from the HSP70, HSP90, and HSP40 families were upregulated in sea cucumber coelomocytes, likely to promote degradation of misfolded proteins and to maintain cell homeostasis. These results provide important evidence of an evolutionarily ancient role for the ER-induced UPR under HS in echinoderms. Further research will be required to confirm this in the future. Additionally, the hypoxia-induced gene, hyou1, was also upregulated by HS and is involved in "protein processing in the ER" pathway. This suggests that thermal and hypoxia tolerance may positively related with each other (Pörtner and Peck, 2010), as has previously been reported for heat-stressed rainbow trout (Huang et al., 2018).

Immune Response
Several previous studies have focused on immune responses within different high temperature stress and/or duration conditions in sea cucumber tissues (intestine and muscle) (Huang et al., 2018;Huo et al., 2019;Li et al., 2019). Here we reported that three enriched pathways were linked to the immune system in response to HS in sea cucumber coelomocytes: "Antigen processing and presentation, " "Systemic lupus erythematosus, " and "Herpes simplex infection." GSEA confirmed the activated pathway "antigen processing and presentation" in response to HS, which was also reported in muscle tissues of sea cucumbers in response to HS . The GO/pathway analysis also indicated that several HSP transcripts were associated with this pathway (e.g., hsp70, hsp90α) in response to HS. It has been reported that HSPs in the cytosol such as HSP70 and HSP90, are involved in antigen presentation and crosspresentation by delivering chaperones, antigenic peptides to MHC-I molecules (Tsan and Gao, 2009). Our results suggest a potential role of HSPs in the immune system of sea cucumber coelomocytes under HS. Similar results have also been reported in some fish species such as rainbow trout Oncorhynchus mykiss (Huang et al., 2018) and Atlantic salmon Salmo salar (Beemelmanns et al., 2021).
The complement system is a major effector of the immune response and triggers many immunoregulatory mechanisms (Carroll, 2008). It plays a fundamental role in innate immune responses in invertebrates, and echinoderms have been confirmed to have a simplified complement system (Smith et al., 2010). In the current study, the complement component C3-2 (AjC3-2) was significantly down-regulated after 48 h HS (30 • C), which is contrary to observations in sea cucumber intestine after 48 h HS at 26 • C, suggesting that 30 • C stress may be a physiological upper limit temperature for sea cucumbers and lead to impaired immune function of coelomocytes under HS.

Signal Transduction
When ectotherms are exposed to HS, mitochondrial respiration rate is increased and this can result in oxidative stress and cellular damage (Heise et al., 2006;Cheng et al., 2015;Banh et al., 2016). In our study, heat stressed sea cucumbers appeared to have an enhanced oxidative stress response. Members of the forkhead box, class O (FoxO) transcription factor family are important regulators of the cellular stress response and promote antioxidant defenses (Klotz et al., 2015). The vascular endothelial growth factor receptor (VEGFR) has been reported to be activated under oxidative stress and trigger downstream cell survival signaling cascades (Wang et al., 2004). The cyclic AMPresponsive element-binding protein 3-like protein 4 (CREB3L4) is a transcription factor that responds to oxidative stress signals (Pregi et al., 2017). These three up-regulated genes that are potentially involved in oxidative stress responses were identified in sea cucumber coelomocytes after HS, suggesting their potentially important functions in the regulation of signaling to activate responses that can deal with injury caused by HS.

Post-transcriptional Regulation
The spliceosome is a large ribonucleoprotein complex involved in splicing the coding regions (exons) of precursor messenger RNAs (pre-mRNAs) with the concomitant excision of intervening non-coding introns. In sea cucumber coelomocytes, 15 genes were significantly enriched in the spliceosome pathway after HS. Similar results were found in other aquatic animals (Qian and Xue, 2016;Huang et al., 2018). Serine-arginine-rich (SR) and heterogeneous nuclear ribonucleoprotein (hnRNP) proteins are best known as the major splicing factors involved in AS regulation (Lee and Rio, 2015). In our study, three hnRNP genes (hnRNPA1-like, hnRNPA3, and hnRNPK-like), and three SR genes (srsf2, srsf3-like, and srsf7) were identified as responding to HS. Most of these two classes of genes varied oppositely in sea cucumber coelomocytes in response to HS. Indeed, previous studies have observed that these two classes of genes often have antagonistic effects on splice-site usage. Exonic splicing enhancers (ESEs) are often bound by SR proteins whereas exonic splicing silencers (ESSs) are typically bound by hnRNP proteins (Matera and Wang, 2014). Interestingly, a negative correlation between hnRNPK expression (a significant decrease) and hsp70 expression (a significant increase) was also observed in response to HS in the present study. The hnRNPK protein has also been reported to play a role in the heat shock response of cells and control expression of hsp70 mRNAs by inhibiting heat shock-induced transcriptional activity (Kim et al., 2017). Our present findings suggest that these splicing factors may organize a regulatory network of various groups of genes through AS as well as splicing inhibition to respond appropriately to a particular stress.
Our differential splicing analysis further indicated that AS may be an important potential mechanism used by sea cucumbers to deal with heat shock. AS has been reported to play a key role in enhancing regulatory capacities and proteomic complexities in eukaryotes (Ben-Dov et al., 2008) and allows organisms to rapidly modulate protein functions in response to physiological changes. Accumulating evidence indicates that AS is a key player in the response to environmental stresses including HS (Fujikake et al., 2005;Liu et al., 2013). In the DAS list, we identified three hnRNP genes (hnRNPA1-like, hnRNPA3, and hnRNPK-like) and four SR genes (srsf2, srsf3-like, srsf4-like, and srsf7). Previous studies have reported that temperature stress can dramatically alter the AS of pre-mRNAs of several SR and hnRNP genes in catfish (Tan et al., 2019).
Increasing the expression of heat shock proteins (HSPs) is an evolutionarily conserved strategy for eukaryotes to deal with HS. This process is mediated by a master regulator: the heat shock transcription factor, HSF1. HSF1 exists as a monomer in unstressed cells but under stress conditions forms an active trimer that can bind to promoter regions of heat shock genes (Anckar and Sistonen, 2011;Hentze et al., 2016). Tan et al. (2019) reported that the inclusion level of exon 9a in catfish hsfs1 decreased significantly from 100 to 2.7% after heat shock. In accordance, we also found that the inclusion level of exon 10 in Ajhsf1 was significantly reduced from 86.3 to 7.3% after HS. Research on mammals demonstrated that HSF1 isoform ratios are a major factor in the regulation of expression of heat shock protein genes and exon splicing. The HSF1 isoform may be produced to facilitate the activation of the homotrimer and promote the expression of HSPs after HS (Neueder et al., 2014).

CONCLUSION
The present study focused on transcriptional and posttranscriptional (AS) regulatory levels, combining PacBio Iso-seq and RNA-seq to reveal acute responses to HS in sea cucumber coelomocytes. Strong numbers of up-regulated DEGs were enriched in the protein processing function of the endoplasmic reticulum, providing key evidence of the involvement of an evolutionarily conserved signal transduction mechanism and stress-detecting system in response to HS in echinoderms. Further analysis identified an abundance of RNA-splicing related genes that displayed AS and confirmed that the inclusion level of splicing out an exon in Ajhsf1 was significantly decreased after HS, suggesting that AS is an important potential mechanism used by sea cucumbers in response to HS. Although satisfactory answers to these findings obtained from high-throughput screening will require further functional confirmation, our study improves understanding and opens up new perspectives on the regulatory mechanisms used by echinoderms to respond to environmental stress.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI [accession nos: PRJNA785124, PRJNA757529, and PRJNA687597].

ETHICS STATEMENT
This animal study was reviewed and approved by the Experimental Animal Ethics Committee of the Ocean University of China.

AUTHOR CONTRIBUTIONS
MC originally conceived the study. YW, YY, and XC conducted the experiments. MC, YW, and YY contributed to the acquisition of data. MC, YW, YY, and KS contributed to the interpretation of the data. MC wrote the first draft of the manuscript. KS revised the manuscript. All authors contributed to manuscript read and approved the submitted version.  Supplementary Table 2 | Detailed information of optimized genes in the reference genome based on PacBio Iso-Seq.
Supplementary Table 3 | Detailed information of differential alternative splicing (DAS) events.