SAETAN Wanida, YE Minghui, LIN Xinghua, LIN Xiaozhan, ZHANG Yulei,HUANG Yang, DU Tao, LI Guangli, and TIAN Changxu

Fisheries College, Guangdong Ocean University, Guangdong Research Center on Reproductive Control and Breeding Technology of Indigenous Valuable Fish Species, Guangdong Provincial Engineering Laboratory for Mariculture Organism Breeding, Southern Marine Science and Engineering Guangdong Laboratory, Zhanjiang 524088, China

Abstract Sillago sihama, commonly known as silver sillago, is considered as an economically important fish species in China. It is sensitive to hypoxia stress in the larval stage, and the mechanism has not been understood thoroughly. In this study, we investigated the transcriptome change in heart tissues under hypoxia stress. The fish were divided into four groups, including 1 h of hypoxia (hypoxia1h, dissolved oxygen (DO) = 1.5 ± 0.1 mg L-1), 4 h of hypoxia (hypoxia4h, DO = 1.5 ± 0.1 mg L-1), 4 h of reoxygen (reoxygen4h,DO = 8.0 ± 0.2 mg L-1) after 4 h of hypoxia (DO = 1.5 mg L-1) and normoxia or control (DO = 8.0 ± 0.2 mg L-1) groups. The results showed that a total of 3068 genes were identified as differentially expressed genes (DEGs) based on the criteria |log2 (Fold change)| ＞ 1.0 and adjusted P-value ＜ 0.05. A total of 7761141 and 1151 DEGs were obtained from hypoxia1h, hypoxia4h and reoxygen4h groups,respectively. The enrichment pathway analysis showed that the DEGs were significantly enriched in ribosome biogenesis in eukaryotes, retinol metabolism, DNA replication and the oxidative phosphorylation (OXPHOS) pathways. Thirteen DEGs from the RNAseq results were validated by quantitative real-time polymerase chain reaction (qRT-PCR). These candidate genes are considered as important regulatory factors involved in the hypoxia stress response in S. sihama.

Key words transcriptomes; heart tissues; hypoxia stress; Sillago sihama; gene expression

1 Introduction

Hypoxia is considered as poor solubility of dissolved oxygen (DO) (＜ 2.0 mg L-1), which frequently occurs in aquaculture (Hanet al., 2018). Due to global warming and eutrophication, hypoxia in the aquatic environment has become one of the most critical factors of aquaculture loss(Laiet al., 2016; Liet al., 2018). The hypoxia stress induces both chronic and acute stress responses in fish, which directly affects fish embryogenesis, immunology and growth physiology. For example, three-spine stickleback (Gasterosteus aculeatus) embryos showed delayed development and increased mortalities under hypoxia stress (Fitzgeraldet al., 2017). Hypoxia stress decreased the fertility of Gulf killifish (Fundulus grandis) and common carp (Cyprinus carpio) (Landryet al., 2007; Wuet al., 2017). The molecular mechanisms of fish response to hypoxia stress have become a research hotspot in recent years. The anaerobic and metabolism system in zebrafish (Danio rerio)were quickly changed after exposure to hypoxia (Martinovicet al., 2009). The expression of immune system-related genes was found to be down-regulated due to acute hypoxia in Nile tilapia (Oreochromis niloticus) (Choiet al.,2007). In addition, differentially expressed genes and their regulatory pathways related to hypoxia stress tolerance were observed in various fish species, such as Nile tilapia (Liet al., 2017), blunt snout bream (Megalobrama amblycephala) (Chenet al., 2017), eelpout (Zoarces viviparus)(Askeret al., 2013), Asian Seabass (Lates calcarifer) (Xiaet al., 2013) and fugu (Takifugu rubripes) (Jianget al.,2017).

