Speckled blue grouper (Epinephelus cyanopodus, body weight: 35.0 ± 1.0 g, 4-month-old females) were acclimated for 3 weeks in 300 L recirculating aquaculture tanks. Fish were acclimated without bacterial or viral infection in a freshwater environment of 22 ± 2°C with the normal oxygen level DO (DO:6.0 ± 0.1 mg/L). The hypoxic treatment followed our previous studies (Wu et al., 2023, 2024). Briefly, we stopped the circulating water, sealed the bucket opening with plastic film and injected nitrogen at a constant rate. Within 1 h, we reduced the normal oxygen level to a low DO level (0.5 ± 0.1 mg/L) (Liang et al., 2024), and then maintained this DO concentration. When we reduced the dissolved oxygen from the normoxic level to 0.5 ± 0.1 mg/L, the floating head frequency, respiratory rate and tail flap frequency initially showed a sharp increase. All these manifestations indicated that the grouper was under hypoxic stress. Under normoxic conditions, 3 fish were randomly selected for head kidney (HK) tissue dissection and named Ctrl (control). When hypoxia persisted for 6 h, HK tissues were randomly collected from 3 fish and named Hyp. To explore the physiological processes of grouper under the dual stress of hypoxia and pathogenic bacteria, we conducted the following experiments. A total of 504 fish were divided into 2 groups (control group and hypoxia group), with 3 replicates in each group, and 84 fish in each replicate. Each replicate was subjected to stress from LPS, Vibrio harveyi (V. harveyi), Poly (I:C), and PBS respectively. Fish cultured under normoxia and hypoxia for 6 h were randomly grouped and intraperitoneally injected with 100 μL LPS (50 μg/mL) (Sandamalika et al., 2022), 100 μL V. harveyi (1 × 10⁸ CFU/mL), 100 μL Poly(I:C) (50 μg/mL) (Thanasaksiri et al., 2014), or 100 μL PBS (control), respectively. The injection dose of V. harveyi was determined by testing 1 × 10⁷, 3 × 10⁷, 9 × 10⁷, 2.7 × 10⁸ and 8.4 × 10⁸ CFU/mL ( Supplementary Figure S1 ) in 150 fish. At 0, 3, 6, 12, 24, 48, and 72 h post-injection (p.i.), liver and head kidney tissues were harvested from 3 individual fish, then stored at - 80°C before use. To study the transcriptional expression of EcHIF-1α in healthy speckled blue grouper, a variety of tissues including head kidney (HK), spleen (SP), peripheral blood (PBL), posterior kidney (PK), heart (HE), thymus (TH), muscle (MU), brain (BRA), gills (GI), liver (LI), and stomach (STO) were carefully dissected and gathered from 3 distinct fish and preserved using liquid nitrogen.

Total RNA was extracted from samples using the TRIzol Kit (Promega, USA) following the manufacturer’s instructions. RNA integrity, quality, and quantity were assessed using an Agilent 2100 bioanalyzer, NanoPhotometer spectrophotometer, Qubit 2.0 Fluorometer, and agarose gel electrophoresis. Poly-A-tailed mRNA was enriched using Oligo(dT)-coated magnetic beads and then fragmented by sonication. First-strand cDNA synthesis was initiated using the fragmented mRNA as a template, random oligonucleotide primers, and the M-MuLV reverse transcriptase system. RNase H was employed to degrade the RNA component, enabling second-strand cDNA synthesis with DNA polymerase I and dNTPs. The cDNA fragments were purified using a QIAquick PCR extraction kit (QIAGEN). The purified cDNA fragments were washed with EB buffer solution for end repair with poly (A), and then attached to the sequencing adapter. After agarose gel electrophoresis, about 300 bp cDNA fragments were purified and enriched by PCR to construct the final cDNA library. The cDNA sequencing libraries subsequently underwent sequencing utilizing the Illumina HiSeq™ 2500 system, carried out by Genedenovo Biotechnology Co., Ltd., located in Guangzhou, China.

