Hypoxia tolerance, longevity and cancer-resistance in the mole rat Spalax – a liver transcriptomics approach

The blind subterranean mole rat Spalax shows a remarkable tolerance to hypoxia, cancer-resistance and longevity. Unravelling the genomic basis of these adaptations will be important for biomedical applications. RNA-Seq gene expression data were obtained from normoxic and hypoxic Spalax and rat liver tissue. Hypoxic Spalax broadly downregulates genes from major liver function pathways. This energy-saving response is likely a crucial adaptation to low oxygen levels. In contrast, the hypoxia-sensitive rat shows massive upregulation of energy metabolism genes. Candidate genes with plausible connections to the mole rat’s phenotype, such as important key genes related to hypoxia-tolerance, DNA damage repair, tumourigenesis and ageing, are substantially higher expressed in Spalax than in rat. Comparative liver transcriptomics highlights the importance of molecular adaptations at the gene regulatory level in Spalax and pinpoints a variety of starting points for subsequent functional studies.


Additional Figure 1 -Venn diagram of the number of expressed genes
Shown are genes with detectable expression in Spalax (upper circle), rat (left circle) and Heterocephalus (right circle), and the overlap between the species.

Additional Note 1: Validation of the RNA-Seq analysis by alternative bioinformatics tools
To reassess the main results of our study, we selected a controlled set of gene orthologs for the interspecies comparison, thereby addressing specific issues in RNA-Seq analysis 1 . We show that this approach corroborates the main results.
Spalax vs. rat putative orthologous genes (1:1 orthologs only) were determined by Orthofinder 2 using rat Ensembl protein sequences (www.ensembl.org) and Spalax protein sequences (ftp://ftp.ncbi.nlm.nih.gov). We then prepared cross-species genomic annotation data for Spalax vs. rat. This was done as follows: (1) pair-wise alignments of i orthologous transcripts (i=13,000) were built using MAFFT v7 3,4 ; (2) every alignment, a i , was divided into j 25 bp sub-alignments; (3) for each sub-alignment, a ij , with > 70% identity and gaps < 3 bp, the matching genomic regions g ij,species1 and g ij,species2 were retrieved and stored in gtf format; (4) for each gene, RNA-Seq coverage levels in all g i,species1 and g i,species2 genomic sub-regions were calculated using FeatureCounts 5 ; (5) we excluded orthologous genes whose coverage along the gene was poorly correlated between the species (r < 0.5), unless the coverage was consistently higher in one species compared to the other (sign test P-value < 0.001, using R binom-test function); (6) for each gene i, we fitted all j sub-alignments in g i,species1 vs. g i,species2 to a linear regression model (RLM module in R), since we observed that this model correctly predicts reads coverage ratios between samples of the same species. RLM outliers were excluded, unless the sign test was significant (see previous step); (7) two final coordinatesfiles in gtf format were produced, one for each species, after excluding the above-mentioned incomparable genes and gene regions. Coverage plots were visually inspected using the IGV genome viewer and visualization scripts. The pipeline yielded 7,184 Spalax-rat 1:1 orthologous genes. We normalized the gene counts of all samples using EdgeR TMM method 6 and DESeq2 default normalization 7 and inferred differentially expressed genes (parameters: log 2 fold-change > 1.0, adj. p-value < 0.05).

Clustering gene orthologs by species and treatment factors:
Based on the normalized readcount data of all comparable Spalax-rat orthologs over the twelve tested samples, multidimensional scaling (MDS) clustering clearly approved that the principal factors that govern the samples' similarity are the species identity, and the level of O 2 ( Fig. A11-1). As the MDS shows, the species and the O 2 level effects explain about 80% and 10% of the variation, respectively. This reflects the numbers and ratios of differentially expressed genes obtained by the analysis as presented in the main manuscript text.

Fig. AN1-2: Genes significantly higher expressed by at least two fold in Spalax compared to rat within the significantly enriched Fanconi Anemia functional group
A large overlap exists between genes with significant > 2 fold expression differences in both hypoxia (right) and normoxia (left).

Functional enrichment among hypoxia responsive genes in Spalax and rat:
Using the controlled gene set, hypoxia-induced differential gene regulation was inferred separately for Spalax and rat, and significant functional categories of GO terms and KEGG pathways were identified (Supplementary dataset 10). In Spalax, these included the liver insulin resistance pathway ( Fig. A11-3), and the Foxo signalling pathway (Fig. A11-4), again confirming results from the main text analysis. Both pathways can indeed be functionally interconnected, since the Foxo signalling pathway can be regulated by insulin. and Sanger-sequenced to verify the correct inserts (StarSEQ). For all Spalax and rat liver samples, we chose a subset of 11 genes to evaluate relative differential expression by qRT-PCR as previously predicted by RNA-Seq analyses. The tested genes covered represented the categories "cancer", "ageing", "DNA repair" and "hypoxia".

Tab. AN2-1: Primer sequences used for RT-PCR and quantitative RT-PCR for
In the vast majority of experiments (23 out of 25), the direction (ratio) of differential expression between Spalax and rat or between normoxic and hypoxic samples of the same species agreed very well with the expression changes observed by RNA-Seq (Tab. AN2-2). In particular, elevated normoxic transcript levels in Spalax versus rat were displayed as expected by six genes, highly important for the inter-species comparison (A2M, ATR, CISD2, WRN, XPA, TOP2B). Unfortunately, due to shortage in RNA/cDNA availability for Spalax samples, not all of the 11 genes could be compared across the two species. Analysing the hypoxia response, the inducibility of FGF21 and A2M in Spalax, and of CITED2, ATR, VEGFA and especially A2M in rat was also confirmed by qRT-PCR ( Fig. AN2-1). Significance of gene expression differences was tested applying two-sided t-tests in Excel (Microsoft).