Skip to main content
  • Research article
  • Open access
  • Published:

Differential responses of Lasiopodomys mandarinus and Lasiopodomys brandtii to chronic hypoxia: a cross-species brain transcriptome analysis



Subterranean rodents have evolved many features to adapt to their hypoxic environment. The brain is an organ that is particularly vulnerable to damage caused by exposure to hypoxic conditions. To investigate the mechanisms of adaption to a hypoxic underground environment, we carried out a cross-species brain transcriptome analysis by RNA sequencing and identified genes that are differentially expressed between the subterranean vole Lasiopodomys mandarinus and the closely related above-ground species Lasiopodomys brandtii under chronic hypoxia [10.0% oxygen (O2)] and normoxia (20.9% O2).


A total of 355 million clean reads were obtained, including 69,611 unigenes in L. mandarinus and 69,360 in L. brandtii. A total of 235 and 92 differentially expressed genes (DEGs) were identified by comparing the hypoxic and control groups of L. mandarinus and L. brandtii, respectively. A Gene Ontology (GO) analysis showed that upregulated DEGs in both species had similar functions in response to hypoxia, whereas downregulated DEGs in L. mandarinus were enriched GO terms related to enzymes involved in aerobic reactions. In the Kyoto Encyclopedia of Genes and Genomes pathway analysis, upregulated DEGs in L. mandarinus were associated with angiogenesis and the increased O2 transport capacity of red blood cells, whereas downregulated DEGs were associated with immune responses. On the other hand, upregulated DEGs in L. brandtii were associated with cell survival, vascular endothelial cell proliferation, and neuroprotection, while downregulated genes were related to the synaptic transmission by neurons.


L. mandarinus actively adapts its physiological functions to hypoxic conditions, for instance by increasing O2 transport capacity and modulating O2 consumption. In contrast, L. brandtii reacts passively to hypoxia by decreasing overall activity in order to reduce O2 consumption. These results provide insight into hypoxia adaptation mechanisms in subterranean rodents that may be applicable to humans living at high altitudes or operating in other O2-poor environments.


Oxygen (O2) is essential to sustain most aerobic organisms, and O2 deprivation creates a significant physiological stress. Hypoxia occurs under natural conditions, for instance in underground tunnels [1], at high altitude [2], in aquatic habitats [3], and in the tumor microenvironment [4, 5]. For many vertebrates, disruption of the O2 supply to the brain for more than a few minutes leads to irreversible neurological damage due to neuronal death [6]. Species that inhabit subterranean environments have developed effective strategies to survive under hypoxia and exhibit a variety of convergent morphological, physiological, behavioral, and genomic adaptations [7,8,9].

Recent studies have revealed the molecular mechanisms underlying the response to hypoxia in subterranean mammals, including the differential expression of genes encoding hemoglobin, hypoxia inducible factor (HIF)-1α, erythropoietin, and vascular endothelial growth factor (VEGF) [10,11,12]. The expression patterns of these genes differ between subterranean mammals and mammals living above ground; the former have evolved a unique cardiovascular system as an adaptation of hypoxia. A large-scale transcriptome sequencing study of the blind mole rat (Spalax galili), a typical subterranean rodent, revealed that apoptosis was suppressed and the expression of angiogenic factors was tightly regulated in the hypoxic environment [13, 14]. However, it is possible that different species of subterranean rodents have evolved distinct adaptive mechanisms in response to chronic hypoxia [15, 16].

The Mandarin vole (Lasiopodomys mandarinus) is a subterranean rodent that is widely distributed throughout northeast and central China and north central Mongolia, as well as in the adjacent areas of Siberia south of Lake Baikal and the southern and central Korean Peninsula. For most of its life, L. mandarinus lives in an underground tunnel system characterized by chronic hypoxia and darkness. As a subterranean species, L. mandarinus exhibits remarkable physiological adaptations to hypoxic stress including a higher capillary density and elevated values of blood parameters such as hematocrit, mean corpuscular volume, mean corpuscular hemoglobin (MCH), and MHC concentration (MCHC) [17]. As the sibling species of L. mandarinus, Brandt’s vole (L. brandtii) is mainly distributed in the grasslands of middle-eastern Inner Mongolia, eastern regions of Mongolia, and some parts of southern Russia [18]. Unlike L. mandarinus, L. brandtii spends most of its life above ground. Their close evolutionary relationship and distinct life histories make L. mandarinus and L. brandtii ideal animal models for a comparative study of mechanisms of adaption to hypoxia in subterranean mammals.

In the present study, we sequenced and assembled the brain transcriptomes of L. mandarinus and L. brandtii under conditions of chronic hypoxia and normoxia. Whole brain RNA was extracted and subjected to Illumina sequencing to identify genes that are differentially expressed under the two conditions as well as between the two species.


Illumina sequencing and de novo transcriptome assembly

A total of 355 million reads with 89.4 billion bases were obtained after stringent quality assessment and data filtering (Additional file 1: Table S2). We obtained 144,789 transcripts (mean length: 1989.39) and 69,611 unigenes (mean length: 944.46) with an N50 of 2214 for L. mandarinus, and 167,002 transcripts (mean length: 2501.43) and 69,360 unigenes (mean length: 974.89) with an N50 of 2306 for L. brandtii. The length distribution of assembled unigenes is shown in Additional file 1: Table S3.

Functional annotation