Silver sillago (Sillago sihama) is one of the common species of Sillaginidae family (Tianet al., 2019). It is one of the main tropical shallow fish species, which is generally distributed in the Indian Ocean and along the coasts of China and Southeast Asia (Saetanet al., 2020). This fish species is nutritious and delicious, with a high economic value (Liet al., 2019). Due to overfishing, the population ofS. sihamahas decreased dramatically in recent years.Also, the poor hypoxic tolerance of this species limits the scale of artificial breeding ofS. sihama(Gunnet al., 1985).Up to date, a number of studies have been conducted inS.sihamaon the morphology (Tongnunuiet al., 2010), reproductive biology and artificial breeding (Yoshioka, 2000),population genetics (GUOet al., 2014; Liet al., 2019),and tissue physiology and ecology (Hakimelahiet al., 2012).We have conducted the transcriptome analysis ofS. sihamagill and liver tissues response to hypoxia stress, finding that the expression patterns of hypoxia-related genes were tissue-specific, and further study is necessary in more tissues (Saetanet al., 2020; Tianet al., 2020).

Transcriptome sequencing technology has been widely used for quantitative and qualitative analysis of transcripts in cells. In recent years, different candidate genes related to hypoxia and their signal transduction pathways have been identified based on RNA-seq technology in many fish species, including the blunt snout bream (Chenet al., 2017),Nile tilapia (Liet al., 2017), schizothoracine fish (Gymnocypris eckloni) (Qiet al., 2018) and crucian carp (Carassius auratus) (Liaoet al., 2013). In a previous study, it was observed the cytochrome P450 (CYP) and glutathione Stransferase (GST) gene families were widely expressed under hypoxia stress (Saetanet al., 2020). In teleosts, the heart is the organ in response to the changes of DO level in water, but the molecular mechanisms are still unclear.

In this study, transcriptome analysis was performed on heart tissue inS. sihama. Furthermore, quantitative realtime polymerase chain reaction (qRT-PCR) was used to verify the expression of selected genes. Our study will provide valuable information to understand the molecular mechanisms ofS. sihamaheart response to hypoxia stress.

2 Materials and Methods

2.1 Fish and Hypoxia Experimental Conditions

Healthy adultS. sihama(13.40 cm ± 1.05 cm of total length and 14.57 g ± 3.17 g of body weight) were obtained from Donghai Island, Guangdong, China. The fish were maintained in fiber tanks with a bio-filtered water circulation system at 25℃ for 1 month. The fish were fed with a commercial diet twice per day. The water quality was checked every day, and the dead animals and particles were removed at once. During the test period, the water temperature was maintained at 25℃ ± 1℃, dissolved oxygen (DO)at (8.0 ± 0.5) mg L-1(normoxia) and the salinity at 29.

The experiment methods were the same as described previously (Saetanet al., 2020). Healthy fish were randomly selected and transferred to four aquarium tanks (50-L) at a density of 50 fish per tank. Each tank contained 40 L of seawater. During hypoxia stress period, the concentration of DO was continuously monitored each hour by JPB-607A dissolved oxygen meter (INESA Scientific Instrument Co. Ltd., Shanghai, China). The fish were randomly divided into four groups, including hypoxia for 0 h (normoxia, DO = 8.0 ± 0.2 mg L-1), hypoxia for 1 h (hypoxia1h,DO = 1.5 ± 0.1 mg L-1), hypoxia for 4 h (hypoxia4h, DO =1.5 ± 0.1 mg L-1) and normal oxygen recovered in 4 h after hypoxia4h (reoxygen4h, DO = 8.0 ± 0.2 mg L-1). At each time point, fish heart sample was collected, immediately frozen in liquid nitrogen, and stored at -80℃ for further analysis.

2.2 RNA Extraction and Sequencing

The total RNA extraction and sequence preparing methods used in this study was described previously by Saetanet al. (2020). Total RNA of heart tissue (n= 3 per group)from four groups was extracted with TRIzol reagent (Life Technologies, Carlsbad, CA, USA) following the manufacturer’s instructions. Purified RNA samples were indicated by A260/A280 ratios ranging from 2.0 to 2.2 with NanoPhotometer spectrophotometer (Nanodrop 2000c,Thermo Scientific, Wilmington, DE, USA). RNA integrity samples were obtained by ethidium bromide staining of 28S and 18S ribosomal bands on a 1.0% agarose gel. The high-quality RNA samples were used to generated cDNA libraries using the NEBNext? Ultra? RNA Library Prep Kit for Illumina? (NEB, USA) following the manufacturer’s instructions. A total of 3 μg RNA was prepared for each Illumina library. The libraries were sequenced on the HiSeq platform with 150 bp sequenced from both ends(paired-end). The RNA-Seq data were uploaded to Sequence Read Archive database (SRA) (Accession no.: SRR 9651325-SRR9651336).

