Transcriptome Analysis of Juvenile Black Rockfish Sebastes schlegelii under Air Exposure Stress

Abstract

The study aimed to uncover the molecular response of juvenile Sebastes schlegelii to air exposure stress by identifying differentially expressed genes (DEGs) that may underlie their anti-stress mechanisms. Juvenile Sebastes schlegelii were subjected to varying durations of air exposure stress. The total RNA was extracted from whole tissues and sequenced using the Illumina NovaSeq 6000 platform. The transcriptome data were analyzed to identify DEGs through pairwise comparisons across a control group and two experimental groups exposed to air for 40 s and 2 min 30 s, respectively. The comparative DEG analysis revealed a significant number of transcripts responding to air exposure stress. Specifically, 5173 DEGs were identified in the 40 s exposure group (BS) compared to the control (BC), 6742 DEGs in the 2 min 30 s exposure group (BD) compared to the control (BC), and 2653 DEGs when comparing the BD to the BS group. Notably, Gene Ontology (GO) analysis showed an enrichment of DEGs associated with peptidase activity and extracellular regions, suggesting a role in the organism’s stress response. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis pointed to the involvement of metabolic pathways, which are crucial for energy management under stress. The upregulation of protein digestion and absorption pathways may indicate a physiological adaptation to nutrient scarcity during stress. Additionally, the identification of antibiotic biosynthesis pathways implies a potential role in combating stress-induced infections or damage. The identified DEGs and enriched pathways provide insights into the complex anti-stress response mechanisms in juvenile Sebastes schlegelii. The enrichment of peptidase activity and extracellular region-related genes may reflect the initiation of tissue repair and immune response following air exposure. The connection between protein digestion and absorption pathways and anti-stress capabilities could be interpreted as a metabolic readjustment to prioritize energy-efficient processes and nutrient assimilation during stress. The role of antibiotic biosynthesis pathways suggests a defensive mechanism against oxidative stress or microbial invasion that might occur with air exposure.

1. Introduction

Sebastes schlegelii, a species within the Scorpaeniformes order, Sebastidae family, and the genus Sebastes, is widely known as the black rockfish [1]. This marine fish is predominantly found in the Yellow Sea, East China Sea, the Korean Peninsula, and the coastal waters of Japan [2]. The black rockfish holds significant commercial value and is a key species for cage aquaculture, particularly in the northern coastal regions of China [3]. However, the recent expansion in the cultivation and seedling production of the black rockfish has highlighted the several challenges within the aquaculture sector. A particularly pressing issue occurs during the seedling rearing phase, where juveniles, averaging 2–3 cm in length, are exposed to air for extended periods during activities such as pond draining and splitting. This exposure can lead to significant stress-induced mortality, a condition commonly referred to as “gill pop**” within the industry. This phenomenon is not unique to the black rockfish; it is also observed in the seedling rearing processes of other marine fish species, including the Atlantic croaker [4]. The high mortality rates associated with these conditions present a critical challenge that requires advancements in breeding technology to reduce these losses.

To mitigate the stress behaviors in juvenile fish subsequent to air exposure, prevalent strategies include the application of anesthetics, vitamin supplementation, and the adjustment of water temperature to ameliorate distress [5]. However, from a molecular breeding perspective, there is an urgent need to investigate the molecular mechanisms underlying the stress response in fish following air exposure. Unraveling these mechanisms is essential for the identification of parental candidates and lineages with inherent stress tolerance. This knowledge could significantly accelerate the breeding process, thereby enhancing the resilience of Sebastes schlegelii to the challenges of air exposure.

