Transcriptome analysis of the liver of Eospalax fontanierii under hypoxia

Sampling and RNA sequencing

As the species (E. fontanierii) is considered an important agricultural pest and is not protected under any local, regional, national, or interational decree, we purchase those individuals from local farmers. Eighteen individuals of E. fontanierii (male and female, 220–280 g) were captured from agricultural land in Yan’an (N 35°09′, E 109°22′), Shaanxi Province, China. The species conservation status is Least Concern (LC), and the population trend is unknown. All animals were captured and treated humanely according to guidelines of the Care and Uses of Laboratory Animals of China, and all the procedures were approved by the Animal Care and Use Committee of Shaanxi Normal University (SNNU-IACUC-EAC-008-2010). Field experiments were approved by the Shaanxi Normal University, College of Life Science (project number: 18.11.20). Each animal was housed in a separate cage [475 L × 350 W × 200 H (mm)] maintained at 21 ± 1 °C under a 12-hour light/12-hour dark cycle. All animals were allowed free access to food (carrots). The adapted E. fontanierii were randomly divided into three groups (n = 6 per group): 6.5% O₂ 6 h (acute hypoxia), 10.5% O₂ 44 h (chronic hypoxia), and 21% O₂ (normoxia). In the normoxia group (21% O₂), animals breathed normal air for 1 week. In the chronic hypoxia group, animals were placed in a hypoxia chamber containing 10.5% O₂ for 44 h. In the acute hypoxia group, animals were placed in 6.5% O₂ hypoxia chamber for 4 h. The chamber was ventilated with nitrogen to maintain a constant oxygen concentration, which was monitored using a JRC-1020 thermo-magnetic analyzer. The animals were anesthetized with an intraperitoneal injection of pentobarbital (45 mg/kg) and sacrificed to collect fresh tissues and frozen into liquid nitrogen immediately. Total RNA was extracted using an RNA Simple Total RNA kit (TaKaRa) according to the instructions. The total RNA integrity was tested through 1% agarose gel electrophoresis. The RNA concentration was verified with a NanoDrop-2000 spectrophotometer (Thermo Fisher, USA). The cDNA library preparation was performed according to the standard protocol for Illumina sample preparation. The mRNA was enriched using magnetic beads with oligo(dT), and then random fragments of mRNA were generated by adding fragmentation buffer. The mRNA was reverse transcribed into the first-strand cDNA using a six-base random primer, and the second-strand cDNA was synthesized by addition of dNTPs, DNA polymerase I, RNase H, and buffer. The cDNA was purified with Ampure XP beads to obtain double-stranded cDNA. After end repair, the products were purified by AMPure XP beads and were amplified by PCR to construct cDNA library. The concentration and insert size of cDNA libraries were evaluated by Qubit2.0 and Agilent 2100, and effective concentration was estimated precisely by qRT-PCR. After complying with the quality control criteria, cDNA libraries were sequenced using the Illumina HiSeq™ 2500 sequencing platform, in which the read length is 150 bp, and 258.4 million pair-end reads were obtained. One schematic workflow for RNA-Seq analysis was given in Fig. S1.

Assembly and annotation

To obtain clean data, raw data containing adapters and primer sequences were removed, and low-quality bases were filtered by cutadapt (version 1.16) (Martin, 2011). Trinity (version: 2.4) was applied for de novo assembly of de Bruijn graphs and full-length transcripts (command: Trinity –seqType fq –left reads.fq –right reads_2.fq –CPU 24 –max_memory 256G) (Haas et al., 2013). Unigenes, a uniquely assembled transcript, were used to represent one or more transcripts in the same cluster. To display its sequence and avoid redundancy, the sequence of the primary mRNA with the highest expression level or the longest length was regarded as the unigene sequence when the whole genomic sequences were unknown. An isoform was first identified as the unigene sequence with the highest expression (>50% total expression value). If this criterion was not satisfied, the isoform with the longest length was considered to be the unigene sequence. The unigenes with at least 10 read in at least two samples were retained by in-house perl script. Software CD-HIT (version 4.8.1) was used to cluster unigene sequences with a sequence identity threshold of 0.9. RepeatMasker (version open-4.0.7) was used to mask the repeat elements in unigene sequences with nhmmscan as engine (version 3.1b1) and Mus musculus as the query species (Tarailo-Graovac & Chen, 2009). RNA sequences from 10 species used for homology searches were downloaded from NCBI RefSeq, and blastn was used to search their homologous genes in our data (Table S1). Protein sequences from 10 species were used to search the sing-copy genes by OrthoMCL, and then they were concatenated and aligned by MAFFT. RAxML was used to plot ML tree for the phylogenetic analysis (Fig. S2). BLAST (version: 2.6.0+) searches for unigene sequences were performed against Nr (blastx, E-value<1e−3, command: blastall -p blastx -i input -d nr -e 1e−3 -m 7 -v 20 -b 20 -o output) and Swiss-Prot (blastx, E-value <1e−10) databases.