2.3 Data Analysis

Data filtering, reads mapping and differential expression analysis were conducted in accordance with the methods of Saetanet al.(2020). The assembledS. sihamagenome(Linet al., 2021) was used as a reference genome for mapping reads. The genome assembly included 521.63 Mb in 551 contigs with a contig N50 of 13559141 bp. An index of the reference genome was built, and the paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5. The gene expression levels were estimated by fragments per kilobase of exon model per million reads mapped (FPKM) (Trapnellet al., 2010). Clean data (clean reads) were picked out by removing reads containing adapter, poly-N and low-quality reads from raw data which were processed through in-house perl scripts. DESeq2 R package (version 1.16.1) was used to identify differentially expressed genes (DEGs) between the normoxia and hypoxia groups (Varetet al., 2016) with the threshold of|log2(Fold change)| ＞ 1.0 and Padj ＜ 0.05 (Anderset al.,2010). The DEGs were further mapped to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/) and Gene Ontology (GO) database(Padj ≤0.05).

2.4 qRT-PCR Validation

A total of 13 DEGs were randomly selected from hypoxia1hvs. normoxia, hypoxia4hvs. normoxia, and reoxygenation4hvs. normoxia to verify the expression of DEGs.Thirteen DEGs included 5 genes of co-expression in hypoxia group, 6 genes of hypoxia-related gene and 2 genes of top up- or down-regulated expressed (Table 1). The primers of all selected genes were designed using Primer Premier software v6.0 and listed in Table 1. qRT-PCR was performed using SYBR Green qPCR Mix (Dongsheng Biotech, Guangzhou, China) on a LightCycler real-time quantitative PCR system (Roche, USA) according to the manufacturer’s instructions. The ribosomal protein L7 (rpl7) gene was used as a reference to standardize gene expressionvalues (Zhanget al., 2018). All PCRs amplification were performed in triplicate. Relative expression levels were measured in terms of threshold cycle value and normalized using the 2-ΔΔCtmethod (Livaket al., 2001).

Table 1 PCR primer sequences

3 Results

3.1 Illumina Sequencing Assembly

Twelve cDNA libraries from four groups with triplicates were sequenced by Illumina technology to investigate the transcriptomes of heart tissues during hypoxia stress (Table 2). A total of 62932232 raw reads were obtained. A total of 159436408, 168800258, 143631922 and 136793080 clean reads were obtained from hypoxia1h, hypoxia4h, reoxygen4h and normoxia group, respectively after removing low-quality reads. All Q20 and Q30 values of the read sequences in the samples exceeded 96.22% and 90.72%,respectively.

3.2 Identification and Annotation of DEGs

In total 3068 DEGs were identified, of which 776, 1141 and 1151 DEGs were obtained from hypoxia1h, hypoxia4h and reoxygen4h groups, respectively. The number of significantly up-regulated genes in hypoxia1h, hypoxia4h and reoxygen4h groups were 387, 478 and 446, respectively. The number of down-regulated genes were 389, 663 and 705, respectively (P＜ 0.05) (Table 3). Further analy-sis showed that only 136 genes were expressed in different hypoxia groups as compared to the normoxia group(Fig.1). The top ten up- and down-regulated genes were presented in Table 4. Interestingly, the expression of heat shock protein 30 (Hsp30) and heat shock protein 70 (Hsp70)were shown to be strongly up-regulated under hypoxia stress. Additionally, cardiac muscle genes including troponin C (TnnC), troponin I (TnnI), tropomyosin (Tpm3) and myosin heavy chain (Myh) were significantly down-regulated.