Previous research on the fry cultivation stage of Brachymystax lenok has identified the phenomenon of “gill pop**,” which has been predominantly associated with elevated levels of ammonia nitrogen and infestations by the ciliate, Ichthyophthirius multifiliis. Despite the avoidance of these specific factors in the rearing of marine fish fry, instances of “gill pop**” mortality induced by air exposure remain prevalent [6]. Studies on marine fish subjected to air exposure stress have shown that when the dissolved oxygen content in the water falls below 4.0 mg/L, Sebastes schlegelii exhibits stress responses, including behaviors such as tilting and floating heads [7]. A significant increase in leukocyte counts is observed, which is indicative of an enhanced immune response to counteract the damage caused by hypoxia [8]. Gao et al. conducted research on the critical oxygen partial pressure and dissolved oxygen content at equilibrium, revealing that Sebastes schlegelii is particularly sensitive to hypoxic conditions, which can lead to substantial physiological, biochemical, and morphological alterations [9]. However, these changes are found to be reversible upon the restoration of oxygen levels following a period of low oxygen stress [10]. In acute low-oxygen stress experiments, Mu et al. observed a notable upregulation in the expression of heat shock protein 70 (HSP70) in Sebastes schlegelii. Additionally, another study, using the Procambarus clarkii (red swamp crayfish) as a model, conducted the first assessment of the histological changes, antioxidant capacity, and molecular mechanisms of the immune response in crustaceans under air exposure stress. The results indicated that air exposure stress could lead to oxidative damage, potentially culminating in disease or death [11]. To date, the specific molecular mechanisms that precipitate “gill pop**” in juvenile Sebastes schlegelii post-air exposure have not been elucidated in the literature.

Transcriptome sequencing methods have been employed to detect the gene expression profiles of aquatic animals under air exposure [12,13,14]. For instance, the transcriptome of the Eriocheir sinensis (Chinese mitten crab) under air exposure is a common stressor for this species, particularly during the periods of harvest and transportation [12]. Another study investigated the effects of air exposure on the hepatopancreatic tissue of clams, including changes in tissue structure, physiological responses (mainly including respiratory metabolism, antioxidant capacity, and immune function), and alterations in the transcriptome [13]. Furthermore, the majority of the literature has focused on the physiological responses, behavioral changes, and survival rates of fish under air exposure [15,16,17,18,19,20].

Utilizing high-throughput sequencing technology, we conducted transcriptome analysis at various time points following air exposure in juvenile Sebastes schlegelii. This approach facilitated the screening and identification of the genes and regulatory pathways that are pertinent to stress responses. The objective of this research was to provide a theoretical framework and empirical data that can inform and enhance the healthy rearing practices for juvenile Sebastes schlegelii. By elucidating the molecular-level effects of varying durations of air exposure on stress responses, this study not only contributes to the fundamental understanding of the underlying mechanisms but also provides valuable reference data for the development of robust cultivation strategies for Sebastes schlegelii.

2. Materials and Methods

2.1. Experimental Materials

The juvenile Sebastes schlegelii utilized in this experiment were sourced from Qingdao Clear Source Marine Biotechnology Co., Ltd. (Qingdao, China). These juveniles exhibited an average body length of 2.5 ± 0.24 cm and an average body weight of 2.6 ± 0.31 g. Upon arrival, the individuals were inspected and confirmed to be in good health, with no signs of mechanical injury or disease. The experimental setup employed a standard flow-through aquaculture system, with a water flow rate maintained above 1000 mL/min to ensure optimal water quality. Before the commencement of the experiment, the fish underwent a one-week acclimation period to stabilize their condition in the new environment.

Throughout the acclimation and experimental phases, the water quality parameters were stringently controlled, with the water temperature regulated at 18 ± 1 °C, the dissolved oxygen content maintained at 8.6 ± 0.3 mg/L, the ammonia nitrogen levels kept below the threshold of <0.1 mg/L, and the salinity levels adjusted to fluctuate within the range of 29 to 31 practical salinity units (PSUs). Only those individuals deemed healthy based on the predefined criteria were randomly selected and subsequently included in the experimental procedures.

2.2. Air Exposure Stress Experiment

The cultivation of Sebastes schlegelii was conducted in three square arc-shaped indoor cement pools, each with a surface area of 25 m² and a water depth ranging from 50 to 80 cm. The experimental design was structured as follows:

Control Group (BC): three fish from each pool were carefully collected using enzyme-free forceps, swiftly transferred into cryotubes, and flash-frozen in liquid nitrogen for subsequent analysis.

