Comparative transcriptomic analysis of the brain in Takifugu rubripes shows its tolerance to acute hypoxia

Hypoxia in water that caused by reduced levels of oxygen occurred frequently, due to the complex aquatic environment. Hypoxia tolerance for fish depends on a complete set of coping mechanisms such as oxygen perception and gene-protein interaction regulation. The present study examined the short-term effects of hypoxia on the brain in Takifugu rubripes. We sequenced the transcriptomes of the brain in T. rubripes to study their response mechanism to acute hypoxia. A total of 167 genes were differentially expressed in the brain of T. rubripes after exposed to acute hypoxia. Gene ontology and KEGG enrichment analysis indicated that hypoxia could cause metabolic and neurological changes, showing the clues of their adaptation to acute hypoxia. As the most complex and important organ, the brain of T. rubripes might be able to create a self-protection mechanism to resist or reduce damage caused by acute hypoxia stress.


Introduction
The fluctuation of living environment such as reduced air pressure and poor gas exchange in nature will cause the occurrence of a hypoxic environment, such as in plateaus (Qiu et al. 2012), aquatic environments (Jackson and Ultsch 2010), and underground tunnels (Avivi et al. 2005). The animals and plants living in a hypoxic environment have corresponding coping mechanisms in physiology, biochemistry, and behavior. Plants use substances other than oxygen as terminal electron acceptors, and the redox potential in the rhizosphere is reduced. Through a series of physiological and biochemical metabolic reactions, plants can adapt to or alleviate the damage caused by hypoxia stress (Schmidt et al. 2018). Animals have different ways of adapting to different levels of hypoxia: in mild or moderate hypoxic environments, they adjust to the hypoxic environment through the adjustment of the overall horizontal compensation Abstract Hypoxia in water that caused by reduced levels of oxygen occurred frequently, due to the complex aquatic environment. Hypoxia tolerance for fish depends on a complete set of coping mechanisms such as oxygen perception and gene-protein interaction regulation. The present study examined the short-term effects of hypoxia on the brain in Takifugu rubripes. We sequenced the transcriptomes of the brain in T. rubripes to study their response mechanism to acute hypoxia. A total of 167 genes were differentially expressed in the brain of T. rubripes after exposed to acute hypoxia. Gene ontology and KEGG enrichment analysis indicated that hypoxia could mechanism, such as faster breathing rates, increased pulmonary ventilation, alveolar-blood, and accelerating blood-tissue gas diffusion and strengthening the permanent expansion of capillaries, etc. (Avivi et al. 1999;Liu et al. 2001;Scott 2011;Wan et al. 2013). In severe hypoxic environments, adaptation through overall compensation alone is far from satisfying the body's energy needs, and more importantly, through the adjustment of cell metabolism and the induction of many anti-hypoxic factors adapt to the hypoxic environment at the molecular level. For example, intracellular free calcium increases during hypoxia, then it can activate calcium-related signaling pathways and regulate the increase in transcription of some hypoxia-sensitive genes (Millhorn et al. 1997). Hif-1 binds to hypoxia-responsive elements (HRE) on its target genes, triggering downstream target genes such as epo, inducible nitric oxide synthase (inos), vegf, and other genes. Transcription affects physiological and pathological processes such as erythropoiesis, angiogenesis, apoptosis, and proliferation, and causes the body to produce a series of adaptive responses to hypoxia (Bruzzi et al. 1997;Erkan et al. 2007). In addition to the hypoxia caused by the reduction of the oxygen content in the environment, many physiological and pathological processes of higher animals also have the phenomenon of hypoxia. For example, in the early stages of mammalian embryonic development, a hypoxic environment is created in the womb. This hypoxic microenvironment is essential for embryo development (Simon and Keith 2008).
Whereas aquatic species (such as fish) are particularly vulnerable to hypoxic conditions, the hypoxic state of aquatic systems usually means that the saturation of dissolved oxygen (DO) is between 0 and 30% (Pelster and Egg 2018;Wu et al. 2003). The DO concentration in fish ponds depends generally on many factors including photosynthesis of phytoplankton, respiration of aquatic organisms, and/or the diffusion of atmospheric O 2 , O 2 partial pressure, water temperature, and salinity. In captivity, fish always face repetitive and chronic stress situations (e.g., confinement, crowding, handling, variable water quality including hypoxia) from which they cannot escape. Across a broad range of fishes, hypoxia can cause direct mortality but more commonly results in various sublethal effects, such as behavioral and physiological stress (Abdel-Tawwab et al. 2019). Hypoxia has been shown to retard growth in bivalves, polychaetes, and fish (Shang and Wu 2004). It also can induce apoptosis (Arend et al. 2011). Fish have developed various adaptation strategies in the long-term evolution process, and their tolerance to hypoxia is very different (Bickler and Buck 2007). Studies have shown that fish could initiate special biological processes and molecular mechanisms in low-oxygen environments, and could more efficiently store and use oxygen through the regulation of key gene expression (Rimoldi et al. 2012). Salmo salar could induce vascular endothelial growth factor expression through hypoxia to promote angiogenesis and increase the body's oxygen supply (Vuori et al. 2004). In zebrafish (Danio rerio) embryos, gene expression changed in response to hypoxia. Hif-1a played an important role in regulating the formation of neural crest and the central nervous system (Ton et al. 2003).
Although there are many studies on molecular mechanisms related to hypoxia, the molecular basis of these adaptations has not been fully understood (Xia et al. 2018). The transcriptome analysis of Carassius auratus showed that the expression of glycolytic pathway-related genes in fish exposed to hypoxia for a long time will be significantly enhanced (Liao et al. 2013). The brain is an oxygen-sensitive organ whose functions are highly susceptible to low oxygen levels (Rahman and Thomas 2015). Also, other studies indicate that the C. auratus may be relying on neurotransmitters and neuromodulators to suppress its CNS energy use under hypoxic conditions (Nilsson and Renshaw 2004). Under hypoxic conditions, the expression of hypoxia-inducible factors in the brain tissue of sturgeon increased significantly (Pelster and Egg 2018). The T. rubripes genome sequencing work was completed in 2002. It was found that the T. rubripes genome is about 400 Mb, which is only one-seventh the size of the human genome. It is the smallest genome of known vertebrates and the number of genes is similar to humans (Brenner et al. 1993;Kai et al. 2011), so it has been extensively studied as a model organism. We have conducted a simple acute hypoxia tolerance experiment for T. rubripes in which 36 T. rubripes were randomly divided into two groups, one kept under hypoxic conditions and the other placed in water containing 6.7 ~ 7.0 mg/L DO as a control group. In the anoxic group, the containers were sealed to prevent air entry and to create hypoxic conditions, and the treatment time was 3 h. The criteria for DEGs in this experiment were adjusted P value < 0.05 and |log2FoldChange|> = 1; a total of 32 DEGs were identified. GO and KEGG enrichment analysis showed that when T. rubripes was exposed to hypoxia, the expression of transcription factors was significantly increased, metabolic processes were delayed and neurological development was disturbed (Jiang et al. 2017). In the present study, we set four O 2 concentrations to detect the transcriptome changes in the brain of T. rubripes and aimed to explore the mechanisms of the brain in T. rubripes responses to acute hypoxic stress.

