Transcriptome Analysis Provides Insights Into the Adaptive Responses to Hypoxia of a Schizothoracine Fish (Gymnocypris eckloni)

The schizothoracine fish endemic to the Qinghai-Tibetan Plateau are comparatively well adapted to aquatic environments with low oxygen partial pressures. However, few studies have used transcriptomic profiling to investigate the adaptive responses of schizothoracine fish tissues to hypoxic stress. This study compared the transcriptomes of Gymnocypris eckloni subjected to 72 h of hypoxia (Dissolved oxygen, DO = 3.0 ± 0.1 mg/L) to those of G. eckloni under normoxia (DO = 8.4 ± 0.1 mg/L). To identify the potential genes and pathways activated in response to hypoxic stress, we collected muscle, liver, brain, heart, and blood samples from normoxic and hypoxic fish for RNA-Seq analysis. We annotated 337,481 gene fragments. Of these, 462 were differentially expressed in the hypoxic fish as compared to the normoxic fish. Under hypoxia, the transcriptomic profiles of the tissues differed, with muscle the most strongly affected by hypoxia. Our data indicated that G. eckloni underwent adaptive changes in gene expression in response to hypoxia. Several strategies used by G. eckloni to cope with hypoxia were similar to those used by other fish, including a switch from aerobic oxidation to anaerobic glycolysis and the suppression of major energy-requiring processes. However, G. eckloni used an additional distinct strategy to survive hypoxic environments: a strengthening of the antioxidant system and minimization of ischemic injury. Here, we identified several pathways and related genes involved in the hypoxic response of the schizothoracine fish. This study provides insights into the mechanisms used by schizothoracine fish to adapt to hypoxic environments.

The schizothoracine fish endemic to the Qinghai-Tibetan Plateau are comparatively well adapted to aquatic environments with low oxygen partial pressures. However, few studies have used transcriptomic profiling to investigate the adaptive responses of schizothoracine fish tissues to hypoxic stress. This study compared the transcriptomes of Gymnocypris eckloni subjected to 72 h of hypoxia (Dissolved oxygen, DO = 3.0 ± 0.1 mg/L) to those of G. eckloni under normoxia (DO = 8.4 ± 0.1 mg/L). To identify the potential genes and pathways activated in response to hypoxic stress, we collected muscle, liver, brain, heart, and blood samples from normoxic and hypoxic fish for RNA-Seq analysis. We annotated 337,481 gene fragments. Of these, 462 were differentially expressed in the hypoxic fish as compared to the normoxic fish. Under hypoxia, the transcriptomic profiles of the tissues differed, with muscle the most strongly affected by hypoxia. Our data indicated that G. eckloni underwent adaptive changes in gene expression in response to hypoxia. Several strategies used by G. eckloni to cope with hypoxia were similar to those used by other fish, including a switch from aerobic oxidation to anaerobic glycolysis and the suppression of major energy-requiring processes. However, G. eckloni used an additional distinct strategy to survive hypoxic environments: a strengthening of the antioxidant system and minimization of ischemic injury. Here, we identified several pathways and related genes involved in the hypoxic response of the schizothoracine fish. This study provides insights into the mechanisms used by schizothoracine fish to adapt to hypoxic environments.
Schizothoracine fish (Teleostei: Cyprinidae) are the largest and most diverse taxon within the Qinghai-Tibetan Plateau (QTP) ichthyofauna (Wu and Wu, 1992;Chen and Cao, 2000). The schizothoracine fish are divided into three groups based on morphological characters: primitive, specialized, and highly specialized (Cao et al., 1981). Previous studies have demonstrated that schizothoracine fish restricted to high-altitude or high-latitude environments are well adapted to cold, hypoxic conditions (Qi et al., 2012;Li et al., 2013;Guan et al., 2014;Xia et al., 2016). Therefore, the schizothoracine fish are an excellent model for studies of the fundamental mechanisms of adaptions to hypoxia and low temperatures. Gymnocypris eckloni, is a representative species of the highly specialized schizothoracine fish, and this species is widely distributed in isolated lakes and the upper reaches of the Yellow River (Wu and Wu, 1992;Chen and Cao, 2000). This species is very well adapted to the dissolved oxygen fluctuations of its aqueous environment (Qi et al., 2012;Xia et al., 2016). During artificial breeding experiments, G. eckloni displayed a substantial tolerance for hypoxia, surviving acute hypoxia (∼0.3 mg O 2 /L) for 12 h and chronic hypoxia (∼3.0 mg O 2 /L) for much longer (Xia et al., 2016). Therefore, we hypothesized that G. eckloni would exhibit transcriptional responses to hypoxia that would differ from less hypoxia-tolerant fish species.
Understanding changes in the gene expression in fish exposed to hypoxia might reveal new mechanisms of hypoxic acclimation, and might shed light on the evolution of this adaptive response in vertebrates. In the present study, we sequenced the transcriptomes of five tissues (the liver, muscle, brain, heart, and blood) from G. eckloni under normoxic or hypoxic conditions. We used this transcriptomic data to examine the gene expression profiles of G. eckloni exposed to hypoxia. Our results increase our understanding of the adaptive mechanisms used by the schizothoracine fish in response to hypoxia.