Forty-second Air Exposure Group (BS): A sample of ten juveniles from each pool was netted and subjected to 40 s of air exposure. Subsequently, three individuals were expeditiously selected using enzyme-free forceps, securely placed in cryotubes, and flash-frozen in liquid nitrogen. The remaining juveniles were then returned to their respective pools for continuous monitoring.

Two-minute 30-s Air Exposure Group (BD): A similar procedure was followed for this group, with ten juveniles from each pool being exposed to air for 2 min and 30 s. Post-exposure, three fish were rapidly chosen, cryopreserved, and flash-frozen in liquid nitrogen. The remaining fish were also returned to their pools for further observation.

Each experimental group was replicated three times to ensure the reliability of the results. After the collection of samples, the cryogenically frozen fish were transported and temporarily stored at −80 °C in an ultra-low temperature freezer for preservation.

2.3. Total RNA Extraction

For each experimental group, including three replicates each for the Control Group (BC1, BC2, BC3), 40 s Air Exposure Group (BS1, BS2, BS3), and 2 min 30 s Air Exposure Group (BD1, BD2, BD3), the total RNA was extracted from a pooled sample of three juvenile Sebastes schlegelii individuals, which were homogenized together. The RNA extraction was executed utilizing the Trizol reagent, a widely recognized method for RNA isolation.

The purity and concentration of the extracted total RNA were determined using a Nano 6000 ultramicro spectrophotometer, with the absorbance ratios at OD260 nm/OD280 nm within the range of 1.8–2.2, indicative of RNA purity, and OD260 nm/OD230 nm ≥ 2.0, suggesting minimal contamination by other biomolecules. The total RNA concentration was required to exceed 250 ng/μL to ensure a sufficient yield for downstream applications. Furthermore, the integrity of the total RNA was ascertained using the Agilent 2100 Bioanalyzer system, which provides a precise measurement of the RNA quality. This system rendered the traditional agarose gel electrophoresis method unnecessary, as it offers a comprehensive analysis that includes the evaluation of the 28S and 18S ribosomal RNA bands, ensuring the absence of RNA degradation or contamination.

2.4. mRNA Library Construction and Sequencing

RNA quantification was performed with high precision using the Agilent 2100 Bioanalyzer system, coupled with a specific reagent kit provided by Agilent Technologies (Santa Clara, CA, USA). To isolate mature mRNA, a purification and enrichment process was undertaken, employing magnetic beads conjugated with oligo(dT) primers to selectively bind and extract polyadenylated mRNA from the total RNA sample.

The construction of the sequencing library was carried out in accordance with the established protocols for library preparation. To enrich for cDNA fragments within the size range of 370–420 base pairs (bp), the library fragments were subjected to size selection using the AMPure XP system from Beckman Coulter (Brea, CA, USA). Following this, PCR amplification was conducted to amplify the selected cDNA fragments. The resulting PCR products were further purified using AMPure XP beads to finalize the library construction.

After the library successfully passed a stringent quality assessment, it was sequenced on the Illumina NovaSeq 6000 platform, which generated 150 bp paired-end reads. The raw sequencing data, in the form of electrical signals, were then processed using specialized computer software to convert these signals into a series of sequencing peaks, thereby yielding the nucleotide sequence information for the targeted cDNA fragments.

2.5. Transcriptome Sequence Quality Control and Assembly

The raw sequencing data underwent initial processing using Perl scripts to systematically eliminate the sequences that contained ambiguous ‘N’ bases, exhibited low-quality base calls, or included read adapter contamination. Concurrently, a suite of quality control metrics were computed, including Q20 and Q30 scores, which denote the percentage of bases with a quality score above 20 and 30, respectively, as well as the GC content, which is the proportion of guanine and cytosine nucleotides within the sequence.

Following the initial data cleaning, the Trinity software suite (version 2.6.6) was deployed for de novo assembly analysis. This step aimed to consolidate the processed sequences into a comprehensive set of transcripts, from which the longest sequence representative of each transcript (termed unigene) was extracted to serve as the reference sequence for the subsequent analyses [21].