Methods and materials
Experimental animals and acute hypoxia exposure Animal experiments were approved by the Animal Care and Use committee at Dalian Ocean University. Healthy T. rubripes with body weights of 461.75 ± 40.66 g were obtained from Dalian Tianzheng Industrial Co., in Dalian, Liaoning, China. Forty fish were randomly selected from the stock and divided equally into four groups kept in 380L tanks with a flow-through seawater supply. The temperature of the water was 19 ± 0.5˚C. Sodium sulfite was used to regulate the DO which was detected by dissolved oxygen meter (Hanna HI9146-04, Romania). The DO of the four tanks (four groups) were maintained at 5.4 ± 0.05 mg O 2 /L (ppm5.4), 4 ± 0.05 mg O 2 /L (ppm4), 2 ± 0.05 mg O 2 /L (ppm2), and 0 ± 0.05 mg O 2 /L (ppm0). The fish treated by the DO of 5.4 ± 0.05 mg O 2 /L (ppm5.4) was considered control because the DO was same as that in the stock. After 6 h of treatment, fish were anaesthetized with 80 mg/L concentration of tricaine methane sulfonate (MS-222, Sigma) and then dissected. Brain was obtained and snap-frozen in liquid nitrogen, and then stored at − 80℃. Then total RNA was extracted from the brain for the construction of RNA library and qPCR.
Transcript profiling and sequencing (RNAseq) Total RNA was quantified by the Qubit® RNA Analysis Kit in Qubit® 2.0 (Life Technologies, CA, USA), and its integrity was checked using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Twenty RNA (cDNA) libraries (5 biological replicates of each group, including ppm5.4, ppm4, ppm2, and ppm0) were constructed for RNA sequencing and bioinformatics analysis. The library was then sequenced by Novaseq 6000 (Illumina, USA).

