Transcriptome Sequencing Reveals the Response to Acute Thermal Stress in the Paci ﬁ c Abalone, Haliotis discus hannai

,


Introduction
Temperature is one of the main physical factors affecting marine biological ecosystems [1].Changes in seawater temperature can directly determine the distribution range and regulate all aspects of marine animals, including their growth and development [2][3][4][5], reproduction [6], metabolism, as well as immune responses [7][8][9].The Pacific abalone Haliotis discus hannai is the main abalone species currently cultivated in China [10] and shows high economic and nutritional value.As a kind of coldwater shellfish, it is naturally distributed in the Korean Peninsula, the northern waters of Japan, the coastal areas of the Shandong Peninsula, and the Liaodong Peninsula of China.
With the vigorous development and continuous southward movement of the abalone aquaculture industry, Fujian Province has become the main region for cultivating Pacific abalone in southern China since 2000.According to statistics, the annual output of cultured abalone in China had reached 203,000 tons in 2020, including 155,000 tons from Fujian Province, accounting for 76.4% of the country's total output (China Fisheries Yearbook).However, the mass-scale culture of Pacific abalone in the south region of China has been facing the great challenge of high mortality in summer [11,12].This is because the temperature in these areas is significantly higher than that in the natural habitats of this species [13].The optimum temperature on the north coast of China was 4-26°C; when the temperature was higher than 26°C, the activity decreased and the growth stopped [14], while when the temperature was higher than 28°C, the mortality rate increased significantly [15].Elevated seawater temperature not only increases the abundance of pathogens such as marine Vibrio spp.but also reduces the innate immunity of marine animals to Vibrio spp., increasing their risk of infection and disease outbreaks [14,[16][17][18].
With the rapid development of next-generation sequencing technologies, low-cost high-throughput sequencing has become possible [19].The transcriptome sequencing (RNAseq) technology, which not only enables rapid generation of transcript data for studying species but also increases the depth and coverage of sequencing, has become a common method for studying nonmodel organisms lacking genomic information [20].Recently, a growing number of economic aquatic organisms have been analyzed using this technique, elucidating the impact of external environmental changes on biological processes and the stress response of target species to various stimuli [21][22][23].There have been several reports on the adaptation of shellfish to temperature changes.For example, transcriptome sequencing technology has been used to study the transcriptional responses of the Manila clam Ruditapes philippinarum [24] and the pearl oyster Pinctada fucata martensii [25] to cold stress to reveal the underlying mechanisms.Besides, the Pacific oyster Crassostrea gigas is often exposed to air and can adapt to severe temperature changes.Zhu et al. [26] have identified some key candidate genes, which encode heat shock proteins (HSPs) and apoptosis-related proteins by analyzing their transcriptome data under low-and high-temperature stress.During the hot summer months, the kuruma shrimp Marsupenaeus japonicus in the aquaculture system often suffers from mass die-offs.Zheng et al. [27] have conducted comparative transcriptome analysis of immune response, HSP, antioxidant system, and metabolic changes in kuruma shrimp under heat stress to explore the potential molecular mechanism of summer death.
Although RNA-seq technology has been widely used for research, the molecular mechanism of Pacific abalone's response to heat stress remains unclear.In recent studies, 26 HSP genes from four tissues (blood cells, gill, mantle, and muscle) of Pacific abalone were differentially expressed under cold and heat stress conditions [28].Yu et al. [29] found that under heat stress conditions, Pacific abalone in Dalian, northern China tends to maintain the synthesis and metabolism of macromolecules such as proteins in cells by activating transcription factors, thereby expressing various regulatory genes and protein folding-related genes.Shellfish, as invertebrates, lack an acquired immune system and only have a nonspecific immune system.Its hemolymph is a unique body fluid, and its functions are similar to human blood in many aspects, in which blood cells play an important functional role in hemolymph, and participate in immune responses as multifunctional hemolytic cells, mainly including phagocytosis, packaging, hemagglutination, blackening, etc. [30].Hence, in this study, the hemolymph of Pacific abalone was selected as the target tissue, and the key candidate genes related to high-temperature stress response were screened by transcriptome sequencing of the hemolymph tissue after heat stress.This study can deepen the understanding of the molecular mechanism underlying the response of Pacific abalone to high temperature changes.