To evaluate the quality and integrity of the Trinity assembly, the Benchmarking Universal Single-Copy Orthologs (BUSCO) software (version 5.7.0) was employed. BUSCO assesses the completeness and accuracy of the assembly by measuring the ratio of complete, single-copy orthologs present in the assembly against a database of conserved orthologs.

Ultimately, the functional annotation of the reference sequences was conducted by aligning them against a panel of well-established databases, including the NCBI non-redundant protein sequences (Nr), NCBI non-redundant nucleotide sequences (Nt), Pfam (a database of protein families), Gene Ontology (GO), Clusters of Orthologous Groups of proteins (KOG/COG), Swiss-Prot (a curated database of protein sequences), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) Ortholog database.

2.6. Differential Gene Expression

Differential gene expression analysis was conducted using the DESeq (version 1.20.0) tool, which is specifically designed to identify differentially expressed genes (DEGs) across various conditions. The analysis was stringent, employing a threshold of a Q-value less than 0.05 and a fold change criterion of |log2(fold change)| greater than 1. This approach facilitated the identification of genes whose expression levels varied significantly between the experimental groups. It is important to note that a smaller adjusted p-value (padj) denotes a higher statistical significance of the observed difference in gene expression [22].

The DEGs identified through this process were further subjected to functional annotation using the Gene Ontology (GO) system, which categorizes genes based on their molecular functions, biological processes, and cellular components (http://www.geneontology.org/, accessed on 13 May 2023). Additionally, the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was performed to identify the biological pathways that these DEGs are involved in (http://www.genome.jp/kegg/, accessed on 14 May 2023).

3. Results

3.1. Observation Results of Air Exposure Stress

The experimental observations detailed in

Table 1. Survival of juvenile Sebastes schlegelii at different times of departure from water.

Item	Control Group (BC)			40 s Air Exposure Group (BS)			2 min 30 s Air Exposure Group (BD)
Group	BC1	BC2	BC3	BS1	BS2	BS3	BD1	BD2	BD3
survival rate	100%	100%	100%	100%	100%	100%	-	-	-
mortality rate	-	-	-	-	-	-	100%	100%	100%

reveal a significant contrast in the survival and recovery of juvenile Sebastes schlegelii subjected to different durations of air exposure. Specifically, when the juveniles were exposed to air for a brief period of 40 s and then promptly returned to their cultivation tanks, no mortality was observed. These individuals demonstrated a remarkable resilience, as they were able to recover fully to their pre-exposure conditions. In stark contrast, when the exposure duration was extended to 2 min and 30 s, the outcome was markedly adverse; all individuals perished and were unable to recover upon being returned to the cultivation tanks.

3.2. Sequencing Results and Transcript Assembly

The raw sequencing data, comprising the initial reads obtained from the sequencing platform, were subjected to preprocessing using Perl scripting. This step was essential to eliminate sequences that contained ambiguous ‘N’ bases, exhibited signs of low quality, or were identified as read adapter contaminants. Following this initial data cleaning, the Trinity software (version 2.15.1) was utilized for de novo assembly of the cleaned reads. The assembly process culminated in the generation of 92,469 transcripts and 58,397 unigenes. The N50 and N90 values for the assembled transcripts were determined to be 2772 base pairs (bp) and 2478 bp, respectively, while the corresponding values for the N90 metric were 631 bp and 534 bp. As presented in

Table 2. Statistics of transcriptome sequencing in Sebastes schlegelii.

Length Transcript	Transcript	Unigene
≥500 bp	66,102	39,618
≥1000 bp	44,384	24,848
Count	92,469	58,397
max length/bp	56,776	56,776
min length/bp	301	301
mean length/bp	1606	1420
median length/bp	936	781
N50	2772	2478
N90	631	534

, these assembly metrics are indicative of a high-quality assembly, which is a prerequisite for accurate gene function annotation and subsequent analysis of differentially expressed genes (DEGs).

The experimental design encompassed three distinct groups: a normal control group (BC), a group subjected to air exposure for 40 s (BS), and a group exposed to air for a more extended period of 2 min and 30 s (BD). Each group was represented by three parallel samples, yielding a total of nine samples for transcriptome sequencing. The quality of the sequencing data was assessed by calculating the Q20 and Q30 base percentages, which were found to be above 97.61% and 93.51%, respectively, for all the samples. These metrics are reflective of a high-quality sequencing output. The average GC content ratio varied slightly between the groups, with values of 48.26% for BC, 50.43% for BS, and 49.24% for BD. The average number of clean reads, which refers to the high-quality reads post-trimming and filtering, was 20.81 × 10⁶ for BC, 22.34 × 10⁶ for BS, and 21.69 × 10⁶ for BD. A comprehensive summary of the sequencing data for each sample is provided in

Table 3. Statistical table of the sequencing data.

Samples	Clean Reads	Clean Bases	Error Rate	Q20 (%)	Q30 (%)	GC (%)
BC1	21,438,555	6.4	0.03	97.67	93.62	47.85
BC2	20,241,671	6.1	0.03	97.61	93.51	47.67
BC3	20,761,241	6.2	0.03	97.81	93.98	49.25
BS1	22,557,563	6.8	0.03	97.77	93.92	50.52
BS2	22,889,460	6.9	0.03	97.73	93.66	50.26
BS3	21,915,676	6.6	0.03	97.83	93.96	50.52
BD1	20,364,363	6.1	0.03	97.68	93.53	49.72
BD2	22,770,664	6.8	0.03	97.7	93.65	49
BD3	21,952,790	6.6	0.03	97.78	93.78	49.01

.

3.3. Gene Function Annotation

The unigene information derived from the transcriptome sequencing was subjected to alignment using a suite of established databases to facilitate gene function annotation. These databases included the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes Ortholog (KO), Clusters of Orthologous Groups (KOG), NCBI non-redundant protein sequences (Nr), NCBI non-redundant nucleotide sequences (Nt), Pfam, and Swiss-Prot. A comprehensive annotation was achieved for a total of 58,397 unigenes, with the quantity and percentage of annotations for each database detailed in

Table 4. Statistics of unigenes annotated.

Annotated Database	Number of Unigenes	Percentage/%
GO	24,555	42.04
KO	18,210	31.18
KOG	12,188	20.87
NR	30,544	52.3
NT	42,859	73.39
Pfam	24,559	42.05
Swiss-Prot	25,755	44.1
Total	58,397	100

. Following homology alignment with the Nr database, a significant number of unigenes—specifically, 30,544—were found to exhibit homology with known genes, thereby providing a valuable reference for further analysis. The annotation information specific to Sebastes schlegelii is graphically represented in Figure 1. Among the species with the highest annotation percentages, Yellow perch (Perca flavescens) was the most prevalent, accounting for 29.8% of the annotations, followed by Asian sea bass (Lates calcarifer) at 8.2%, Large yellow croaker (Larimichthys crocea) at 7.7%, Greater amberjack (Seriola dumerili) at 5.2%, and Yellowtail amberjack (Seriola lalandi) at 4.5%.

3.4. GO Functional Analysis of Unigenes

The unigene information obtained from the transcriptome sequencing was further analyzed using Gene Ontology (GO) enrichment analysis to classify the genes based on their associated biological roles. This analysis successfully annotated 24,555 unigenes, categorizing them into three principal GO domains: biological processes (BP), molecular functions (MF), and cellular components (CC). The distribution of annotated unigenes across these domains was examined, with a notable predominance observed in the biological processes category (BP), as depicted in Figure 2.

3.5. KOG Functional Analysis of Unigenes

The unigene dataset, derived from the transcriptome sequencing of Sebastes schlegelii, was further analyzed and annotated using the Clusters of Orthologous Groups of proteins (KOG) functional classification system. This classification process successfully annotated 12,188 unigenes, providing insights into their evolutionary relationships and conserved functions. Among the annotated unigenes, the largest proportion was associated with signal transduction mechanisms, which constituted 17.6% of the total annotated genes. This finding underscores the importance of signal transduction in the cellular response to environmental stimuli. The second-largest category was posttranslational modification, protein turnover, and chaperones, which made up 14.9% of the annotations, highlighting the significance of protein homeostasis in cellular processes. Transcription-related unigenes were also prominently represented, accounting for 11.7% of the annotations, which is reflective of the critical role of gene regulation in cellular function.

In contrast, the category with the smallest proportion was cell motility, which accounted for less than 1% of the annotated unigenes. This suggests that, within the context of the studied transcriptome, cell motility may be less prominent or that the genes involved in this process are less conserved across species, leading to lower annotation rates. The distribution of annotated unigenes across the various KOG categories is graphically represented in Figure 3, offering a visual summary of the functional diversity of the Sebastes schlegelii transcriptome.

3.6. KEGG Pathway Analysis of Unigenes

Utilizing the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, an in-depth functional annotation of the Sebastes schlegelii unigene dataset was performed. This comprehensive annotation effort resulted in the classification of 20,214 unigenes across five major KEGG-defined categories, as illustrated in Figure 4: Cellular Processes (A), Environmental Information Processing (B), Genetic Information Processing (C), Metabolism (D), and Organismal Systems (E). Of particular note, the category with the most substantial number of annotated unigenes was signal transduction, which falls under the Environmental Information Processing category (B). This category alone accounted for 2929 unigenes, reflecting the complexity and diversity of the signal transduction pathways in Sebastes schlegelii and their potential significance in the organism’s adaptation and response to environmental stimuli.

3.7. Statistical Analysis of Differentially Expressed Genes (DEGs)

The DESeq (version 1.20.0) software package was employed to identify differentially expressed genes (DEGs) with a stringent threshold set at a Q-value of 0.05 and a fold change of 2. By integrating the data from three biological replicates, the comparative analysis between the experimental groups yielded the following findings: the BS (40 s air exposure) versus BC (control) group comparison identified 5173 DEGs, comprising 2609 upregulated and 2564 downregulated genes; the BD (2 min and 30 s air exposure) versus BC group comparison revealed 6742 DEGs, with 2459 upregulated and 4283 downregulated genes; and the BD versus BS group comparison disclosed 2653 DEGs, including 497 upregulated and 2156 downregulated genes (

Table 5. Statistical table of the number of differentially expressed genes.

Compare	DEG Number	Upregulated Number	Downregulated Number
BS vs. BC	5173	2609	2564
BD vs. BC	6742	2459	4283
BD vs. BS	2653	497	2156

and Tables S1–S3).

To visualize the overlap of DEGs between the comparison groups, Venn diagrams were constructed based on the specified criteria for each group. This analysis identified a total of 9149 DEGs across all the groups. Notably, the BS versus BC and BD versus BC comparisons shared 2752 common DEGs. The BD versus BC and BD versus BS comparisons had 764 DEGs in common, while the BD versus BS and BS versus BC comparisons shared 809 DEGs. Crucially, 547 DEGs were found to be consistently differentially expressed across all three comparison groups (Figure 5).

3.8. GO Functional Enrichment Analysis of Differentially Expressed Genes (DEGs)

In the comparative analysis of the BS (40 s air exposure) versus BC (control) group, which encompassed 5173 DEGs with 2609 upregulated and 2564 downregulated genes, the GO functional analysis indicated significant enrichment across 134 subcategories within the three main ontologies: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). Notably, within the Cellular Component (CC) ontology, enrichment was observed for terms such as Extracellular matrix (GO:0031012), Chromosome (GO:0005694), and Extracellular region (GO:0005576). Within the Molecular Function (MF) ontology, the term Peptidase activity (GO:0008233) was significantly enriched (Figure 6A).

