To identify HSP90 and HSP70 genes in L. maculatus, the reference genome (BioProject: PRJNA408177) and transcriptomic databases (SRR4409341, SRR4409397) were searched using TBLASTN based on the query sequences of HSP90 and HSP70 genes from human, mouse (Mus musculus), chicken (Gallus gallus), tropical clawed frog (Xenopus tropicalis), zebrafish, spotted gar (Lepisosteus oculatus), Mexican tetra (Astyanax mexicanus), tongue sole (Cynoglossus semilaevis), Atlantic salmon (Salmo salar), and channel catfish (Ictalurus punctatus) retrieved from the National Center for Biotechnology Information (NCBI)¹ and previous research (Song et al., 2014) (The query sequence identifiers are listed in Supplementary Table 1), with a cutoff E-value of 1e^––5. The open reading frames (ORFs) were predicted, and the retrieved sequences were translated by ORF Finder.² Predicted ORFs were validated by BLASTP against NCBI non-redundant protein database.

The gene copy numbers of HSP90 and HSP70 family were compared based on the genome databases of spotted seabass and several selected vertebrates, such as human, mouse, chicken, zebrafish, channel catfish, Japanese medaka (Oryzias latipes), fugu (Takifugu rubripes), spotted gar, large yellow croaker, and Nile tilapia (Oreochromis niloticus). The chromosomal location of each HSP gene was shown according to the coordinates of each HSP gene on the spotted seabass genome.

Phylogenetic analysis was conducted using the amino acid sequences of HSP90 and HSP70 genes from several representative vertebrates retrieved from NCBI and Ensembl³ databases, such as human, mouse, chicken, zebrafish, spotted gar, Atlantic salmon, channel catfish, Nile tilapia, large yellow croaker, fugu, and Japanese medaka (Supplementary Table 1). After conducting multiple alignments of HSP sequences by MUSCLE with default parameters, the phylogenetic tree was constructed using MEGA 7 with the neighbor-joining method. Jones-Taylor-Thornton (JTT) + invariant sites (I) + gamma distribution for modeling rate heterogeneity (G) model was selected and bootstrapping with 1,000 replications was conducted to evaluate the phylogenetic tree. The tree was further modified using Interactive Tree of Life 6 (iTOL 6⁴).

The syntenic analysis was performed to provide additional evidence for the annotation of the duplicated HSP90 and HSP70 genes. The neighboring genes of hsp90aa1, hspa8a, hyou1, hsp70.1, and hsp70.2 were identified from spotted seabass reference genome and further compared with the genome regions of selected representative vertebrates, which were obtained from NCBI and Genomicus databases 100.01⁵.

Protein characteristics (MW and isoelectric points, pIs) of HSP genes in spotted seabass were predicted by the Prot Param tool⁶ using the deduced amino acids. The homologous domain was surveyed by the SMART 7.0 program⁷ and the NCBI conserved domain database.⁸ The gene structures of HSPs were obtained from the L. maculatus reference genome database and were visualized using the GSDS 2.0 software.⁹ The conserved motifs of HSPs were observed with the MEME online tool¹⁰ and were displayed using TBtools 1.068.

For the prediction of heat shock elements (HSEs) and hypoxia response elements (HREs), 2-kb upstream regions from the transcription start site (TSS) of each HSP gene were extracted from the spotted seabass reference genome and analyzed using JASPAR 2018 server.¹¹

The Protein-Protein Interaction (PPI) relationships of the HSPs were predicted by constructing a zebrafish association model using STRING 11.0 software¹² with deduced amino acid sequences of spotted seabass.

The heat stress experiment was conducted at Shuangying Aquaculture Company, Dongying, Shandong, China. Sixty spotted seabass juveniles (body length: 13.33 ± 0.24 cm, body weight: 38.96 ± 2.01 g) were collected and acclimated for 2 weeks in a tank [5 m × 5 m × 1.5 m (L × W × H)]. During acclimation period, the experimental conditions were kept constant, such as water temperature (25.0 ± 1.0^oC), pH (7.5 ± 0.4), salinity (31 ± 1.0 ppt), and dissolved oxygen (DO) (7.0 ± 0.5 mg/L).

