Transcriptomic Response of the Liver Tissue in Trachinotus ovatus to Acute Heat Stress

Simple Summary In this study, the transcriptomic response of Trachinotus ovatus liver was observed against high-temperature stress. Through differential expression and short time-series expression miner (STEM) analyses, some high-temperature-related genes and biological pathways were screened, which were mainly related to protein balance, hypoxia adaptation, and energy metabolism. Our results suggest that protein dynamic balance and function, hypoxia adaptation, and energy metabolism transformation are crucial in response to acute high-temperature stress. These results contribute to understanding the molecular response mechanism of T. ovatus under acute heat stress and provide novel insights into the selection and breeding of heat-tolerant cultivars and the high-quality development of aquaculture. Abstract Trachinotus ovatus is a major economically important cultured marine fish in the South China Sea. However, extreme weather and increased culture density result in uncontrollable problems, such as increases in water temperature and a decline in dissolved oxygen (DO), hindering the high-quality development of aquaculture. In this study, liver transcriptional profiles of T. ovatus were investigated under acute high-temperature stress (31 °C and 34 °C) and normal water temperature (27 °C) using RNA sequencing (RNA-Seq) technology. Differential expression analysis and STEM analysis showed that 1347 differentially expressed genes (DEGs) and four significant profiles (profiles 0, 3, 4, and 7) were screened, respectively. Of these DEGs, some genes involved in heat shock protein (HSPs), hypoxic adaptation, and glycolysis were up-regulated, while some genes involved in the ubiquitin-proteasome system (UPS) and fatty acid metabolism were down-regulated. Our results suggest that protein dynamic balance and function, hypoxia adaptation, and energy metabolism transformation are crucial in response to acute high-temperature stress. Our findings contribute to understanding the molecular response mechanism of T. ovatus under acute heat stress, which may provide some reference for studying the molecular mechanisms of other fish in response to heat stress.


Introduction
Aquaculture is one of the fastest-growing sectors of the food industry worldwide and provides about 50% of animal protein for humans [1], which is predicted to grow to 53% by 2030 [2]. In the past few years, the aquaculture industry has developed rapidly worldwide, especially in China. However, there still are problems hindering the healthy development of the aquaculture industry, such as heat stress. The earth's climate is rapidly warming, and spells of high temperature are becoming more frequent and severe, resulting in an overall higher rate of high-temperature stress for fish [3]. High temperatures can also accelerate fed twice a day with commercial pelleted feed, and 30-40% of the seawater in each tank was regularly changed every day. Three parallel groups were set up in the experiment, and three temperature gradients, including 27 °C (control group, CT), 31 °C (high-temperature group, HT31), and 34 °C (high-temperature group, HT34), were set up in each parallel group. Before the heat stress experiment, all experimental fish were fasted for 24 h and were not fed during the experiment. After the beginning of the experiment, water temperature in three parallel groups was gradually increased from 27 °C to 34 °C at 1 °C/h and maintained at 34 °C for 2 h using heating rods (WN-B15). Three fish were selected randomly from each parallel group at 27 °C, 31 °C, and 34 °C (kept for 2 h) and anesthetized with phenolic. After measuring body length and weight, fish were dissected, and their liver tissues were removed and put into liquid nitrogen for cryopreservation. Nine liver samples were collected at each temperature (3 parallel × 3 fish/parallel = 9 fish). A total of 27 samples were collected at three temperature groups (CT, HT31, and HT34).

RNA Extraction, Library Construction, and Sequencing
For each temperature group, three liver samples were randomly selected for library construction and sequencing. Nine samples (3/group × 3 groups) at three temperature groups were used for all experiments and subsequent analyses. According to the manufacturer's instructions, total RNA was extracted from nine liver samples using a Trizol reagent (Invitrogen, Carlsbad, CA, USA). Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and RNase-free agarose gel electrophoresis were used to detect the integrity of total RNA and the presence of contamination. Qubit 2.0 Fluorometer (Invitrogen) was used to quantify RNA concentration accurately, and a NanoPhotometer spectrophotometer (Implen, Westlake Village, CA, USA) was used to detect the purity of total RNA. Only RNA samples with the required integrity, concentration, and purity were used for subsequent library construction and Illumina RNA-Seq.
The messenger RNAs (mRNAs) with poly (A) tails were enriched from qualified total RNA using magnetic beads with oligo (dT), and the enriched mRNA was fragmented into short fragments by ultrasound. The short mRNA fragments were reverse-transcribed into first-strand complementary DNA (cDNA) using random primers, and M-MuLV reverse transcriptase, followed by second-strand cDNA synthesis using RNase H and DNA polymerase I. The cDNA fragments were purified with a QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands), end-repaired, poly (A) added, and ligated to Illumina sequencing adapters. The cDNA of about 200 bp was screened using AMPure XP beads. PCR amplification was then performed, and the PCR product was purified using AMPure XP beads to obtain a cDNA library. To obtain high-quality sequencing data, the quality of Three parallel groups were set up in the experiment, and three temperature gradients, including 27 • C (control group, CT), 31 • C (high-temperature group, HT31), and 34 • C (high-temperature group, HT34), were set up in each parallel group. Before the heat stress experiment, all experimental fish were fasted for 24 h and were not fed during the experiment. After the beginning of the experiment, water temperature in three parallel groups was gradually increased from 27 • C to 34 • C at 1 • C/h and maintained at 34 • C for 2 h using heating rods (WN-B15). Three fish were selected randomly from each parallel group at 27 • C, 31 • C, and 34 • C (kept for 2 h) and anesthetized with phenolic. After measuring body length and weight, fish were dissected, and their liver tissues were removed and put into liquid nitrogen for cryopreservation. Nine liver samples were collected at each temperature (3 parallel × 3 fish/parallel = 9 fish). A total of 27 samples were collected at three temperature groups (CT, HT31, and HT34).

