RNA sequencing and de novo transcriptome assembly
In E. fontanierii, we performed RNA sequencing of liver organ, and 258.4 million pair-end raw reads were generated. After removal of adapters, primer sequences, and low-quality reads, 239.2 million pair-end clean reads were retained for further analysis. The guanine and cytosine (GC) contents ranged from 49.73% to 52.48%, and Q30 (99.9% base accuracy) scores were ≥85.00% (Table 1). An average of 26.6 million pair-end reads was generated for each sample. All reads were pooled for the de novo assembly. There were 39,439 unigenes, in which the N50 length was 2,556 nt, and the median length was 1,412 nt. The lengths of the unigenes ranged from 200~21,949, with the majority distributed in the range of 600~800 nt (Fig. S1A).
In total, 11,667,577 nt (15.6%) of unigene sequences were masked by RepeatMasker, the majority of which were retrotransposons (12.54%) and DNA transposons (0.96%). Alu/B1, B2, and B4 were the three most abundant short interspersed nuclear elements (SINEs). By calculating the best matches with genes/RNAs from nine different rodent species and Homo sapiens, homologous genes were identified for a total of 37,789 unigenes (ranges: 30,574~36,539) (blastn, E-value < 1e-5, Table S1).
Gene function annotation and CDS prediction of unigenes
Several databases, such as GO (Fig S2), Nr, Swiss-Prot, Pfam, and KEGG, were searched for functional annotation of unigenes. In total, 56.78% (22,395) unigenes were annotated at least once in the databases searched (Table 2). Furthermore, 16,754 unigenes ≥1,000 nt in length were annotated by the different databases.
A total of 33,003 CDSs were predicted from 30,893 unigenes, of which 1,888 unigenes possessed two or more CDSs. The mean and N50 lengths of CDSs were 814.9 nt and 1,554 nt, respectively. Many short CDSs were in the range of 200-300 nt in length (Fig. S1B). Among the CDSs identified, 58.7% (19,378) were complete. In addition, 46.4% (15,304) and 41.1% (13,561) of the CDSs were annotated by the UniRef90 and Pfam databases (E-value < 1e-5), respectively.
Read-mapping and differentially expressed gene identification
Based on RNA-Seq data, gene expression levels can be quantified by counting reads mapped to transcripts. This process is often influenced by changes in gene length and sequencing depth, alternative splicing, and gene duplication. RSEM enables accurate transcript quantification for species without sequenced genomes. To compare gene expression levels in multiple samples, FPKM values were used to measure the expression levels of unigenes. Overall, expression levels of 87.7% (34,591 of 39,439) of the unigenes were available with FPKM values ≥1 in at least one sample. A majority of the most highly expressed genes in the live are very important for liver functions (Table S2). For example, apolipoproteins (APOE, APOAI, APOA2, and APOC1) are involved in the metabolism of lipoproteins and their uptake in tissues [21]. These genes showed dominant or tissue-specific expression in the liver, with no significant changes in expression under hypoxia, indicating that these proteins facilitate the maintenance of liver function integrity. Hemopexin (HPX) prevents the pro-oxidant and pro-inflammatory effects of heme and also promotes its detoxification [22]. Fructose-bisphosphate aldolase B (ALDOB) is essential in fructose metabolism [22]. Glutathione S-transferase A1 (GSTA1) protects the cells from reactive oxygen species and the products of peroxidation [23]. HPX, ALDOB, and GSTA1 under normoxia and hypoxia could protect liver cells from oxidative damage and energy deficiency.
To identify hypoxia-induced DEGs, DESeq2 was used for comparisons between the treatment groups and control groups: 10.5% O2 vs. 21% O2; 6.5% O2 vs. 21% O2; 6.5% O2 vs. 10.5% O2, Fig. S3)[24]. A total of 725 unigenes were considered to be DEGs in response to different hypoxia conditions (Table 3, Table S3). In total, 71.2% (516) of the DEGs were identified in the 6.5% O2 vs. 10.5% O2 comparison, while only 4.83% of the DEGs (35 genes) were identified in the 10.5% O2 and 21% O2 comparison, suggesting a weak response to chronic hypoxia. Sixty- three upregulated DEGs and 62 downregulated DEGs were found in both the 6.5% vs. 10.5% and 6.5% vs. 21% comparisons. In all the DEG sets, the numbers of downregulated genes exceeded the number of upregulated genes under hypoxia when compared with the numbe. It can be speculated that reduced transcriptional activities may contribute to the tolerance to hypoxia in E. fontanierii.
GO enrichment of DEGs
To identify the functions of DEGs induced by hypoxic stress, GO terms were assigned and enriched. The GO term annotations of DEGs were mainly related to metabolic process, biological regulation, response to stimulus, catalytic activity, transporter activity, and molecular function regulator (Fig 1A, Table S4). We found the most enriched GO terms were commonly associated with lipid metabolism (58 genes), heme binding (18 genes), oxygen binding (7 genes), and response to oxygen-containing compounds (48 genes) (Fig. 1B&C, Fig. S4, Table S5). Many GO terms were associated with energy metabolism, such as fatty acid transporter activity, glucose metabolism, threonine metabolism, and acylglycerol metabolism. Along with the metabolism processes under hypoxia, we also identified some GO terms with response functions, such as response to carbohydrate, response to ethanol, response to zinc ion, response to hormone, response to hexose, and response to monosaccarides. We also identified GO terms associated with regulation of these metabolic and response processes, including regulation of hormone levels, negative regulation of gluconeogenesis, regulation of insulin secretion, regulation of hormone secretion, regulation of lipid metabolism, and regulation of peptide transport, of which the former one is statistically significant and the latter five are not significant. These results showed that hypoxic environments alter metabolism in the liver in E. fontanierii at numerous molecular levels, and the corresponding responses and regulations accommodate the changes in metabolism to maintain liver cell function.
KEGG enrichment of DEGs
To further explore the interactions among DEGs, the pathways containing DEGs were annotated and enriched using the KEGG database. The DEGs were mainly involved in lipid, amino acid, and carbohydrate metabolism, signal transduction, the digestive and immune system, as well as infectious diseases, endocrine and metabolic diseases (Fig. 2A). As crucial signaling pathways, PPAR and AMPK signaling pathways were found to be enriched among the DEG sets (Fig. 2B, Fig. S5), and both were associated with energy metabolism. In the AMPK signaling pathway, fatty acid oxidation and gluconeogenesis are enhanced under hypoxia stress, while protein synthesis and fatty acid synthesis are repressed (Fig 3A). In the PPAR signaling pathway, fatty acid oxidation, gluconeogenesis, and cholesterol metabolism were enhanced by upregulating the expression levels of related genes under hypoxic stress (Fig. 3B). Other enriched pathways playing key roles in liver functions were also detected, including steroid hormone biosynthesis, cholesterol metabolism, and bile secretion (Fig. 2B). Energy metabolism-associated pathways, such as glycolysis/gluconeogenesis, tyrosine metabolism, and linoleic acid metabolism, and pathways related to the immune system or disease, such as chemical carcinogenesis, antigen-processing and presentations, and metabolism of xenobiotics by cytochrome P450, were also observed (Fig. 2B). Those pathways enriched for DEGs were associated with basic liver functions, which may contribute to the adaptation of the liver under hypoxic stress. Other important pathways and key genes involved in liver cell survival and functional integrity were explored.
FoxO signaling pathway
The growth arrest and DNA damage-inducible (GADD45beta) gene, which encodes a component of FoxO signaling pathway (Fig. 4A), is involved in oxidative stress resistance and DNA repair. In this study, we detected 6.86-fold (padj<0.05) and 4.45-fold (padj<0.05) higher expression in acute hypoxia (6.5% O2) compared with normoxia (21% O2) and chronic hypoxia (10.5% O2), respectively. GADD45beta could prevent autophagy and apoptosis in rat [25]. Therefore, it can be speculated that upregulated GADD45beta expression protects E. fontanierii liver cells from apoptosis under acute hypoxia.
Phosphatidylinositol 3,4,5-trisphosphate 3-phosphatase and dual-specificity protein phosphatase PTEN (Pten) expression was downregulated by 2.48-fold (padj < 0.05) and 2.71-fold (padj < 0.05) down-regulated expression in acute hypoxia (6.5% O2) compared with normoxia (21% O2) and chronic hypoxia (10.5% O2), respectively. PTEN catalyzes dephosphorylation of protein acts as a phosphatase to dephosphorylate phosphatidylinositol 3,4,5-trisphosphate (PIP3), which functions as a tumor suppressor by negatively regulating the PI3K-Akt signaling pathway (Fig. 4A). The phosphatase activity of Pten may be involved in regulation of the cell cycle, preventing cells from growing, and dividing too rapidly [26]. Furthermore, Pten deletion allows nerve regeneration in mice [27]. It can be speculated that Pten downregulation inhibits liver cell apoptosis and promotes cell regeneration, thus contributing to the tolerance to acute hypoxia in E. fontanierii liver.
BCL2/adenovirus E1B 19 kDa protein-interacting protein 3 (BNIP3) displayed 3.47-fold (padj < 0.05) and 3.83-fold (padj < 0.05) higher expression under acute hypoxia (6.5% O2) compared with normoxia (21% O2) and chronic hypoxia (10.5% O2), respectively. BNIP3 is a member of the apoptotic Bcl-2 family that induces autophagy, apoptosis, and necrosis [28, 29]. Nuclear BNIP3 has been shown to acts as a transcriptional repressor to reduce apoptosis-inducing factor expression and increase resistance to apoptosis in human malignant gliomas [30]. Thus, combined with the other genes described, BNIP3 upregulation under acute hypoxia may help to control apoptosis to avoid irreversible liver injury.
Glycolysis /gluconeogenesis pathway
The glycolysis-related genes Glycolysis pathway, hexokinase (HK) (4.01 fold-change and padj < 0.05, 2.85 fold-change and padj = 0.15) and glucokinase (GCK) (5.05 fold-change and padj < 0.05, 3.19 fold-change and padj < 0.1) were downregulated under acute hypoxia compared with normoxia and chronic hypoxia, respectively, which impedes the conversion of glucose to glucose-6-phosphate (G6P) under acute hypoxia (Fig. 4B). G6PC encodes glucose-6-phosphatase, which catalyzes the conversion of G6P to a phosphate group and free glucose, which is the last step in gluconeogenesis [31]. In the liver of E. fontanierii, the upregulated expression of G6PC (4.17 fold-change, padj < 0.05) under acute hypoxia compared with chronic hypoxia will facilitate gluconeogenesis (Fig. 4B). These results showed that under hypoxia, glucose consumption in the liver is reduced, while production supply for other tissues, such as the brain, is increased.
Validation of six DEGs
We selected six genes that were found to be upregulated in the liver under hypoxia by RNA-seq for validation by quantitative real-time PCR (qRT-PCR) (Fig. 5A-F). Three of the genes [very long-chain acyl-CoA synthetase (Slc27a2), carnitine O-palmitoyltransferase 1 (Cpt1a), and cholesterol 7-alpha-monooxygenase (Cyp7a1)] were included in PPAR pathway, while the other three genes [Metallothionein (MT), C4b-binding protein beta chain (C4bpb), and Hyaluronidase-2 (Hyal2)] were key genes in liver functions (Fig. 3B). The qRT-PCR results were generally consistent with RNA-seq data, which showed that all of genes were significantly upregulated under acute hypoxia compared with normoxia or chronic hypoxia. In the qRT-PCR results, MT, Cpt1a, and Cyp7a1 were significantly upregulated under acute hypoxia compared with normoxia and chronic hypoxia, simultaneously. Slc27a2, C4bpb, and Hyal2 were significantly up-regulated in acute hypoxia and chronic hypoxia compared with normoxia. As these genes are highly important in maintaining liver function, we investigated their functions under hypoxic stress.
Very long-chain acyl-CoA synthetase (Slc27a2), which converts free long-chain fatty acids into fatty acyl-CoA esters, and plays key role in lipid biosynthesis and fatty acid degradation [32]. Carnitine O-palmitoyltransferase 1 (Cpt1a) catalyzes the transfer of the acyl group from CoA to carnitine to form palmitoylcarnitine as the first and rate-limiting step in the carnitine palmitoyltransferase system. The acyl carnitine is then shuttled across the inner mitochondrial membrane by a translocase [33]. Slc27a2 and Cpt1a are crucial for the activation and oxidation of fatty acids. In the E. fontanierii liver transcriptome, Slc27a2 and Cpt1a were upregulated by 2.14-fold (padj < 0.05) and 3.19-fold (padj < 0.05) under acute hypoxia compared with normoxia, respectively. The results suggested that the upregulated expression of Slc27a2 and Cpt1a in the liver of E. fontanierii under acute hypoxia could increase the local energy supply through oxidation of fatty acids, especially when the levels of glucose as the energy supply in the liver were decreased.
Cholesterol 7-alpha-monooxygenase (Cyp7a1), which plays an important role in cholesterol metabolism, converts cholesterol to 7-alpha-hydroxycholesterol, in the first and rate-limiting step in bile acid synthesis [34]. Our results showed that Cyp7a1 in the liver was upregulated by 3.29-fold (padj < 0.05) under acute hypoxia compared with normoxia, which could contribute to bile acid biosynthesis and regulation of cholesterol levels.
Metallothionein (MT) has been reported to function as a negative regulator of apoptosis [35]. In E. fontanierii liver, MT contains 61 residues, including 20 cysteine residues. As a cysteine-rich protein, MT performs important functions in the control of oxidative stress, by capturing damaging oxidant radicals, such as superoxide, and hydroxyl radicals via the cysteine residues [36]. In our study, MT was upregulated by 3.26-fold (padj < 0.05) and 3.22-fold (padj < 0.1) under acute hypoxia compared with normoxia and chronic hypoxia, respectively. These results indicate that MT protects E. fontanierii liver tissue against the harmful influences of acute hypoxia.
C4b-binding protein beta chain (C4bpb), which is the main inhibitor of the classical complement activation pathway, accelerates the decay of C3-convertase and hydrolyzes the C4b complement fragment [37]. It also interacts with anticoagulant protein S, and binds apoptotic and necrotic cells as well as DNA to clean up after injury and limit the inflammatory potential of necrotic cells [38, 39]. In this study, C4bpb was upregulated by 2.78-fold (padj < 0.05) under acute hypoxia compared with normoxia. The finding indicated that upregulated C4bpb in acute hypoxia plays a role in preventing the inflammation and blood coagulation induced by acute hypoxia. C4bpb negatively regulates complement activation, and its upregulated expression may helps to reduce coagulation under hypoxia, which is consistent with previous study in E. fontanierii heart tissue [18].
Hyaluronidase-2 (Hyal2) is thought to be involved in cell proliferation, migration, and differentiation [40, 41]. Various functions have been described for the gene, such as response to reactive oxygen species, positive regulation of inflammatory response, negative regulation of protein kinase, cellular response to tumor necrosis factor, and homeostatic processes [42]. In our study, Hyal2 was upregulated by 2.85-fold (padj < 0.05) and 3.04-fold (padj < 0.05) under acute hypoxia compared with normoxia and chronic hypoxia, respectively, suggesting that upregulated Hyal2 is an important gene in regulating multiple responses to hypoxia.
DEGs with the GO term “response to hypoxia”
Apart from the DEGs validated by qRT-PCR, we identified some DEGs assigned the GO term “response to hypoxia”, which may play key roles in hypoxia adaptation in subterranean animals and have important additional functions. Compared with normoxic conditions, cystathionine beta-synthase (Cbs) and heme oxygenase 1 (Hmox1) were upregulated by 2.11-fold (padj < 0.05) and 3.29-fold (padj < 0.05) under acute hypoxia, respectively. Cbs catalyzes the first step in the trans-sulfuration pathway, from homocysteine to cystathionine, which is the precursor of cysteine. In mammals, Cbs is a highly regulated enzyme, which contains a heme cofactor that functions as a redox sensor [43]. In our study, Cbs was assigned to multiple GO terms that may contribute to liver hypoxia tolerance, such as oxygen binding, carbon monoxide binding, cysteine synthase activity, nitrite reductase (NO-forming) activity, selenocystathionine beta-synthase activity, superoxide metabolic process and negative regulation of apoptotic. Hmox1 catabolizes free heme, produces carbon monoxide (CO), and induces the upregulation of interleukin 10 (IL-10) and interleukin 1 receptor antagonist (IL-1RA), which form the basis of its anti-inflammatory properties [44]. Certain important GO terms associated with hypoxia adaptation, such as liver regeneration, negative regulation of extrinsic apoptotic signaling pathway via death domain receptors, positive regulation of angiogenesis, erythrocyte homeostasis, regulation of blood pressure and cellular iron ion homeostasis, were assigned to Hmox1. Under hypoxia, hypoxia-inducible factor 1-alpha (HIF-1-alpha) is often significantly upregulated [45], and is considered to be the master transcriptional regulator of cellular and developmental responses to hypoxia [45, 46]. As a component of the HIF signaling pathway, HIF-1 could induce upregulation of Hmox1 expression under acute hypoxia to promote protection of the liver.
Compared with chronic hypoxia, heme oxygenase 2 (Hmox2) and eukaryotic translation initiation factor 4E-binding protein (EIF4EBP1 or 4E-BP1) were upregulated by 2.02-fold (padj < 0.05) and 3.20-fold (padj < 0.05) under acute hypoxia, respectively. As a modifier in the regulation of hemoglobin metabolism, Hmox2 has been reported to contribute to high-altitude adaptation in Tibetans [47]. EIF4EBP1 is a repressor of translation initiation that regulates eIF4E activity by preventing its assembly into the eIF4E complex [48]. In the AMPK signaling pathway, upregulated EIF4EBP1 inhibits protein synthesis and reduces cell activity to save energy (Fig. 3A).