According to the BLASTX results, 20,011 (28.75%) unigenes of L. mandarinus and 19,120 (27.70%) unigenes of L. brandtii had homologous proteins in the NCBI non-redundant (Nr) database (Additional file 1: Table S4). Based on annotated unigenes in the database, 10,593 and 9838 unigenes of L. mandarinus and L. brandtii, respectively, were assigned to one or more Gene Ontology (GO) terms, with 33.5%/32.5% in cellular components, 21.1%/21.7% in molecular functions, and 45.4%/45.8% in biological processes (Fig. 1). To identify biological pathways that are differentially regulated between L. mandarinus and L. brandtii, the annotated unigenes were mapped to reference pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. L. mandarinus and L. brandtii unigenes were mapped to 368 pathways.

Fig. 1
figure 1

GO annotations for L. mandarinus and L. brandtii transcriptomes

Gene expression pattern analysis

To detect genes that are differentially expressed under hypoxia and normoxia, RSEM and edgeR were used with a false discovery rate (FDR) threshold of ≤0.05 and fold change of ≥2. For L. mandarinus, 695 differentially expressed genes (DEGs) were identified from a total of 69,611 unigenes, of which 425 were upregulated and 270 were downregulated in the hypoxic brain relative to the normoxic brain (Additional file 1: Figure S1). Of the 695 DEGs, 289 were annotated in at least one of the following databases: Nr (n = 270), Nt (n = 208), Swissprot (n = 239), and KEGG (n = 170). Of the 270 DEGs annotated in the Nr database, 47 were described as “hypothetical” in the GenBank database and were therefore excluded from further analysis. Ultimately, 223 DEGs that were annotated in Nr were screened, of which 110 were upregulated and 113 were downregulated (Additional file 2: Table S5).

For L. brandtii, there were only 158 DEGs among a total 69,360 unigenes, with 104 upregulated and 54 downregulated in the hypoxic brain relative to the normoxic brain (Additional file 1: Figure S1). Of the 158 DEGs, 110 were annotated in at least one of the following databases: Nr (n = 103), Nt (n = 76), Swissprot (n = 93), and KEGG (n = 73). Of the 103 DEGs annotated in the Nr database, 14 were described as “hypothetical” in the GenBank database and were therefore excluded from further analysis. Ultimately, 89 DEGs that were annotated in Nr were screened, of which 49 were upregulated and 40 were downregulated (Additional file 2: Table S5).

DEG functional enrichment analysis

The GO enrichment analysis indicated that the upregulated DEGs of L. mandarinus were enriched in seven GO terms in three categories—i.e., biological process (n = 2), cellular component (n = 3), and molecular function (n = 2). There were also five enriched GO terms among the downregulated DEGs in L. mandarinus. For L. brandtii, all the DEGs were upregulated and were enriched in 10 GO terms. (Fig. 2a and Additional file 1: Table S6).

Fig. 2
figure 2

a GO terms and b KEGG pathways significantly enriched for up- and downregulated DEGs in L. mandarinus and L. brandtii

The GO categories enriched for the upregulated DEGs of L. mandarinus represented biological functions such as endothelial cell proliferation (Extracellular region and growth), cell migration, cell differentiation (Biological adhesion), gene expression (Helicase activity), and energy harvesting (Peptidase activity); whereas enriched GO terms for downregulated DEGs included coenzymes (Generation of precursor metabolites and energy), oxidoreductases (Oxidoreductase activity), and proteases (Peptidase activity) related to aerobic reactions. Notably, peptidase activity was found to be enriched among DEGs that were up- and downregulated in L. mandarinus (Additional file 1: Table S7).

Three GO terms were shared by L. brandtii and L. mandarinus, whereas six differed between the two species that were related to vascular endothelial cell proliferation (Cell proliferation), cell differentiation (Cell differentiation), and cell migration (Anatomical structure formation involved in morphogenesis).

DEG pathway analysis

To clarify the relationships between the DEGs, we mapped the genes in the KEGG pathway database and performed enrichment analysis with Fisher’s exact test. Pathways with fewer than three genes were discarded, and those with both P value < 0.05 and FDR < 0.05 were selected as enriched pathways. We identified eight enriched pathways for upregulated and five for downregulated DEGs in L. mandarinus, and three pathways for upregulated and six for downregulated DEGs in L. brandtii (Fig. 2b and Additional file 1: Table S8).

Among the eight enriched pathways for upregulated DEGs in L. mandarinus, the major functions were angiogenesis, cell proliferation, and apoptosis (AGE-RAGE signaling pathway in diabetic complications, Bladder cancer, and HIF-1 signaling). Some were involved in the regulation of cell transendothelial migration (Leukocyte transendothelial migration and Cell adhesion molecule), inhibition of angiogenesis, or induction of cell cycle arrest, cell repair, and apoptosis (p53 signaling). The DEGs in the malaria pathway included two subunits of heme (HBA and HBB). An increase in heme can significantly increase the O2-carrying capacity of red blood cells. Of the five enriched pathways for downregulated DEGs, most were related to the immune response and shared a common feature of O2 consumption.

In L. brandtii, the pathways associated with the three upregulated DEGs were linked to cell survival (Maturity onset diabetes of the young), extracellular matrix (ECM)-related vascular endothelial cell proliferation (ECM-receptor interaction), and neuroprotection (Focal adhesion). However, downregulated DEGs were enriched in six KEGG pathways that were related to synaptic transmission in neurons.

Comparative analysis of unique differentially expressed genes