Experimental Fish and Hypoxia Treatment
Healthy G. eckloni (average body weight: ∼170 g) were obtained from the Native Fish Artificial Proliferation and Release Station, Xunhua, Qinghai Province, China. Fish were transported to our laboratory in aerated buckets. We used six adult individuals for our experiments. All of the fish were allowed to acclimate for 3 days before experimentation, in 10-13 • C freshwater almost identical in composition to the river water at the capture site. All of the research involving animals followed the guidelines of, and was conducted under the approval of, the Animal Care and Use Committee, Qinghai University, China.
Three fish were randomly assigned to the control group (normoxia); the remaining three fish were assigned to the hypoxia group. The two groups were kept in separate, identical 20 L translucent plastic tanks. The freshwater in the tanks was bubbled with nitrogen and with the ambient air in our laboratory. The dissolved oxygen (DO) level was continuously monitored using an AZ8402 dissolved oxygen meter (AZ Instrument Corp., Taiwan) (Xia et al., 2016). In the normoxia (control) tank, DO was 8.4 ± 0.1 mg/L. In the hypoxia tank, DO was reduced from 8.4 ± 0.1 to 3.0 ± 0.1 mg/L over 1 h, and then maintained at 3.0 ± 0.1 mg/L for 72 h.
After the experiment, fish were stunned with a blow to the head and blood was quickly collected from a caudal puncture in an EDTA-treated tube, then transferred to a microfuge tube. The liver, brain, and heart were removed from each fish, and a sample was taken of the muscle tissue. Blood and tissue samples were stored in liquid nitrogen until RNA isolation.

RNA Extraction, Illumina Library Preparation, and Sequencing
An Ambion MagMAX-96 total RNA isolation kit (Life Sciences, United States) was used to isolate total RNA from each sample, following the manufacturer's instructions. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was analyzed with a NanoPhotometer Spectrophotometer (Implen, United States). The RNA concentration was quantified with a Qubit RNA Assay Kit in a Qubit 2.0 Fluorometer (Life Technologies, United States). RNA integrity was analyzed using a RNA Nano 6000 Assay Kit and an Agilent Bioanalyzer 2100 (Agilent Technologies, United States). Thirty separate Illumina sequencing libraries were prepared using 3.0 µg RNA per sample.
Sequencing libraries were created with the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, United States), following the manufacturer's instructions. Index codes were added to assign sequences for each sample. In brief, mRNA was purified from the total RNA with poly-T oligo-attached magnetic beads. Fragmentation was done with divalent cations under an elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (NEB, United States). First strand cDNA was synthesized with M-MuLV Reverse Transcriptase (RNase H − ) and random hexamer primers. Second strand cDNA synthesis was subsequently carried out with RNase H and DNA Polymerase I. Exonuclease/polymerase activity converted the remaining overhangs into blunt ends. After adenylation of the 3 ends of the DNA fragments, NEBNext Adaptors with hairpin loop structures were ligated to prepare for hybridization. Library fragments were purified with the AMPure XP system (Beckman Coulter, United States) to preferentially select 150-200 bp cDNA fragments. The adaptor-ligated, size-selected cDNA was then incubated with 3.0 µl USER Enzyme (NEB, United States) for 15 min at 37 • C followed by 5 min at 95 • C before PCR. PCRs were performed using universal PCR primers, Phusion High-Fidelity DNA polymerase, and the Index (X) Primer. Finally, PCR products were purified (with the AMPure XP system). Library quality was assessed with an Agilent Bioanalyzer 2100.
Clustering of the index-coded samples was carried out on a cBot Cluster Generation System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina), following the manufacturer's instructions. Library preparations were sequenced on an Illumina HiSeq 2500 platform after cluster generation, which generated 125 bp paired-end reads. All data sets have been submitted to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database (accession number SRP150490).