In the BD (2 min and 30 s air exposure) versus BC comparison, which identified 6742 DEGs comprising 2459 upregulated and 4283 downregulated genes, the Molecular Function (MF) ontology again highlighted Peptidase activity (GO:0008233). Similarly, within the Cellular Component (CC) ontology, the terms Extracellular matrix (GO:0031012) and Chromosome (GO:0005694) were significantly enriched (Figure 6B).

For the BD versus BS comparison, which revealed 2653 DEGs with 497 upregulated and 2156 downregulated genes, the Cellular Component (CC) ontology showed enrichment for the term Extracellular region (GO:0005576). Within the Biological Process (BP) ontology, the terms Cell adhesion (GO:0007155) and Nervous system process (GO:0050877) were identified as enriched. Concurrently, the Molecular Function (MF) ontology once more highlighted Peptidase activity (GO:0008233) (Figure 6C).

3.9. KEGG Pathway Enrichment Analysis of Differentially Expressed Genes (DEGs)

In the BS (40 s air exposure) versus BC (control) group comparison, the enriched KEGG pathways identified were Metabolic pathways (ko01100), Protein digestion and absorption (ko04974), and Biosynthesis of antibiotics (ko01130) (Figure 7A). The notable enrichment in metabolic pathways suggests that there may be significant alterations in the cellular metabolic processes in response to air exposure stress.