After acclimation, 60 individuals were transferred to three tanks at the density of 20 per tank. The temperature was settled as 25^oC, and fish were acclimated for 48 h. After that, the water temperature was increased at a constant rate of 1^oC/h until it reached 32°C, and thereafter, the temperature was maintained. Three individuals per tank were sampled at each time point, such as 0 h (H0), 10 h (3 h after heat stress, H3), 19 h (12 h after heat stress, H12), and 31 h (24 h after heat stress, H24) after heat stress (Supplementary Figure 1). Sampled fish were anesthetized with tricaine methanesulfonate (MS-222, 200 mg/L), and gill, liver, and muscle tissues were dissected and flash-frozen in liquid nitrogen for RNA extraction.

Total RNA of gill, liver, and muscle tissues in heat stress experiment was extracted using TRIzol reagent (Invitrogen, Waltham, CA, United States) according to the instructions of manufacturer. RNA concentration and integrity were measured by using the Biodropsis BD-1000 nucleic acid analyzer (Beijing Oriental Science & Technology Development Ltd., Beijing, China) and 1.5% agarose gel electrophoresis, respectively. cDNA was synthesized using PrimeScriptTM RT reagent Kit (Takara, Otsu, Japan) following the protocol of the manufacturer. qPCR was performed using the StepOne Plus Real-Time polymerase chain reaction (PCR) system (Applied Biosystems, Foster City, CA, United States) to detect the expression patterns of HSP90 and HSP70 mRNA. Primers of HSP90 and HSP70 genes were designed by Primer 6 software (Table 1). 18S rRNA was used as the internal control to correct the qPCR veracity, and the samples were run in triplicate (Wang et al., 2018). The 20 μl qPCR reaction volume included 2 μl template cDNA, 0.4 μl each forward/reverse primer, 10 μl 2 × ChamQ SYBR Color qPCR Master Mix, and 7.2 μl of nuclease-free water. The qPCR amplification was carried out using the following procedures: 95^oC for 30 s and 40 cycles of 95^oC for 5 s, 55^oC for 30 s, and 72^oC for 30 s. The relative mRNA expression levels of HSP90 and HSP70 genes were calculated according to the 2^–ΔΔCT method.

ANOVA and Duncan’s multiple tests were applied to assess the means of the relative mRNA expression level using SPSS 26.0 software. The differences were considered as statistically significant when the P-value < 0.05. The graphs were depicted by the software of GraphPad Prism.8.0.2.

To investigate the responsiveness of HSP90 and HSP70 genes to different abiotic stress conditions, RNA-Seq datasets generated by our previous challenge experiments of alkalinity and hypoxia were used to determine the gene expression pattern.

Briefly, for the alkalinity challenge experiment, 1-year-old spotted seabass individuals (body weight: 140.32 ± 2.56 g) were firstly acclimated in fresh water (pH: 7.8 ± 0.4) for 30 days. Carbonate-alkalinity solution was prepared by adding NaHCO₃ (12.8 mmol/L) and Na₂CO₃ (2.6 mmol/L) to fresh water and aerated for 24 h before the experiment. The carbonate alkalinity (mmol/L) was monitored every day during the exposure period using acidimetric titrations. After acclimation, 45 spotted seabass individuals were immediately transferred to three replicated 100 L square tanks with alkaline water (carbonate alkalinity: 18 ± 0.2 mmol/L), and the dissolved oxygen concentration, temperature, and pH were maintained at 7.1 ± 0.4 mg/L, 22 ± 1^oC and pH: 9.0 ± 0.2, respectively. The fish were not fed during the stress experiment. Gill tissues of three fish individuals in each tank were sampled at several time points, such as 0, 12, 24, and 72 h after alkalinity stress. The samples were frozen in liquid nitrogen and stored at − 80^oC for RNA extraction. RNA of the three individuals from the same tank was pooled as one sample, and 12 sequencing libraries (3 replicated samples × 4 time points) were generated. Total 150 bp paired-end reads (PRJNA611641) were obtained by the Illumina HiSeq X Ten platform.