Processing of RNAseq data
Raw data were firstly qualified by FastQC (www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc/). Clean data (clean reads) were obtained after removing read with containing adapter, poly-N and lowquality bases by Trimmomatic v0.38 (Bolger et al. 2014). All the downstream analyses were based on clean data. Reference genome (assembly fTakRub1.2, www. ncbi. nlm. nih. gov/ genome/ 63) was downloaded from the genome database of NCBI (The National Center for Biotechnology Information). Paired-end clean reads were mapped to the indexed reference genome using Hisat2 v2.1.0 (Kim et al. 2015). The results of mapping were also qualified by FastQC. HTSeq v0.11.2 was used to count the reads mapped to each gene (Anders et al. 2015). Differential gene expression analysis was performed by an R package DESeq2 v3.10 (Love et al. 2014). Principal component analysis (PCA) was performed using an in-house python script (Supplementary Information File S12). Obtained DEGs were annotated by KOBAS v3.0 (Xie et al. 2011), and subsequently subjected to GO functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis by an R package clusterProfiler v3.18.1 (Yu et al. 2012). The software Cytoscape v3.7.1 was used to map the enriched gene and KEGG signaling pathway network interactions (Shannon et al. 2003). Several R packages (ggplot2, pheatmap, venn) were used to draw graphics for the visualization of analyzed data (cran.rproject.org/).

Validation of DEGs by q RT-PCR
Total RNA was extracted from the brain samples using RNAprep pure tissue kit (Tiangen, China). First-strand cDNA was synthesized from 1 μg of total RNA using a PrimeScript™ RT Master Mix (Takara, Japan). qPCR was performed using Trans-Start® Top Green qPCR SuperMix (Transgen, China) on an ABI StepOnePlus™ Real-Time PCR System (Life Technologies, USA). The primers for each gene are listed in Supplementary Information Table S7. The gene expression levels were evaluated relative to the expression level of β-actin using the 2 −ΔΔCT method (Livak and Schmittgen 2001). The correlation between normalized counts from RNAseq and the relative expression from qPCR was calculated.

Sampling and sequencing
No death was found in all the experimental groups. But during hypoxia, especially lower DO, individuals were manifested by reddening and bleeding of the fins. Total RNA with high quality was extracted from brain for library construction. Total 20 libraries were constructed and sequenced from 20 samples in four treatment groups. The average size of insert fragments in these libraries was 275-315 bp (Supplementary Information Table S1). After quality control, total 102.75 Gb clean nucleotides and 685.06 million clean reads were obtained finally (Supplementary Information Table S2). All Q20 and Q30 values of the read sequences in the samples exceeded 95% and 89%, respectively. All clean data was submitted to SRA databases in NCBI (BioProject: PRJNA645780).