Table 2 Summary of heart transcriptome sequencing data of S. sihama

Table 3 Summary of differentially expressed genes (DEGs)in S. sihama based on the criteria |log2 (Fold change)| ＞ 1.0 and Padj ＜ 0.05

Fig.1 Analysis of differentially expressed genes (DEGs).Venn diagram of corresponding significantly up-regulated or down-regulated genes in hypoxia1h, hypoxia4h and reoxygen4h groups compared to the normoxia group(log2 (Fold chang) ＞ 1.0 and Padj ＜ 0.05).

Table 4 Top 20 differentially expressed genes (DEGs) in the heart of S. Sihama

3.3 GO and KEGG Enrichment Analyses of DEGs

The DEGs were classified into three major functions, including biological process (BP), cellular component (CC)and molecular function (MF) according to GO enrichment analysis. The GO terms of DEGs in each group were shown in Table 5. Among significantly top 30 BP terms (P＜ 0.05),the BP terms were mostly enriched to ribosome biogenesis (GO:0042254) and DNA replication initiation (GO:000 6270). The majority of DEGs in the CC terms were related to the cytoskeleton (GO:0005856), actin cytoskeleton (GO:0015629) and cytoskeletal part (GO:0044430). In the MF terms, the DEGs were significantly enriched to heme binding (GO:0020037) and tetrapyrrole binding (GO:0046 906) (Table 5).

KEGG pathway analysis was annotated to obtain significantly enriched pathways. There were 3 KEGG pathways,including ribosome biogenesis in eukaryotes, retinol metabolism and DNA replication pathways, which were significantly enriched in the hypoxia1h group (Table 6). Inaddition, the oxidative phosphorylation pathway (Fig.2) was significantly enriched in the reoxygen4h group. There was no KEGG pathway significantly enriched in the hypoxia-4h group (Table 6).

Table 5 Gene Ontology (GO) enrichment of differentially expressed genes (DEGs)

Table 6 The significant enrichment of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of treatment groups

Fig.2 Effect of hypoxia on the oxidative phosphorylation pathway. The green frames represent the genes were up-regulated, while the red frames represent that the genes were down-regulated, respectively.

3.4 Validation of Gene Expression Levels

A total of 13 DEGs were selected and analyzed using qPCR (Fig.3). Our results demonstrated that the changing trends of those genes from q-PCR were similar to the results from RNA-seq expression analysis, which supported the reliability of the transcriptome data.

4 Discussion

Hypoxia is a common phenomenon that frequently occurred in aquatic environment. It severely affects various physiological functions in fish, such as metabolism and cardiovascular regulation (Abdel-Tawwabet al., 2019). The fishes have evolved various adaptation methods to hypoxia stress through a complex suite of molecular regulation (Qiet al., 2018). It was showed that the gill and liver tissues ofS. sihamaresponded to rapid changes of the DO level in water, while the expression patterns of hypoxia-related genes were tissue-specific (Saetanet al., 2020; Tianet al.,2020). Thus it is necessary to carry out research on hypoxia in more tissues. The heart is one of the major organs for fish to sense changes of DO level (Muet al., 2020). Hypoxia stress usually reduces the rate of pumping oxygenrich blood to various organs in the body (Nemtsaset al.,2010; Incardonaet al., 2016). Therefore, we conducted a comparative transcriptome analysis ofS. sihamaheart tissue under hypoxia stress to understand the molecular mechanisms of heart tissue’s response to hypoxia stress.

Fig.3 Comparison of gene expression data between RNA-seq and quantitative real-time PCR (qRT-PCR) after hypoxia acclimation compared to the normoxia. The x-axis presents the gene name and the y-axis presents fold change in gene expression. All data represent the mean value of three biological replicates. Error bars represent the standard errors of three replicates. Statistically significant differences from control are presented, with * P ＜ 0.05.