For the BD (2 min and 30 s air exposure) versus BC comparison, the enriched KEGG pathways included Biosynthesis of secondary metabolites (ko01110), Protein digestion and absorption (ko04974), Biosynthesis of antibiotics (ko01130), and Adrenergic signaling in cardiomyocytes (ko04261) (Figure 7B). The emphasis on the biosynthesis of secondary metabolites implies that there may be changes in the cellular response mechanisms, potentially involving the production of compounds that help the organism cope with stress.

In the BD versus BS comparison, the enriched KEGG pathways were Metabolic pathways (ko01100), Biosynthesis of secondary metabolites (ko01110), Neuroactive ligand–receptor interaction (ko04080), Microbial metabolism in diverse environments (ko01120), and Protein digestion and absorption (ko04974) (Figure 7C). The focus on metabolic pathways in this comparison implies that there are likely to be changes in the organism’s energy and substance metabolism, which could be crucial for its survival under stress conditions.

4. Discussion

During the cultivation or breeding processes of marine fish fry, such as Sebastes schlegelii, the stress response induced during water transfer or relocation can profoundly impact their survival rates [23]. The prolonged exposure to air can result in severe consequences, including the phenomenon known as “gill pop**”, which poses significant challenges to the viability of aquaculture operations. To address these challenges, high-throughput transcriptome sequencing has emerged as an indispensable tool. It enables the rapid and comprehensive capture of gene expression profiles under defined conditions, effectively bridging the gap between genetic information and the intricate details of protein transcription [24,25]. To elucidate the stress response mechanisms in Sebastes schlegelii during the critical phase of water transfer, this study undertook transcriptome sequencing on samples that were subjected to varying durations of air exposure, categorized into the BC (control), BS (40 s exposure), and BD (2 min and 30 s exposure) groups.