Quality Control and de novo Assembly
Raw data (raw reads) in FASTQ format were sorted by individual fish and processed using self-written Perl scripts. Our scripts generated cleaned data (cleaned reads) by removing reads containing adapter sequences, reads containing poly-N, and low quality reads. We also calculated Q20, Q30, GC-content, and sequence duplication level for the cleaned data. All of the downstream analyses were based on high-quality cleaned data.
De novo transcriptome assembly was performed using the short read assembly program Trinity, with min_kmer_cov set to 2 by default and all of the other parameters set default (Grabherr et al., 2011). Briefly, reads of a certain length with overlapping areas were joined initially in order to form longer fragments (contigs) without gaps. The paired-end reads were then mapped back onto the contigs, which were subsequently connected. These steps were repeated until the sequences could no longer be extended (Stewart et al., 2015). The final sequences were considered to be unigenes, which still were, strictly speaking, gene fragments.

Transcriptome Annotation
We annotated the gene fragments against the NCBI nonredundant (Nr) protein database, the Swiss-Prot database, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa et al., 2004), the Clusters of Orthologous Groups of proteins (COG) database, the Pfam database, and the nucleotide (Nt) database using BLASTX, with an E-value 1.0 × 10 −5 . All of the annotated gene fragments were translated into amino acid sequences using ESTScan (Iseli et al., 1999), and then screened against the Gene Ontology (GO) database using Blast2GO (Conesa et al., 2005). The clusters of orthologous groups (COG) database was used to identify the putative functions of each gene fragment product based on known orthologous gene products (Tatusov et al., 2003). KEGG pathways were identified using the KEGG Automatic Annotation Sever (KAAS 1 ) with the bi-directional best-hit method and an E-value of 1.0 × 10 −10 .

Differential Expression Analysis and Functional Enrichment
To test for differentially expressed genes (DEGs), individual sequence reads from each sample were mapped back to the assembled transcriptome using the alignment program Bowtie (Langmead et al., 2009). Differential expression analysis was performed using the DESeq R package (Anders and Huber, 2010). This package includes statistical routines that determine differential expressions in digital gene expression data based on a model of a negative binomial distribution. Mean expression levels were generated using three tissue samples per treatmen. The resulting P-values were adjusted using the Benjamini-Hochberg approach for controlling the false discovery rate. Genes identified by DEGSeq with an adjusted P-value < 0.05 were considered DEGs. GO and KEGG pathway analyses were then performed on the DEGs.

Quantitative Real-Time PCR Validation (qRT-PCR)
RNA samples from three biological replicates were used for qRT-PCR validation of the transcriptome data. RNA used for transcriptomic sequencing was converted into cDNA using the RNAprep Pure Tissue Kit (TIANGEN Biotech Co., Ltd., China), and treated with RNase-free DNase I (Takara, Japan). Gene-specific primers were designed using Primer Premier 5.0 (Supplementary Table 1). qRT-PCR was performed in an iQ5 Multicolor Real-Time PCR Detection System (Bio-Rad Laboratories, Inc., United States) using SuperReal PreMix Plus (SYBR Green) (TIANGEN Biotech Co., Ltd.). All of the samples were analyzed in triplicate (as technical replicates) and fold changes in gene expression were calculated using the 2 − Ct method (Livak and Schmittgen, 2001). GAPDH served as the internal control.