In the present study, numbers of down-regulated DEGs were increased with the increase of exposure time to hypoxia, which was in agreement with a previous study (Genget al., 2014). A series of down-regulated DEGs, such asTnnC,TnnI,Tpm3andMyh, are associated with cardiac muscle function, demonstrating that the fish heart responds to external environmental stress by down-regulating the gene expression levels related to energy metabolism.Tnnwas a complex of skeletal and cardiac muscle thin filaments, which consists of three subunits, includingTnnI,TnnTandTnnC.Tnnplays an important role in muscle activity and changing intracellular Ca2+concentration (Katrukhaet al., 2013).Tpm3is a component of thin muscle that is associated with cardiac muscle activation (Marston, 2008;Baiet al., 2013).Myhconverts the chemical energy of Adenosine triphosphate (ATP) to mechanical energy in eukaryotic cells (Lianget al., 2007). All of these genes are associated with stress response in cardiac muscle (Mizbaniet al., 2016; Stelzleet al., 2018), and play a pivotal role in predicting the expression of heart failure in patients. The heart ofS. sihamacan also respond to hypoxia stress through changes in the gene expression levels related to energy metabolism and oxygen consumption.

The detoxification proteins,heat shock protein 30(Hsp30)and heat shock protein 70(Hsp70) were significantly upregulated under hypoxia stress. Heat shock proteins (Hsps)are important to protect cells and prevent aggregation of proteins (Junprunget al., 2019). Currieet al.(2000) reported thatHsp30was induced by different environmental stressors, which was involved in the inhibition of apoptosis of cells.Hsp70is mainly involved in protecting the cells from extra stress to improve the cell survival. It is used as a bioindicator of cellular stress in animals (Zhouet al.,2019). In this experiment, the expression levels ofHsp30andHsp70were significantly up-regulated (P＜ 0.05)inS. sihamahearts under hypoxic stress, which was in line with the results reported in many fish species (Qian and Xue, 2016; Liuet al., 2019; Gaoet al., 2020). Our results suggested that Hsps 30/70 might play an essential role in protecting the heart tissues ofS. sihamafrom hypoxia stress.

In the present study, the functional classification of DEGs was carried out by GO enrichment and KEGG pathway analyses. As showed in Table 6, DEGs were mainly enriched in metabolic (e.g., oxidative phosphorylation and retinol metabolism) and genetic information processingrelated pathways (e.g., ribosome biogenesis in eukaryotes and DNA replication), which were also observed in schizothoracine fish (Qiet al., 2018), blue tilapia (Oreochromis aureus) (Nitzanet al., 2019) and blunt snout bream(Chenet al., 2017) under hypoxia stress, suggesting that these categories pathways may play an essential role under hypoxia stress inS. sihamaas well as other fishes.

ATP is the main source of energy within a cell. ATP is generated in mitochondria by the oxidative phosphorylation (OXPHOS) pathway (Wanget al., 2020). OXPHOS is a metabolic process, in which electrons produced by the citric acid cycle are transferred to the mitochondrial respiratory complexes (Silva-Marreroet al., 2017; Luoet al.,2019; Wanget al., 2020). This pathway is also involved in multiple cellular processes, such as calcium homeostasis, cyclic adenosine monophosphate (cAMP)/protein kinase A (PKA) signaling, inflammation, reactive oxygen species(ROS) production and apoptosis (Bergman and Ben-Shachar, 2016). Our study showed that most of the DEGs in OXPHOS pathway were down-regulated (Fig.2), suggesting that the OXPHOS pathway in heart tissue was suppressed by hypoxia stress. Fourteen DEGs related to Nicotinamide adenine dinucleotide (NADH) dehydrogenase were down-regulated in the reoxygen4h group, such as NADH dehydrogenase (ubiquinone) subunit 3 (Nduf3), Fe-S protein 1 (Ndufs1) and Fe-S protein 2 (Ndufs2). The decrease in NADH dehydrogenase expression can promote a stress-adaptive response in different aquatic animals under different stress conditions (Olsviket al., 2013; Chakrapaniet al., 2016; Mohapatraet al., 2018), demonstrating that the OXPHOS pathway plays a vital role in fish adapting to hypoxic environments.