Materials and Methods
2.1.Pacific Abalone Samples.The Pacific abalone samples were collected on May 14, 2021 at Changqing Ocean Technology Co., Ltd. in Weihai City (Shandong Province, China).The average shell length of Pacific abalone was 49.07 AE 3.81 mm.After being transported to the laboratory, all abalones were temporarily cultured in a water tank, with the temperature set at 17°C and supply of oxygen and the salinity of the seawater is 35.They were fed with an appropriate amount of salted kelp every afternoon.In each next morning, a sewage pump was used to clean up the residual bait and feces and replace two-thirds of the water body with fresh one.Subsequent experiment was carried out after a week of temporary feeding.

Determination of Feeding Rate and Oxygen Consumption
Rate of Pacific Abalone at Different Temperatures.To determine the optimal temperature for heat stress experiment, we measured the changes in feeding rate (FR) and oxygen consumption rate (OR) of abalone at 7, 14, 21, and 27°C (Figure 1).It was found that as the temperature exceeded the optimum growth temperature of the abalone, its FR and OR decreased.

Aquaculture Research
Correlation fitting analysis found that FR, OR, and temperature (T) fit a cubic function relationship: FR = (−1.2× 10−5) T 3 + 0.00047T 2 −0.0038T + 0.046 (R 2 = 0.709); OR = (−2.8× 10 −4)T 3 + 0.01166T 2 −0.111T + 0.464 (R 2 = 0.850).In summary, the final decision was made to set the temperature of the thermal stress group at 28°C and the control group at 17°C.The formulas for calculating feeding rate and oxygen consumption rate were as follows [31]: In the formula, F represents the food intake (g/(ind•d)), F 1 is the dry weight of the kelp fed (gram, g), F 2 is the dry weight of the remaining kelp (gram, g), t 1 is the feeding time (day, d), and N is the number of abalone in the experimental group; FR represents the feeding rate (%) and G is the average wet weight of the abalone in the experimental group (gram, g); OR stands for oxygen consumption rate per body weight (mg/(g)), DO 0 is the initial dissolved oxygen in water (mg/L), DO t is the dissolved oxygen in the water after the experiment (mg/L), V is the water volume (liter, L) when measuring oxygen consumption, W is the dry weight of the abalone soft body (gram, g), and t 2 is the duration of the experiment (hour, hr).

Thermal Stress
Treatment and RNA Extraction.After temporary culture, abalones were grouped and transferred to new breeding tanks with the tank temperatures set to 17°C (control group) and 28°C (thermal stress group), respectively.At 6, 24, and 48 hr of the experiment, a total of nine live abalones were randomly selected from each of two groups, respectively.The foot of each abalone was cut with a sterilized scalpel blade, and the hemolymph was extracted using a 1 ml syringe.The hemolymph of three abalone individuals was mixed in a 10 ml cryogenic tube and taken as one sample, and three biological replicates were set for each group.According to the sample processing requirements of Shanghai Biozeron Co., Ltd., the hemolymph and Trizol reagent (Invitrogen, Waltham, USA) were mixed at a ratio of 1 : 3, shaken with a vortex oscillator until the white flocs disappeared, and then immediately placed in liquid nitrogen for freezing.The frozen sample was then stored at −80°C for subsequent RNA extraction.Total RNAs were extracted from the hemolymph using a traditional Trizol reagent (Invitrogen, Waltham, USA) according to the manufacturer's instructions, and potential genomic DNA contamination was removed using DNase I (Takara, San Jose, USA) reagent.RNA degradation and contamination were detected by 1% electrophoresis gel.Then, the integrity of RNA was assessed by Agilent 2100 bioanalyzer, and the purity and concentration of RNA were detected by NanoDrop-2000 spectrophotometer.