To identify and remove non-coding RNAs, the clean reads were aligned with data from the GenBank and Rfam databases. The DEGs among the 3 groups were determined by the DESeq2 algorithm. Genes with more than 2-fold the Change (|log2 (Fold Change) |>1) and a false detection rate (FDR) < 0.05 were marked as significant DEGs. The DEGs were subjected to GO and KEGG enrichment analyses. In review enriched GO terms or pathways were determined employing an FDR-adjusted p value threshold of 0.05. Further GO and KO enrichment analysis (p < 0.05) was conducted using a hypergeometric distribution test.

The cDNA sequence of EcHIF-1α gene was subsequently identified using transcriptome analysis. The EcHIF-1α-specific primers were designed using Primer Premier 5.0 and presented in Table 1 . The PCR product was sequenced at BGI TECH SOLUTIONS CO., LIMITED (Guangzhou, China) using an ABI 3730XL automated sequencer. The predicted amino acid (aa) sequence of EcHIF-1α was assessed via the Expert Protein Analysis System on http://www.expasy.org/. The molecular weight (MW) and theoretical isoelectric point (PI) were calculated using the ProtParam tool at https://web.expasy.org/protparam/?tdsourcetag=s_pcqq_aiomsg. The Simple Modular Architecture Research Tool (SMART) from http://smart.emblheidelberg.de/ was applied for protein domain prediction. Multiple sequence alignments were performed using DNAMAN software (LynnonBiosoft, USA), and the phylogenetic tree for HIF-1α sequences was generated using the neighbor-joining (N-J) method in MEGA 6.0, with reliability evaluated through 1000 bootstrap replicates.

To verify the RNA-seq findings, 12 key DEGs were randomly chosen for qRT-PCR validation. DEGs-specific primers were designed with Primer Premier 5.0, with β-actin and Ef1α serving as reference genes ( Supplementary Table S1 ). The qRT-PCR assay was conducted on the ABI 7500 apparatus (Applied Biosystems, USA) utilizing Yeasen’s SYBR qPCR Master Mix. The reaction mixture with a volume of 20 μL was prepared as follows: 2 μL of each primer (2 μM), 10 μL of 2× Hieff^® qPCR SYBR Green Master Mix (Low Rox Plus) (Yeasen Biotechnology, China), 3 μL of diluted cDNA, and 3 μL of DEPC-treated water were combined. The thermal cycling protocol consisted of an initial denaturation at 95°C for 30 seconds, followed by 40 cycles of 95°C for 10 seconds, 60°C for 30 seconds, and a final extension. Ultimately, a melting curve analysis was performed following these steps: 95°C for 15 s, 60°C for 60 s, and 95°C for 15 s. This step was conducted to identify any non-specific amplification or primer dimers (Wu et al., 2024). All DEG expression levels were calculated by 2^−ΔΔCT (Schmittgen and Livak, 2008).

In healthy fish, the tissue distribution of EcHIF-1α was analyzed across various tissues, including HK, PBL, SP, MU, IN, PK, TY, HE, STO, LI, BR, and GI, by qRT-PCR. Relative expression analyses of tissue distribution were normalized against β-actin and Ef1α by the 2^−ΔCT method (Gao et al., 2020). To further investigate the expression of EcHIF-1α against the LPS, V. harveyi, Poly(I:C), and PBS, the differential expression levels of EcHIF-1α were computed using β-actin and ef1α as internal standards, and the findings were benchmarked against the corresponding control group to assess the variation in gene expression. The EcHIF-1α expression was calculated by 2^−ΔΔCT method.

Data handling was performed using Microsoft Office Excel 2021. All statistical analyses and figures were generated with GraphPad Prism 9.3.1. Data analyzed by one-way ANOVA are presented as the mean ± SD, and *p< 0.05 and **p< 0.01 were used to defined the statistical significance.

RNA-Seq of 6 head kidney tissue samples yielded 33,586,082 to 55,062,978 raw reads. After removing low-quality sequences, we obtained 33,482,278 to 54,827,026 clean reads for further analysis. The parameter statistics of the clean reads in each group were Q20, 96.28 − 98.77%; Q30, 88.69 – 96.26%; GC content, 47.97 − 49.25%, and total mapped rate, 81.13 – 84.16% ( Table 2 ). The above results suggest that the RNA-seq data are reliable and suitable for subsequent analysis.