Mapping and counting reads
Total and unique mapping rates of all samples were more than 92% and 89% respectively (Supplementary  Information Table S3). According to the structure of the reference genome, mapped reads were mainly positioned in exons (77-81.5%) (Supplementary Information Table S4). Analysis of mapping results by FastQC showed that the level of duplication in mapped reads was about 15-21% from all samples (Supplementary Information Table S5). These results showed a good mapping quality for further analysis.

Analysis of differentially expressed genes (DEGs)
The counted reads by HTSeq were loaded by DESeq2 for differential expression analysis (Supplementary  Information Table S6). The results showed that 167 genes were differentially expressed among the different treated groups with adjusted P value < 0.05 (Supplementary Information Table S8). Of those 167 DEGs, only 19 genes had a difference greater than twofold (Fig. 1). In this experiment, we defined the criteria for DEGs as adjusted P value < 0.05. The comparison group with the highest number of DEGs was ppm0Vsppm5.4, and the groups with a wide range of DEGs fold change were ppm0Vsppm2, ppm-2Vsppm5.4 and ppm2Vsppm4, respectively (Fig. 2). The numbers of DEGs in different comparison groups were showed in (Fig. 3). Interestingly, the expression of vegfa and flt1, genes associated with angiogenesis, was upregulated, and the expression of cry1 and bhlhe40, genes associated with circadian rhythm, was downregulated. In contrast, hif1 (and hif1an) showed no significant changes.
Transformed count data (Supplementary Information Table S13) was extracted by DESeq2 using the method of variance stabilizing transformation in order to visualization, including clustering, heatmap and principal components analysis (Figs. 3, 4 and 5). All 20 samples were clustered into two main branches by the 167 DEGs (Figs. 4 and 5). The group ppm0 and ppm2 were clustered into one clade; ppm5.4 and ppm4 were clustered into another. These results indicated that the level of DO might make more serious effects on the brain of T. rubripes when the concentration of DO was less than 4 mg O 2 /L. In heatmap, 167 DEGs were mainly clustered into 3 clades based on their expression levels (Fig. 4). About a quarter of the DEGs exhibited high expression levels, suggesting that they could play important roles when T. rubripes was under the condition of low DO.

Annotation and function analysis of DEGs
The GO classification system grouped the DEGs into three main categories: biological process, cellular component, and molecular function (Fig. 6). The possible roles of DEGs were investigated by GO enrichment analysis. The enrichment analysis showed that most of the DEGs were significantly enriched to the category of biological processes (Fig. 6a). For example, the three terms ("response to steroid hormone," "response to hypoxia," and "response to decreased oxygen levels") were enriched significantly enriched with a high number of DEGs. The next significant enrichment was in the category of molecular function (Fig. 6b), and the top two enriched terms with a high number of DEGs were "dioxygenase activity" and "RNA polymerase II − specific DNA − binding transcription factor binding." However, for the category of cellular components, we found that the number of DEGs enriched was small and insignificant (Fig. 6c). Interestingly, we clustered all these enriched terms and found that they could be well clustered into eight major functional categories (Fig. 7). Among them, most of the significantly enriched terms were clustered into two functional categories, decreased oxygen levels hypoxia and corticosteroid corticosterone death glucocorticoid.
KEGG is a database resource containing many metabolic pathways and their relationships (Kanehisa and Goto 2000). In our study, 167 DEGs were grouped into 187 known pathways and the largest group was MAPK signaling pathway, containing 10 DEGs, followed by PI3K-Akt signaling pathway (9), Axon guidance (6), IL-17 signaling pathway (5), HIF-1 signaling pathway (5). Among them, the MAPK signaling pathway enrichment was most significant (Fig. 8). Interestingly, through the network diagram we found that MAPK signaling pathway and PI3K-Akt signaling pathway had the highest degree and connected the most genes (Fig. 9). For genes, fos, jun, rxra, vegfa, and flt1 had the highest degree and connected the most network pathways. Metabolism-related pathways and genes also clustered well together. In addition, five genes were found to be enriched in HIF-1 signaling pathways related to hypoxia stress: LOC101071669 (egln2), flt1, LOC101079462 (hk2), LOC101063282 (slc2a1), and vegfa. Among those HIF-1 signaling pathway-related genes, as the DO concentration decreased, all five enriched genes were upregulated, while hif1a did not change significantly.