Raw Sequencing Data and de novo Transcriptome Assembly
Here, 30 separate Illumina sequencing libraries produced close to 1,837.7 million paired-end reads (2 × 125 bp). After trimming and quality filtering, 4.5% of all of the reads were discarded, leaving over 1,753.7 million reads (219.12 Gb) for downstream analysis. The number of qualified Illumina reads per sample ranged from 50.9 million (6.36 Gb) to 69.8 million (8.72 Gb), with a mean of 58.5 ± 4.7 million. The number of Illumina reads in each treatment group was well balanced, with 871.97 million reads (109.04 Gb) in the control group and 881.70 million reads (110.08 Gb) in the hypoxia group (Table 1).

Functional Annotation
To annotate the G. eckloni transcriptome, all 337,481 gene fragments were annotated against six public databases (Nr protein, Nt, Swiss-Prot, KEGG, GO, and COG) using BLASTX and BLASTN, both with an E-value cut-off of <10 −5 . We annotated 46,418 gene fragments against the Nr protein database (13.75%; Table 2). The E-value distribution of the top hits in the Nr database ranged from 0 to 1.0E −15 (Supplementary Figure 1A). The similarity distribution of the top BLAST hits for each sequence ranged from 18 to 100% (Supplementary Figure 1B). The most common species with similar sequences in the Nr database, were Danio rerio (62.1%), followed by Astyanax mexicanus (4.7%), Oncorhynchus mykiss (3.9%), Oreochromis niloticus (1.8%), and Oryzias latipes (1.5%) (Supplementary Figure 1C). The previously published genome annotations for these species are comprehensive and largely accepted, suggesting that the G. eckloni gene fragments were correctly assembled and annotated.
Gene Ontology (GO) assignments were used to classify the functions of the predicted G. eckloni gene fragments. Based on sequence homology, 50,492 sequences were categorized into 51 functional groups in three functional divisions (Figure 1). Many of the gene fragments were classified as cellular processes, metabolic processes, cell, cell part, single-organism processes, binding, and catalytic activity.
To test the effectiveness of the annotation process and evaluate the integrity of the transcriptome libraries, gene fragments were classified based on KOG categories. Of the 46,418 Nr hits, 21,106 sequences were assigned to 26 KOG categories (Figure 2). Of these, "General function prediction only" was assigned the most sequences (4,956), followed by "Signal transduction mechanisms" (4,412 sequences), "Posttranslational modification, protein turnover, chaperones" (1,917 sequences), "Transcription" (1,303 sequences), "Intracellular trafficking, secretion, and vesicular transport" (1,211 sequences), and "Translation, ribosomal structure and biogenesis" (929 sequences). The KOG categories "Cell motility" and "Unnamed protein" had the fewest sequences (64 and 3 sequences, respectively).
The 46,418 annotated sequences were mapped to the canonical reference pathways in the KEGG database. Of these, 26,660 gene fragments were successfully assigned to 278 KEGG pathways in five categories (Figure 3). Most gene fragments fell into the KEGG pathway category organismal systems (8,025 sequences), followed by metabolic pathways (5,216 sequences), environmental information processing (4,868 sequences), cellular   processes (3,548 sequences), and genetic information processing (2,555 sequences).

DEGs
Differential expression analysis of gene fragments indicated that 462 sequences were differentially expressed in fish exposed to hypoxia. Of these, 28 were differentially expressed in liver (20 upregulated and 8 downregulated), 30 were differentially expressed in the heart (18 upregulated and 12 downregulated), 74 were differentially expressed in the blood (59 upregulated and 15 downregulated), 60 were differentially expressed in the brain (48 upregulated and 12 downregulated), and 282 were differentially expressed in the muscle (247 upregulated and 35 downregulated). To validate our Illumina sequencing results, 21 DEGs from liver, heart, brain, and muscle were selected for qRT-PCR analysis. The expression patterns of 20 of these genes were consistent with those indicated by our RNA-Seq analysis. The remaining DEG, polyprotein (Gene ID: c227508_g1), was upregulated according to the RNA-Seq analysis, but downregulated according to the qRT-PCR analysis (Supplementary Figure 2). GO functional enrichment analysis using GO slim terms showed that only the DEGs in the muscle, liver, and brain were significantly enriched in GO terms (Supplementary Table 2). In the muscle, 231 of the 247 upregulated DEGs were significantly enriched in 22 biological process categories, one cellular component category, and six molecular function categories (Figure 4A), while 20 of the 35 downregulated DEGs were significantly enriched in five biological process categories, two cellular component categories, and three molecular function categories (Figure 4B). In the liver, the 20 upregulated DEGs were significantly enriched in five molecular function categories ( Figure 4C). In the brain, 18 of the 48 upregulated DEGs were significantly enriched in four molecular function categories ( Figure 4D).
Only DEGs in the muscle, liver, brain, and blood were enriched in KEGG pathways. For example, 60 significantly enriched KEGG pathways, including 119 DEGs, were identified in the muscle (see Figure 5A for the top 20 most enriched pathways; other enriched pathways are shown in Supplementary  Table 3). The ascorbate and aldarate metabolism pathway, which is associated with antioxidant defense, was the most enriched pathway, followed by pentose and glucuronate interconversion, the pentose phosphate pathway, and galactose metabolism. In the blood, 11 DEGs were enriched in nine KEGG pathways ( Figure 5B). In the brain, four DEGs were enriched in four KEGG pathways ( Figure 5C). In the liver, seven DEGs were enriched in seven KEGG pathways ( Figure 5D).