For hypoxic treatment experiment, spotted seabass individuals (body weight: 178.25 ± 18.56 g) were randomly divided into two groups: normoxic group (6.89 ± 0.25mg/L) and hypoxic group (1.1 ± 0.14 mg O₂/L) in triplicate tanks at the density of 20 fish per tank. The threshold DO level was set as 1.1 ± 0.14 mg O₂/L, which was determined by preliminary experiments. The oxygen level of the hypoxia group was reduced to 1.1 ± 0.14 mg/L by bubbling nitrogen gas for 30 min. During the experiment period, the dissolved oxygen levels were maintained constant by a mixture of air and nitrogen gas. Gill tissues of three fish per tank (a total of nine individuals for each time point) were sampled at 0, 3, 6, and 12 h after the hypoxia challenge. A total of 12 sequencing libraries (3 replicated samples X 4 time points) were generated. Total 150 bp paired-end reads (unpublished data) were generated using the Illumina HiSeq 4000 platform.

For each RNA-seq project, the high-quality clean reads from each library were generated using Trimmomatic 0.36 and then mapped to the reference genome of spotted seabass (PRJNA408177) using HISAT 2 (-p 4 –dta -t –phred33) with no mismatch. The mapped reads from alignments were counted and then normalized to determine expected number of Fragment Per Kilobase of transcript sequence per Million base pairs sequenced (FPKM). Differential expression in fold change of each HSP gene was determined based on the ratio of the expression value. Differential expression statistical analysis was performed using the DESeq 2 R package with P-value < 0.05.

A total of 5 HSP90 (hsp90aa1.1, hsp90aa1.2, hsp90ab1, hsp90b1, and trap1) and 17 HSP70 (hspa8a.1, hspa8a.2, hspa8b, hsc70, hspa1b, hsp70.1, hsp70.2, hspa5, hspa9, hspa13, hyou1, hspa4a, hspa4b, hspa4l, hspa14, hspa12a, and hspa12b) genes were identified in spotted seabass. Their detailed information is summarized in Table 2. The ORF of HSP90 genes was ranged from 1,944 to 2,241 bp, encoding protein lengths different from 647 to 744 amino acid (aa). Their predicted MW was varied from 74.22 to 85.31 kDa, and pI were from 4.65 to 5.99. HSP70 genes contained ORFs varying from 1,329 to 2,982 bp in length and encoded proteins with 442 to 993 aa. MWs of HSP70 genes were distinct from 48.10 to 110.39 kDa, and pIs were from 4.95 to 7.22. All HSPs sequences of spotted seabass had been submitted to GenBank databases, and their accession numbers are listed in Table 2.

As shown in Figure 1A, HSP90 genes of spotted seabass were clustered with respective counterparts and four clades were generated, namely, HSP90AA1, HSP90AB1, HSP90B1, and TRAP1, which were consistent with their annotation. Two copies of hsp90aa1 (hsp90aa1.1 and hsp90aa1.2) in spotted seabass were orthologous to teleost hsp90aa1.1 and hsp90aa1.2, respectively.