Validation of DEGs by q RT-PCR
To validate our Illumina sequencing results, 6 genes were randomly selected for qRT-PCR analysis. As shown in Fig. 10, although the relative expression levels were not completely consistent, the expression patterns of these genes identified by qRT-PCR were similar to those obtained in RNA-Seq analysis, which indicated the expression data from RNAseq was reliable.

Discussion
Maintaining homeostasis is a key function of the brain, and homeostasis includes blood oxygen levels (van der Velpen et al. 2017). Hypoxia-inducible factor 1α (hif-1α) and vascular endothelial growth factor (vegf) are activated by hypoxia, which then leads to a disruption of the total blood-brain barrier (BBB),

Fig. 2
Boxplot and density of the foldchange of differentially expressed genes in each comparison group. The horizontal coordinate indicates the comparison groups between two different oxygen concentrations, and the vertical coordinate indicates the log 2 (foldchange). log 2 FoldChange > 0 indicates the genes were upregulated, and log 2 FoldChange < 0 indicates the genes were downregulated resulting in brain edema (Lafuente et al. 2016;Mohaddes et al. 2017). In mammals, hypoxia is a trigger stimulus for vascular remodeling, altered cell permeability, and angiogenesis (Zimna and Kurpisz 2015). Vegf has been proven to be an important determinant of angiogenesis and plays an important role in hypoxic-ischemic brain injury, functioning to promote angiogenesis and neuroprotection (Cao et al. 2018). Vegf could stimulate axonal growth and improve neuronal cell survival. Under pathological Fig. 4 Heatmap of differentially expressed genes. DEGs with adjusted P value < 0.05 were selected for heatmap. Heatmap was drawn based on the normalized expression data (Supplementary Information Table S13). 20 samples were marked under the heatmap conditions, vegf has a protective effect on the central nervous system (Plaschke et al. 2008). Vascular endothelial growth factor-α (vegfa) is a pro-angiogenic member of the vascular endothelial growth factor (vegf) family (Ferrara et al. 2003). Here, we showed that vegfa was upregulated with decreasing DO concentration. It participates in the AGE-RAGE signaling pathway. In addition, egr1, a gene of the  EGR family involved in the AGE-RAGE signaling pathway, showed a decreasing trend. Early growth response 1 (egr1) is considered a transcription factor sensitive to ischemia and hypoxia. In the brain tissue after ischemic stroke, the mRNA level of egr1 was significantly increased, while egr1 was overexpressed induced post-stroke inflammatory response and led to secondary brain damage after cerebral infarction (Tureyen et al. 2008). The main mechanism of the AGE-RAGE signaling pathway is that the interaction between AGE and RAGE causes an oxidative stress response, which in turn promotes the increase in diacylglycerol (DAG) synthesis and activates PKC. Activation of PKC increases the expression of vascular endothelial growth factor (vegf), transforming growth factor-β (tgf-β), endothelin-1, and prostaglandin, which proliferates the extracellular matrix, leading to vasomotor dysfunction and capillaries changes in permeability (Ishii et al. 1998;Toth et al. 2008). In our study, under hypoxic stress, vegfa expression was upregulated and egr1 was downregulated. It is shown that the brain of T. rubripes might create hypoxia tolerance to resist brain damage caused by hypoxia stress.
HIF-1 signaling pathway is a critical pathway during hypoxia. vegfa not only mediates the AGE-RAGE signaling pathway but also mediates the HIF-1 signaling pathway. Hypoxia-inducible factor 1 (hif1) is a transcription factor expressed in all metazoans and consists of hif-1alpha and hif-1beta subunits. Under hypoxic conditions, hif1 regulates the transcription of hundreds of genes in a celltype-specific manner and was a major regulator in the HIF-1 signaling pathway (Semenza 2007). In our study, genes enriched for the HIF-1 signaling pathway: egln2, flt1, sl2a1, vegfa, and hk2 were all upregulated. However, hif1 (and hif1an) did not show a significant change during hypoxia in Fig. 7 Functional grouping network diagram for GO enrichment analysis. The annotated GO terms were plotted in a network diagram, and the ellipses of different colors in the diagram represented clustered functional groups. Each point represented a GO term our results. In Eurasian perch, hif1a mRNA levels were upregulated after acute severe hypoxia exposure in the brain and liver, but not in muscle tissue, whereas significant changes were detected in muscle, but not in the brain and liver after chronic moderate hypoxia exposure (Rimoldi et al. 2012). Ndubuizu et al. demonstrated that despite the lack of hif1 activation, the relative mRNA levels of vegf in the aged cortex significantly increased (Ndubuizu et al. 2010). Fong et al. demonstrated that the knockout of hif1α in colon cancer cells had no effect on angiogenesis (Fong 2008). Liu et al. demonstrated that the large yellow croaker (Larimichthys crocea) has low hypoxia tolerance compared with other fish species, and the mRNA levels of hif1α in its brain did not change markedly under hypoxic conditions (Liu et al. 2018), which was consistent with our results. However, Shan et al. found that the upregulated expression of hif-1a and subsequent production of vegfa under hypoxic conditions play an important role in initiating pituitary tumor neovascularization (Shan et al. 2012). Hif1α is specific for the hypoxia response, and its degradation mediated by three enzymes egln1, egln2, and egln3 (Zhang et al. 2019). Tcf7l2 positively regulated aerobic glycolysis by suppressing Egl-9 family hypoxia inducible factor 2 (egln2), leading to the upregulation of hif1α (Xiang et al. 2018). Egln2 can hydroxylate foxo3a on two specific prolyl residues in vitro and in vivo. Hydroxylation of these sites prevents the binding of USP9x deubiquitinase, thereby promoting the proteasomal degradation of foxo3a (Zheng et al. 2014). Flt-1 (vegfr-1) and kdr (vegfr-2) were two highly homologous tyrosine kinase receptors of vegf, which can synergize with vegf to promote angiogenesis (Gille et al. 2000). Solute Carrier Family 2 Member 1 (slc2a1) was the most important energy carrier of the brain: present at the blood-brain barrier and assures the energyindependent, facilitative transport of glucose into the brain (Klepper et al. 1999). Sattiraju et al. found that the expression of the hypoxia-responsive gene slc2a1 was upregulated in hypoxic UnaG + GBM cells, which is consistent with our findings (Sattiraju Fig. 9 Network of enriched KEGG pathways and related genes. Red points represented the DEGs with upregulated trend, green points represented the DEGs with upregulated trend, and gray points represented enriched KEGG signaling pathway. Enriched metabolism-related pathways were shaded by a gray ellipse. Annotations of KEGG pathway was in Supplementary Information Table S9 et al. 2020). However, it has also been reported that slc2a1 was downregulated after hypoxia exposure and upregulated after reoxygenation (Wang et al. 2020). The hyperpolarization has been proposed to occur via stimulation of Na + /K + ATPase pumps caused by a glucokinase (gck) induced rise of ATP levels within neurons, leading to inhibition of neuronal activity (De Backer et al. 2016). In summary, hif1-mediated gene expression may be related to hypoxia-induced tolerance (Jones and Bergeron 2001), or it may be related to the different stresses of different tissues on hypoxia stimulation. Hif1 was inhibited in the T. rubripes brain to prevent brain damage and activate related genes of angiogenesis, promotes blood vessel growth, maintains normal blood vessel density, and was protected from damage caused by ischemia and hypoxia. At the same time, hypoxic stress also promoted the brain to reduce energy consumption.
In our research, we observed that when T. rubripes faced with hypoxia, they stopped swimming, laid on the bottom of the tank, and maintained their balance, with only slight swings in the fins. We speculated that the T. rubripes had a certain ability to tolerate hypoxia, and it can be made tolerant to hypoxia by changing its activity mode. The central nervous system of vertebrates controls the body's cognitive functions and autonomous motor activities (Paridaen and Fig. 10 Quantitative real-time PCR verification. The gene expression levels were evaluated relative to the expression level of β-actin using the 2 −ΔΔCT method. The correlation between normalized counts from RNAseq and the relative expression from qPCR was calculated Huttner 2014). Axon guidance represents a key stage in the formation of neuronal networks (Negishi et al. 2005). Axons were guided by a variety of guidance factors and these guidance cues are read by growth cone receptors, and signal transduction pathways downstream of these receptors converge onto the Rho GTPases to elicit changes in the cytoskeletal organization that determine which way the growth cone will turn (Govek et al. 2005). Recent work in the D. rerio has shown that developmental hypoxic injury disrupts pathfinding of forebrain neurons in D. rerio, leading to errors in which commissural axons fail to cross the midline. EphrinB2a acts as a ligand for one of the receptor tyrosine kinases (RTK) of the epha3, epha4, or ephb4 families, which in turn sets off an intracellular signaling cascade in the RTK-expressing cell (Stevenson et al. 2012). Hypoxia stimulated the uptake of 5-bromo-2′-deoxyuridine (BrdU) and reduced cell death Coincident with these proliferative changes, both hif1-α and phosphor (p)-AKT were increased while ephb3 expression was decreased (Baumann et al. 2013). These reports are consistent with our results. Our results showed that the expression of ephb3, ntng1 and rnd1 was downregulated with decreasing DO concentration, while epha4, sema5b and nck2 appeared to be upregulated. Among the downregulated genes, rnd1 was a small signal transduction G protein and a member of the rnd subgroup of the Rho family of GTPases (Ridley 2006). It contributes to the regulation of the actin cytoskeleton in response to extracellular growth factors (Nobes et al. 1998). Ntng1 serves as an axonal guidance cue during vertebrate nervous system development (Nakashiba et al. 2000). Among the upregulated genes, sema5b regulates the development and maintenance of synapse size and number in hippocampal neurons (O'Connor et al. 2009). Nck2 adaptor proteins were involved in signaling pathways mediating proliferation, cytoskeleton organization, and integrated stress response (Labelle-Cote et al. 2011). From this, we can speculate that the brain of T. rubripes may be affected by acute hypoxic stress. However, the brain responded promptly by inhibiting the expression of hif1, reducing the damage caused by hypoxic stress and repairing the damaged nerves in time. This may also be the reason why the brain of T. rubripes is able to tolerate hypoxia.
Hypoxia induces bhlhe40 expression independent of hif1α but through a novel p53-dependent signaling pathway, and inhibition of bhlhe40 or p53 may facilitate muscle regeneration after ischemic injuries (Wang et al. 2015). Studies have shown that prenatal hypoxia in mice can cause continuous changes in circadian rhythms in born mice (Joseph et al. 2002). Hypoxia significantly reduced clock, cry2, and per3 in GF and cry1, cry2, and per3 in PDLF (Janjic et al. 2017). Hif-2α increased the expression levels of clock, bmal1, per1, cry1, cry2, and ckiɛ, and decreased the expression levels of per2 and per3 (Yu et al. 2015). The negative circadian regulator cry1 was a negative regulator of hif1a (Dimova et al. 2019). Per3 was one of the primary components of circadian clock system. It was found to play a pivotal role in corticogenesis via regulation of excitatory neuron migration and synaptic network formation (Noda et al. 2019). In our study, the expression levels of both genes cry1 and bhlhe40 showed a decreasing trend due to acute hypoxic stress. Hypoxia has caused brain damage in T. rubripes to some extent, thus causing circadian rhythm disturbance in T. rubripes.
In conclusion, the brain of T. rubripes has a certain tolerance to hypoxia, but hypoxia also causes damage to it. Under hypoxic stress, the brain of T. rubripes struggles to maintain its central system from damage. Choose to maintain only basic life activities to resist the effects of hypoxic stress, such as stopping swimming to reduce oxygen consumption.
In addition, through the KEGG network interaction map we found that metabolism-related pathways and genes could be well clustered together, including fatty acid metabolism, phosphate metabolism, nitrogen metabolism, bile secretion, and steroid hormone synthesis. The results of Ding et al. on L. crocea showed that under acute hypoxic stress rhododendron stabilized energy supply by promoting glycogenolysis, glycolysis and gluconeogenesis and inhibiting gluconeogenesis (Ding et al. 2020). Cui et al. found that the pathways with the highest enrichment of DEGs in Callosobruchus chinensis under hypoxic conditions were carbohydrate, energy, lipid and amino acid metabolism, respectively (Cui et al. 2019). This remains consistent with our findings. Acl3 encodes a protein that is an isozyme of the long-chain fatty acid coenzyme A ligase family and plays a key role in lipid biosynthesis and fatty acid degradation. And this isozyme is highly expressed in the brain. Furthermore, acadl, the gene encoding LCAD (acyl coenzyme A dehydrogenase), has been shown to consume less energy in LCAD-deficient mice and also suffers from hypothermia, which can be explained by the fact that the reduced rate of fatty acid oxidation correlates with a reduced ability to generate heat (Diekman et al. 2014). From this, we can speculate that regulation of energy metabolism may be another effective way for T. rubripes to cope with hypoxic stress.