Library Preparation and Illumina Hiseq Sequencing.
RNA-seq transcriptome libraries were created using 1 g of total RNA and the TruSeq TM RNA Sample Preparation Kit from Illumina (San Diego, CA).Shortly, oligo (dT) beads were used to isolate messenger RNA, which was then fragmented with fragmentation buffer.Illumina's methodology was followed for cDNA synthesis, end repair, A-base addition, and ligation of the Illumina-indexed adaptors.Libraries were amplified for 15 cycles using Phusion DNA Polymerase (New England Biolabs, Ipswich, USA), and then screened 200-300 bp cDNA target fragment with 2% super-agarose gel.Paired-end libraries were sequenced by Illumina NovaSeq 6000 sequencing (150 bp * 2, Shanghai Biozeron Co., Ltd.) after being quantified by TBS 380.

Differential Expression Analysis and Functional
Enrichment.The expression level of each gene was estimated by calculating the fragments per kilobase of exon per million mapped reads (FRKM).The R statistical package edgeR (Empirical Analysis of Digital Gene Expression in R, version 3.40.2;https://bioconductor.org/packages/release/bioc/html/ edgeR.html)was used to screen the differentially expressed genes (DEGs) and analyze their expression according to the fold change and P-value of the expression change.The selection standard is |log2 (fold change)| > 2, P-value < 0.05 and false discovery rate (FDR) < 0.05.The analysis steps of Gene Ontology (GO) included mapping all DEGs to each term in the GO database, calculating the number of genes in each term, and finding the terms significantly enriched.Then, the DEGs contained in different levels of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway were counted to determine the major metabolic pathways and signal transduction pathways involved.The DEGs between two samples were chosen according to the following criteria: logarithmic fold change of greater than 2 and FDR lower than 0.05.Goatools (https://github.com/tanghaibao/Goatools)and KOBAS (http://bioinfo.org/kobas)were used to perform GO functional enrichment and KEGG pathway analyses to better understand the functionalities of DEGs.When the Bonferroni-corrected P-value for DEGs was lower than Aquaculture Research 0.05, they were considered highly enriched in GO keywords and metabolic pathways.

Validation of RNA-Seq
Data by Quantitative Real-Time PCR.Several interesting genes in the thermal stress group (6 hr), including HSPs, calreticulin (CALR), and immunerelated genes, were selected for validation by quantitative real-time PCR (qRT-PCR), using the same samples used for RNA-seq.The cDNA of samples was synthesized using HiScript III RT Super Mix for qPCR (Vazyme, China) with 1 μg of total RNA, and the cDNAs were diluted 10 times for using as amplification template.qRT-PCR was performed on the Applied Biosystems StepOnePlus™ Real-Time PCR System (Applied Biosystems, USA) with ChamQ SYBR Color qPCR Master Mix (Vazyme, China) according to the manufacturer's instructions.Three biological and three technical replicates were set, with β-actin [32] serving as the internal standardization reference.The specific primers for qRT-PCR are listed in Table 1.The PCR cycling conditions were as follows: 95°C for 30 s, then 40 cycles of 95°C for 10 s, and 60°C for 30 s. Finally, the 2 −ΔΔCt method was used to calculate the relative expression of each gene.

Quality Control Evaluation and Reference Genome
Alignment.A total of 595,226,452 raw reads were obtained by transcriptome sequencing (Table 2), and 515,533,974  2(d)).Many DEGs of interest were used to draw heat maps to observe the expression changes of these genes in different heat stress stages (Figure 3).