RNA Extraction, Library Construction, and Sequencing
For each temperature group, three liver samples were randomly selected for library construction and sequencing. Nine samples (3/group × 3 groups) at three temperature groups were used for all experiments and subsequent analyses. According to the manufacturer's instructions, total RNA was extracted from nine liver samples using a Trizol reagent (Invitrogen, Carlsbad, CA, USA). Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and RNase-free agarose gel electrophoresis were used to detect the integrity of total RNA and the presence of contamination. Qubit 2.0 Fluorometer (Invitrogen) was used to quantify RNA concentration accurately, and a NanoPhotometer spectrophotometer (Implen, Westlake Village, CA, USA) was used to detect the purity of total RNA. Only RNA samples with the required integrity, concentration, and purity were used for subsequent library construction and Illumina RNA-Seq.
The messenger RNAs (mRNAs) with poly (A) tails were enriched from qualified total RNA using magnetic beads with oligo (dT), and the enriched mRNA was fragmented into short fragments by ultrasound. The short mRNA fragments were reverse-transcribed into first-strand complementary DNA (cDNA) using random primers, and M-MuLV reverse transcriptase, followed by second-strand cDNA synthesis using RNase H and DNA polymerase I. The cDNA fragments were purified with a QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands), end-repaired, poly (A) added, and ligated to Illumina sequencing adapters. The cDNA of about 200 bp was screened using AMPure XP beads. PCR amplification was then performed, and the PCR product was purified using AMPure XP beads to obtain a cDNA library. To obtain high-quality sequencing data, the quality of the cDNA library was strictly checked using the same method described for total RNA. The qualified cDNA library was sequenced using Illumina HiSeqTM 4000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Sequencing Quality Detection and Comparison
To ensure the quality of data and downstream analyses, the raw reads obtained from sequencing were filtered by fastp version 0.18.0 [45] to remove low-quality reads containing adapter sequences, reads with more than 10% of unknown nucleotide (N), all A base, and low-quality reads (more than 50% of the bases with a quality score Q ≤ 20). The quality of clean reads was evaluated by Q20, Q30, GC content, and error rate. Clean reads were mapped to the ribosome RNA (rRNA) database using the short-reads alignment tool bowtie2 (version 2.2.8) [46] to remove rRNA-mapped reads. The remaining high-quality clean reads were used for subsequent analyses.
The gene annotation file of the reference genome and gene model of T. ovatus was downloaded from Figshare [47], and the index of the reference genome was constructed. Highquality clean reads were mapped to the reference genome of T. ovatus by HISAT 2.2.4 [48], obtaining the position and annotation information on the reference genome or gene.

Differential Gene Expression Analysis
The mapped reads were normalized as fragments per kilobase of transcript per million mapped reads (FPKM), and further data analyses were based on normalized data. Pearson's correlation coefficient was calculated to evaluate the repetition correlation between parallel samples within every temperature group. Differentially expressed gene (DEG) analysis was performed using DESeq2 [49], and the genes with false discovery rate (FDR) < 0.05 and |log2 (Fold Change)| > 1 were considered significant DEGs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed for all significant DEGs. The GO terms and KEGG pathways with p < 0.05 were considered significantly enriched.

STEM Analysis
Based on the transformed FPKM values, all the genes having similar expression patterns were classified by STEM [50]. The parameters used for STEM analysis were as follows: number of modules 20, correlation coefficient > 0.7, and significance threshold p < 0.05. GO term and KEGG pathway enrichment analyses were performed on each significant gene expression profile, understanding the specific biological processes in response to heat stress.
Based on the protein interaction data from T. ovatus, protein-protein interaction (PPI) network analysis of DEGs in each significant profile was performed using the String v10 [51] database to find the module network and prioritize genes. The first 300 genes were used to construct gene co-expression networks, and the network graph was visualized using Cytoscape v3.8.2 [52]. Finally, candidate hub genes were screened according to their highest connectivity with other genes.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR) Validation of mRNA Expression Patterns
For confirming the reliability of transcriptome data, 11 DEGs were randomly selected for qRT-PCR validation, and the β-actin gene was used as an internal reference gene. The gene-specific primers were designed using Primer-BLAST in National Center for Biotechnology Information (NCBI) based on the target gene sequences (Table 1). Total RNA was extracted from nine liver samples of CT, HT31, and HT34 groups using TRIzol Reagent, and RNA concentration and integrity were detected using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis. About 100 ng RNA was used for reverse transcription using TransScript ® Uni All-in-One First-Strand cDNA Synthesis SuperMix for qPCR Kit (TransGen Biotech, Beijing, China) following the manufacturer's instructions.
The final volume of qRT-PCR reaction was 20 µL, including 8.2 µL RNase-free water, 10 µL 2× Perfect Start ® Green qPCR Super Mix (Beijing TransGen Biotech Co., Ltd., Beijing, China), 1 µL cDNA, 0.4 µL forward primer (10 mmol), and 0.4 µL reverse primer (10 mmol). Light Cyeler 96 (Roche, Germany) was applied for qRT-PCR, and the PCR reaction cycle program was as follows: pre-denaturation at 94 • C for 30 s, followed by 40 cycles of 94 • C for 5 s and 60 • C for 30 s. Each sample was prepared in three technical replicates, and each experiment contained a negative control without template DNA. After qRT-PCR amplification, the amplification curve and melting curve were examined to evaluate specific amplification and the relative expression levels of target genes were calculated by the 2 − CT method. The data of 11 selected genes were compared and obtained from qRT-PCR and transcriptome sequencing to validate the reliability of transcriptome sequencing data.