To further analyze the similarities and differences between the L. mandarinus and L. brandtii hypoxia adaptations, we removed DEGs shared by the two species and constructed protein interaction networks for their unique DEGs (Fig. 3). Among the chronic hypoxia-specific differentially expressed genes in L. mandarinus, the most prominent cores were MMP2, THBS1, SERPINE1, and VEGFR-2, which are responsible for angiogenesis regulation /inhibition. It is worth noting that MBP is a marker protein of cranial nerve damage that is regulated by multiple proteins. The down-regulation of MBP as well as the expression of other interacting proteins indicates that neurons in L. mandarinus are less prone to hypoxia-induced injury. The most significant cores of the chronic hypoxia-specific differentially expressed gene protein interaction network in L. brandtii are dopamine synthesis and transport regulation-related proteins. The reduction of these proteins leads to decreased neurogenic excitability and decreased activity in L. brandtii.

Fig. 3
figure 3

Protein interaction network for specific DEGs in L. mandarinus (a) and L. brandtii (b) brain under chronic hypoxia

Validation of DEGs by reverse-transcription quantitative (RT-q)PCR analysis

To validate the expression data obtained from RNA sequencing (RNA-Seq), seven genes with different expression patterns between L. mandarinus and L. brandtii were randomly selected to perform qPCR. Among the seven selected genes, five genes were up-regulated and two were down-regulated in L. mandarinus, compared with L. brandtii. The results showed a strong correlation between the data of RNA-Seq and those of qPCR (R = 0.887, P = 0.004 for L. mandarinus and R = 0.838, P = 0.019 for L. brandtii; Fig. 4a), indicating that the data from our transcriptome analysis were reliable. The expression profile of the seven genes validated via RT-qRCR is shown in Fig. 4b.

Fig. 4
figure 4

Validation of RNA sequencing results by RT-qPCR. a Correlations between gene expression levels measured by RT-qPCR and RNA-Seq methods. b Comparison of RNA-Seq log2FoldChange read counts with log2FoldChange RT-qPCR copy numbers. The upper panel shows the RNA-Seq read counts (log2FoldChange) for seven genes, of which five are upregulated in chronic hypoxia vs. normoxia in L. mandarinus and L. brandtii brain and two are downregulated. The lower panel shows log2FoldChange by RT-qPCR copy numbers for the same genes in chronic hypoxia vs. normoxia in L. mandarinus and L. brandtii brain samples. (LM: L. mandarinus; LB: L. brandtii)


In this study, we carried out a comparative transcriptome analysis between L. mandarinus and L. brandtii whole brain tissues exposed to chronic hypoxia vs. normoxia. Our results reveal that similar GO terms were enriched for upregulated DEGs in L. mandarinus and L. brandtii under chronic hypoxia; these were associated with endothelial cell proliferation, cell migration, gene expression, angiogenesis, angiogenesis inhibition, energy acquisition, O2 transport, neuroprotection, and protection from carbon dioxide [13, 19] (Fig. 5 and Additional file 1: Table S6). This suggests that both L. mandarinus and L. brandtii respond to chronic hypoxic stress by increasing the number of blood vessels, improving blood flow, and enhancing O2 transport capacity. Additionally, in L. mandarinus, enriched GO terms for downregulated DEGs were related to cofactors and enzymes such as NADH dehydrogenase, NAD-dependent malic enzyme, prenylcysteine oxidase, and peroxisomal sarcosine oxidase (Fig. 2a and Additional file 2:Table S5), which are also enriched in the Upper Galilee Mountains blind mole-rat Spalax galili [20]. A previous study of S. galili indicated that tight control of angiogenesis may be a novel mechanism of hypoxia tolerance in subterranean rodents [14]. Our results demonstrate that as O2 supply decreases, L. mandarinus and L. brandtii reduce O2 consumption—the former, by suppressing O2-consuming reactions such as redox reactions, proteolysis, and coenzyme metabolism—and/or redirect O2 usage. Besides, hypoxia-tolerance studies of heterothermic Heterocephalus glaber (the naked mole-rat) found that it can endure tissue hypoxic conditions by actively reducing brain oxygen consumption, and it was suggested that its protective strategies include modulation of immune response, thrombolysis, antioxidant defense, and activation or inactivation of pre-existing proteins [21, 22]. These results suggest that reducing brain oxygen consumption may be a common way for subterranean rodents to adapt to hypoxia conditions.

Fig. 5
figure 5

Presumed responses of L. mandarinus and L. brandtii to chronic hypoxia. Red and green lines indicate up- and downregulated gene enrichment, respectively. The images of voles were photographed in the animal laboratory of the school of life sciences, Zhengzhou University

Peptidase activity was enriched among up- and downregulated DEGs in L. mandarinus (Fig. 2a and Additional file 1: Table S7). The upregulated group showed the activation of genes such as A disintegrin and metalloproteinase with thrombospondin motifs (ADAMTS)1 and ADAMTS9, which are involved in angiogenesis inhibition. In addition, complement factor B (CFB) and MMP2 are important for angiogenesis [23]; plasminogen activator tissue type (PLAT) contributes to thrombolysis [24], and autophagy 4C (ATG4C) is related to autophagy [25, 26]. Among the downregulated DEGs, transmembrane protease serine 5 (TMPRSS5) and chymotrypsin-like elastase family member 1 (CELA1) are serine proteases that activate the complement system and induce blood coagulation; fibrinogen-like protein 2 (FGL2) and coagulation factor IX (F9) are highly expressed in vascular endothelial cells and accelerate blood coagulation; and ADAMTS4 and ADAMTS2 are collagen precursors that are involved in coagulation and cell support.