GO and KEGG Enrichment Analysis of DEGs. GO and KEGG enrichment analysis of DEGs in Pacific abalone
hemolymph after heat stress were carried out, respectively (Table 3).The GO enrichment analysis results of DEGs were consistent after 6, 24, and 48 hr of heat stress.The DEGs were mainly involved in cellular process (GO: 0009987), biological regulation (GO: 0065007), metabolic process (GO: 0008152), multicellular organismal process (GO: 0032501), response to stimulus (GO: 0050896), developmental process (GO: 0032502), and other biological processes that maintain cell homeostasis.The results of KEGG enrichment analysis showed that DEGs after 6 hr of thermal stress were significantly enriched in antigen presentation and processing, female signaling pathway, and endoplasmic reticulum protein processing pathway in genetic information processing in organic systems.The involved genes were heat shock protein 70 gene (HSP70), HSP90, glucose-regulated protein 78 gene (GRP78), endoplasmic reticulum protein 29 gene (ERP29), CALR, etc.In terms of metabolism, DEGs after 6 hr of thermal stress were significantly enriched in cysteine and methionine metabolic pathways.The KEGG enrichment analysis results of TH 48 hr group and TH 24 hr group were similar to those of TH 6 hr group.Many metabolism-related DEGs were significantly enriched in tyrosine metabolic pathways after 24 hr of heat stress and were significantly enriched in cysteine and methionine metabolic pathways after 48 hr of heat stress.Additionally, immune-related Toll-like receptor 6 (TLR6) and myeloid differentiation factor 88 (MyD88) genes of Toll and Imd signaling pathways in organic systems were also observed to be significantly upregulated after 6 hr of thermal stress.
3.4.qRT-PCR Validation Analysis.The 13 candidate genes that were significantly upregulated or downregulated and involved in different biological processes were selected for qRT-PCR analysis to verify the accuracy of DEGs as revealed by RNA-seq (Figure 4).The qRT-PCR results showed that although the relative expression levels of the 13 candidate genes were slightly different from the transcriptome results, the overall upregulated and downregulated trend was consistent, which verified the accuracy of the transcriptome sequencing data.

Discussion
As a key abiotic factor for the high mortality of Pacific abalone in summer, high temperature restricts the development of abalone aquaculture industry in southern China.Therefore, it is necessary to understand the molecular regulatory mechanism of heat stress response in Pacific abalone and further analyze its thermal adaptation capacity.In this study, the transcriptome changes of Pacific abalone after acute heat stress at 28°C were explored by RNA sequencing in order to screen genes related to thermal stress response and the underlying regulatory mechanism.Here, 2,213, 2,337, and 1,420 DEGs were identified in the thermal stress group after 6, 24, and 48 hr, respectively, as compared with the control group without heat stimulation.Further, GO and KEGG enrichment analysis of DEGs suggested that Pacific abalone may maintain cellular homeostasis through the upregulated expression of molecular chaperones associated with the endoplasmic reticulum-associated degradation (ERAD) pathway.
In addition, the expression of genes of TLR signaling pathways such as TLR6 and MyD88 was also induced.Thus, Pacific abalone may activate the immune defense system by upregulating immune-related genes to defend the invasion of pathogenic bacteria after heat stress.

Effects of Thermal Stress on Protein Folding and
Degradation in Pacific Abalone.Maintenance of protein homeostasis is necessary for cells to respond to environmental stresses, and HSPs play an important role in this process.HSPs are molecular chaperones, which are highly conserved in terms of sequence among different species.When organisms are subjected to environmental stress, HSPs can promote the folding of cellular proteins and the degradation of misfolded proteins to prevent the accumulation of unfolded proteins and maintain cellular homeostasis, thus protecting the body from damage [33][34][35].In recent years, there have been studies on the significant expression of HSPs in shellfish under heat stress.The induction of HSPs seems to be an important indicator for monitoring the state of intracellular 6 Aquaculture Research heat stress.The studies on greenlip abalone Haliotis laevigata [36], pearl oyster [37], Pacific oyster [38], Hong Kong oyster Crassostrea hongkongensis [39], Manila clam [40], Yesso scallop Patinopecten yessoensis [41], razor clam Sinonovacula constricta [42], Pacific abalone [28,43], and other related studies have strongly proved that the induction of HSPs is a biological self-protection mechanism in response to heat stress, and the upregulation of HSPs can help cells resist damage and apoptosis caused [44].The present study found that HSPs such as HSP70A, HSP70B, HSP90A, HSP90B, and HSC70-4 (heat shock protein 70 cognate 4) were significantly upregulated after 6, 24, and 48 hr of thermal stress, indicating that HSPs are also the thermal stress markers of Pacific abalone.
Moreover, in the present study, the expression of genes related to endoplasmic reticulum protein-processing pathway (GRP78, ERP29, CALR, and PDIA6) in Pacific abalone hemolymph was significantly upregulated after heat stress.These genes play important roles in the unfolding protein response (UPR) and ERAD pathways and can perform