DISCUSSION
Limited oxygen availability (hypoxia) is an important environmental stressor that has a profound impact on the survival, reproduction, growth, and development of aquatic life, and thus the aquatic ecosystem as a whole (Shang and Wu, 2004;Storz et al., 2010). Many fish species are able to cope with or adapt to hypoxia via a complex suite of molecular strategies that operate at the transcriptional level (Gracey et al., 2001;van der Meer et al., 2005;Zhong et al., 2009;Kang et al., 2017;Liu et al., 2018). Here, we describe the transcriptomic response of five G. eckloni tissues to hypoxia, affording new insights into the adaptive mechanisms employed by fish subjected to hypoxic conditions. Under hypoxia, fish may either reduce their metabolic rate to match reduced energy supply, or may maintain their metabolic rate, compensating for the ATP deficit with anaerobic metabolism (glycolysis) (Val et al., 1995;Gracey et al., 2001;Nikinmaa, 2002;Ton et al., 2003;Zhong et al., 2009). Here, several GO slim terms related to the aerobic metabolism of sugar, fat, and protein (including acetyl-CoA transport, fatty-acyl-CoA, coenzyme transport, acetyl-CoA transporter activity, and coenzyme transporter activity) were downregulated in the muscle of the hypoxic fish, as compared to the muscle of the normoxic fish. In addition, several GO slim terms related to the glycolysis pathway were upregulated in the muscle of the hypoxic fish as compared to the muscle of the normoxic fish, including the 6-phosphofructokinase complex, phosphofructokinase activity, and 6-phosphofructokinase activity. These enriched GO terms FIGURE 5 | An overview of the KEGG pathways significantly enriched in differentially expressed genes in the (A) muscle, (B) blood, (C) brain, and (D) liver. The specific pathways are plotted along the y-axis, and the x-axis indicates the enrichment factor. The size of the colored dots indicates the number of significantly differentially expressed genes associated with each corresponding pathway: pathways with larger-sized dots contain a higher numbers of genes. The color of the each dot indicates the corrected P-value for the corresponding pathway.
were supported by our KEGG pathway analysis. These results were consistent with the expression patterns of these genes in several other fish species (Gracey et al., 2001;Nikinmaa and Rees, 2005;Zhong et al., 2009).
The downregulation of genes associated with the tricarboxylic acid cycle (TCA) concurrent with the upregulation of genes associated with the glycolysis pathway in the muscle of hypoxic fish indicated a shift from aerobic oxidation to anaerobic glycolysis. However, glycolysis pathway genes were not upregulated in other tissues, possibly because muscle tissue has a higher glycolytic capacity than other tissues, and is thus better suited to supply the reduced metabolic needs of hypoxic G. eckloni. In addition, three genes involved in the pentose phosphate pathway were induced in the muscle of hypoxic fish, suggesting that muscle tissue may use pyrophosphate (PPi) as an energy source for hypoxic metabolism in G. eckloni.
Catabolism of gluconeogenic amino acids produces either pyruvate or an intermediate of the TCA cycle to generate carbon skeletons for gluconeogenesis (Gracey et al., 2001). Here, two amino acid pathways ("cysteine and methionine metabolism" and "arginine and proline metabolism") were enriched in the muscle. Several genes involved these two pathways were upregulated under hypoxia, including S-adenosylmethionine synthase, S-adenosylmethionine decarboxylase, brachypodium distachyon adenosylhomocysteinase, homocysteine methyltransferase, and prolyl 4-hydroxylase.
Glucose-6-phosphatase (G-6-Pase) was strongly expressed in the liver, and it was accompanied by hypoxia-induced gluconeogenesis. G-6-Pase catalyzes the dephosphorylation of glucose-6-phosphate to glucose, which is then circulated to other tissues to fuel glycolysis (Gracey et al., 2001;Zhong et al., 2009).
It is notable that the expression of 6-phosphofructo-2kinase/fructose-2, 6-biphosphatase (PFK-2) was downregulated in the blood under hypoxia, which is not consistent with observations in the muscles and livers of other fish (Gracey et al., 2001;Zhong et al., 2009). PFK-2 may contribute to the regulation of glycolytic flux under hypoxia (Yalcin et al., 2009;Zhong et al., 2009). Thus, in G. eckloni, amino acid catabolism coupled with gluconeogenesis in the liver, concurrent with a decrease in the transcriptional activity of PFK-2, may represent an alternative mechanism for the maintenance of blood glucose levels during hypoxia.
In the muscle, several genes involved in amino acid catabolism (e.g., S-adenosylmethionine synthase and S-adenosylmethionine decarboxylase) were upregulated under hypoxia. S-adenosylmethionine synthase and S-adenosylmethionine decarboxylase catalyze stages of methionine degradation. Consistent with the upregulation of these two genes, glutathione metabolism, which plays an important role in ammonia detoxification in the muscle (Gracey et al., 2001), increased significantly under hypoxia.
Studies have shown that toxic reactive oxygen species (ROS) are overproduced in fish muscle under hypoxic conditions (Lushchak et al., 2001;Bagnyukova et al., 2007;Clanton, 2007). Indeed, the activity levels of antioxidant enzymes, such as catalase, glutathione peroxidase, and superoxide dismutase, increase significantly to protect fish exposed to hypoxic stress (Lushchak et al., 2001;Zhong et al., 2009;Huang et al., 2015;Zhang et al., 2016;Wang et al., 2017). It is interesting that the key pathway related to the antioxidant defense system, ascorbate and aldarate metabolism, was significantly enriched in the muscle of G. eckloni, and that muscle was most affected by hypoxic stress. Three genes involved this pathway were upregulated under hypoxia: L-ascorbate oxidase, L-ascorbate peroxidase, and GDP-L-galactose phosphorylase. Both L-ascorbate oxidase and L-ascorbate peroxidase are oxidoreductases. The former acts on diphenols and related substances as the donor, with oxygen as the acceptor, while the latter acts on peroxides as the acceptor (peroxidases). They are the most important enzymes in ascorbic acid-glutathione detoxification, playing an important role in scavenging ROS and protecting cells from the destructive effects of ROS.
Consistent with this increase in ascorbate and aldarate metabolism, GDP-L-galactose phosphorylase, the key enzyme for ascorbic acid synthesis, was upregulated in the hypoxic G. eckloni. In addition, DyP-type peroxidase, a member of the haem peroxidase family, was upregulated in muscle, indicating the detoxification of the muscle exposed to hypoxic stress. DyPtype peroxidase degrades typical peroxidase substrates, as well as hydroxyl-free anthraquinone. Hypoxic conditions induce the expression of genes involved in the cellular stress response, such as heat shock proteins (HSPs). Ischemia-induced HSP70 expression has been demonstrated in the tissues of intact animals (Dillmann et al., 1986), and in cultured mammalian cells (Benjamin et al., 1990). HSP70 functions as a chaperone, protects cells from apoptosis, and plays important roles in reducing intracellular oxygen free radicals and maintaining the intracellular redox state (Fang et al., 2012). During hypoxia, regulation of the antioxidant system may be a biochemical mechanism that minimizes the destructive effects of free radicals (Zhong et al., 2009). While the enzymes involved in the antioxidant system differ somewhat among fish species, the activation of these enzymes also varies dependent on hypoxia regime.
There was significantly more insulin-like growth factor binding protein 1 (IGFBP-1) mRNA in the blood of G. eckloni under hypoxia as compared to the blood under normoxia. IGFBP-1 is a negative regulator of insulin and insulin-like growth factor in circulating blood (Gracey et al., 2001). Studies have shown that IGFBP-1 upregulation is correlated with intrauterine growth restriction in fetuses subjected to long-term chronic hypoxia (Chard, 1994); IGFBP-1 is also involved in the suppression of cell growth and proliferation in the euryoxic fish Gillichthys mirabilis after exposure to hypoxia (Gracey et al., 2001).
The ras-related and estrogen-regulated growth inhibitorlike protein (RERG), which negatively regulates the growth of cells involved in cardiac muscle development, was strongly upregulated in the heart of G. eckloni under hypoxia. In contrast, eukaryotic elongation factor 2 kinase (eEF-2K) was downregulated in G. eckloni muscle under hypoxia. eEF-2K is a dedicated kinase that regulates the activity of eukaryote extension factor 2 (eEF-2), its only known substrate, via phosphorylation at Thr-56 (Ryazanov et al., 1988;Dumont-Miscopein et al., 1994;Tavares et al., 2014). This decreases the affinity of phospho-eEF-2 for the ribosome, which impedes the elongation phase of protein translation, inhibiting the activity of eEF-2 (Carlberg et al., 1990;Dumont-Miscopein et al., 1994).
Interestingly, translationally-controlled tumor protein (TCTP) was also upregulated in the muscle. TCTP has similar properties to tubulin binding proteins, and associates with microtubules in a cell cycle-dependent manner (Bommer and Thiele, 2004). The overexpression of TCTP results in cell growth retardation (Gachet et al., 1999). Therefore, under hypoxia, the expression patterns of genes involved in the suppression of cell growth and proliferation changed, and protein synthesis rates decreased. These mechanisms may represent an especially important energy-saving strategy for the survival of G. eckloni in hypoxic environments.
Hypoxia-inducible factors (HIFs) have been recognized as the master regulators of the hypoxia signaling pathway; HIFs play important roles in maintaining dynamic oxygen balance and regulating changes in oxygen concentration under hypoxia (Xiao, 2015). An increase in Hif-1α/2α mRNA expression in response to hypoxia has been observed in mouse skeletal muscle (Stroka et al., 2001), in various tissues of Bos grunniens (Wang et al., 2006), and the liver and gill of Megalobrama amblycephala (Li et al., 2015). However, we detected no changes in the expression of Hif-1α/2α in G. eckloni subjected to hypoxia. Instead, egl nine homolog 1 (EGLN1) was significantly upregulated in the brain of G. eckloni. EGLN1 is considered homologous to prolylhydroxylase domain-containing protein 2 (PHD2), which plays an important role in the stabilization of Hif-1α. It is thus clear that the response of the HIF signaling pathway to hypoxia differs between mammals and fish. The HIF response also appears to vary among teleosts depending on species, hypoxia regime, and tissue. Further investigation is required to fully understand these variations.

CONCLUSION
We sequenced, assembled, and characterized the transcriptome of the Tibetan schizothoracine fish G. eckloni. Hypoxia had a greater effect on muscle than on any other tissue. Under hypoxia, G. eckloni underwent changes in gene expression, suggesting several key adaptive strategies. First, a reorganization of metabolic pathways and a transition from aerobic oxidation to anaerobic glycolysis. Second, a suppression of major energyrequiring processes, such as cell growth/proliferation and protein synthesis. In addition to these two strategies, which have also been observed in other teleosts, G. eckloni strengthens its antioxidant system and minimizes ischemic injury in response to hypoxia. This study is the first to analyze gene expression patterns in the Tibetan schizothoracine fish G. eckloni in response to hypoxic stress. Nevertheless, further studies using additional physiological and molecular approaches are necessary to verify the transcriptomic response of schizothoracine fish to hypoxia.

AUTHOR CONTRIBUTIONS
DQ designed the study and wrote the manuscript. YC and ZZ conducted the hypoxia treatment and collected the materials. RW, MX, and QC performed the experiments and analyzed the data. All authors have read and approved the final manuscript.