Conclusions
The purpose of this study was to explore the mechanisms of the brain in T. rubripes responses to acute hypoxic stress. The transcriptomes of their brain were sequenced and showed small changes under acute hypoxia. We also found that the circadian rhythm, neurodevelopment, and energy metabolism were affected to some extent under acute hypoxic stress. In addition, acute hypoxia also affected the HIF-1 signaling pathway and AGE-RAGE signaling pathway, probably promoting increased cerebral blood flow. Finally, our results indicated that the brain of T. rubripes was able to adapt to acute hypoxia and showed high tolerance to acute hypoxia.
Author contribution M. Bao performed the experiments and sampling, initially analyzed the results, and drafted the manuscript. F. Shang further analyzed the experimental results, made graphs, and revised the paper. M. Bao and F. Shang contributed equally to this work. F. Liu, Z. Hu, S. Wang, X. Yang, Y. Yu, H. Zhang, C. Jiang, and J. Jiang participated in the experiment and sampling. Y. Liu and X. Wang wrote and reviewed the manuscript. All authors reviewed and approved the final manuscript.
Funding This work was supported by the National Key R&D Program of China (2018YFD0900301-10) and the China Agriculture Research System (CARS-47).
Data availability All sequencing data were submitted to the Sequence Read Archive (SRA) public database in NCBI (www. ncbi. nlm. nih. gov/), under the accession code PRJNA645780. All other data included in this study are available upon request by contact with the corresponding author (Yang Liu).
Code availability Python code for PCA was included in Supplementary Information File S12. Other codes used in this study are available upon request by contact with the corresponding author (Yang Liu).

Declarations
Ethics approval and consent to participate Animal experiments were approved by the Animal Care and Use committee at Dalian Ocean University. All names in the author list have been involved in various stages of experimentation or writing.

Consent for publication
All authors agreed to submit the paper for publication in the Journal of Fish Physiology and Biochemistry.

Conflict of interest
The authors declare no competing interests.