To annotate the transcriptome, blastx searches against the Nr database were performed for all unigenes (E-value <1e−3). Functional annotation with gene ontology (GO) terms was conducted using Blast2GO software, which is designed for the high-throughput and automatic functional annotation of DNA or protein sequences based on the gene ontology vocabulary. Blast2GO Command Line (version:1.3.2) (downloaded from http://www.blast2go.com/blast2go-pro/blast2go-command-line) used the BLAST output to map and annotate unigenes (E-value <1e−6, command: Blast2GO_HOME/blast2go_cli.run -properties cli.prop -loadfasta input.fasta -loadblast blastResult.xml -mapping -annotation -saveb2g -savedat -annex -useobo go-basic.obo) (Conesa et al., 2005). Unigene sequences were uploaded into KOBAS3.0 to search for functional annotation in the Kyoto Encylopedia of Genes and Genomes (KEGG) (Wu et al., 2006). The predicted peptide sequences of unigenes were searched against Pfam (E-value ≤ 1e−5, command: phmmer -E 1e−5 -cpu 8 -pfamtblout output input.pep Pfam-A.hmm) using HMMER (version: 3.1b2) (Mistry et al., 2013). TransDecoder (version: 3.0.1) (Punta et al., 2012; Haas et al., 2013) was used to predict potential confident coding sequences (CDS, minimum length: 150 nucleotides) in unigenes and corresponding amino acid sequences based on open reading frame, log-likelihood score, and Pfam alignment information (command 1: TransDecoder.LongOrfs -t input.fasta; command 2: TransDecoder.predict -t input.fasta –retain_pfam_hits) (Haas et al., 2013). The code for data processing and analysis can be accessed in GitHub (https://github.com/Zhiqiang-hao/Scripts-for-transcriptome-analysis).

Expression analysis and DEG identification

We used a script from Trinity toolkit to quantitate transcript abundance (command: TRINITY_HOME/util/align_and_estimate_abundance.pl –transcript Trinity.fasta –seqType –left reads_1.fq –right reads_2.fq –est_method RSEM –aln_method bowtie –trinity_mode –prep_reference –out_dir rsem_outdir). To quantify gene expression levels, we used bowtie (version: 1.2.2) align reads from all samples against transcript set, and RSEM (version: 1.3.0) was used to estimate the expression abundance. We plotted the CPCoA plot using the expression matrix on webset: http://www.ehbio.com/ImageGP/index.php/Home/Index/CPCoAplot.html# (Fig. S3).

DESeq2 from R packages was used to identify DEGs (Love, Huber & Anders, 2014). Unigenes with a false discovery rate (FDR) <0.05 and fold change in expression >2 or <0.5 were considered to be DGEs. We used a perl script (run_DE_analysis.pl) in Trinity2.4 toolkits for DEG detection (Command: TRINITY_HOME/Analysis/Differential Expression/run_DE_analysis.pl –matrix counts.matrix –method DESeq2 –samples_file samples_ described.txt).

Functional enrichment analysis for GO terms was performed using the binGO app in Cytoscape3.4 (assesses over- or under-representation: overrepresentation, statistical test: hypergeometric test, multiple testing correction: Benjamini & Hochberg False Discovery Rate (FDR) correction, significance level: 0.05, reference set: Use the whole annotation as the reference set, ontology file: GO_Cellular_Component/Molecular_Function/Biological_Process) (Maere, Heymans & Kuiper, 2005). GO annotation was conducted with Blast2GO. KEGG enrichment for DEGs was completed in KOBAS3.0 (hypergeometric test, FDR correction <0.05) (Wu et al., 2006).

Real-time PCR validation of RNA-seq data

Six DEGs identified in response to hypoxia in the liver were selected to validate the reliability of the RNA-seq data. Gene expression was measured using Step One Real-Time System (ABI) with SYBR Premix ExTaq. Relative gene expression levels were normalized against that of an internal reference gene (β actin) and calculated using the Δ ΔCt method. Primers were designed using Primer-BLAST on the NCBI website. The expression of each gene was analyzed using three biological replications for each condition. Data were presented as the mean ± standard deviation (SD). SPSS 17.0 statistical software was used for statistical analysis of the data. The the statistical significance of the differences between the groups were evaluated using student’s t-test. P-values of <0.05 were considered to indicate statistical significance.

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. S4A).

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 S2).

Gene function annotation and CDS prediction of unigenes