The comparison of Ctrl and Hyp revealed 2,887 DEGs, with 926 up-regulated and 1,961 down-regulated DEGs ( Supplementary Figures S2A, B ). The 2,887 DEGs were analyzed for GO and KEGG enrichment. As shown in Figure 1A , the comparison of Ctrl and Hyp showed GO enrichment results. GO enrichment between Ctrl and Hyp was mainly concentrated in biological processes (BP), molecular functions (MF), and cellular components (CC).

In the BP category, the top three terms were cellular process, metabolic process, and biological regulation. In the CC category, the top three terms were cellular anatomical entity, protein-containing complex, and virion component. In the MF category, binding, catalytic activity, and transcription regulator activity ranked as the top three terms. The top 20 enriched GO terms were primarily concentrated in the MF category. The top 3 GO terms, ranked by q-values, were transmembrane signaling receptor activity, G protein-coupled receptor activity, and molecular transducer activity ( Figure 1B ).

After acute hypoxia, KEGG pathway analysis revealed that multiple pathways were enriched, including six categories: metabolism, human diseases, organic systems, genetic information processing, cellular processes, and environmental information processing. Among them, the pathways related to metabolism, human diseases and organic systems were 12, 11 and 10 respectively. These were mainly associated with lipid and aa metabolism, infectious diseases, and immune regulation ( Figure 2A ).

To explore the biological pathways involved in hypoxic adaptation within HK tissues, the top 20 enriched KEGG pathways were analyzed. Key pathways related to hypoxia included cytokine-cytokine receptor interaction, PPAR signaling pathway, and fatty acid metabolism, HIF-1 signaling pathway, PI3K-Akt signaling pathway, and AMPK signaling pathway ( Figure 2B ).

In order to investigate the hypoxia adaptation mechanism of HIF-1 signaling pathway in grouper, the DEGs of HIF-1 signaling pathway enriched in the HK tissue were identified. After hypoxia, a total of 26 DEGs were enriched compared with the Ctrl group, and the transcription levels were shown in Supplementary Figure S3 . There were 21 DEGs showed significant up-regulation and 5 DEGs showed significant down-regulation in the Hyp group.

DEGs with significantly up-regulated expression levels in the Hyp group included in glycolytic related genes (enolase 1 (ENO1), 6-phosphofructo-2-kinase/fructose-2,6-bisphosphatase 3 (PFKFB3), ATP-dependent 6-phosphofructokinase, liver type (PFKL), hexokinase-1 (HK1), LDHA, fructose-bisphosphate aldolase C-B (aldocb), ATP-dependent 6-phosphofructokinase, platelet type-like (PFKP), glyceraldehyde-3-phosphate dehydrogenase 2, cytosolic (GAPC)), and hypoxia responses (HIF1A, facilitated glucose transporter member 1 (SLC2A1), transferrin receptor 1b (TFRC)), cell proliferation and apoptosis (BCL2), and other unnamed genes. DEGs with markedly significantly down-regulated expression levels in the Hyp group included hexokinase HKDC1-like (Hkdc1), dual specificity mitogen-activated protein kinase kinase 2b (map2k2), eukaryotic translation initiation factor 4E family member 1c (EIF4E), cyclin-dependent kinase inhibitor 1 (CDKN1B), and other unnamed genes. These DEGs were mainly involved in human diseases, metabolism, organic systems, environmental information processing, and cellular processes.

To verify the accuracy of the RNA-seq data obtained, HIF-1α, TFRC, NFKBIA, HSP30, HSP70, IL21R, TLR9, ddit4, Egln1, CXCR1, CXCR2, CXCR4a, CXCR4b, CCR9a, CCR9b, ccnb1, ccnd1, egr1, CD274, and MRC1 were selected in Ctrl and Hyp groups for qRT-PCR validation ( Figure 3 ). log2FC was used for correlation analysis. Meanwhile, qRT-PCR validation was performed for BCL2, mknk2, PIK3R1, vegfaa, SLC2A1, PFKP, HK1, LDHA, Hkdc1, map2k2, EIF4E, and CDKN1B in the HIF1 pathway ( Figure 4 ). The qRT-PCR results confirmed the RNA-Seq data, showing consistent expression patterns.