Enriched pathways for upregulated DEGs in L. mandarinus were associated with angiogenesis (HIF-1 signaling, bladder cancer, and AGE-RAGE signaling pathway in diabetic complications), transendothelial migration (leukocyte transendothelial migration and cell adhesion molecules), angiogenesis inhibition, cell cycle arrest for cell repair, and apoptosis (p53 signaling pathway). These pathways may be involved in increasing blood vessels, enhancing red blood cell O2 transport capacity, and repairing damaged cells. Enriched pathways for downregulated DEGs in L. mandarinus were mainly related to aerobic metabolism and immune response, which could supply additional O2 for cell survival (Fig. 5).

The metabolic pathway network in L. mandarinus can be roughly divided into three categories (Fig. 6a). The first includes pathways and genes related to hypoxia stress, angiogenesis and its inhibition, apoptosis, and cell repair; and the second involves genes that are associated with the ECM and regulate blood flow. Most of these genes were upregulated. Interestingly, the results for the first category were consistent with research findings for S. galili compared with Spalax judaei. Avivi et al. uncovered that S. galili exhibits species-specific responses to hypoxic stress via numerous genes involved in angiogenesis, apoptosis, and oxidative stress management [27]. In contrast, the third category comprises downregulated genes and pathways involved in amino acid synthesis and energy metabolism. Other studies that compared transcript abundance in Spalax vs. rat whole brain tissues showed that down-regulated genes in Spalax were significantly associated with carbohydrate metabolism, lipid metabolism, redox metabolism, mitochondria inner membrane activity, and oxidative phosphorylation [20]. These findings on the hypoxic adaptation of Spalax are similar to those of our third category in L. mandarinus. Genes with more connections in the gene-metabolism network such as VEGFA, matrix metallopeptidase (MMP)2 [28], serpin peptidase inhibitor [29], and thrombospondin 1 are associated with angiogenesis and its inhibition; their upregulation indicates that angiogenesis in L. mandarinus is in a state of dynamic equilibrium under chronic hypoxia.

Fig. 6
figure 6

Gene-pathway networks for DEGs in (a) L. mandarinus and (b) L. brandtii under chronic hypoxia. Nodes (circles) represent DEGs and enriched pathways for DEGs, respectively. Genes indicated in red and blue are up- and downregulated, respectively. Lines between nodes represent connections between genes and pathways in the network

Upregulated pathways unique to L. brandtii hypoxic brain were involved in ECM-related endothelium cell proliferation and cell survival (ECM-receptor interaction and focal adhesion) as well as neuroprotection (maturity onset diabetes of the young) [30]. It is worth noting that all six downregulated pathways in L. brandtii were related to synaptic transmission, suggesting that L. brandtii reduces its activity level in order to reduce O2 consumption under hypoxic conditions (Fig. 5). The DEG-metabolic pathway network of L. brandtii can be divided into two parts according to gene function (Fig. 6b). The first part mainly includes ECM and adhesion-related genes such as Vwf (Von Willebrand factor) [31], COL4A1 (collagen, type IV, alpha 1), ITGA4 (integrin alpha 4), and VEGFA. The increased expression of these genes is linked to vascular endothelial cell proliferation and angiogenesis [32]. The second part of the metabolic network is related to synaptic transmission. Key genes in the network including TH (tyrosine hydroxylase), GRIN2D (glutamate ionotropic receptor NMDA type subunit 2D), solute carrier family (SLC)6A3 (neurotransmitter transporter, dopamine), and SLC18A2 were downregulated. These four genes are involved in the production, transport, and synaptic transmission of the neurotransmitter dopamine; their downregulation decreases dopamine release, thereby reducing neuronal excitability [33].

Differences in gene expression between the two species of vole in this study under chronic hypoxia are mainly determined by their distinct life histories. L. mandarinus lives in an underground tunnel system for most of its life; during summer when there is abundant rainfall, tunnels can collapse as the earth becomes wet. To repair the tunnel under hypoxic conditions, L. mandarinus maintains a high blood pressure and neuronal excitability and has evolved an enhanced O2-transport capacity as well as mechanisms to reduce O2 consumption and protect the brain from hypoxic injury.

In contrast, L. brandtii also burrows in tunnels but this is mainly used for rest; O2-consuming activities such as mating, feeding, and excavation are performed above ground. A hypoxic environment develops when snowfall covers the ground in winter. During this time, L. brandtii performs few activities in the tunnel and there is little risk of tunnel collapse. This explains the relatively passive response of L. brandtii to chronic hypoxia, including reductions in blood pressure and nerve excitability.

Presently, research on activities related to hypoxia, such as human life on the plateau or diving, has received much attention. Our findings may shed new light on human adaptation to hypoxia conditions. In fact, many people, including Tibetans, Andeans, and Ethiopians are currently living in long-term hypoxic plateau environments [34]. However, the mechanisms of high-altitude adaptation in different human populations are distinct and complex. Mounting evidence suggests that the observed physiological adaptations are controlled by interactions among multiple genes, especially those that are part of the HIF pathway. Tibetans inhale more air with each breath and breathe more rapidly than either sea-level inhabitants or Andeans [35]. In addition, they have high levels of NO (nitric oxide) in their blood when compared to low land dwellers, and this probably aids the dilation of their blood vessels to enhance blood circulation [36]. The patterns of genetic adaptation among Andeans are largely distinct from those of Tibetans. For instance, the lack of significant associations between EPAS1 or EGLN1 SNP genotypes and hemoglobin concentration is characteristic of Tibetans [37]. For Ethiopians, adaptation to hypoxia involves several candidate genes including CBARA1, VAV3, ARNT2, and THRB [38]. THRB and ARNT2 are known to play a role in the HIF-1 pathway, a pathway implicated in a previous study conducted in Tibetan and Andean populations [39]. In addition to facing hypoxic conditions in the living environment, some special human activities, such as diving, frequently expose certain individuals to hypoxic conditions. Early studies of diving populations suggested that they were adapted to hypoxic conditions via bradycardia and peripheral vasoconstriction, which lower oxygen consumption and selectively redistribute blood flow to the organs most sensitive to hypoxia [40, 41]. A recent comparative genomic study revealed natural selection of genetic variants in PDE10A in the “Sea Nomads” (The indigenous Bajau people) and suggested that mutations in this gene led to an increase in spleen size, thereby providing the people with a larger reservoir of oxygenated red blood cells [42].