Several databases, such as GO (Fig. S5), 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. S4B). 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. To cluster the expression pattern, co-expression cluster was performed for all genes (Fig. S6). In subcluster_1, 10,114 genes have the pattern with higher expression in high oxygen level in turn. In subcluster_8, 6,092 genes have higher expression in 10.5% and 21% O₂ than in 6.5%. In subcluster_2, 867 genes have higher expression pattern in 6.5% O₂ than in 10.5% and 21% O₂. In subcluster_4 and subcluster_9, there are 2,391 and 498 genes with higher and lower expression than in 6.5% and 21% O₂, separately. In subcluster_3, 1,242 genes have higher gene expression than in 6.5% and 10.5% O₂. A majority of the most highly expressed genes in the liver are very important for liver functions (Fig. S7, Table S3), in which multiple genes had higher expression in acute hypoxia.

To identify hypoxia-induced DEGs, DESeq2 was used for comparisons between the treatment groups and control groups (Love, Huber & Anders, 2014): 10.5% O₂ vs. 21% O₂; 6.5% O₂ vs. 21% O₂; 6.5% O₂ vs. 10.5% O₂, (Fig. S8). A total of 725 unigenes were considered to be DEGs in response to different hypoxia conditions (Table 3, Table S4). In total, 71.2% (516) of the DEGs were identified in the 6.5% O₂ vs. 10.5% O₂ comparison, while only 4.83% of the DEGs (35 genes) were identified in the 10.5% O₂ and 21% O₂ 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 (Fig. 1). In all the DEG sets, the numbers (176, 281, 29) of downregulated genes exceeded the corresponding numbers (134, 235, 6) of upregulated genes under low oxygen concentration when compared with high oxygen concentration (Table 3). It can be speculated that reduced transcriptional activities may contribute to the tolerance to hypoxia in E. fontanierii.

Protein-protein interaction for all DEGs

We got the protein-protein interaction from String (https://string-db.org) for all DEGs, and then Cytoscape was used to visualize the network (Fig. S9). We observed that up-regulated genes often were clustered closely, and down-regulated genes also have a center part. Some up-regulated gene nodes were connected with down-regulated gene nodes. As gene play functions by protein level finally and the post-transcriptional regulation act roles after gene transcription, the pattern of RNA level is not always consistent with the protein, which should be noticed in inferred phenotypic changes by RNA approach.

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. 2A, Table S5). 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) (Figs. 2B & 2C, Fig. S10, Table S6). 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. Other important DEGs with GO term “response to hypoxia” were also explored.

DEGs with the GO term “response to hypoxia”

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 (Kabil, Zhou & Banerjee, 2006). 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 (Ozono, 2006). 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 (Minet et al., 1999), and is considered to be the master transcriptional regulator of cellular and developmental responses to hypoxia (Iyer et al., 1998; Minet et al., 1999). 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 (Yang et al., 2016). EIF4EBP1 is a repressor of translation initiation that regulates eIF4E activity by preventing its assembly into the eIF4E complex (Yanagiya et al., 2012). In the AMPK signaling pathway, upregulated EIF4EBP1 inhibits protein synthesis and reduces cell activity to save energy.

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. 3A). As crucial signaling pathways, PPAR and AMPK signaling pathways were found to be enriched among the DEG sets (Fig. 3B, Fig. S11), 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. 4A). 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. 4B). Other enriched pathways playing key roles in liver functions were also detected, including steroid hormone biosynthesis, cholesterol metabolism, and bile secretion (Fig. 3B). 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. 3B). 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, such as FoxO signaling pathway and Glycolysis/gluconeogenesis pathway involved in liver cell survival and functional integrity were also explored.

FoxO signaling pathway

The growth arrest and DNA damage-inducible (GADD45beta) gene, which encodes a component of FoxO signaling pathway (Fig. 5A), 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% O₂) compared with normoxia (21% O₂) and chronic hypoxia (10.5% O₂), respectively. GADD45beta could prevent autophagy and apoptosis in rat (He et al., 2016). 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) was down-regulated by 2.48-fold (padj <0.05) and 2.71-fold (padj <0.05) in acute hypoxia (6.5% O₂) compared with normoxia (21% O₂) and chronic hypoxia (10.5% O₂), respectively. It can be speculated that PTEN downregulation contributes 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% O₂) compared with normoxia (21% O₂) and chronic hypoxia (10.5% O₂), respectively. It can be speculated that BNIP3 upregulation contributes to the tolerance to acute hypoxia in E. fontanierii liver.

Glycolysis /gluconeogenesis pathway

In 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 down-regulated under acute hypoxia compared with normoxia and chronic hypoxia, respectively. The two genes have the functions with the conversion of glucose to glucose-6-phosphate (G6P) (Fig. 5B). 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 (Ghosh et al., 2002). 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. 5B). These results showed that under hypoxia, glucose consumption in the liver is reduced.

Validation of six DEGs

In order to validate DEGs found in our result, we selected six genes that were found to be upregulated in the liver under hypoxia by RNA-seq by quantitative real-time PCR (qRT-PCR) (Figs. 6A–6F, Table 4, Table S7 showed raw data). 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.