All the HSP70 genes of spotted seabass were well distributed into distinct clades and grouped with homologous genes of selected species, supported by strong bootstrap values (Figure 1B). Eleven clades were generated, such as HSPA8, HSPA1, HSPA5, HSPA9, HSPA13, HYOU1, HSPA4, HSPA14, HSPA12, HSPA2/6/7, and HSPH1, and HSP70 genes of spotted seabass were distributed in nine of them (HSPA8, HSPA1, HSPA5, HSPA9, HSPA13, HYOU1, HSPA4, HSPA14, and HSPA12). HSPH1 genes were only discovered for selected tetrapods and a few teleosts, meanwhile hspa2/6/7 was specific to tetrapods. The annotation of five genes (hspa5, hspa9, hspa13, hyou1, and hspa14) with a single copy and duplicated genes, such as hspa12 (hspa12a and hspa12b) and hspa4 (hspa4a, hspa4b, and hspa4l), could be well supported by the phylogenetic relationships, but it is hard to distinguish and name the remaining HSP70 genes with duplicated copies solely by phylogeny.

Thus, syntenic analysis was performed to provide additional evidence for the annotation of three duplicated HSP genes in spotted seabass (Figure 2). For hsp90aa1, the hsp90aa1.1 and hsp90aa1.2 of zebrafish and spotted seabass were tandem located, with highly similar upstream neighbor genes like rps29 and mgat2 and downstream genes as ppp2r5cb, dio3b, slc25a47a and slc25a29. A unique gene copy of HSP90AA1 was found for human, and its neighbor genes, such as WDR20, PPP2R5C, and DIO3, were conserved with tested teleosts (Figure 2A). For all tested three species, hspa8/hspa8a and hyou1 were located at the same chromosome. Duplicated gene copies of hspa8a (hspa8a.1 and hspa8a.2) were identified for spotted seabass, alternatively, human and zebrafish harbored single hspa8 and hspa8a genes, respectively. Neighboring regions of hspa8a were more conserved between zebrafish and spotted seabass in comparison with human (Figure 2B). The synteny analyses of hsp70.1 and hsp70.2 were only conducted in zebrafish and spotted seabass due to the absence of hsp70.1 and hsp70.2 genes in human. As shown in Figure 2C, hsp70.1 and hsp70.2 were tandem duplicated in the zebrafish genome, while hsp70.1 and hsp70.2 in spotted seabass were separated by nlk2, cdipt, and TAO2. However, they shared a relatively conserved genomic neighborhood.

Taken together, the phylogenetic and syntenic analysis results not only well support the annotation and nomenclature of the HSP90 and HSP70 genes in spotted seabass but also indicate that these HSP gene families were highly conserved in the evolution.

For the chromosomal location of target HSP genes, 5 HSP90 genes were distributed on three chromosomes, and 17 HSP70 genes were located on 10 chromosomes and one un-anchored scaffold, respectively (Figure 3). In detail, for the HSP90 family, hsp90ab1 and tandem duplicated hsp90aa1.1 and hsp90aa1.2 were located on chr18, meanwhile, hsp90b1 and trap1 were located on chr6 and chr15, respectively. For the HSP70 family, chr12 possessed the largest gene number of HSP70s, that four genes, such as the tandem arranged hspa8a.2 and hspa8a.1, and hyou1 and hspa4a, were located on it. The tandem duplicated hsp70.1 and hsp70.2, and hsc70 were positioned on chr15. Additionally, two HSP70 genes, hspa13 and hspa8b, were located on chr5. For the rest nine HSP70 genes (hspa1b, hspa12a, hspa4b, hspa5, hspa4l, hspa14, hspa12b, and hspa9), they were distributed on the different chromosomes and scaffold fragment (Chr9, Chr13, Chr14, Chr17, Chr19, Chr22, Chr23, and scaffold_52), respectively.

Copy numbers of the HSP90 and HSP70 genes in representative vertebrates are summarized in Table 3. In general, the number of HSP90 genes was relatively conserved across a broad spectrum of vertebrate species from mammals to fishes, that the total gene numbers of HSP90 were slightly varied from 4 to 5 among different species. Only one copy of each gene was present except for that hsp90aa1. hsp90aa1 was duplicated in several detected teleosts, such as zebrafish, Japanese medaka, fugu, Nile tilapia, and spotted seabass, but a single copy was detected in tetrapods.