To investigate the role of EcHIF-1α in E. cyanopodus under hypoxia stress, we cloned and identified HIF-1α. The EcHIF-1α gene, containing 2262 bp and encoding 753 aa, had a predicted MW of approximately 84.38 kDa and a PI of 5.11 ( Figure 5A ). The EcHIF-1α amino acid sequence contains key domains, including a HLH domain (residues 22–77), two PAS domains (residues 88–154 and 230–296), a PAC domain (residues 302–345), a HIF-1 domain (residues 551–582), and a HIF-1_CTAD domain (residues 716–752) via the SMART program ( Figures 5A, B ). The aa sequences of HIF-1α from speckled blue grouper and other species were compared. The results revealed that the gene has homology to other known and predicted HIF-1α ( Supplementary Figure S5 ). Moreover, the EcHIF-1α shared high sequence identity with giant grouper-HIF-1α (98.01%), zebrafish-HIF-1α (62.16%), african clawed frog-HIF-1α (51.05%), mouse-HIF-1α (51.05%), and human-HIF-1α (52.85%). To analyze the evolutionary relationships of EcHIF-1α, a phylogenetic tree was constructed using the NJ method ( Supplementary Figure S6 ). The tree demonstrated that EcHIF-1α is highly homologous to HIF-1α sequences from other fish species, with the closest relationship to that of the giant grouper (E. lanceolatus).

qRT-PCR was utilized to assess the tissue-specific expression pattern of EcHIF-1α in healthy grouper. As depicted in Figure 6 , EcHIF-1α was detected across all examined tissues. The highest expression of EcHIF-1α was observed in the HE, with subsequent notable expression in the TH, LI, PBL, and so on.

To investigate how bacterial infection affects EcHIF-1α expression in HK and LI, qRT-PCR was employed to investigate alterations with the V. harveyi, Poly(I:C) and LPS. In head kidney, the EcHIF-1α was expressed from 3 h to 72 h and reached the highest fold change at 6 h (25.53-fold) when stimulated with LPS ( Figure 7A ). Upon Poly(I:C) infection, EcHIF-1α exhibited sustained and significant up-regulation, peaking at 6 h with a 21.17-fold increase ( Figure 7B ). After V. harveyi challenged, it also showed significant up-regulation, reaching a maximum 16.05-fold increase at 12 h ( Figure 7C ).

In the liver, EcHIF-1α increased strongly after the challenge, although the increase was not as high as in HK. Challenged with LPS and V. harveyi, the EcHIF-1α expression was markedly up-regulated, peaking at 12 h with an 8.98-fold and 5.27-fold increase ( Figures 7D, E ), respectively. EcHIF-1α was prominently up-regulated after Poly (I:C) challenged and reached the peak at 6 h ( Figure 7F ).

To further investigate the effects of hypoxia and pathogens on EcHIF-1α expression, we cultured healthy speckled blue grouper in a bucket with low-xygen conditions 0.5 ± 0.1 mg/L) for 6 h, and then transferred it to a bucket with normal oxygen for challenge. As shown in Figure 8 , EcHIF-1α was in a differential pattern after treatment of different immune challenges.

In HK tissue, EcHIF-1α was prominently up-regulated at 3 h after LPS and Poly(I:C) stimulation (26.05-fold and 30.57-fold, respectively) ( Figures 8A, B ). Post stimulated with V. harveyi, EcHIF-1α was exhibited significant up-regulation, peaking at 6 h post-infection with a 25.77-fold increase ( Figure 8C ). In LI tissue, from 3 to 48 h after LPS infection, EcHIF-1α was significantly up-regulated and reached the peak at 3 h (14.86 -fold) ( Figure 8D ). The tendency of EcHIF-1α was similar to that after both Poly(I:C) and V. harveyi infections, EcHIF-1α exhibited significant up-regulation, reaching its peak at 6 hours p.i., which were 9.45-fold and 11.44-fold, respectively ( Figures 8E, F ).