Compared with human hypoxia adaptation studies, although some physiological or pathway changes, such as adjustment of angiogenesis and contraction or interaction of some genes in the HIF signaling pathway or other pathways (VEGF signaling pathway), were similar to the results of our studies on L. mandarinus, the genes involved in these pathways differed between humans and L. mandarinus. These may be related to the large differences in hypoxia conditions between humans and subterranean rodents.


Long-term existence in underground tunnels has enabled L. mandarinus to better adapt to hypoxia than the closely related species L. brandtii. At 10% O2, L. mandarinus actively adapts its physiological functions by increasing O2 transport capacity and reducing O2 consumption, whereas L. brandtii reacts passively by decreasing its activity. These results provide insight into mechanisms of hypoxia adaptation in subterranean rodents that may be applicable to humans living in high-altitude or other O2-poor environments.


Animals and hypoxia treatment

L. mandarinus was trapped live from croplands in Xinzheng, Henan, China (N 34°52′, E 113°85′), and L. brandtii was imported from the Chinese Academy of Agricultural Science. The animals were maintained individually in polycarbonate cages (37 × 26 × 17 cm) in the laboratory on a 14:10-h light/dark cycle at a temperature of 20 °C–24 °C.

To mimic chronic hypoxic stress, 12 healthy adult male voles (3 months of age; n = 6 of each species) were randomly divided into the hypoxia (10% O2 for 48) and normoxia (20.9% O2 for 48 h) groups. A DS-II hyperbaric cabin (Huaxin Hyperbaric Cabin, Weifang, China) was used to simulate chronic hypoxia. O2 level in the cabin was maintained at a constant level by balancing the flow rate of O2 and N2, and was monitored with an Oxygen analyzer (Talantek, Beijing, China). A bottle of sodium hydroxide was placed in the cabin to absorb the CO2 released by the animals. Once the treatment was completed, the animals were immediately sacrificed with an overdose of pentobarbital sodium. The brain was removed and immediately frozen in liquid nitrogen and stored at − 80 °C until use.

RNA extraction, cDNA library preparation, and RNA-Seq

Experimental procedures including sample preparation and RNA-Seq were performed according to standard protocols. Total RNA was extracted from each sample using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The RNA was treated with RNase-free DNase I (Takara Bio, Dalian, China) to remove residual DNA. RNA integrity was verified by agarose gel electrophoresis (1.2%), and RNA concentration was measured with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

High-quality RNA samples were sent to Biomarker Technologies Corp. (Beijing, China) for cDNA library construction and sequencing, with mRNAs purified through interaction of the poly(A) tails and magnetic oligo(dT) beads. RNA-Seq libraries were generated using the TruSeq RNA Sample Prep kit (Illumina, San Diego, CA, USA) and multiplexing primers according to the manufacturer’s protocol. The cDNA library was constructed with average inserts of 250 bp (150–250 bp) by non-stranded library preparation. The cDNA was purified using a QiaQuick PCR extraction kit (Qiagen, Hilden, Germany). The short cDNA fragments were subjected to end repair, adapter ligation, and agarose gel electrophoresis filtration, and suitable fragments were selected as templates for PCR amplification. Sequencing was performed via a paired-end 125-cycle rapid run on the Illumina HiSeq2500 system.

Read filtering and sequence assembly

High-quality clean reads were obtained by removing adaptor sequences, duplicated sequences, and ambiguous (‘N’) and low-quality reads. Transcriptomes of the two species were separately assembled de novo using the short-reads assembly program Trinity [43]. After assembly, the TIGR Gene Indices clustering tools were used to cluster and remove redundant transcripts [44]. The longest transcripts were considered as unigenes after removing redundancies, and these were subjected to downstream functional annotation and coding sequence (CDS) prediction [45].

Functional annotation