FIGURE 3:
The heat map of the expression of many differentially expressed genes (DEGs) in Haliotis discus hannai hemolymph after thermal stress (6, 24, and 48 hr).The expression of these genes reached the highest value after thermal stress (6 hr) and then gradually decreased with the extension of thermal stress time.
protein quality control and maintain ER protein homeostasis [45].High-temperature treatment results in ER stress that activates the UPR pathway.UPR enhances protein folding ability by increasing the synthesis of endoplasmic reticulum protein chaperones [46].ERAD pathway can identify, locate, and eliminate unfolded or misfolded proteins [47], maintaining the dynamic balance of protein folding in the ER.As the dominant player of UPR in normal intracellular environment, GRP78 could protect cells from damage caused by harmful chemicals, oxidative stress, calcium depletion, programmed cell death, etc. by forcing unfolded proteins to refold or degrade through cellular degradation mechanisms [48].Besides, GRP78 may play an important role in the defense of the giant tiger prawn to environmental stress [49], as well as in response of the giant river prawn to high temperature and pathogen stress [50].ERP29 has been speculated to be a chaperone protein, which is mainly involved in the processing and secretion of secreted proteins [51].CALR is the binding partner of calcium ions in the ER, mainly involved in calcium storage and signal transduction [52].PDIA6 is a member of the protein disulfide isomerase (PDI) family, which acts as a molecular chaperone in the folding and repair of faulty proteins during stress responses [53].It has been observed that the expression of PDIA6 is significantly upregulated after heat stress in the Pacific oyster [54] and razor clam [42].In conclusion, maintaining cellular homeostasis through the upregulated expression of molecular chaperones associated with ERAD pathways may be central to the defense of Pacific abalone against thermal stress.4.2.Inducing Immune Defense System after Thermal Stress in Pacific Abalone.Rising seawater temperature not only leads to increase in the abundance of pathogens such as marine Vibrio spp., but also reduces the natural immunity of marine animals, resulting in easier bacterial infection.It has been reported that thermal stress negatively affects the immune  [56], the Taiwan abalone (Haliotis diversicolor supertexta) [57], and Pacific abalone [58].In this study, DEGs were found significantly enriched in immune-related signaling pathway of antigen processing and presentation in hemolymph tissue after high-temperature challenge.HSP family genes, including HSP90A, HSP90B, HSP70A, and HSP70B, were involved in this pathway.In addition to acting as molecular chaperones, HSP90 proteins are also necessary for the assembly, activation, and maintenance of nucleotide-binding oligomerization domain-like receptor (NLR) protein complexes [59].These genes can participate in the turnover control of NLRs to maintain proper NLR protein levels and avoid immune deficiency caused by excessive accumulation of immune receptors [60].HSP70 proteins have been confirmed to be induced by pyrogen and Vibrio and play a crucial role in immune response and regulation of TLR signaling pathway [61].The results were consistent with those reported in the Pacific oyster [26] and razor clam [42], in which the upregulation of HSP70 and HSP90 family genes has been also observed after heat treatment.
In addition, the expression of other immune-related genes such as TLR6 and MyD88 gene of TLR signaling pathway was also induced in the Pacific abalone under thermal stress.Since TLRs were first discovered to play an important role in the immune response against fungal infection in Drosophila [62], a large number of TLRs have been identified as important genes for recognizing pathogen structures [63].TLRs play a crucial role in activating the immune system by sensing pathogen-associated molecular patterns of viruses and Vibrio spp.[64] and their family members are conservative during evolution [65].As a member of TLRs, TLR6 plays an important role in the immune defense of mollusks.Studies on the Pacific oyster [66] and Pearl oyster [67] have shown that TLR6 expression is markedly upregulated by bacterial infection.MyD88, as a key connector protein in the TLR signal transduction pathway, can stimulate downstream signal transduction and induce the expression of proinflammatory cytokines or antiviral genes, thereby participating in the innate immune process.MyD88 and TLRs can interact each other [68].Research on miiuy croaker (Miichthys miiuy) has found that TLR2 plays an indispensable role in the transcription of MyD88 [69].MyD88 has been considered to represent the intracytoplasm partner of TLRs, since the expression of MyD88 and TLR1 is both upregulated in Mediterranean mussel (Mytilus galloprovincialis) after injection with filamentous fungus (Fusarium oxysporum) [70].In the present study, HSP70, TLR6, and MyD88 were highly expressed in the hemolymph of Pacific abalone after thermal stress, indicating that the Pacific abalone might upregulate some TLR signaling pathway genes to resist pathogen infection after temperature stress.In the study of Ren et al. [61], the expression of MyD88 and TLR13 is significantly reduced after silencing HSP70 by RNA interference, confirming that HSP70 plays an important role in regulating the TLR signaling pathway.These results may support that thermal stress could activate the immune defense system in the Pacific abalone.