Retinol metabolism plays an important role in cell signal transduction in embryonic development and adult physiology (Perlmann, 2002; Kamet al., 2012). Excessive retinol can lead to hypoxia and pathological endosteum mineralization in rats (Lindet al., 2011). In this study, several genes, such as retinol dehydrogenase 10 (Rdh10),retinal dehydrogenase (Aldh1a) and cytochrome P450 27C1(Cyp27c1) were up-regulated. In contrast, the retinol dehydrogenase 8 (Rdh8), aldehyde oxidase (Aox) and cytochrome P450 26B1 (Cyp26b1) genes were down-regulated.The cytochrome P450 (Cyp) gene was one of the environmental stress-induced genes, and hypoxia exposure influenced these genes in teleost fishes (Escobar-Camachoet al., 2019). The different expression patterns ofCypgene responding to stress were also reported inS. sihama(Saetanet al., 2020), zebrafish (Ben-Mosheet al., 2014) and Nile tilapia (Fenget al., 2015). The Aldh1a2protein, belonging to the aldehyde dehydrogenase (ALDH) family,can catalyze the synthesis of retinoic acid (RA) from retinaldehyde (Liet al., 2015). Hypoxia stress can induce the up-regulation ofAldh1a2gene expression, which has also been observed in multiple species, including zebrafish (D’Anielloet al., 2015), Nile tilapia (Fenget al., 2015)and rabbits (Oryctolagus cuniculus) (Jacksonet al., 2011).The genes related to retinol metabolism play an important role in response to hypoxic stress inS. sihama, and the mechanism remains to be studied.

Ribosomes are large ribonucleoproteins responsible for translating mRNA into protein complex in cells. The ribosome biogenesis played a crucial role in biological processes, such as cell growth and proliferation (Chaillouet al.,2014). Under the hypoxia1h group, many DEGs were significantly up-regulated in ribosome biogenesis of the eukaryotes pathway, including fibrillarin (Nop1), NOP58 ribonucleoprotein (Nop58) and small nuclear ribonucleoprotein (Snu13). These genes were involved in the box C/D snoRNA binding protein and responsible for rRNA modification process (Makimotoet al., 2006), demonstrating that the protein synthesis in heart tissues increased under hypoxia stress. Besides, DNA replication pathway was enriched with several up-regulated genes in the hypoxia1h group, such as mini-chromosome maintenance complexes 3, 4, 5 and 6 (Mcm3,Mcm4,Mcm5andMcm6). It had been reported that these genes families were the crucial components for the formation of the pre-replication complex(Wuet al., 2012). Consistent with the previous study, these gene families were differently expressed in response to stress in medaka fish (Oryzias latipes) (Chataniet al., 2016) and Nile tilapia (Kwoket al., 2015; Majerskaet al., 2018).The up-regulation of genes belonging to genetic information processing-related pathway, such as ribosome biogenesis in eukaryotes and DNA replication pathway, indicate that heart tissue needs sufficient energy for blood circulation under hypoxia stress.

5 Conclusions

In our study, heart transcriptome response to hypoxia stress was examined using RNA-Seq technology inS. sihama. A total of 3068 DEGs were identified, which represented the strongly down-regulated DEGs involved in the cardiac muscle function. Furthermore, the up-regulatedHsp30andHsp70genes were related to hypoxia stress. In addition, several DEGs were enriched in the OXPHOS pathway during hypoxia exposure. Our data revealed that candidate genes are important regulatory factors involved in the hypoxia stress response inS. sihama.

Acknowledgements

This study was supported by grants from the National Natural Science Foundation of China (Nos. 41706174 and 31702326), the Natural Science Foundation of Guangdong Province (No. 2019A1515110619), the Department of Education of Guangdong Province (Nos. 2018KQNCX111 and 2019KTSCX060), the College Students’ Innovation and Entrepreneurship Project of Guangdong Province (No. CX XL2019138) and the Program for Scientific Research Startup Funds of Guangdong Ocean University (No. R19026).