The experimental observations have delineated two pivotal time points in the survival response of Sebastes schlegelii to air exposure: the species can fully recover upon return to the culture tank after a 40 s exposure, whereas a duration of 2 min and 30 s results in complete mortality. These findings underscore the sensitivity of this species to air exposure and the need for careful handling during the transfer processes. Previous research on the impact of MS-222 anesthesia on the survival of Nile tilapia post-air exposure has demonstrated extended survival times, exceeding 7 h at a temperature of 15 °C with an anesthetic concentration of 100 mg/L [26]. These results highlight the potential of anesthetics in mitigating stress during fish handling, which could have significant implications for aquaculture practices.

Given the severe mortality associated with “gill pop**” in Sebastes schlegelii following air exposure, this study provides a foundational understanding of the differential gene expression that occurs post-exposure. This knowledge can serve as a springboard for future research aimed at identifying the specific anesthetics or anti-stress agents that could potentially extend the survival time of Sebastes schlegelii after air exposure. By doing so, it may be possible to develop the strategies to reduce the risks inherent in handling operations such as tank inversion and fry separation, thereby improving the overall survival rate and welfare of the species in aquaculture settings.

This study conducted a comparative analysis across different time groups—BS (40 s air exposure) versus BC (control), BD (2 min and 30 s air exposure) versus BC, and BD versus BS—and identified a substantial number of differentially expressed genes (DEGs). Within the BS vs. BC group, the analysis revealed 2609 upregulated genes and 2564 downregulated genes. The BD vs. BC group exhibited the highest number of DEGs, with downregulated genes outnumbering upregulated genes by approximately twofold. In contrast, the BD vs. BS group displayed a relatively smaller number of DEGs, with downregulated genes being roughly four times the number of upregulated genes. Notably, as the duration of air exposure increased, a progressive suppression of a multitude of genes was observed, which may underlie the eventual mortality observed under water stress conditions [27].