Conclusions
In this study, the transcriptome of Pacific abalone under thermal stress was sequenced and analyzed.A total of Aquaculture Research 2,213, 2,337, and 1,420 DEGs were identified in Pacific abalone after 6, 24, and 48 hr of thermal stress in 28°C, respectively, as compared with the control group.We found HSPs such as HSP70A, HSP70B, HSC70-4, HSP90A, and HSP90B were significantly upregulated after thermal stress, while the expression of genes related to endoplasmic reticulum protein-processing pathway (GRP78, ERP29, CALR, and PDIA6) was significantly upregulated.It indicated that Pacific abalone against thermal stress by maintaining cellular homeostasis through the upregulated expression of molecular chaperones associated with ERAD pathways.In addition, DEGs were found significantly enriched the immune-related signaling pathway of antigen processing and presentation in hemolymph tissue after high-temperature challenges.The immune-related genes like HSP70, TLR6, and MyD88 were highly expressed in the hemolymph of Pacific abalone after thermal stress.This revealed that the Pacific abalone might activate the immune defense system to resist the invasion of pathogenic bacteria by upregulating immune-related genes after thermal challenges.This study preliminarily explored the molecular mechanism of thermal stress adaptation in Pacific abalone and provided a reference for further study of the response patterns of related genes and pathways in Pacific abalone under high-temperature stress.

FIGURE 1 :
FIGURE 1: Ingestion rate (a) and oxygen consumption rate (b) of Haliotis discus hannai at different temperatures.Different letters among groups indicate significant differences (P <005).

FIGURE 2 :
FIGURE 2: The diagram and Venn diagram of DEGs in hemolymph of Haliotis discus hannai after thermal stress (FDR ≤ 0.05, |log2FC| > 1).(a) The volcano diagram of DEGs after 6 hr of thermal stress.(b) The volcano diagram of DEGs after 24 hr of thermal stress.(c) The volcano diagram of DEGs after 48 hr of thermal stress.(d) The Venn diagram of DEGs in different stages after thermal stress.

FIGURE 4 :
FIGURE 4:  The relative expression levels of 13 genes as revealed by qRT-PCR to confirm the accuracy of transcriptome data.

TABLE 1 :
The primer information of genes detected by qRT-PCR.

TABLE 2 :
Analysis results of transcriptome sequencing data of the hemolymph of Haliotis discus hannai.obtained after filtering the raw data.The total sequencing depth of this transcriptome was 276x.Among them, 17.32 G clean bases and 127,129,600 clean reads, 23.22 G clean bases and 168,733,690 clean reads, and 15.17 G clean bases and 110,991,322 clean reads were identified in the Pacific abalone in the three different stages of thermal stress (6, 24, 48 hr), respectively.Note that 14.81 G clean bases and 108,679,362 clean reads were obtained for the control group.The contents of Q20 and Q30 bases of clean reads of each group were calculated.It was found that the lowest contents of Q20 and Q30 bases were 97.29% and 90.08%, respectively.The clean reads after quality control were mapped to the genome of Haliotis discus hannai, and the mapping rate was in a range of 85.99%-88.17%.

TABLE 3 :
Significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in the hemolymph of Haliotis discus hannai under thermal stress.