Gene copy number of the HSP70 family ranged from 13 to 19 among these detected species. In accordance with the phylogenetic analysis results, hspa6 and hspa7 genes were only discovered in human, and hspa2 was only present in higher vertebrates, such as human, mouse, and chicken. Four HSP70 genes were single-copy in all selected species, such as hspa9, hspa13, hspa14, and hyou1. Multiple gene copies were found in five HSP70 genes of teleosts, such as hspa1, hspa4, hspa5, hspa8, and hspa12. Among them, hspa1 harbored the highest copy numbers in most detected species, and hspa5 and hspa8 were only duplicated in teleost species (Table 3).

Exon-intron structures and motif analysis were employed to investigate the conservatism and diversity of gene structures in HSP90 and HSP70 genes (Figure 4). The intron numbers of the HSP90 genes varied from 8 to 18. Among them, hsp90aa1.1 and hsp90aa1.2 harbored highly similar exon-intron structures (Figure 4A). Contrastingly, the intron numbers of the HSP70 genes were remarkably varied from 0 to 22 (Figure 4B). According to the number of introns, HSP70 genes could be divided into two modes: mode 1 contained no intron, such as hspa1b, hsp70.1, and hsp70.2, and the rest genes belonged to mode 2, which possessed multiple introns, with intron numbers varied from 4 to 22. In general, the paralogous gene pairs shared similar exon-intron structures (Figure 4B).

For the motif analysis results, the motif numbers among HSP90 genes varied from 7 to 12, whereas all HSP90 members harbored seven conservative motifs, such as motif 1, motif 2, motif 4, motif 5, motif 6, motif 7, and motif 9 (Figure 4A). For the HSP70 family, most genes shared five conservative motifs, namely, motif 1, motif 2, motif 6, motif 10, and motif 11. Among them, hsc70, hspa1b, hsp70.1, hsp70.2, and hspa5 shared highly similar motif pattern, although belonging to two distinct exon-intron structure modes (with or without intron). Moreover, the paralogous gene groups, such as hspa4 genes (hspa4a, hspa4b, and hspa4l) and hspa12 genes (hspa12a and hspa12b), shared conserved motif pattern (Figure 4B). The length of each HSP90 and HSP70 motifs was varied from 21 to 50 and 29 to 50 amino acids, respectively.

As cis-regulatory elements, HSEs and HREs are the master regulators which bind to HSFs and hypoxia-inducible factors (HIFs) to mediate the transcriptional response of their target genes, respectively. Therefore, the distributions of these elements in the promoter regions of HSP90 and HSP70 genes were predicted for spotted seabass. As shown in Figure 5, 28 HSEs and 43 HREs were predicted in the selected promoter regions of HSPs genes. The HSE numbers were various among different HSP genes ranging from 0 to 3. In the 2 kb upstream region from the TSS, three HSEs were detected for hsp90ab1, hspa8a.1, hsp70.1, and hspa4a, two HSEs were identified for hsp90aa1.2, hsp90b1, hspa8b, hspa1b, and hspa12a genes, meanwhile only one HSE was existed in hsp90aa1.1, hsp70.2, hspa5, hspa13, hspa4l, and hspa12b (Figure 5).

For the prediction results of HREs, 0–6 HRE numbers were discovered in the 2 kb upstream of TSS of target HSP genes. Hyou1 harbored six HREs, meanwhile, hspa14 and hspa12b harbored five HREs. Four HREs were identified in hspa4b, three HREs were detected for trap1, hspa1b, hsp70.2, hspa13, and hspa4a. Hsc70 and hspa4l contained two HREs. Unique HRE position has existed in 2 kb upstream from the TSS of hsp90aa1.2, hspa8a.2, hspa8b, and hsp70.1 (Figure 5).