Gene Ontology (GO) enrichment analysis of the DEGs uncovered sustained alterations in the genes associated with peptidase activity, as well as those involved in the extracellular region and extracellular matrix, in response to water stress in juvenile Sebastes schlegelii. These findings are consistent with previous studies that documented damage to gill filaments, epithelial cell flattening, and mitochondrial changes following acute water stress [28]. In the BD vs. BC group, the most significantly enriched DEGs were linked to peptidase activity and the extracellular region, suggesting that the genes related to enzyme activities are among the primary responders to water stress [29]. When fish are exposed to air, they experience a sudden and severe stress due to the lack of water, which is essential for their respiratory process [17]. This immediate stress triggers a rapid response from the body to meet the energy demands for survival. Peptidases, as enzymes that break down peptides into amino acids, play a critical role in this process by accelerating protein catabolism to quickly release energy [30,31]. As previously reported, dipeptidyl peptidase 3 (DPP3), a zinc-dependent aminopeptidase, moonlights as a regulator of the cellular oxidative stress response pathway, e.g., the Keap1-Nrf2 mediated antioxidative response [32]. Furthermore, the genes associated with the nervous system were notably affected in the BD vs. BS group, implying that prolonged air exposure can lead to severe neural system disruption, potentially resulting in irreversible damage and a diminished capacity for recovery after 2 min and 30 s of air exposure.

The KEGG pathway enrichment analysis revealed a notable concentration of metabolic pathways and protein digestion and absorption pathways across all three comparison groups. This suggests that Sebastes schlegelii may modulate its substance or energy metabolism in response to stress, as a strategy to mitigate stress-induced damage and regulate the pathways associated with apoptosis, thereby alleviating stress-induced cell apoptosis [33]. In the BD vs. BS group, there was a significant enrichment in the pathways related to the biosynthesis of secondary metabolites, indicating a rapid accumulation of stress-related metabolites and the onset of apoptosis, which may culminate in irreversible death [34].

5. Conclusions

In summary, the study’s findings underscore the significant impact of air exposure duration on the gene expression profiles of juvenile Sebastes schlegelii, which is crucial for understanding the species’ survival mechanisms and stress response. Through transcriptome sequencing analysis, we have compared three groups: BC (control group with no air exposure), BS (group exposed to air for 40 s), and BD (group exposed to air for 2 min and 30 s). This comparative analysis has unveiled the potential molecular pathways that are affected by water stress and influence gene expression. The discovery of specific molecular signatures in response to stress, such as peptidase activity and changes in protein digestion, is intriguing. While these findings could be indicative of a secondary stress response, they also prompt the hypothesis that primary stress responses, possibly involving the hypothalamic–pituitary–adrenal (HPA) axis and cortisol production, may trigger a cascade of metabolic changes. It is plausible that cortisol, a hormone known to mediate the body’s response to stress, could alter metabolic pathways and subsequently affect processes like peptidase activity and protein digestion.

Our study not only deepens the understanding of how Sebastes schlegelii adapts to and copes with air exposure stress but also sets a foundation for future research. The identification of differentially expressed genes (DEGs) and their potential roles in the stress response provides a scientific basis for develo** targeted interventions. These interventions could aim to mitigate the adverse effects of water stress by modulating the expression or activity of these DEGs, thereby enhancing the resilience of juvenile Sebastes schlegelii in aquaculture.