Transcriptome Sequencing
A total of 386,879,166 and 58,031,874,900 bp raw data were obtained from the transcriptome sequencing of 9 samples in three temperature groups (  (Table 3), of which that of eight samples was higher than 90%, which indicates that the transcriptome sequencing data have a good alignment result with the genome. The results of sequencing, quality control, and genome comparison showed that the transcriptome sequencing data were reliable, which could ensure the accuracy of subsequent analysis.

Differential Gene Expression Analysis
As shown in Figure 2A, the Pearson correlation coefficients between parallel samples within the CT, HT31, and HT34 groups were 0.98-0.99, 0.93-1.00, and 0.77-0.96, respectively, indicating good reproducibility between biological replicates. Based on FPKM values, DEGs between temperature groups were screened. A total of 297 (188 up-regulated and 109 downregulated), 842 (413 up-regulated and 429 down-regulated), and 665 (334 up-regulated and 331 down-regulated) DEGs were detected in CT vs. HT31, CT vs. HT34, and HT31 vs. HT34 comparison groups, respectively ( Figure 2B). DEGs increased significantly (from 297 to 842) with the increase of water temperature (from 31 • C to 34 • C), indicating that the increasing temperature could induce changes in the expression levels of more related genes. A total of 1347 DEGs were identified in three comparison groups by temperature ( Figure 2C), of which nine were shared by three comparison groups; 439 were shared by two comparison groups, and the other 899 were unique for the single comparison group. A clustering heat map of DEGs ( Figure 2D) showed that 1347 DEGs had different expression patterns at three temperature points. A total of 55 level 2 GO terms were obtained by GO enrichment analysis of 1347 DEGs, of which 24 belonged to biological processes (BP), 12 belonged to molecular functions (MF), and 19 belonged to cellular components (CC). Separate GO enrichment analysis of each comparison group ( Figure 3) showed that CT vs. HT31, CT vs. HT34, and HT31 vs. HT34 had 52, 54, and 54 level 2 GO terms, respectively, of which 51 were shared by three comparison groups, and only four (cell killing of biological processes, structural molecule activity, electron carrier activity, and translation regulator activity in molecular function) were different. These results indicated that the physiological and biochemical processes related to high temperature in the liver of T. ovatus were activated when water temperature gradually rose to 31 °C. Although the three comparison groups had similar GO terms, the DEGs of the main GO terms increased significantly with the increase in water temperature. For example, cellular processes in biological processes contained 124 and 410 DEGs in CT vs. HT31 and CT vs. HT34, respectively. The above indicated that more genes participated in the biological processes activated at 31 °C along with a further increase in water temperature (from 31 to 34 °C). A total of 55 level 2 GO terms were obtained by GO enrichment analysis of 1347 DEGs, of which 24 belonged to biological processes (BP), 12 belonged to molecular functions (MF), and 19 belonged to cellular components (CC). Separate GO enrichment analysis of each comparison group ( Figure 3) showed that CT vs. HT31, CT vs. HT34, and HT31 vs. HT34 had 52, 54, and 54 level 2 GO terms, respectively, of which 51 were shared by three comparison groups, and only four (cell killing of biological processes, structural molecule activity, electron carrier activity, and translation regulator activity in molecular function) were different. These results indicated that the physiological and biochemical processes related to high temperature in the liver of T. ovatus were activated when water temperature gradually rose to 31 • C. Although the three comparison groups had similar GO terms, the DEGs of the main GO terms increased significantly with the increase in water temperature. For example, cellular processes in biological processes contained 124 and 410 DEGs in CT vs. HT31 and CT vs. HT34, respectively. The above indicated that more genes participated in the biological processes activated at 31 • C along with a further increase in water temperature (from 31 to 34 • C).
Level 3 GO enrichment analysis showed (Table S1) that a total of 2395, 3454, and 3109 GO terms were detected in CT vs. HT31, CT vs. HT34, and HT31 vs. HT34, respectively, of which 395, 495, and 343 were significant (p < 0.05). The significantly enriched GO terms included angiogenesis, antigen processing and presentation of peptide antigen via MHC class I and apoptotic process, carbohydrate metabolic process, cellular response to stress, immune system process, misfolded or incompletely synthesized protein catabolic process, regulation of protein folding, response to oxidative stress, tRNA metabolic process, proteasome-mediated ubiquitin-dependent protein catabolic process, regulation of ubiquitin-protein transferase activity, and so on. Level 3 GO enrichment analysis showed (Table S1) that a total of 2395, 3454, and 3109 GO terms were detected in CT vs. HT31, CT vs. HT34, and HT31 vs. HT34, respectively, of which 395, 495, and 343 were significant (p < 0.05). The significantly enriched GO terms included angiogenesis, antigen processing and presentation of peptide antigen via MHC class I and apoptotic process, carbohydrate metabolic process, cellular response to stress, immune system process, misfolded or incompletely synthesized protein catabolic process, regulation of protein folding, response to oxidative stress, tRNA metabolic process, proteasome-mediated ubiquitin-dependent protein catabolic process, regulation of ubiquitinprotein transferase activity, and so on.
The big categories of KEGG pathway enrichment analysis of DEGs ( Figure 4A-C) showed that the significantly enriched pathways (p < 0.05) in three comparison groups were classified into six major categories. Metabolism, genetic information processing, and organismal systems were common in three comparison groups. Cellular processes and human diseases were common to CT vs. HT31, and HT31 vs. HT34, and environmental information processing was unique to CT vs. HT31. Environmental information processing processes, such as signaling, were primarily activated when water temperature rose from 27 to 31 °C, further activating more physiological processes in the liver of T. ovatus. The big categories of KEGG pathway enrichment analysis of DEGs ( Figure 4A-C) showed that the significantly enriched pathways (p < 0.05) in three comparison groups were classified into six major categories. Metabolism, genetic information processing, and organismal systems were common in three comparison groups. Cellular processes and human diseases were common to CT vs. HT31, and HT31 vs. HT34, and environmental information processing was unique to CT vs. HT31. Environmental information processing processes, such as signaling, were primarily activated when water temperature rose from 27 to 31 • C, further activating more physiological processes in the liver of T. ovatus.

Expression Trend and Enrichment Analysis of DEGs
A total of 1347 DEGs were clustered into eight profiles ( Figure 5A), of which four profiles (4, 3, 0, and 7) were significant (p = 4.5 × 10 −34 , 4.0 × 10 −31 , 1.4 × 10 −10 , and 2.6 × 10 −8 ). Profiles 4, 3, 0, and 7 had 315, 342, 180, and 170 genes, respectively, which were significantly higher than those of the other four non-significant profiles (49-108 genes). Four significant expression profiles showed different expression trends ( Figure 5B). The expression levels of most DEGs in profile 3 and profile 4 remained unchanged as the water temperature was raised from 27 • C to 31 • C. When the water temperature reached 34 • C, the gene expression decreased in Profile 3 but increased in Profile 4. Different from Profile 3 and Profile 4, the DEGs in Profile 0 and Profile 7 were down-regulated and up-regulated at two high-temperature points, respectively. Overall, Profile 4 and Profile 7 showed up-regulation, and Profile 3 and Profile 0 showed down-regulation with the increase in water temperature.
GO and KEGG enrichment analyses were performed for DEGs in profiles 4, 3, 0, and 7 to screen key biological processes related to high temperature. The GO enrichment analysis results showed ( Figure 6) that profiles 4, 3, 0, and 7 were enriched to 51, 51, 47, and 46 levels 2 GO terms, respectively. Four profiles had similar GO terms. The first four terms in biological processes were cellular processes, single-organism processes, metabolic processes pathway, and biological regulation; the first two terms in molecular function were binding and catalytic activity; the first three terms in cellular component were cell, cell part, and organelle. gene expression decreased in Profile 3 but increased in Profile 4. Different from Profile 3 and Profile 4, the DEGs in Profile 0 and Profile 7 were down-regulated and up-regulated at two high-temperature points, respectively. Overall, Profile 4 and Profile 7 showed upregulation, and Profile 3 and Profile 0 showed down-regulation with the increase in water temperature. GO and KEGG enrichment analyses were performed for DEGs in profiles 4, 3, 0, and 7 to screen key biological processes related to high temperature. The GO enrichment analysis results showed ( Figure 6) that profiles 4, 3, 0, and 7 were enriched to 51, 51, 47, and 46 levels 2 GO terms, respectively. Four profiles had similar GO terms. The first four terms in biological processes were cellular processes, single-organism processes, metabolic processes pathway, and biological regulation; the first two terms in molecular function were binding and catalytic activity; the first three terms in cellular component were cell, cell part, and organelle.
GO enrichment analysis at level 3 (Table S2) showed that DEGs in profiles 4, 3, 0, and 7 were enriched to 2453, 2171, 1760, and 1771 GO terms, respectively, in which 358, 277, 230, and 231 were significant (p < 0.05). The significantly enriched GO terms of DEGs in profile four mainly included carbohydrate catabolic process, angiogenesis, cellular carbohydrate metabolic process, single-organism carbohydrate metabolic process, regulation of carbohydrate metabolic process, regulation of cellular carbohydrate metabolic process, and so on. The significantly enriched GO terms for DEGs in profile 3 included antigen processing and presentation of exogenous peptide antigen via MHC class I, positive regulation of ubiquitin-protein transferase activity, positive regulation of protein ubiquitination, proteasome complex, proteasome-mediated ubiquitin-dependent protein catabolic process, and so on. The predominant enriched GO terms of DEGs in profile 0 included RNA binding, arginine N-methyltransferase activity, rRNA methyltransferase activity, tRNA methyltransferase activity, ubiquitin-like protein transferase activity, ubiquitinprotein transferase activity, and so on. The significantly enriched GO terms for DEGs in Profile 7 mainly include sperm-egg recognition, protein folding, 'de novo' protein folding, small molecule binding, regulation of ATPase activity, complex protein assembly, etc. GO enrichment analysis at level 3 (Table S2) showed that DEGs in profiles 4, 3, 0, and 7 were enriched to 2453, 2171, 1760, and 1771 GO terms, respectively, in which 358, 277, 230, and 231 were significant (p < 0.05). The significantly enriched GO terms of DEGs in profile four mainly included carbohydrate catabolic process, angiogenesis, cellular carbohydrate metabolic process, single-organism carbohydrate metabolic process, regulation of carbohydrate metabolic process, regulation of cellular carbohydrate metabolic process, and so on. The significantly enriched GO terms for DEGs in profile 3 included antigen processing and presentation of exogenous peptide antigen via MHC class I, positive regulation of ubiquitin-protein transferase activity, positive regulation of protein ubiquitination, proteasome complex, proteasome-mediated ubiquitin-dependent protein catabolic process, and so on. The predominant enriched GO terms of DEGs in profile 0 included RNA binding, arginine N-methyltransferase activity, rRNA methyltransferase activity, tRNA methyltransferase activity, ubiquitin-like protein transferase activity, ubiquitin-protein transferase activity, and so on. The significantly enriched GO terms for DEGs in Profile 7 mainly include sperm-egg recognition, protein folding, 'de novo' protein folding, small molecule binding, regulation of ATPase activity, complex protein assembly, etc.
In KEGG pathway enrichment analysis (Figure 7), profiles 4, 3, 0, and 7 were enriched in 218, 237, 132, and 182 pathways, respectively, among which 27, 19, 6, and 33 pathways were significantly enriched (p < 0.05). The significantly enriched KEGG pathways in profile 4 included primarily amino acid metabolism, cancer, signal transduction, carbohydrate metabolism, and other related pathways, which were refined to the biosynthesis of amino acids, phenylalanine metabolism, microRNAs (miRNAs) in cancer, HIF-1 signaling pathway, glycolysis/gluconeogenesis, and so on. The significantly enriched KEGG pathways in profile 3 included folding, sorting and degradation, translation, carbohydrate metabolism, and other related pathways, and more specific pathways included proteasome, protein export, ribosome biogenesis in eukaryotes, and so on. The significantly enriched KEGG pathways in profile 0 included primarily translation, transcription, folding, sorting and degradation, and other related pathways, which were refined to ribosome biogenesis in eukaryotes, spliceosome, ubiquitin-mediated proteolysis, and so on. The significantly enriched KEGG pathway in profile 7 included carbohydrate metabolism, amino acid metabolism, immune system, folding, sorting and degradation, and other related pathways, which could be refined to more specific pathways, such as glycolysis/gluconeogenesis, citrate cycle (TCA cycle), cysteine and methionine metabolism, antigen processing and presentation, protein processing in the endoplasmic reticulum, and so on. In KEGG pathway enrichment analysis (Figure 7), profiles 4, 3, 0, and 7 were enriche in 218, 237, 132, and 182 pathways, respectively, among which 27, 19, 6, and 33 pathway were significantly enriched (p < 0.05). The significantly enriched KEGG pathways in pro file 4 included primarily amino acid metabolism, cancer, signal transduction, carbohy drate metabolism, and other related pathways, which were refined to the biosynthesis o amino acids, phenylalanine metabolism, microRNAs (miRNAs) in cancer, HIF-1 signalin pathway, glycolysis/gluconeogenesis, and so on. The significantly enriched KEGG path ways in profile 3 included folding, sorting and degradation, translation, carbohydrate me tabolism, and other related pathways, and more specific pathways included proteasom protein export, ribosome biogenesis in eukaryotes, and so on. The significantly enriche The hub genes were screened from four significant profiles (4, 3, 0, and 7) by protein interaction network analysis. As shown in Figure 8, five (PRDM10, MYC, HIF1α, JUNB, and GAPDH), three (PSMD8, PSMD3 and SEC61AL1), four (DDX5, DKC1, FBL, and PRMT1), and four (HSP90AB1, HSPD1, TCP1, and CCT5) hub genes were screened from profiles 4, 3, 0, and 7, respectively. The 16 hub genes had relatively high connectivity with other DEGs in the profiles, which suggests that these genes may play a relatively important role in response to acute high-temperature stress. which could be refined to more specific pathways, such as glycolysis/gluconeogenesis, citrate cycle (TCA cycle), cysteine and methionine metabolism, antigen processing and presentation, protein processing in the endoplasmic reticulum, and so on. The hub genes were screened from four significant profiles (4, 3, 0, and 7) by protein interaction network analysis. As shown in Figure 8, five (PRDM10, MYC, HIF1α, JUNB, and GAPDH), three (PSMD8, PSMD3 and SEC61AL1), four (DDX5, DKC1, FBL, and PRMT1), and four (HSP90AB1, HSPD1, TCP1, and CCT5) hub genes were screened from profiles 4, 3, 0, and 7, respectively. The 16 hub genes had relatively high connectivity with other DEGs in the profiles, which suggests that these genes may play a relatively important role in response to acute high-temperature stress.

Validation of Gene Expression Patterns by qRT-PCR
To verify the DEGs identified by RNA-seq, we randomly selected 11 genes involved in HSPs (HSP90AB1), Ubiquitination-Proteasome (PSMD1, PSMC5, and PSMB3), and HIF1A signaling pathway (HIF1α, PFKL, GAPDH, ALDOCB, ENO1, EPO, and MKNK2) for qPCR. As shown in Figure 9, the relative expression levels of 11 genes were consistent with the RNA-Seq results, illustrating the accuracy and validity of our RNA-seq data and subsequent analyses.

Validation of Gene Expression Patterns by qRT-PCR
To verify the DEGs identified by RNA-seq, we randomly selected 11 genes involved in HSPs (HSP90AB1), Ubiquitination-Proteasome (PSMD1, PSMC5, and PSMB3), and HIF1A signaling pathway (HIF1α, PFKL, GAPDH, ALDOCB, ENO1, EPO, and MKNK2) for qPCR. As shown in Figure 9, the relative expression levels of 11 genes were consistent with the RNA-Seq results, illustrating the accuracy and validity of our RNA-seq data and subsequent analyses.

Discussion
Based on RNA-Seq technology, 1,347 DEGs and four significant profiles were screened from the transcriptomic data of 9 liver samples of T. ovatus. The enrichment analysis of both DEGs and significant profiles showed that T. ovatus liver cells initiated biological processes, such as protein and amino acid metabolism, glucose metabolism, transcription, translation, folding, sorting and degradation to cope with acute high-temperature stress and reduce liver damage. Among these processes (Figure 10), transcriptional changes of genes related to heat shock proteins (HSPs), the ubiquitin-proteasome system (ubiquitinmediated proteolysis and proteasome), hypoxia adaptation (HIF-1 signaling pathway), energy metabolism (glycolysis/gluconeogenesis and fatty acid biosynthesis), etc. play very important roles in high-temperature adaptation.

Discussion
Based on RNA-Seq technology, 1,347 DEGs and four significant profiles were screened from the transcriptomic data of 9 liver samples of T. ovatus. The enrichment analysis of both DEGs and significant profiles showed that T. ovatus liver cells initiated biological processes, such as protein and amino acid metabolism, glucose metabolism, transcription, translation, folding, sorting and degradation to cope with acute high-temperature stress and reduce liver damage. Among these processes (Figure 10), transcriptional changes of genes related to heat shock proteins (HSPs), the ubiquitin-proteasome system (ubiquitin-mediated proteolysis and proteasome), hypoxia adaptation (HIF-1 signaling pathway), energy metabolism (glycolysis/gluconeogenesis and fatty acid biosynthesis), etc. play very important roles in high-temperature adaptation.

HSPs and Correct Folding of Proteins
Acute high-temperature stress can cause a cellular stress response in fish, leading to misfolding of intracellular proteins. If misfolded proteins accumulate in the body, it harms cells and interferes with their regular physiological functions [53]. HSPs are a highly conserved and ubiquitous protein family whose members participated in many important biological processes under abiotic (e.g., high-temperature, hypoxia, and heavy metal ions) and biotic stress, such as proper protein folding and processing, classical antigen presentation and cross-presentation [54,55]. By coordinating various functions of cells, HSPs protect the body from harm caused by high-temperature stress [56].
In this study, under acute heat stress, the expression of HSP70 and HSP90AB1 in the T. ovatus liver significantly increased, and the increase was positively correlated with temperature. HSP70 is the most important member of the HSP family, with a molecular weight of around 70 KD. It is not expressed or expressed only in small amounts in cells under normal conditions, while its expression increases rapidly under stress conditions

HSPs and Correct Folding of Proteins
Acute high-temperature stress can cause a cellular stress response in fish, leading to misfolding of intracellular proteins. If misfolded proteins accumulate in the body, it harms cells and interferes with their regular physiological functions [53]. HSPs are a highly conserved and ubiquitous protein family whose members participated in many important biological processes under abiotic (e.g., high-temperature, hypoxia, and heavy metal ions) and biotic stress, such as proper protein folding and processing, classical antigen presentation and cross-presentation [54,55]. By coordinating various functions of cells, HSPs protect the body from harm caused by high-temperature stress [56].
In this study, under acute heat stress, the expression of HSP70 and HSP90AB1 in the T. ovatus liver significantly increased, and the increase was positively correlated with temperature. HSP70 is the most important member of the HSP family, with a molecular weight of around 70 KD. It is not expressed or expressed only in small amounts in cells under normal conditions, while its expression increases rapidly under stress conditions [57]. This protein is involved in folding nascent polypeptides, repairing misfolded proteins, and enhancing the heat tolerance of cells or the body, and it can speed up the restoration of protein synthesis [58]. HSP90 is one of the molecular chaperones found in cells, and its family members include HSP90α, HSP90β, and HSP83, etc. Of these, HSP90β is encoded by the HSP90AB1 gene and participates in the regulation of cell signal transduction, protein folding, and long-term cell adaptability [59]. The upregulated expression levels of HSP70 and HSP90AB1 in the liver of T. ovatus are helpful for the refolding of misfolded proteins in cells, maintaining protein homeostasis and regular physiological functions, and minimizing damage to the fish body during heat stress. Consistent with this result, as the water temperature increased, the expression levels of HSP70 and HSP90 were also enhanced in brown trout (Salmo trutta fario) skin tissue [60], grass carp (Ctenopharyngodon idellus) liver tissue [61], spotted sea bass (Lateolabrax maculatus) liver and muscle tissue [62], O. mykiss head kidney and brain tissue [63,64], starry flounder (Platichthys stellatus) kidney tissue [65], and lake whitefish (Coregonus clupeaformis) liver tissue [66]. This regulation helps maintain intracellular protein homeostasis, minimizing the damage caused by high-temperature stress to the body. In summary, under high-temperature stress, fish can regulate the homeostasis of proteins by regulating the expression level of HSPs to maintain the regular physiological functions of the body and adapt to the high-temperature environment.

Ubiquitin-Proteasome System (UPS) and Protein Degradation
In addition to HSPs, another important means to regulate protein homeostasis and function in eukaryotic cells is UPS, and the degradation of more than 80% of proteins in cells depends on this system [67]. Usually, the UPS can recognize and degrade unstable, denatured, or misfolded proteins to ensure that proteins in cells are in a homeostatic state, which is important to maintain homeostasis in cells [68]. The degradation of proteins by the UPS is an energy-consuming process consisting of a stepwise enzymatic cascade reaction. During this process, through the sequential actions of ubiquitin-activating enzyme (E1), ubiquitin-conjugating enzyme (E2), and ubiquitin ligase (E3), ubiquitin (Ub) is bound to the specific target protein, and the ubiquitinated target protein is recognized and degraded by the 26S proteasome.
The recognition and degradation of many proteins in cells depend on ubiquitination. This study found four significantly downregulated ubiquitination-related genes in T. ovatus liver, namely, ubiquitin-conjugating enzyme E2D4 (UBE2D4), ubiquitin-protein ligase E3B (UBE3B), photomorphogenic regulatory factor 1 (COP1), and ring finger and CHY zinc finger domain containing 1 (RCHY1) after acute heat stress treatment. The UBE2D4 gene belongs to the ubiquitin-conjugating enzyme E2 gene family and encodes the ubiquitinconjugating enzyme (E2), which acts as an intermediate to deliver activated ubiquitin molecules to E3 connected to the target protein [69,70]. This intermediate step is essential in the whole ubiquitination process. The UBE3B, COP1, and RCHY1 genes encode ubiquitin ligase E3, a key enzyme in the ubiquitination reaction that binds the correct E2 and the substrate to increase the rate of ubiquitin transfer [71,72]. The gene encoding E2 in the liver tissue of antarctic eelpout (Pachycara brachycephalum) was suppressed at high temperatures [73], consistent with the results of this study. Under CO 2 and high-temperature stress, the gill tissues of blue-green damselfish (Chromis viridis) [74] and O. mykiss [75] showed an enhanced ubiquitination ability by upregulating the gene encoding E2, accelerating the degradation of damaged proteins in vivo. In addition, the degradation of many proteins (e.g., eukaryotic low-density lipoprotein receptor (LDLR) and mouse p53) was inhibited, along with the reduction of E2 and E3 activities [70,76]. E2 and E3 play important roles in the degradation process of eukaryotic proteins. The downregulation of E2 and E3-related genes in the liver of T. ovatus suggests that the activities of E2 and E3 may be reduced due to acute high-temperature stress, which slows down the ubiquitination process of target proteins (denatured or misfolded ones).
Different Ub molecules specifically mark the target protein, and the target protein marked by Ub will be degraded into small peptides or amino acids by the 26S proteasome [68]. The 26S proteasome is a highly conserved and ATP-dependent complex enzyme composed of 19S Regulatory Particle (RP) and a 20S Core Particle (CP) [77]. The ubiquitinated protein first binds to 19S RP and is degraded in 20S CP to regulate cellular protein levels. 19S RP can be separated biochemically into the lid and base subcomplexes, among which the lid is composed of nine protein subunits, and its main function is to cut off the ubiquitin chain from the protein substrate before the protein substrate is degraded for recycling (de-ubiquitination) [78]. The Base is composed of 10 protein subunits, and its main function is to use the energy provided by ATP hydrolysis to unfold the protein substrate that has separated from the ubiquitin chain and transport it to 20S CP for degradation [79]. In this study, after high-temperature stress, the expression levels of proteasome 26S subunit non-ATPase 3 (PSMD3), non-ATPase 8 (PSMD8), non-ATPase 1 (PSMD1), ATPase 4 (PSMC4), and ATPase 5 (PSMC5) in T. ovatus liver were significantly downregulated at 34 • C. PSMD3 and PSMD8 encode two subunits in the lid; PSMD1, PSMC4, and PSMC5 encode three subunits in the base. Some studies have indicated that the decreased expression of proteasome subunits leads to the loss of proteasome activity [80]. Therefore, in T. ovatus, downregulation of PSMD3, PSMD8, PSMD1, PSMC4, and PSMC5 may reduce the de-ubiquitination and protein unfolding abilities of 19S RP and its protein transfer function of presenting deubiquitinated and unfold proteins to 20SCP, slowing down the subsequent protein degradation process. 20S CP is the catalytically active part of the 26S proteasome and is responsible for recognizing, stretching, and degrading deubiquitinated proteins [81]. It mainly comprises two α-rings and two β-rings, and each ring contains seven structurally similar subunits. In this study, after high-temperature stress treatment, the expression levels of proteasome 20S subunit alpha 4 (PSMA4), alpha 7 (PSMA7), beta 3 (PSMB3), and beta 7 (PSMB7) in T. ovatus were significantly downregulated at 34 • C, which may cause the reduction or loss of 20S CP catalytic activity. That is, the degradation of denatured or misfolding proteins was likely to slow down. In C. idellus intestinal mucosa and O. mykiss gill tissue under the conditions of oxidized fish oil and high CO 2 concentration, respectively, PSMA4, PSMB3, and PSMB7 were significantly upregulated, increasing the activity of 26S proteasome and helping to clear the damaged or incorrect proteins to protect tissues and cells from impairment [75,82]. The expression levels of proteasome subunit-encoding genes were increased in the muscles of rats with long-term high-intensity exercise and C. idellus infected with saprolegniasis, resulting in enhanced activity of 26S proteasome [83,84]. Genes encoding proteasome subunits in C. idellus head kidney tissue responded to high-density and saline-alkali stress mainly in negative regulation [85], similar to the results found in this study. The above results indicate that under acute heat stress, the significant downregulation of proteasome-related genes in T. ovatus may lead to lower levels of proteasome subunits, leading to the reduction or loss of 26S proteasome activity.
In summary, under acute high-temperature stress, the ubiquitination-proteasome system in the liver tissue of T. ovatus is inhibited to a certain extent. The degradation rate of misfolded proteins is reduced, allowing sufficient time for the refolding of misfolded proteins and the quality control of candidate degrading proteins. Under high-temperature stress, the concentration of DO in water decreases, and energy metabolism undergoes compensatory changes (see Section 4.2 below), yet protein synthesis and degradation consume large amounts of energy [86]. Therefore, the UPS is inhibited, representing an energy consumption reduction strategy [73]. As an important protein modification pathway, the UPS participates in multiple biological processes. However, there are few studies on UPS in fish, and the regulatory mechanisms of the genes are intricate. Thus, further in-depth exploration is required.

Adaptation to Hypoxia and Transition of Energy Metabolism
Water temperature is a crucial factor affecting aquatic organisms' metabolic rate and water's DO content. With the increase in water temperature, the metabolism of aquatic organisms is enhanced, and the oxygen consumption rate is also significantly increased. However, the temperature increase may reduce hemoglobin's binding ability to oxygen [87]. In addition, the amount of DO in water is correlated with water temperature, i.e., an in-crease in water temperature leads to a decrease in DO in water, aggravating the hypoxia of aquatic organisms [88]. To adapt to the low-oxygen environment under high temperatures, some physiological activities in the metabolism of aquatic organisms may undergo adaptive changes to maintain homeostasis and normal physiological functions. In this study, with the increase in water temperature, the expression levels of genes related to hypoxia response (hypoxia-inducible factor 1, HIF-1), oxygen transport (oxygen transportation erythropoietin, EPO), and energy metabolism (phosphofructokinase L, PFKL; fructose-bisphosphate aldolase B, ALDOB; glyceraldehyde phosphate dehydrogenase, GAPDH; and enolase-1, ENO1) in the liver of T. ovatus showed significant changes. These changes may alter the body's oxygen-carrying capacity and energy metabolism, resulting in adaptation to the hypoxic environment generated by acute heat stress.
HIF-1 is a key molecule regulating the hypoxic response of tissue cells and is composed of HIF-1α and HIF-1β subunits [89]. Of these, HIF-1α is the active subunit of HIF-1. This subunit is widely present in most tissues, which are regulated by hypoxic signals and can regulate the activity of HIF-1 [90]. The effect of HIF-1α depends on the oxygen concentration in the cell. Under normoxic conditions, HIF-1α is rapidly degraded through UPS, but under hypoxic conditions, HIF-1α will continue to accumulate and bind to HIF-1β inside the nucleus, forming active HIF-1 [91]. This study found that the expression of HIF-1α in the liver of T. ovatus decreased at 31 • C but was significantly upregulated at 34 • C. Most tissues of channel catfish (Ictalurus punctatus) [92] and the embryonic tissue of zebrafish (Danio rerio) [93] showed a similar expression trend of HIF-1α (decreased first and then increased) during hypoxic stress. During the early stage of hypoxia, the expression of HIF-1α was not upregulated but downregulated, which may be because HIF-1 gradually accumulated during the early stage, and this accumulation then negatively regulated gene transcription [93]. Other research indicates that hypoxic shock may also cause downregulated expression of HIF-1α [94,95]. Under hypoxic stress, the expression of HIF-1α was significantly upregulated in D. rerio embryos [94], oscar (Astronotus ocellatus) livers [96], yellow catfish (Pelteobagrus fulvidraco) livers [97], O. mykiss kidney [98], etc. It can be seen that the continuous hypoxic environment under high-temperature stress can significantly increase the expression of HIF-1α in T. ovatus. HIF-1α further forms different signaling pathways with various upstream and downstream proteins, mediating hypoxic signals and regulating a series of compensatory responses to hypoxia in cells.
HIF-1α is an important factor in regulating EPO expression [99]. This study showed that the expression of EPO, a downstream gene of HIF-1α, was slightly downregulated at 31 • C but significantly upregulated at 34 • C, which was consistent with the expression trend of HIF-1α. EPO is a glycoprotein hormone that can stimulate the differentiation and proliferation of erythroid precursor cells [100]. It is a key factor in regulating erythropoiesis and ensuring tissue oxygen supply [101]. The upregulation of EPO may improve the differentiation of red blood cells in the liver, effectively increase the number of red blood cells and hemoglobin, and then increase the body's O 2 -carrying capacity and utilization of oxygen to compensate for the body's hypoxia state in a high-temperature hypoxic environment [102,103]. Therefore, under high temperatures and a hypoxic environment, the upregulation of EPO's expression level helps T. ovatus effectively use low oxygen concentrations and maintain normal physiological processes.
In addition to the EPO gene, the upregulated HIF-1α can also regulate the transcription of downstream glycolysis-related genes [104,105], which may lead to compensatory changes in energy metabolism to adapt to the hypoxic environment caused by high temperature. This study showed that when the temperature increased from 27 • C to 34 • C, the expression of genes encoding several key enzymes of glycolysis (PFKL, ALDOB, GAPDH, and ENO1) in the liver of T. ovatus showed significant upregulation trends. Among them, PFKL (a glycolytic rate-limiting enzyme) plays a decisive role in regulating the direction and rate of the glycolytic metabolic pathway and can catalyze the reaction of fructose-6-phosphate (F6P) to generate fructose-1,6-phosphate (F1,6P). Both ALDOB and GAPDH are key enzymes of glycolysis, and ENO1 is the rate-limiting enzyme of glycolysis. F1,6P is gradually converted into phosphoenolpyruvate (PEP) with enzymes, such as ALDOB, GAPDH, and ENO1, which promotes the generation of more pyruvate to participate in the tricarboxylic acid cycle or other metabolic pathways. In this study, the significant upregulation of PFKL, ALDOB, GAPDH, and ENO1 expressions indicates that in a high-temperature and low-oxygen environment, anaerobic glycolysis is activated in the liver of T. ovatus, becoming one of the primary energy supply pathways for the body. Similar findings have been reported in previous fish studies. In O. Niloticus [106], largemouth bass (Micropterus salmoides) [107], and large yellow croaker (Larimichthys crocea) [108], the livers provide energy for the body by activating the anaerobic glycolysis pathway to adapt to the hypoxic environment.
With the enhanced glycolysis, transcriptome sequencing results showed that the fatty acid metabolism-related genes long-chain acyl-CoA synthetases 1 (ACSL1), long-chain acyl-CoA synthetase 3 (ACSL3), fatty acid desaturase 2 (FADS2), and elongase of very long chain fatty acid 6 (ELOVL6) were significantly downregulated in the liver of T. ovatus at 34 • C, implying that fatty acid absorption and unsaturated fatty acids synthesis may be reduced. The long-chain fatty acyl-CoA synthetases (ACSLs) family is a crucial enzyme in the synthesis and catabolism of fatty acids and can catalyze the formation of acyl-CoA from free fatty acids (an energy-consuming process) [109]. ACSL1 and ACSL3 are important members of the ACSLs family, and they can accelerate the synthesis and catabolism of fatty acids. ACSL1 can transport saturated fatty acids and unsaturated fatty acids into cells and catalyze them, and this protein is highly expressed in tissues related to energy metabolism, such as fat and liver tissues [110,111]. In the liver tissue of C. idellus, ACSL1 catalyzes the acylation of endogenous fatty acids, activating the β-oxidation pathway and providing energy for fish [112]. ACSL3 is a lipid droplet-associated protein that participates in the absorption of fatty acids and helps maintain lipid homeostasis. Under an acute hyperthermic hypoxic environment, the significant downregulation of ACSL1 and ACSL3 in the liver of T. ovatus may hinder fatty acid absorption and β-oxidation, thereby reducing the liver tissue's demand for O 2 .
In addition, the formation of unsaturated fatty acids is an oxygen-consuming process, which uses fatty acid chains as substrates to synthesize unsaturated fatty acids through a series of alternating reactions of dehydrogenation and carbon chain elongation under the catalysis of FADS and ELOVLs [113,114]. FADS2 and ELOVL6, belonging to the FADS and ELOVLs families, respectively, are rate-limiting enzymes in synthesizing unsaturated fatty acids and catalyzing fatty acid dehydrogenation and carbon chain elongation reactions respectively. In this study, the significant downregulation of expression levels of FADS2 and ELOVL6 indicates that fatty acid dehydrogenation and carbon chain elongation were inhibited in the liver of T. ovatus. That is, the body's ability to synthesize unsaturated fatty acids was reduced. Under high-temperature stress, the expression levels of genes related to unsaturated fatty acid syntheses, such as FADS2, ELOVL2, and ELOVL5, were also significantly decreased in the muscle of fathead minnows (Pimephales promelas) [115]. It can be seen that under high-temperature stress, fatty acid metabolism in fish is inhibited to a certain extent, thus reducing the body's consumption of O 2 and energy to adapt to the hypoxic environment and compensatory changes in energy metabolism caused by high temperature.
In summary, as the water temperature increases, the DO in the water decreases, and the continuous hypoxic environment activates the HIF-1 signaling pathway in T. ovatus, which in turn regulates the cells to produce a series of compensatory responses to hypoxia, such as enhancing O 2 -carrying ability and O 2 -utilization rate, activating the anaerobic glycolysis pathway, and inhibiting fatty acid metabolism. These responses are conducive to the body's effectively using low-concentration oxygen to maintain normal energy metabolism and other physiological processes to adapt to hypoxia caused by a high-temperature environment.

Conclusions
Through liver transcriptomic analysis, our studies found that genes and biological processes related to HSPs, UPS, HIF-1 signaling pathway, and energy metabolism (gly-Animals 2023, 13, 2053 20 of 25 colysis and fatty acid metabolism) play a critical role in response to heat stress in the liver of T. ovatus. Under acute heat stress, the expression of HSPs was increased, and the expression of UPS was inhibited, which helped to reduce the aggregation of denatured and misfolded proteins and further maintain normal protein and liver physiological functions. Moreover, the HIF-1 signaling pathway was activated by a hypoxic environment caused by acute heat stress, further inducing a series of compensatory responses, such as the increase in oxygen-carrying ability and oxygen utilization, activation of the anaerobic glycolytic pathway, inhibition of fatty acid metabolism, and so on, which helps effective utilization of low concentration of oxygen at the organismal level to maintain normal energy metabolism and other physiological processes. The above findings indicated that maintaining normal protein function, energy metabolism, and other physiological functions contribute to heat stress adaption. The present study preliminary explores the molecular response mechanism of T. ovatus under acute heat stress. These genes, pathways, and biological processes related to high temperature screened in this study would be of great value for selecting and breeding heat-tolerant cultivars and the high-quality development of T. ovatus farming.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/ani13132053/s1, Table S1: GO enrichment analysis of three comparison groups at Level 3; Table S2: GO enrichment analysis of four profiles at Level 3.