It can be deduced from the above result that six HSP genes (hsp90aa1.1, hspab1, hsp90b1, hspa5, hspa8a.1, and hspa12a) contained only HSE in their promoter regions, meanwhile six genes (trap1, hspa4b, hspa8a.2, hspa14, hsc70, and hyou1) only possessed HREs. Alternatively, nine genes, namely, hsp90aa1.2, hspa8b, hspa1b, hsp70.1, hsp70.2, hspa13, hspa4a, hspa4l, and hspa12b, harbored both HSE and HRE in their promoter regions (Figure 5).

As the transcription of HSPs is mainly regulated by heat shock factor 1 [HSF1 (Buckley and Hofmann, 2001], the PPI networks of HSPs and HSF1 in spotted seabass were constructed using STRING database based on the orthologs in zebrafish. The results indicated that 22 HSPs of spotted seabass were closely associated with 21 known zebrafish HSPs participated in the interaction network (Supplementary Figure 2 and Supplementary Table 2). As expected, extensive interactions were observed between the HSF1 and HSPs, which was revealed by a direct link representing “transcriptional regulation” in Supplementary Figure 2.

Additionally, HSP70s interacted with HSP70s, and HSP90s also be related to HSP90s, respectively. For example, HSPA4A, HSPA5, and HSPA9 showed a wide range of interactions (reaction, binding, and catalysis) with other HSP70 family proteins. Meanwhile, HSP90AA1.1, HSP90AA1.2, and HSP90AB1 proteins were connected by blue lines, suggesting their “binding” relationship (Supplementary Figure 2). Furthermore, strong interaction relationships were commonly found between HSP70 and HSP90 proteins, for instance, HSP90AA1.1, HSP90AA1.2, and HSP90AB1, were closely connected to HSPA4A, implying that the former may play a negative role in transcriptional regulation of the latter (Supplementary Figure 2).

The expression patterns of HSP90 and HSP70 genes in gill, liver, and muscle tissues of spotted seabass were systematically examined at 0, 3, 12, and 24 h after heat stress. Results showed that except for trap1 in the muscle, all HSP90 genes were significantly induced after heat stress in tissue- and time-dependent manner (Figure 6A). Among them, hsp90aa1.1 and hsp90aa1.2 exhibited the biggest expression changes in all three detected tissues, with their expression levels remarkably increased to at least 14.8-fold. The highest expression level of hsp90aa1.1 and hsp90aa1.2 appeared at 24 h in the muscle, with a 2,330-fold and 285-fold increase compared with 0 h, respectively. In gill, three genes (hsp90aa1.1, hsp90aa1.2, and hsp90b1) exhibited similar expression patterns after heat stress, which were significantly upregulated at 3 and 12 h and then reduced their expression after 24 h. The results were paralleled with the observation of their expression patterns in the liver (Figure 6A). In contrast, in muscle, the induced expressions of hsp90aa1.1, hsp90aa1.2, and hsp90ab1 genes were appeared for 3 h, maintaining the high expression levels until the end of the challenge experiment. A significant heat stress-responsive expressions were also detected for hsp90ab1 and trap1 in the gill and liver and hsp90b1 in the muscle, with their upregulated expression levels appeared since 12 h after heat stress treatment (Figure 6A).

Similarity, significant expression variations were commonly identified for several HSP70 genes in three target tissues after heat stress (Figure 6B). In gill tissue, excerpt for hspa8a.1, all HSP70 genes displayed significant differential expression patterns after heat treatment. Among them, highly heat stress-induced expressions in gill were discovered for hspa8b, hspa1b, hsp70.1, hsp70.2, and hspa4a genes, with the maximum expression values exceeded 10-fold compared with those under normal conditions (0 h). Their expression was dramatically induced at the initial tested phase of post-heat stress (3 h) and gradually scaled down during subsequent time points (Figure 6B). In the liver, remarkable upregulations were detected for all HSP70 genes with one exception that the expression change of hspa12b was not statistically significant (Figure 6B). The highest expression variations were discovered for hspa8b, hsc70, hspa1b, hsp70.1, hsp70.2, hspa4l, and hspa4a, which showed similar expression profile after heat stress treatment, that their expressions were dramatically induction at 3 h, remaining highly until 12 h, and then decreased to their normal expression level at 24 h (Figure 6B). The other heat stress-responsive HSP70 genes, namely, hspa8a.1, hspa8a.2, hspa9, hspa4b, and hspa14, exhibited dramatically induced expressions at 12 and 24 h. The expression of hyou1 was highly upregulated at all three tested time points after heat treatment (Figure 6B). For muscle tissues, the expressions of 13 HSP70 genes, such as hspa8a.1, hspa8a.2, hspa8b, hsc70, hspa1b, hsp70.1, hsp70.2, hspa5, hspa9, hyou1, hspa4a, hspa4b, and hspa4l, were significantly increased after heat stress. Among them, the upregulated expression levels of hspa1b, hsp70.1, hsp70.2, and hspa4l were extremely significant, exceeding hundreds of times than their normal expression values. Unlike the expression pattern in gill and liver, most differential expressed HSP70s, such as hspa8b, hspa1b, hsp70.1, hsp70.2, hspa5, hspa9, hspa4a, and hspa4l, displayed upregulation at all time points after heat treatment (Figure 6B). It was worth noting that hsp70.1 and hsp70.2 were the most intense responsive genes in all three tissues after heat stress treatment, suggesting the two genes may play significant roles in response to heat stress in spotted seabass.

As revealed by RNA-Seq analysis, after alkalinity stress, the HSP genes in gills of spotted seabass exhibited dynamic time-dependent expression pattern, and the number of differential expressed genes increased as the alkalinity treating time was prolonged (Figure 7A and Supplementary Table 3). A total of 11 genes, namely, hsp90aa1.2, hsp90b1, trap1, hsp70.2, hspa8a.1, hsc70, hspa13, hspa4a, hspa14, hspa5, and hyou1, showed downregulated expression profile (fold change: − 1.23 to − 3.06, P < 0.05). In detail, downregulated expressions of hspa8a.1 and hspa13 were found since 12 h after alkalinity treatment, a trend that continues until the final experimental time point (72 h). The other genes showed suppressed expression levels starting from the later time points (24 or 72 h). The highest expression changes of HSP genes during the alkalinity challenge experimental period were detected for hspa13, with fold change values nearly double in comparison with 0 h. Conversely, only two HSP genes (hspa12a and hspa12b) showed upregulated expressions after alkalinity stress (P < 0.05).

After the hypoxia stress, a time-dependent expression profile was also detected for HSP genes in gill, however, the expression trend was distinct from those under alkalinity stress. As shown in Figure 7B and Supplementary Table 4, a total of 12 genes (hsp90aa1.2, hspab1, hsp90b1, trap1, hsp70.1, hsp70.2, hspa8b, hspa13, hspa4a, hspa14, hspa12b, and hyou1) in the gill of spotted seabass were changed significantly after the hypoxia challenge (P < 0.05). In general, five genes, namely, hsp90ab1, trap1, hyou1, hspa13, and hspa14, exhibited downregulated trends (fold change: − 1.23 to − 1.47, P < 0.05), conversely, more genes, such as hsp90aa1.2, hsp90b1, hsp70.1, hsp70.2, hspa8b, hspa4a, and hspa12b, showed remarkably induced expression at 12 h during low-oxygen treatment experiment (fold change: 1.50–3.97, P < 0.05). Among the hypoxia-responsive expressed genes, hsp90aa1.2, hsp70.1, and hsp70.2 had more significant upregulated changes than other HSPs, with expression changes at 12 h exceeding threefold (Figure 7B and Supplementary Table 4).