The unigenes of L. mandarinus and L. brandtii were compared using BLASTX against the Nr, KEGG [46], GO [47], the Eukaryotic Orthologous Groups (KOG) [48], and Swiss-Prot [49] databases (E-value ≤1e− 5) to retrieve protein functional annotations based on sequence similarity. Gene names were assigned based on the best BLAST hit. High-priority databases (followed by Nr, Swiss-Prot, and KEGG) were selected to determine the direction of unigene sequences. Sequences showing the best alignment were used to predict the CDSs, and TransDecoder (Find Coding Regions Within Transcripts, was used to identify candidate coding regions within transcript sequences. CDSs were translated into amino sequences using the standard codon Table. GO terms including molecular function, biological process, and cellular component categories were assigned to each sequence using Blast2GO software with an E-value threshold of 1e− 6 for further functional categorization [50]. The distribution of the GO functional classifications of unigenes was plotted using the OmicShare tool ( The unigenes were also aligned to the Cluster of Orthologous Groups/KOG database to predict and classify possible functions. KOBAS 2.0 software [51] ( was used to assign unigenes to KEGG pathway annotations and analyze metabolic pathways.

Identification of DEGs

The clean reads of the L. mandarinus and L. brandtii sequences were mapped with their respective unigenes by RSEM package [52]. The Fragments Per Kilobase of exon per Million mapped fragments (FPKM) was used to eliminate the influence of different gene lengths and sequencing levels on the calculation of gene expression, and FPKM values were directly used to compare gene expression differences between samples. The edgeR package [53] ( was used to obtain the base mean value for identifying DEGs. To correct for multiple testing, the FDR was calculated to adjust the P value threshold. Transcripts with an FDR ≤ 0.05 and a minimum of 2-fold difference in expression (|log2 ratio| ≥ 1) were considered as thresholds for the significance of gene expression differences between two groups. In addition, information for DEGs was collected from unigene annotations, and these genes were subjected to GO and KEGG significant enrichment analyses to identify biological functions and metabolic pathways involving these genes.

Validation of RNA-Seq results by RT-qPCR

To validate the reliability of DEGs identified by RNA-Seq, the mRNA expression levels of seven selected genes were measured by RT-qPCR using the same samples. Primers were designed using Primer-BLAST; the sequences are shown in Additional file 1: Table S1. All primer sets yielded a single peak in the dissociation curves with an amplification efficiency of approximately 1.0. Three technical replicates were prepared for each gene in 96-well plates and amplification was performed on a LightCycler® 480 Instrument II (Roche Diagnostics, Mannheim, Germany). Relative gene expression levels were normalized to that of the internal reference gene β-actin and were calculated with the 2−ΔΔCt method. Correlation analysis in SPSS v.19.0 (SPSS Inc., Chicago, IL, USA) was used to evaluate the concordance between RT-qPCR results and RNA-Seq data. Differences were defined as significant at P < 0.05 and highly significant at P < 0.01 (Welch’s t-test).



Differentially expressed genes


Gene Ontology


Kyoto encyclopedia of genes and genomes

L. brandtii :

Lasiopodomys brandtii

L. mandarinus :

Lasiopodomys mandarinus


Reverse-transcription quantitative PCR


  1. Roper TJ, Bennett NC, Conradt L, Molteno AJ. Environmental conditions in burrows of two species of African mole-rat, Georychus capensis and Cryptomys damarensis. Proceed Zool Soc London. 2010;254(1):101–7.

    Article  Google Scholar 

  2. Qiu Q, Zhang G, Ma T, Qian W, Wang J, Ye Z, Cao C, Hu Q, Kim J, Larkin DM. The yak genome and adaptation to life at high altitude. Nat Genet. 2012;44(8):946–9.

    Article  CAS  PubMed  Google Scholar 

  3. Jackson DC, Ultsch GR. Physiology of hibernation under the ice by turtles and frogs. J Exp Zool Part A Ecol Genet Physiol. 2010;313A(6):311–27.

    Article  CAS  Google Scholar 

  4. Spill F, Reynolds DS, Kamm RD, Zaman MH. Impact of the physical microenvironment on tumor progression and metastasis. Curr Opin Biotechnol. 2016;40:41–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Gilkes DM, Semenza GL, Wirtz D. Hypoxia and the extracellular matrix: drivers of tumour metastasis. Nat Rev Cancer. 2014;14(6):430.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Larson J, Drew KL, Folkow LP, Milton SL, Park TJ. No oxygen? No problem! Intrinsic brain tolerance to hypoxia in vertebrates. J Exp Biol. 2014;217(Pt 7):1024–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Nevo E. Stress, adaptation, and speciation in the evolution of the blind mole rat, Spalax, in Israel. Mol Phylogenet Evol. 2013;66(2):515–25.

    Article  PubMed  Google Scholar 

  8. Šumbera R, Chitaukali WN, Elichov M, Kubov J, Burda H. Microclimatic stability in burrows of an Afrotropical solitary bathyergid rodent, the silvery mole-rat (Heliophobius argenteocinereus). Proc Zool Soc London. 2010;263(4):409–16.

    Article  Google Scholar 

  9. Kim EB, Fang X, Fushan AA, Huang Z, Lobanov AV, Han L, Marino SM, Sun X, Turanov AA, Yang P. Genome sequencing reveals insights into physiology and longevity of the naked mole rat. Nature. 2011;479(7372):223.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Shams I, Avivi A, Nevo E. Hypoxic stress tolerance of the blind subterranean mole rat: expression of erythropoietin and hypoxia-inducible factor 1 alpha. Proc Natl Acad Sci U S A. 2004;101(26):9698–703.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Shams I, Nevo E, Avivi A. Ontogenetic expression of erythropoietin and hypoxia-inducible factor-1 alpha genes in subterranean blind mole rats. FASEB J. 2005;19(2):307–9.

    Article  CAS  PubMed  Google Scholar 

  12. Schmidt HMA, Bicker A, Poetzsch G, Avivi A, Shams I. Hypoxia tolerance, longevity and cancer-resistance in the mole rat Spalax – a liver transcriptomics approach. Sci Rep. 2017;7(1):14348.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Malik A, Korol A, Huebner S, Hernandez AG, Thimmapuram J, Ali S, Glaser F, Paz A, Avivi A, Band M. Transcriptome sequencing of the blind subterranean mole rat, Spalax galili: utility and potential for the discovery of novel evolutionary patterns. PLoS One. 2011;6(8).

  14. Assaf M, Abraham K, Mathias W, Thomas H, Aaron A, Mark B. Transcriptome analysis of the spalax hypoxia survival response includes suppression of apoptosis and tight control of angiogenesis. BMC Genomics. 2012;13(1):615.

    Article  Google Scholar 

  15. Fang X, Nevo E, Han L, Levanon EY, Zhao J, Avivi A, Larkin D, Jiang X, Feranchuk S, Zhu Y. Genome-wide adaptive complexes to underground stresses in blind mole rats Spalax. Nat Commun. 2014;5(3966):3966.

    Article  CAS  PubMed  Google Scholar 

  16. Xiao BLL, Xu C, Zhao S, Lin L, Cheng J, et al. Transcriptome sequencing of the naked mole rat (heterocephalus glaber) and identification of hypoxia tolerance genes. Biol Open. 2017;6(12):1904–12.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Liu B, Wang Z, Lu J. Response to chronic intermittent hypoxia in blood system of mandarin vole (Lasiopodomys mandarinus). Comp Biochem Physiol A Mol Integr Physiol. 2010;156(4):469–74.

    Article  PubMed  Google Scholar 

  18. Zhang Z, Pech R, Davis S, Shi D, Wan X, Zhong W. Extrinsic and intrinsic factors determine the eruptive dynamics of Brandt's voles Microtus brandti in Inner Mongolia, China. Oikos. 2003;100(2):299–310.

    Article  Google Scholar 

  19. Wang HC, Thampatty BP, Lin JS, Im HJ. Mechanoregulation of gene expression in fibroblasts. Gene. 2007;391:1):1–15.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Malik ADV, Han L, Fang XD, Korol A, Avivi A, et al. Genome maintenance and bioenergetics of the long-lived hypoxia-tolerant and cancer-resistant blind mole rat, Spalax: a cross-species analysis of brain transcriptome. Sci Rep. 2016;6:38624.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Drew KL, Buck CL, Barnes BM, Christian SL, Rasley BT, Harris MB. Central nervous system regulation of mammalian hibernation: implications for metabolic suppression and ischemia tolerance. J Neurochem. 2007;102(6):1713–26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Drew KL, Harris MB, LaManna JC, Smith MA, Zhu XW, Ma YL. Hypoxia tolerance in mammalian heterotherms. J Exp Biol. 2004;207(Pt 18):3155–62.

    Article  CAS  PubMed  Google Scholar 

  23. Schnabolk G, Coughlin B, Joseph K, Kunchithapautham K, Bandyopadhyay M, O'Quinn EC, Nowling T, Rohrer B. Local production of the alternative pathway component factor B is sufficient to promote laser-induced choroidal neovascularization. Invest Ophthalmol Vis Sci. 2015;56(3):1850.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Wang J, Li J, Liu Q. Association between platelet activation and fibrinolysis in acute stroke patients. Neurosci Lett. 2005;384(3):305–9.

    Article  CAS  PubMed  Google Scholar 

  25. Wu C, Wen Y, Guo X, Yang T, Shen H, Chen X, Tian Q, Tan L, Deng HW, Zhang F. Genetic association, mRNA and protein expression analysis identify ATG4C as a susceptibility gene for Kashin-Beck disease. Osteoarthritis Cartil. 2017;25(2):281.

    Article  CAS  Google Scholar 

  26. Farfariello V, Amantini C, Santoni G. Transient receptor potential vanilloid 1 activation induces autophagy in thymocytes through ROS-regulated AMPK and Atg4C pathways. J Leukoc Biol. 2012;92(3):421–31.

    Article  CAS  PubMed  Google Scholar 

  27. Avivi A, Brodsky L, Nevo E, Band MR. Differential expression profiling of the blind subterranean mole rat Spalax ehrenbergi superspecies: bioprospecting for hypoxia tolerance. Physiol Genomics. 2006;27(1):54–64.

    Article  CAS  PubMed  Google Scholar 

  28. Mook ORF, Frederiks WM, Van Noorden CJF. The role of gelatinases in colorectal cancer progression and metastasis. Biochimica Biophys Acta-Rev Cancer. 2004;1705(2):69–89.

    Article  CAS  Google Scholar 

  29. Vaughan DE. PAI-1 and atherothrombosis. J Thromb Haemost. 2005;3(8):1879–83.

    Article  CAS  PubMed  Google Scholar 

  30. Cherry TJ, Wang S, Bormuth I, Schwab M, Olson J, Cepko CL. NeuroD factors regulate cell fate and neurite stratification in the developing retina. J Neurosci Official J Soc Neurosci. 2011;31(20):7365–79.

    Article  CAS  Google Scholar 

  31. Randi AM, Laffan MA. Von Willebrand factor and angiogenesis: basic and applied issues. J Thromb Haemost. 2017;15(1):13–20.

    Article  CAS  PubMed  Google Scholar 

  32. Luissint AC, Lutz PG, Calderwood DA, Couraud PO, Bourdoulous S. JAM-L-mediated leukocyte adhesion to endothelial cells is regulated in cis by alpha4beta1 integrin activation. J Cell Biol. 2008;183(6):1159–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Brew N, Azhan A, den Heijer I, Boomgardt M, Davies GI, Nitsos I, Miller SL, Walker AM, Walker DW, Wong FY. Dopamine treatment during acute hypoxia is neuroprotective in the developing sheep brain. Neuroscience. 2016;316:82–93.

    Article  CAS  PubMed  Google Scholar 

  34. Bigham AW. Genetics of human origin and evolution: high-altitude adaptations. Curr Opin Genet Dev. 2016;41:8–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Moore LG, Niermeyer S, Zamudio S. Human adaptation to high altitude: regional and life-cycle perspectives. Am J Phys Anthropol. 1998;(Suppl 27):25–64.

  36. Beall CM, Laskowski D, Erzurum SC. Nitric oxide in adaptation to altitude. Free Radic Biol Med. 2012;52(7):1123–34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Bigham AW, Wilson MJ, Julian CG, Kiyamu M, Vargas E, Leon-Velarde F, Rivera-Chira M, Rodriquez C, Browne VA, Parra E, et al. Andean and Tibetan patterns of adaptation to high altitude. Am J Hum Biol. 2013;25(2):190–7.

    Article  PubMed  Google Scholar 

  38. Alkorta-Aranburu G, Beall CM, Witonsky DB, Gebremedhin A, Pritchard JK, Di Rienzo A. The genetic architecture of adaptations to high altitude in Ethiopia. PLoS Genet. 2012;8(12):e1003110.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Scheinfeldt LB, Soi S, Thompson S, Ranciaro A, Woldemeskel D, Beggs W, Lambert C, Jarvis JP, Abate D, Belay G, et al. Genetic adaptation to high altitude in the Ethiopian highlands. Genome Biol. 2012;13(1):R1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Zapol WM, Liggins GC, Schneider RC, Qvist J, Snider MT, Creasy RK, Hochachka PW. Regional blood flow during simulated diving in the conscious Weddell seal. J Appl Physiol Respir Environ Exerc Physiol. 1979;47(5):968–73.

    CAS  PubMed  Google Scholar 

  41. Kooyman GL, Campbell WB. Heart rates in freely diving Weddell seals, Leptonychotes weddelli. Comp Biochem Physiol A Comp Physiol. 1972;43(1):31–6.

    Article  CAS  PubMed  Google Scholar 

  42. Ilardo MA, Moltke I, Korneliussen TS, Cheng J, Stern AJ, Racimo F, Damgaard PB, Sikora M, Seguin-Orlando A, Rasmussen S, et al. Physiological and genetic adaptations to diving in sea nomads. Cell. 2018;173(3):569.

    Article  CAS  PubMed  Google Scholar 

  43. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Pertea G, Huang X, Liang F, Antonescu V, Sultana R, Karamycheva S, Lee Y, White J, Cheung F, Parvizi B. TIGR gene indices clustering tools (TGICL): a software system for fast clustering of large EST datasets. Bioinformatics. 2003;19(5):651–2.

    Article  CAS  PubMed  Google Scholar 

  45. Wu QBX, Zhao W, Xiang D, Wan Y, Yan J, et al. De novo assembly and analysis of tartary buckwheat (fagopyrum tataricum garetn.) transcriptome discloses key regulators involved in salt-stress response. Genes. 2017;8(10):255.

    Article  PubMed Central  Google Scholar 

  46. Kanehisa M, Goto S, Kawashima S, Okuno Y, Hattori M. The KEGG resource for deciphering the genome. Nucleic Acids Res. 2004;32(Database issue):D277.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25(1):25–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Koonin EV, Fedorova ND, Jackson JD, Jacobs AR, Krylov DM, Makarova KS, Mazumder R, Mekhedov SL, Nikolskaya AN, Rao BS. A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biol. 2004;5(2):R7.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Apweiler R, Bairoch A, Wu CH, Barker WC, Boeckmann B, Ferro S, Gasteiger E, Huang H, Lopez R, Magrane M. UniProt: the universal protein knowledgebase. Nucleic Acids Res. 2004;32(Database issue):D115.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Conesa A, Götz S, Garcíagómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.

    Article  CAS  PubMed  Google Scholar 

  51. Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, Kong L, Gao G, Li CY, Wei L. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011;39(Web Server issue):W316–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12.

  53. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  CAS  PubMed  Google Scholar 

Download references


Aside from funding support, we also thank Pengcheng He, Xiujuan Li and Yimeng Du for their help in data analysis and the feeding of experimental animals.


This work was supported by the National Natural Science Foundation of China (grant no. 31372193).

Availability of data and materials

Raw Illumina sequences were deposited in the National Center for Biotechnology Information (NCBI) database and our sequence read archive (SRA) records will be accessible via the following links after the indicated release date:, SRA accession: SRP156832; Temporary Submission ID: SUB4380898.

Author information

Authors and Affiliations



QD and LS contributed equally to this research. ZW conceived and designed the research. QD and LS wrote the manuscript and analyzed data. QD and YL participated in the design of the work and performed the experiments. MJ, HS, BW, HC, YZ, TS, and YS responsible for animal collection and husbandry. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Zhenlong Wang.

Ethics declarations

Ethics approval and consent to participate

The experimental protocol was approved by the Animal Care and Use Committee of Zhengzhou University and was in accordance with the Guide for the Care and Use of Laboratory Animals of China.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. RT-qPCR primers for validation of RNA-Seq data. Table S2. Illumina sequencing data for analyzed samples. Table S3. Length distribution of assembled unigenes. Table S4. Functional annotation results for L. mandarinus and L. brandtii transcriptomes. Table S6. GO terms significant enriched for up- and downregulated DEGs in L. mandarinus and L. brandtii. Table S7. Genes associated with the GO term “peptidase activity” among up- and downregulated DEGs in L. mandarinus. Table S8. KEGG pathways enriched for up- and downregulated DEGs in L. mandarinus and L. brandtii under acute hypoxia. Figure S1. DEGs in the brain of L. mandarinus and L. brandtii under chronic hypoxia vs. normoxia. FC, fold change; FDR, false discovery rate. Red, blue, and green dots represent up- and downregulated and unchanged genes, respectively. (DOCX 481 kb)

Additional file 2:

Table S5. DEGs for L. mandarinus and L. brandtii. (XLSX 48 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Dong, Q., Shi, L., Li, Y. et al. Differential responses of Lasiopodomys mandarinus and Lasiopodomys brandtii to chronic hypoxia: a cross-species brain transcriptome analysis. BMC Genomics 19, 901 (2018).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: