Increased longevity due to sexual activity in mole-rats is associated with transcriptional changes in the HPA stress axis
Abstract
Sexual activity and/or reproduction are associated with a doubling of life expectancy in the long-lived rodent genus Fukomys. To investigate the molecular mechanisms underlying this phenomenon, we analyzed 636 RNA-seq samples across 15 tissues. This analysis suggests that changes in the regulation of the hypothalamic–pituitary–adrenal stress axis play a key role regarding the extended life expectancy of reproductive vs. non-reproductive mole-rats. This is substantiated by a corpus of independent evidence. In accordance with previous studies, the up-regulation of the proteasome and so-called ‘anti-aging molecules’, for example, dehydroepiandrosterone, is linked with enhanced lifespan. On the other hand, several of our results are not consistent with knowledge about aging of short-lived model organisms. For example, we found the up-regulation of the insulin-like growth factor 1/growth hormone axis and several other anabolic processes to be compatible with a considerable lifespan prolongation. These contradictions question the extent to which findings from short-lived species can be transferred to longer-lived ones.
Introduction
Most of our current understanding of the underlying mechanisms of aging comes from short-lived model species. It is, however, still largely unclear to what extent insights obtained from short-lived organisms can be transferred to long-lived species, such as humans (Parker et al., 2004; Keller and Jemielity, 2006). Comparative approaches, involving species with particularly long healthy lives and seeking the causative mechanisms that distinguish them from shorter-lived relatives, try to overcome this limitation (Austad, 2009). Many studies that involved organisms with particularly long lifespans, for example, queens in social hymenoptera, birds, bats, African mole-rats, and primates, have produced findings that were not always congruent with established aging theories (Keller and Jemielity, 2006; Austad, 2009; Salmon et al., 2009; Austad, 2011; Dammann, 2017; Bens et al., 2018). Species comparisons, however, also have their limitations. Many observed differences between species with differing lifespans are influenced by phylogenetic constraints, ecophysiological differences, or both, rather than being causal for the species-specific differences in aging and longevity.
Bimodal aging occurs naturally in the genus Fukomys from the rodent family Bathyergidae (African mole-rats). These animals live in families (often called colonies) of usually consisting of 9–16 individuals (Sichilima et al., 2008; Šklíba et al., 2012), although single families may occasionally grow considerably larger in some species (Jarvis and Bennett, 1993; Scharff et al., 2001). Regardless of group size, an established family typically consists of only one breeding pair (the founders of the family, often called king/queen) and their progeny from multiple litters (often called workers). Because of strict avoidance of incest (Burda, 1995), the progeny do not engage in sexual activity in the confines of their natal family, even after reaching sexual maturity. Hence, grown Fukomys families are characterized by a subdivision into breeders (the founder pair) and non-breeders (all other family members). Interestingly, breeders reach the age of 20 years or more in captivity, whereas non-breeders usually die before their tenth birth date (Figure 1A). This divergence of survival probabilities between breeders and non-breeders is found in all Fukomys species studied so far, irrespective of sex. Because no difference in diet or workload has been observed between breeders and non-breeders in captivity, status-specific changes of gene expression after the transition from non-breeder to breeder are considered the most likely explanation of the differing lifespans (Dammann and Burda, 2006; Dammann et al., 2011).

Motivation (A) and principle of the experimental setup (B).
(A) For both species of the Fukomys genus that were examined in this study – Fukomys mechowii and Fukomys micklemi – it was shown that, in captivity, breeders live significantly longer than non-breeders. Lifespan data were redrawn from Dammann and Burda, 2006 and Dammann et al., 2011. (B) Schematic overview of animal treatments. Non-breeders (open shapes) are offspring of the breeder pair of their family (filled shapes) and do not mate with other members of the same family because of incest avoidance in Fukomys (Burda, 1995). Therefore, non-breeders of opposite sexes were taken from different families – labeled as ‘Family A’ (blue) and ‘Family B’ (red) – and permanently housed in a separate terrarium. The two unrelated animals mated with each other, thus producing offspring and becoming breeders of the new ‘Family C’ (green). In addition to the animals that were promoted to be slower-aging breeders, age-matched controls that remained in the faster-aging non-breeders of ‘Family A’ and ‘Family B’ were included in our study – in most cases full siblings (ideally litter mates) of the respective new breeders. After at least 2 years and two pregnancies in ‘Family C’, breeders from ‘Family C’ and their controls from Colonies A and B were put to death, and tissues were sampled for later analysis, which included identification of differentially expressed genes. The shown experimental scheme was conducted with 5–7 (median 6) specimens per cohort (defined by breeding status, sex, species) and 12–15 tissues (median 14) per specimen: in total, 46 animals and 636 samples.
In the wild, non-breeders must meet a member of another family by chance to ascend to breeder status; in captivity, the establishment of new breeder pairs is subject to human control. Allowing an animal to breed in captivity can be regarded as a simple experimental intervention that results in an extension of life expectancy of approximately 100%. This extension is far more than most experimental interventions in vertebrates can achieve, for example, by caloric restriction (e.g., Carmona and Michan, 2016) or diets containing resveratrol or rapamycin (e.g., Valenzano et al., 2006; Johnson et al., 2013). Furthermore, this relative lifespan extension starts from a non-breeder lifespan that is already more than twice as long as that of the mammalian model organisms most widely used in aging research, such as mice or rats.
Until now, relatively few studies have addressed the potential mechanisms behind this natural status dependency of aging in Fukomys sp. Contrary to the predictions of both the advanced glycation formation theory (Morimoto and Cuervo, 2009) and the oxidative stress theory of aging (Harman, 1956; Harman, 2001), markers of protein cross-linking and -oxidation were surprisingly higher in breeders of Ansell's mole-rats (Fukomys anselli) than in age-matched non-breeders (Dammann et al., 2012). On the other hand, in the Damaraland mole-rat (Fukomys damarensis), oxidative damage to proteins and lipids was significantly lower in breeding females than in their non-reproductive counterparts (Schmidt et al., 2014), a finding that is compatible with the oxidative stress theory of aging. In good agreement with their overall longevity irrespective of social status, Ansell's mole-rats produce less thyroxine (T4) and recruit smaller proportions of their total T4 resources into the active unbound form than do euthyroid mammals. Still, nonetheless the levels of unbound T4 (fT4) do not explain the intraspecific longevity differences between F. anselli breeders and non-breeders because the levels of this hormone did not differ between the two cohorts (Henning et al., 2014). Closely connected to the topic of this paper is the finding that non-breeding giant mole-rats (Fukomys mechowii) maintain fairly stable gene expression into relative old ages, quite in contrast to the shorter-lived Norway rat (Sahm et al., 2018a). It is, however, still unclear what happens on the gene expression level when an individual attains breeding status.
Interestingly, a recent study by Bens et al., 2018 found that the longest-lived rodent, the naked mole-rat Heterocephalus glaber, tended globally to show opposite changes in the transition from non-breeders to breeders compared to shorter-lived guinea pigs. Heterocephalus and Fukomys are similar in their mating and social behavior, but differences appear to exist regarding the effect of breeding on longevity: until very recently, naked mole-rat non-breeders have been reported to be as long-lived as breeders (Sherman and Jarvis, 2002; Buffenstein, 2008). In 2018, a lifespan advantage of breeders over non-breeders was reported in females (but not males), yet the divergence of the two groups appeared to be considerably smaller than in Fukomys (Ruby et al., 2018) and underlying data is being debated (Dammann et al., 2019). In summary, status-dependent aging is either absent in Heterocephalus or less pronounced than in Fukomys. Nevertheless, both Fukomys and naked mole-rats can be seen as counterexamples to the disposable soma theory of aging, which assumes that the aging process is linked to an evolutionary tradeoff in which the preservation of the organism is essentially sacrificed to reproductive success (Kirkwood, 1977).
To gain insights into the molecular basis of the bimodal aging pattern in Fukomys, it would be of course interesting to compare sufficiently large cohorts of old breeders and non-breeders both at similar chronological and biological ages. But as old animals, particularly of a long-lived species, are extremely valuable due to obvious logistical, timely, and financial reasons, such an approach was impossible for us. As an alternative, we decided to study young breeders and non-breeders that are much easier to obtain and hypothesized that the status change marks the beginning of a slowdown in the aging process. In this paper, we studied two Fukomys species (F. mechowii and Fukomys micklemi, Figure 1A) by comparing the gene expression profiles of breeders (n = 24) and age-matched non-breeders (n = 22) in 16 organs or their substructures (hereinafter referred to as tissues, Figure 1—figure supplement 1, Figure 1B). Our main aim was to identify genes and pathways whose transcript levels are linked to the status-dependent longevities and to relate these patterns to insights into aging research obtained in shorter-lived species.
Results
We measured gene expression differences between breeders and non-breeders in two African mole-rat species, F. mechowii and F. micklemi. Altogether, we performed RNA-seq for 636 tissue samples covering 16 tissue types from both species, sexes, and reproductive states (breeders and non-breeders). Each of the four groups (male/female breeders/non-breeders) of each species consisted of 5–7 animals (see Supplementary file 1a–e for sample sizes, animal data, and pairing schemes). For each tissue, we conducted a multifactorial analysis of differentially expressed genes (DEGs): the analysis was based on the variables reproductive state, sex, and species. During this exercise, we focused on the differences between slower-aging breeders and faster-aging non-breeders. This approach increases our statistical power by giving us a fourfold increase of sample size in comparison to species- and sex-specific breeder vs. non-breeder analyses. At the same time, we can additionally reduce the number of false-positive DEGs by restricting the analysis to those breeding status-related genes that show the same direction in both sexes and species. We deliberately focused on those genes to concentrate our study on universal mechanisms that hold for both sexes and species. The number of samples collected and processed for this work required splitting them into different batches for sequencing. To analyze the degree this procedure potentially biased our results, we systematically investigated potential batch effects. Our results suggest that the sequencing strategy has little effect on the outcomes reported in the following (see Control analyses).
To globally quantify the transcriptomic differences between the reproductive states, we performed three analyses: clustering of the samples based on pairwise correlation, principal variant component analysis, and an overview of the number of DEGs between reproductive states in comparison to DEGs between species and sex. Clustering of the samples based on pairwise correlations showed a full separation of the two species at the highest cluster level (Figure 2—figure supplement 1). Below that level, an almost complete separation according to tissues was observed. Within the tissue clusters, the samples did not show a clear-cut separation between sex or breeder/non-breeder status. Accordingly, a principal variance component analysis showed that species, tissue, and the combination of both variables accounted for 98.4% of the total variance in the data set; individual differences explained 1.4% of the variance, and only 0.004% was explained by breeder/non-breeder status (Figure 2A). Regarding the numbers of DEGs, we found – unsurprisingly considering the aforementioned facts – by far the highest number of DEGs in the species comparison (Figure 2B). Although in almost every examined tissue the numbers of detected DEGs were also high between sexes, most tissues exhibited very few DEGs due to breeder/non-breeder status. Exceptions were liver, spleen, ovary, and, especially, tissues of the endocrine system (adrenal gland, pituitary gland, thyroid), in which the number of DEGs between breeders and non-breeders ranged from more than 60 to several thousand.

Total variance distribution (A) and numbers of differentially expressed genes (B).
(A) Relative contribution of the model factors (breeding status, sex, species, tissue) and their combinations (:) to the total variance in the examined data set. The relative contributions were determined by principal variance component analysis. (B) Numbers of identified differentially expressed genes per tissue and model factor (first column, species; second, sex; third, status). Stacked bars indicate the proportions of up- and down-regulated genes (red and green, respectively; directions: F. mechowii vs. F. micklemi, female vs. male, breeder vs. non-breeder).
Since the change from non-breeder to breeder status apparently marks the beginning of a slowdown in the aging process, we first wanted to find out whether and where there are intersections of reproductive status DEGs with those whose expression level is known to change during aging. Therefore, we determined overlaps by using the Digital Aging Atlas (DAA) – a database of genes that show aging-related changes in humans (Craig et al., 2015). Across species and sexes, significant overlaps (false discovery rate [FDR] < 0.05, Fisher’s exact test) with the DAA were found in three tissues: adrenal gland, ovary, and pituitary gland (FDR=0.005, each; Figure 3A). Among these three endocrine tissues, the DEGs of the ovaries overlapped significantly with those from adrenal (p=2.8*10−27) and pituitary glands (p=0.005), but there was no significant overlap between the two glands (Figure 3A). Thus, together, we found indications for aging-relevant expression changes after the transition from non-breeders to breeders in three tissues of the endocrine system, which presumably affect separate aspects of aging in adrenal and pituitary glands.

Assessment of the aging relevance of genes that are differentially expressed between breeders and non-breeders.
(A) For each tissue, we separately tested whether the identified differentially expressed genes between status groups significantly overlapped with the genes within the Digital Aging Atlas database (Fisher’s exact test, false discovery rate [FDR] < 0.05). Significant overlaps were found for three tissues: adrenal gland, ovary, and pituitary gland. The Venn diagram depicts the overlaps of these three tissues with the Digital Aging Atlas and with each other (**FDR < 0.01; ***FDR < 0.001). (B) A similar experiment comparing the transcriptomes of breeders vs. non-breeders was recently conducted in naked mole-rats (NMRs) and guinea pig (GPs) (Bens et al., 2018). For NMR, there is also evidence that breeders have a (slightly) longer lifespan than non-breeders, whereas for GP the opposite is assumed (Bens et al., 2018; Ruby et al., 2018). Across 10 tissues that were examined in both studies, the analysis determined whether status-dependent differentially expressed genes identified in the current study were regulated in the same or opposite direction in NMR and GP (Supplementary file 1f). The listed p-values (two-sided binomial exact test; hypothesized probability, 0.5) describe how extremely the ratio of genes expressed in the same and opposite directions deviates from a 50:50 ratio. Green and orange indicate the majority and minority of genes within a comparison, respectively. Figure created with BioRender.com.
Moreover, we compared the DEGs with respect to the reproductive status that we identified in Fukomys with regard to their direction to transcript-level changes observed in similar experiments using naked mole-rats and guinea pigs (Bens et al., 2018). The direction of the status-dependent DEGs regulation in Fukomys, as found in this study, was significantly more often the same rather than opposite as in the naked mole-rat (females, 60%, p=5.2*10−58; males, 62%, p=10−44 for females, Figure 3B, Supplementary file 1f). In the guinea pig, on the contrary, the Fukomys reproductive status DEGs were significantly more often regulated in the opposite direction (females, 57%, p=9*10−25; males, 59%, p=4.7*10−23, Figure 3B, Supplementary file 1f). Thus, at the single-gene level, the expression changes linked to reproductive status may affect lifespan differently in long-lived African mole-rats than in shorter-lived guinea pigs.
Beyond the single-gene level, we aimed to identify metabolic pathways and biological functions whose gene expression significantly depends on reproductive status. For this, we used Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Kanehisa et al., 2017) and Molecular Signatures Database (MSigDB) hallmarks (Liberzon et al., 2015) as concise knowledge bases. As standard approach, we used an established method combining all p-values of genes in a given pathway in a threshold-free manner (Figure 4). The advantage of this approach is that it bundles the p-values from test results of individual gene expression differences at the level of pathways (see Materials and methods for details). In addition, we applied a second enrichment method that aggregates fold-changes instead (Figure 4—figure supplement 11; for a comparison of the results of both approaches, see Control analyses). Altogether, the gene expression of 55 KEGG pathways and 41 MSigDB hallmarks was significantly affected by reproductive status in at least one tissue (Figure 4—figure supplements 1 and 2). Because the individual interpretation of each of these pathways/hallmarks would go beyond the scope of this study, we focus here on those 14 pathways and 13 hallmarks that were significantly different between breeders and non-breeders (FDR < 0.1) in a global analysis across all tissues (Figure 4). Because many pathways are driven mainly by gene expression in subsets of tissues, we weighted gene-wise the differential expression signals from the various tissues by the respective expression levels in the tissues. For instance, the expression level of the growth hormone (GH) gene GH1 is known to be almost exclusively expressed in the pituitary glands. In our data set, the GH1 level of the pituitary gland accounted for 99.96% of the total GH1 across all tissues. Accordingly, in pathways that contain GH, our weighted cross-tissue differential expression signal for this gene is almost exclusively determined by the pituitary gland. On the contrary, a differential expression signal of this gene in another tissue with a very low fraction of the gene’s total expression would have almost no impact on the weighted cross-tissue level – even if that signal were very strong (see Materials and methods for details).

Pathways and metabolic functions enriched for status-dependent differential gene expression.
Shown are all Kyoto Encyclopedia of Genes and Genomes pathways (A) and Molecular Signatures Database hallmarks (B) that are enriched for differential gene expression between breeders and non-breeders at the weighted cross-tissue level (false discovery rate [FDR], <0.1). This enrichment analysis was carried out threshold-free, which is why tissues without differentially expressed individual genes (see Figure 2B) can also show differentially expressed pathways. The numbers within the cells give the FDR, that is, the multiple testing corrected p-value, for the respective pathway/hallmark and tissue. As indicated by the color key, red and green indicate up- or down-regulated in breeders, respectively. White indicates a pathway/hallmark that is significantly affected by differential expression and whose signals for up- and down-regulation are approximately balanced. Dark colors up to black mean that there is little or no evidence that the corresponding pathway/hallmark is affected by differential gene expression. Figure 4—figure supplements 1 and 2 provide detailed overviews of all pathways/hallmarks that are enriched in at least one tissue for status-dependent differential expression signals.
We found strong indications for increase in the activity of certain anabolic functions in breeders: ribosomal protein expression (hsa03010 Ribosome, hsa03008 Ribosome biogenesis in eukaryotes, Figure 4A) was elevated in most tissues (and accordingly also in the weighted cross-tissue analysis). In MSigDB hallmarks, the strongest enrichment signal came from MYC targets (HALLMARK_MYC_TARGETS_V1), which can largely be considered a reflection of enhanced ribosomal protein expression and the fact that MYC is a basal transcription factor up-regulating genes involved in protein translation (Hofmann et al., 2015; Figure 4B). In functional correspondence, we observed an increase in the expression of mitochondrial respiratory chain components (hsa00190 Oxidative phosphorylation, HALLMARK_OXIDATIVE_PHOSPHORYLATION, Figure 4A, B). We also found strong indications for increased protein degradation (hsa03050 PROTEASOME, Figure 4A). This weighted cross-tissue signal was, in contrast to the situation regarding ribosomes, driven mainly by two tissue types: the adrenal gland and the gonads.
To examine whether the simultaneous up-regulation of the ribosome, proteasome, and oxidative phosphorylation is a coordinated regulation, we performed a weighted gene co-expression network analysis (WGCNA) (Langfelder and Horvath, 2008) from our gene count data and examined the connectivity between pairs of those KEGG pathways flagged in the weighted cross-tissue analysis. We found that the expression of ribosomal genes (hsa03010) was significantly linked to those of ribosome biogenesis (hsa03008, FDR = 4.59*10−3), oxidative phosphorylation (hsa00190, FDR = 4.05*10−4), and proteasome (hsa03050, FDR = 4.59*10−3), whereas no other examined pathway pair exhibited a significant connectivity (Figure 4—figure supplement 3). Interestingly, ribosome, proteasome, and oxidative phosphorylation pathways also shared other characteristics of their differential expression signals: subtle fold-changes, that is, up-regulation of 3–9% on average. Thus, statistically significant signals at the pathway level resulted from relatively small shifts in all genes of these pathways in a seemingly coordinated manner and across multiple tissues (Source data 1a). In addition, ribosome (hsa03010, in ovary), proteasome (hsa03050, in ovary and adrenal gland), and RNA-transport (hsa03013, in adrenal gland) are enriched in those Fukomys status-dependent DEGs that show, in a similar experimental setting (Bens et al., 2018), the same direction in both naked mole-rat sexes and the opposite direction in both guinea pig sexes (FDR < 0.05, Fisher’s exact test).
The myogenesis hallmark (Figure 4B) was also found to be up-regulated in breeders. Expectedly, this weighted cross-tissue result was driven mainly by differential expression signals from muscle tissue: muscle from all tissues exhibited the lowest p-value (Figure 4A), and 15 of 20 up-regulated genes that contributed most to the weighted cross-tissue differential myogenesis signal exhibited their highest expression in muscle. These genes were involved mainly in calcium transport or part of the fast-skeletal muscle-troponin complex (Source data 1a). A clear exception of this muscle-dominated expression is found in the gene that exhibited the highest relative contribution to the differential myogenesis signal, insulin-like growth factor 1 (IGF1). This gene was found to be expressed most strongly in ovary and liver and was strongly up-regulated in the breeders' ovaries and adrenal glands (Table 1). IGF1 codes for a well-known key regulator of anabolic effects such as cell proliferation, myogenesis, and protein synthesis (Schiaffino and Mammucari, 2011; Jung and Suh, 2014) and has a tight functional relation to GH (gene: GH1) another key anabolic regulator upstream of IGF1; together, these factors form the so-called GH/IGF1 axis (Cannata et al., 2010; Junnila et al., 2013; Bodart et al., 2015; Raisingani et al., 2017; Carotti et al., 2018; Lozier et al., 2018). Also, GH1 was strongly up-regulated in breeders in its known principal place of synthesis, the pituitary gland (Table 1).
Top 10 genes regarding weighted cross-tissue differential expression signal.
Weighted cross-tissue | Tissue with highest expression | ||||||
---|---|---|---|---|---|---|---|
p-Value | FDR | log2-foldchange* | Tissue | % of cross-tissue expression | FDR | log2-foldchange* | |
SULT2A1 | 0.00E+00 | 0.000 | −3.19 | Liver | 97.42 | 5.27E-06 | -3.26 |
MC2R | 6.00E-06 | 0.046 | −0.53 | Adrenal gland | 89.53 | 8.78E-05 | -0.56 |
INHA | 1.60E-05 | 0.081 | 1.57 | Ovary | 67.91 | 1.54E-05 | 2.22 |
INHA – tissue with second highest expression | Testis | 27.61 | 0.56 | 0.21 | |||
CYP11A1 | 5.90E-05 | 0.224 | 0.61 | Ovary | 45.55 | 5.49E-05 | 0.91 |
CYP11A1 – tissue with second highest expression | Adrenal gland | 43.63 | 1.48E-03 | 0.47 | |||
NLRP14 | 1.17E-04 | 0.355 | −0.80 | Ovary | 68.98 | 2.84E-05 | −1.25 |
NLRP14 – tissue with second highest expression | Testis | 18.67 | 0.54 | 0.19 | |||
GH | 2.58E-04 | 0.618 | 0.45 | Pituitary gland | 99.96 | 1.99E-02 | 0.45 |
IGF1 | 3.62E-04 | 0.618 | 0.59 | Ovary | 51.73 | 1.20E-04 | 0.76 |
IGF1 – tissue with second highest expression | Liver | 36.28 | 0.39 | 0.24 | |||
ZP4 | 3.74E-04 | 0.618 | −1.16 | Ovary | 95.30 | 1.99E-03 | −1.23 |
TCL1A | 5.33E-04 | 0.618 | −1.05 | Ovary | 90.12 | 2.14E-03 | −1.18 |
PNLIPRP2 | 5.69E-04 | 0.618 | 1.36 | Ovary | 76.08 | 2.23E-03 | 1.65 |
-
*Direction: breeder/non-breeder.
Tissues are listed if they contribute at least 10% to cross-tissue expression.
-
FDR: false discovery rate.
With xenobiotic metabolism and TNF-α signaling, two defense hallmarks were also found to be up-regulated in breeders by the weighted cross-tissue analysis (HALLMARK_XENOBIOTIC_METABLISM and HALLMARK_TNFA_SIGNALING_VIA_NFKB, Figure 4B). The up-regulation of the reactive oxygen species (ROS) hallmark comprising genes coding for proteins that detoxify ROS (HALLMARK_REACTIVE_OXYGEN_SPECIES_PATHWAY, Figure 4B) falls into a similar category.
Another interesting aspect that was found to be significantly altered in breeders is steroid hormone biosynthesis (hsa00140 Steroid hormone biosynthesis, Figure 4A). In this case, both up- and down-regulated genes were in the pathway, and their absolute fold-changes were roughly balanced. Steroid hormones, on the one hand, comprise sex steroids – because these hormones are important players in sexual reproduction, such differences should be expected given the experimental setup. On the other hand, the class of steroid hormones – corticosteroids – has regulatory functions in metabolism, growth, and the cardiovascular system, as well as in the calibration of the immune system and response to stress (Liu et al., 2013). The most influential contributor to the differential pathway signal by far on the weighted cross-tissue level was CYP11A1, which codes for the (single) enzyme that converts cholesterol to pregnenolone. This is the first and rate-limiting step in steroid hormone synthesis (Miller and Auchus, 2011). Because CYP11A1 was found to be up-regulated in breeders in its main places of synthesis – the gonads and the adrenal gland (Table 1) – it can be assumed that the total output of steroid hormone biosynthesis in breeders is increased. The pattern of up- and down-regulation on the KEGG pathway, however, suggests that sex steroids especially are produced at a higher rate in breeders, whereas circulating levels of glucocorticoids – such as cortisol – should be lower than in non-breeders (Figure 4—figure supplement 4; see also our discussion on ACTH-R below).
Finally, several pathways flagged by the weighted cross-tissue analysis seem to be derivatives of the abovementioned differentially expressed pathways instead of representing altered functions on their own. For example, Huntington’s (hsa05016), Parkinson’s (hsa05012), and Alzheimer’s (hsa05010) diseases could, in principle, be interpreted as highly relevant for aging and lifespan. A closer inspection of these pathways reveals, however, that the genes of the mitochondrial respiratory chain – which is the core of the oxidative phosphorylation pathway – are in all three cases the main contributors to the respective differential expression signals (Figure 4—figure supplements 5–7, Source data 1a). Similarly, we see in the case of fat digestion (hsa04975 Fat digestion and absorption) that two of the three largest contributors to the differential expression signal of that pathway – ABCG8 and SCARB1 – are directly involved in the transport of cholesterol (Liu et al., 2008; Wang et al., 2015). Therefore, it seems likely that this signal is an expression of the altered steroid hormone biosynthesis rather than indicating altered fat digestion.
Interestingly, three genes, which we had already mentioned as potential regulators during the pathway analysis, also appeared among the 10 most clearly altered genes on the weighted cross-tissue level: GH1, IGF1, and CYP11A1 (Table 1). The top 2 among these 10 are sulfotransferase family 2A member 1 (SULT2A1) and melanocortin two receptor (MC2R). SULT2A1 is the main catalyzer of the sulfonation of the steroid hormone dehydroepiandrosterone (DHEA) to its non-active form DHEA-S (Hammer et al., 2005). DHEA has repeatedly been proposed to be an ‘anti-aging hormone’ because its levels are negatively associated with chronological aging (Baulieu, 1996; Celec and Stárka, 2003; Rutkowski et al., 2014). We found that SULT2A1 is strongly down-regulated in breeders’ liver, which is also the main location of its enzymatic action. The second candidate, MC2R, encodes the adrenocorticotropin hormone (ACTH) receptor, which is the main inducer of glucocorticoid synthesis and a crucial component of the hypothalamic–pituitary–adrenal (HPA) axis (Walker et al., 2015). In humans and many other mammals, prolonged glucocorticoid excess leads to Cushing’s syndrome. Affected individuals exhibit muscle weakness, immune suppression, impairment of the GH/IGF1 axis, higher risk of diabetes, cardiovascular disease (hypertension), osteoporosis, decreased fertility, depression, and weight gain (Chabre, 2014; Ferraù and Korbonits, 2015). The large overlap of these symptoms with those of aging could explain to some extent that Cushing’s syndrome patients exhibit considerably higher mortality rates (Etxabe and Vazquez, 1994). We hypothesized that the increased expression of the ACTH receptor in Fukomys non-breeders can cause similar expression patterns and consequences (Figure 5).

Model of the stress axis as a key mechanism for status-dependent lifespan differences in Fukomys.
From a wide range of mammals, including humans (Ferraù and Korbonits, 2015), dogs (Kooistra and Galac, 2012), horses (McCue, 2002), cats (Meijer et al., 1978), and guinea pigs (Zeugswetter et al., 2007), it is known that chronic glucocorticoid excess leads to a number of pathological symptoms that largely overlap with those of aging and result in considerably higher mortality rates for affected individuals (Etxabe and Vazquez, 1994, #808). The most common cause of chronic glucocorticoid excess is excessive secretion of the adrenocorticotropic hormone (ACTH) by the pituitary gland. ACTH is transported via the blood to the adrenal cortex where it binds to the ACTH-receptor (ACTHR; encoded by the gene MC2R), which induces the production of glucocorticoids, especially cortisol. Glucocorticoids are transported to the various tissues, where they exert their effect by activating the glucocorticoid receptor (NR3C1) that acts as a transcription factor and regulates hundreds of genes. The constant overuse of this transcriptional pattern eventually leads to the listed symptoms. Our hypothesis is that the permanent, high expression of the ACTH-receptor in Fukomys non-breeders causes effects similar to those known from overproduction of the hormone. In line with this hypothesis, (i) cortisol levels are increased in non-breeders and (ii) target genes of the glucocorticoid receptor are highly enriched for status-dependent differential gene expression. Furthermore, the animals were examined for common symptoms of chronic glucocorticoid excess: (iii) non-breeders gained more weight during the experiment than breeders, (iv) exhibited lower bone density at the end of the experiment, and (v) lower gene expression in the growth hormone/insulin-like growth factor 1 axis than breeders.
We subjected this hypothesis (Figure 5) to a first test by matching the global fold-changes of the breeder/non-breeder comparison with those of a gene expression comparison of human controls/Cushing patients (Hochberg et al., 2015) and found a significant correlation (R = 0.11, p=6.6*10−36). Since the target tissue of Hochberg et al., subcutaneous adipose tissue was unfortunately not included in our study, we used our cross-tissue data for comparison. If we use instead our data from skin as adjacent tissue, a similar correlation results (R = 0.11, p=1.1*10−40).
We tested the hypothesis further by checking five of its key predictions. Altered MC2R expression (Figure 6A) coincides with higher cortisol levels in hair samples from non-breeding F. mechowii than in those from breeders of the same species (Begall et al., 2021). Furthermore, glucocorticoids such as cortisol exert their effect by binding to the glucocorticoid receptor that, in turn, acts as a transcription factor for many genes (Gjerstad et al., 2018). We tested whether the expression of targets of the glucocorticoid receptor (NR3C1) was significantly altered throughout our data using two gene lists (Phuc Le et al., 2005): about 300 direct target genes of the receptor that were identified by chromatin immunoprecipitation (i) and about 1300 genes that were found to be differentially expressed depending on the presence or absence of exogenous glucocorticoid (ii). Both gene lists were found to be significantly affected by differential expression at the weighted cross-tissue level (i, p=0.001; ii, p<10−9) as well as in five (i) and eight (ii) single tissues (Supplementary file 1g, h). In line with our hypothesis, we observed that the weight gain in non-breeders was, on average, twice as strong compared to the weight gain in breeders during the experiment (p=7.49*10−3, type II ANOVA, Figure 6B). In addition, we found a subtle but significant influence of reproductive status on the density of the vertebrae: the vertebrae of breeders were slightly denser than those of age-matched non-breeders (p=0.03 for vertebra T12 only, and p=0.01 across all examined vertebrae L1, L2, and T12; ANOVA, Figure 6C).

MC2R gene expression and physiological measurements.
(A) Gene expression of MC2R, coding for the adrenocorticotropic hormone receptor, in breeders and non-breeders of Fukomys mechowii and Fukomys micklemi. (B) Weight gain of the animals during the experiment. (C) Measured optical densities of vertebra T12 of F. micklemi breeders and non-breeders. Red, breeders; blue, non-breeders; filled, F. mechowii; unfilled, F. micklemi; circles, females; squares, males; dashed line, median. Statistically significant differences between breeders and non-breeders were determined with (A) DESeq2 (Love et al., 2014) and (B, C) analysis of variance with status, species, sex, and age as independent variables (see Materials and methods).
Discussion
The vast lifespan differences between Fukomys breeders and non-breeders are, according to our RNA-seq data, associated with only subtle global pattern shifts in transcript levels. Concerning the tested explanatory variables (Fukomys species, sex, breeding status), we found by far the highest number of DEGs at the level of the species comparison. Although the number of DEGs between the sexes was comparably high in almost each examined tissue, only very few DEGs were found comparing breeders and non-breeders. One exception is the ovary, whose high number of DEGs corresponds well with the disparity in reproductive activity. Other exceptions are liver, spleen, and, especially, the tissues of the endocrine system (adrenal gland, pituitary gland, thyroid), in which the number of DEGs between breeders and non-breeders ranged from more than 100 to more than 2500.
Changes in the gene expression of the endocrine system are well known to play an important role both in sexual maturation and in aging and the development of aging-associated diseases, for example, diabetes and cardiovascular diseases (Tatar et al., 2003; Chahal and Drake, 2007; Jones and Boelaert, 2015). This finding fits well with the observation of substantial changes in the endocrine system after the transition from non-breeders to breeders in the related naked mole-rat – one of the key results of a recent study in that species (Bens et al., 2018). The dominance of differential expression in endocrine tissue is also plausible insofar as these tissues exert a strong control function for other tissues via hormone release.
Steroid hormone biosynthesis exhibits a bipartite pattern in breeders, with up-regulated sex steroid genes and a simultaneous down-regulation of corticosteroid synthesis genes. The former could be expected as a consequence of sexual activity in breeders. For aging, it is, however, interesting that SULT2A1, a gene that codes for the specialized sulfotransferase converting the sex steroid DHEA to its non-active form DHEA-S, was found to be heavily down-regulated in breeders. DHEA is the most abundant steroid hormone; it serves as a precursor for sex steroid biosynthesis but also has various metabolic functions on its own (Allolio and Arlt, 2002; Webb et al., 2006). DHEA levels decrease continuously during the human aging process to an extent that favors it as an aging biomarker (Maggio, 2007; Traish et al., 2011). As a result, an aggressive advertising of DHEA as ‘anti-aging hormone’ in the form of dietary supplements could be observed in recent years (Baulieu, 1996; Celec and Stárka, 2003; Webb et al., 2006; Rutkowski et al., 2014). Despite conflicting experimental data from various animal studies and clinical trials, a positive effect of DHEA on human health is frequently considered to be likely in the literature. Large-scale and long-term studies, however, are still pending (Traish et al., 2011; Rutkowski et al., 2014; Samaras et al., 2015). Interestingly, DHEA is not only described as an aging marker but also as a marker for clinically relevant glucocorticoid excess (Burkhardt et al., 2013).
Regarding corticosteroid synthesis, we hypothesize that the regulation of adrenal gland MC2R, coding for the ACTH receptor, causes a considerable proportion of the overall observed expression patterns and of the lifespan extension in breeders (Figure 5). As a critical component of the HPA axis, ACTH is a stress hormone that is produced by the pituitary gland and transported by the blood to the adrenal cortex, where it binds to the ACTH receptor (Fridmanis et al., 2017). Subsequently, the ACTH receptor induces the synthesis of glucocorticoids (Walker et al., 2015), for example, cortisol, which in turn cause immunosuppressive and various metabolic effects throughout the organism (Becker, 2013). In many mammals (such as humans, dogs, and guinea pigs) glucocorticoid excess disorder, Cushing’s syndrome, is caused by overproduction of the hormone ACTH. Our results for Fukomys suggest that the increased levels of the ACTH receptor may lead to symptoms and expression patterns that could resemble some of molecular and phenotypic aspects of this pathological condition (Figure 5).
We looked for evidence of a regulation upstream of the ACTH receptor. Given the known positive feedback loop between the ACTH receptor and its own ligand, decreased ACTH synthesis in breeders would have been an obvious explanation (Imai et al., 2001). We found that the expression of the ACTH polypeptide precursor gene (POMC) is unchanged. However, also other post-transcriptional or post-translational mechanisms such as cleaving may influence the ACTH levels in breeders and non-breeders. Of those genes known to be involved in the regulation of MC2R (Beuschlein et al., 2001; Lin et al., 2007), we found that only PRKAR1B was differentially regulated in the adrenal gland. However, this gene codes for only one of seven subunits of the involved protein kinase A. Alternative explanations could be that epigenetic modifications, other still unidentified regulators of the transcript level, or both are responsible for the differential expression of MC2R.
While the opportunities to pursue the hypothesis upstream of the ACTH receptor were limited with our data set, many of the predicted downstream effects (Figure 5) could be verified. One of the clearest symptoms of Cushing's syndrome known in many animal species is an increase in body fat and a consequent significant weight gain (Ferraù and Korbonits, 2015). Our recording of body weights at the beginning and end of the experiment showed that non-breeders gained twice as much weight on average as breeders (p=7.49*10−3, type II ANOVA, Figure 6B). However, it should be noted that our experimental design does not allow to distinguish between increased body fat content and other possible reasons (i.e., overall growth, muscle gain) for these differences in weight gain. Another known complication associated with glucocorticoid excess is osteoporosis (Lodish et al., 2010). A computed tomography measurement of three vertebrae (L1, L2, and T12), in accordance with this, showed significantly lower bone densities of non-breeders (p=0.01, ANOVA, Figure 6C). The absolute differences, however, are probably too small to assume that the non-breeders are already physically impaired by that. It seems possible that a longer experiment duration than the 1–2 years we conducted would have aggravated this difference, but this remains speculative. Directly from our gene expression analysis we know that the GH/IGF1 axis of the non-breeders is impaired, which is another well-known symptom of chronic glucocorticoid excess (Ferraù and Korbonits, 2015).
In addition, we have systematically compared our observed and expected – based on the Cushing’s syndrome hypothesis – expression patterns. A key element in this context is the glucocorticoid receptor (NR3C1), which, when activated by glucocorticoids such as cortisol, is considered to mediate their effects, including those of Cushing's syndrome in the case of chronic activation (Gjerstad et al., 2018). According to our hypothesis, the known, experimentally determined NR3C1 target genes (Phuc Le et al., 2005) should therefore show a clear pattern of expression differences between breeders and non-breeders. This was the case for both the 1300 indirect (p<10−9) and the 300 direct target genes (p=0.001). Furthermore, according to our hypothesis, the expression changes in the comparison of control individuals against patients with Cushing's syndrome should correspond to those of the comparison of breeders vs. non-breeders. A global comparison of fold-changes from human subcutaneous adipose tissue (Hochberg et al., 2015) against our cross-tissue data showed a significant correlation (R = 0.11, p=6.6*10−36). A very similar result was obtained by comparing the human data with our data from skin as the tissue closest to subcutaneous adipose tissue that was examined in our study (R = 0.11, p=1.1*10−40). Despite the clear significance of these correlations, the reported values of the correlation coefficients could be considered to be only modest. It should, however, be taken into account that two evolutionarily widely separated species and not exactly matching tissues were compared.
For further validation, we applied the combination of limma (Smyth, 2004) and voom (Law et al., 2014) as an alternative DEG detection method to DESeq2 to account for potential dispersion problems in our data. All single DEGs mentioned in the paper were confirmed by this validation. This specifically includes all genes supporting our HPA hypothesis, for example, SULT2A1 (liver), MC2R (adrenal gland), CYP11A1 (adrenal gland and ovary), GH (pituitary gland), and IGF1 (ovary). In a similar way, we checked the main findings obtained with our standard bioinformatic p-value-based enrichment approach (Fridley et al., 2010; Evangelou et al., 2012; Poole et al., 2016), using a second method that aggregates fold-changes instead. In particular, this step confirmed our initial finding of an enrichment of genes involved in the steroid hormone biosynthesis and ribosome metabolism as a central part of up-regulated anabolism in breeders, consistent with our HPA hypothesis (Figure 4—figure supplement 11) (see Control analyses).
In view of the different evidence presented, we consider it justified to assume that the transcriptional, metabolic and physiological condition of Fukomys non-breeders resembles that of chronic glucocorticoid excess described from other species. Since this condition shortens life in other species (Etxabe and Vazquez, 1994) and several of its symptoms also occur during aging, we hypothesize that it may shorten the lifespans in Fukomys non-breeders as well. The ultimate proof and the quantification of its impact on Fukomys lifespan, however, will be difficult in such a long-lived species, which cannot yet be genetically manipulated. Nevertheless, a recent study confirmed one of the most important predictions of the hypothesis: hair cortisol levels of non-breeders are on average about twice as high as those of breeders. However, a subgroup of non-breeders, defined by social environmental factor (one or both parents had died), exhibited cortisol levels nearly as low as in breeders (Begall et al., 2021). The drop of cortisol levels in orphaned non-breeders is accompanied by a substantial drop of mortality rates in this subgroup approaching the values of breeders (Figure 6—figure supplement 1). This means that across status groups (breeders, orphaned non-breeders, and non-breeders from intact families) mortality rates are clearly correlated with long-term cortisol levels. Taken together, all evidence suggests that the presence of breeders causes distress in non-breeders, resulting in constantly activating the HPA axis and contributing to reduced life expectancy.
In the circulatory system, represented by blood and spleen, down-regulation of coagulation factors was observed in breeders. Because coagulation factors are known to be up-regulated during aging in humans mice, rats, and even fish, this can be interpreted as a sign of a more juvenile breeders’ transcriptome (Ochi et al., 2016; Benayoun et al., 2019). Furthermore, the activity of coagulation factors is associated with a higher risk of coronary heart disease (Lowe and Rumley, 2014). However, we found no obvious histopathological lesions in the hearts or other organs (spleen, kidney, liver, lung) of non-breeders. Therefore, if down-regulation of coagulation attenuates the aging process in Fukomys, it seems to exert its effect only, if at all, at a latter age not investigated here.
Two defense mechanisms were also found to be up-regulated in breeders: xenobiotic metabolism and TNF-α signaling. Increased TNF-α signaling often leads to the induction of apoptosis (Annibaldi and Meier, 2018). In line with this, we found that apoptosis and P53 signaling were also up-regulated in breeders. Apoptosis is considered an important anticancer mechanism (Baig et al., 2016; Pistritto et al., 2016). We hypothesize this to compensate cancerogenic effects of the anabolic alterations described above (especially the up-regulation of the GH/IGF1 axis). In line with this, our more than 30 years’ breeding history with several Fukomys species in Germany and the Czech Republic suggests that breeders are as ‘cancer-proof’ as non-breeders despite their much longer lifespan (own unpublished data and R. Šumbera, personal communication).
Many anabolic pathways are up-regulated in breeders across tissues: protein biosynthesis, myogenesis, and the GH/IGF1 axis. In line with this finding and with the fact that protein synthesis consumes 30–40% of a cell’s ATP budget (Hands et al., 2009), we observed increased expression of mitochondrial respiratory chain components, implying an increase in the capacity for cellular ATP production. On the other hand, protein degradation and clearance by the proteasome are also up-regulated in breeders. The fact that the expression of proteasomal genes is significantly linked to the genes of ribosome biogenesis and oxidative phosphorylation indicates that those processes influence, or even trigger, each other and hence are regulated in a coordinated manner in Fukomys mole-rats.
An important positive regulator of ribosome biogenesis and protein synthesis is mTOR (Johnson et al., 2013). Based on this, one could expect an up-regulation of the corresponding gene in our scenario. At the same time, however, the inhibition of mTOR is one of the best-documented life-prolonging interventions from invertebrates to mammals (Table 2). Given these conflicting premises, we find mTOR in almost all Fukomys tissues as not significantly altered; the exception is the adrenal gland, where mTOR is significantly down-regulated in breeders. It is obvious that, during evolution of Fukomys mole-rats, both the extension of breeders’ lifespan and an increase in their anabolic processes provided fitness benefits. Consequently, in the underlying organismic framework mTOR may have acquired expression patterns and functions, different from those in organisms studied before. Thus, in the future it may be worth to study Fukomys mTOR biology in more detail, particularly in respect to its potential role in naturally evolved ways towards lifespan extension and healthy aging.
Behavior of important aging-relevant genes and pathways in this study.
Gene/pathway | Regulation in indicated direction and species can reduce lifespan | Regulation in indicated direction and species can extend lifespan† | Differentially expressed in this study in indicated tissues and direction* | Gene/pathway is expressed mainly in the following tissues |
---|---|---|---|---|
GH1 | Mouse ↑1 | Mouse ↓2, rat ↓2 | Pituitary gland ↑ | Pituitary gland |
IGF1 | – | Mouse ↓2 | Adrenal gland ↑, ovary ↑ | Liver, ovary |
IGF1R | – | Mouse ↓1, worm ↓3, fly ↓3 | Adrenal gland ↓, ovary ↓ | Many |
KL | Mouse ↓3 | Mouse ↑3, worm ↑3 | Ovary ↓ | Endocrine tissue, kidney |
SIRT1 | – | Mouse ↑6, worm ↑7, fly ↑7 | – | All |
MYC | – | Mouse ↓8 | Thyroid ↑ | All |
mTOR | – | Mouse ↓9, worm↓9, fly ↓9 | Adrenal gland ↓ | All |
PRKAA2 (AMPK) | – | Worm ↑9 | – | Many |
TP53 | Mouse ↑11, Fly ↑11 | Worm ↓11, mouse ↑11, fly ↑11 | – | All |
SOD2 | – | Worm ↓10, fly ↑10 | – | All |
FOXO3 | Fly ↑12 | – | Ovary ↓, adrenal gland ↓ | All |
Protein synthesis | – | Mouse ↓8, worm ↓9, fly ↓9 | Many ↑ | All |
Proteasome | – | Worm ↑14, fly ↑14 | Gonads ↑, adrenal gland ↑ | All |
Lysosome | – | Worm ↑13 | – | All |
Respiratory chain | – | Worm↓15, fly ↓16, killifish ↓17 | Gonads ↑, Adrenal gland ↑, Blood ↑, Spleen ↑ | All |
Apoptosis | Mouse ↑11, fly ↑11 | Worm ↓11, mouse ↑11, fly ↑11 | Skin ↑, pituitary gland ↑ | All |
-
*Direction: breeder/non-breeder.
†Direction: old/young.
-
1Allolio and Arlt, 2002, 2Bartke, 2019, 3van Heemst, 2010, 4Elibol and Kilic, 2018, 5Kwon et al., 2015, 6Satoh et al., 2013, 7Burnett et al., 2011, 8Hofmann et al., 2015, 9Johnson et al., 2013, 10Edrey and Salmon, 2014, 11Feng et al., 2011, 12Giannakou et al., 2004, 13Carmona-Gutierrez et al., 2016, 14Saez and Vilchez, 2014, 15Dillin et al., 2002, 16Copeland et al., 2009, 17Baumgart et al., 2016.
–: either not affecting lifespan or not known to the best of our knowledge (columns 1 and 2); no change (column 3).
The results of differential expression of anabolic components such as the GH/IGF1 axis are surprising. They fall within a debate in aging research that has been highly controversial over time: based on the well-known fact that the expression and secretion of GH and IGF1 decline with age in humans and other mammals (Bartke, 2019), Rudman et al., 1990 administered synthetic GH to elderly subjects in 1990, thereby reversing a number of aging-associated effects such as expansion of adipose mass. This led to GH being celebrated as an anti-aging drug (Junnila et al., 2013), including dubious commercial offers. Today's aging research, on the contrary, strongly assumes that the enhanced activity of the GH/IGF1 axis accelerates aging and that its suppression could extend lifespan even in humans (López-Otín et al., 2013; Longo et al., 2015; Pitt and Kaeberlein, 2015). In addition to several studies of synthetic GH in humans yielding less convincing results than those of Rudman et al., the main reasons for this turn are the results of studies on short-lived model organisms. From worms to mice, the impairment of the GH/IGF1 axis by genetic intervention consistently led to longer lifespans (Table 2), for example, the up-regulation of Klotho – an IGF1 inhibitor – extended the mouse lifespan by as much as 30% (Kurosu et al., 2005). As with the impairment of the GH/IGF1 axis, reducing of protein synthesis by decreasing the expression of MYC, a basal transcription factor, extended the mouse lifespan by as much as 20% (Hofmann et al., 2015), whereas the impairment of the respiratory chain by rotenone resulted in prolongation of the killifish lifespan by 15% (Baumgart et al., 2016; Table 2).
Therefore, it is astonishing that massive up-regulation of these anabolic key components accompanies a lifespan extension of approximately 100% in long-lived mammals and potentially even contributes to it. Several points could help to resolve this apparent contradiction: first, the up-regulation of anabolic pathways and key genes is at least partially accompanied by the regulation of other mechanisms that could plausibly compensate for deleterious effects. For example, it is, widely assumed that the negative impact of enhanced protein synthesis on lifespan is to a large extent caused by the accumulation of damaged or misfolded proteins, which is also known to contribute to aging-associated neurodegenerative diseases (Hipkiss, 2007; Saez and Vilchez, 2014; Carmona and Michan, 2016). Up-regulation of the proteasome, as we observed in breeders in a weighted cross-tissue approach and especially in endocrine tissues, is known to counteract these effects by clearing damaged proteins, leading to lifespan extension in worm and fly (Saez and Vilchez, 2014). Enhanced proteasome activity has also been linked with higher longevity of the naked mole-rat compared to the laboratory mouse (Pérez et al., 2009; Rodriguez et al., 2012) and in comparative approaches across several mammalian lineages (Pride et al., 2015). We hypothesize that the simultaneously high anabolic synthesis and catabolic degradation of proteins will lead to a higher protein turnover rate in breeders and, accompanied with that, a reduced accumulation of damaged and misfolded proteins. Similarly, it seems plausible that the up-regulation of the mitochondrial respiratory chain (oxidative phosphorylation) in breeders is compensated for by simultaneous up-regulation of the reactive oxygen hallmark: the mitochondrial respiratory chain is the main source of cellular ROS that can damage DNA, proteins, and other cellular components (Balaban et al., 2005; Starkov, 2008); the reactive oxygen hallmark consists by definition of genes that are known to be up-regulated in response to ROS treatments. Unsurprisingly, at least 25% of these genes code for typical antioxidant enzymes such as thioredoxin, superoxide dismutase, peroxiredoxin, or catalase that can detoxify ROS. Furthermore, the known cancer-promoting effects of an enhanced GH/IGF1 axis (Junnila et al., 2013) could, to some extent, be compensated for by up-regulation of apoptosis and p53 signaling because these are the major mechanisms of cancer suppression (Bieging et al., 2014). More generally, potential lifespan-extending effects of moderate up-regulation of both the GH/IGF1 axis and ROS production can also be viewed in the light of the hormesis hypothesis (Ristow, 2014), which postulates that mild stressors can induce overall beneficial adaptive stress responses. In line with these arguments, we found higher resting metabolic rates in breeders compared to non-breeders in F. anselli (Schielke et al., 2017), a species closely related to F. micklemi.
A second point that could help to resolve this apparent contradiction concerns the time of intervention. The transition from non-breeder to breeder takes place in adulthood when by far the largest portion of the growth process has already been completed. In contrast, genetic interventions aimed at prolonging the lifespan by inhibiting the GH/IGF1 axis (Table 2) often affect the organisms throughout their entire lifespan, including infancy and youth. However, there are some exceptions: for example, late treatments with rapamycin prolong the life of mice almost as much as those that are started early (Johnson et al., 2013). Therefore, it is still under debate whether the up-regulation of translation and anabolic processes by the GH/IGF1 axis independently enhances growth and aging or enhances aging only secondarily as a consequence of accelerated growth (Bartke, 2017). Our results are an argument for the latter.
A third point is the question of the transferability of knowledge obtained in one species to other species. Most insights into current aging research originate from very short-lived model organisms (Table 2). It is clear, on the other hand, that the observed effects of lifespan-prolonging interventions listed in Table 2 are by far the smallest in the model organisms with the relatively longest lifespans: mice and rats. For example, knockout or knockdown of IGF1R in worms that normally live only a few days can result in a lifespan extensions of more than 500% (Houthoofd et al., 2004). For mice that have a maximum lifespan of 4 years (Tacutu et al., 2013), comparable interventions only lead to an extension of the lifespan by about 25% (Holzenberger et al., 2003).It seems reasonable to hypothesize that the effect size of such interventions, regardless of their direction, would diminish further in mole-rats or humans that can live more than 20 or 100 years, respectively, especially if the magnitude of gene expression change is small as in the comparison of breeders and non-breeders. Therefore, it may also be possible that those gene expression patterns caused by our lifespan-extending intervention in Fukomys mole-rats highlight particularly the differences in aging mechanisms between short-lived and long-lived species.
As a fourth perspective, one could interpret the down-regulation of the GH/IGF1 axis in non-breeders as a by-product of the apparent up-regulation of the HPA axis in non-breeders, which may well be adaptive in itself. In the wild, non-breeding Fukomys mole-rats can maximize their inclusive fitness by either supporting their kin in their natal family or by founding a new family elsewhere. It has therefore been suggested that the shorter lives of non-breeders could be adaptive on the ultimate level if longevity were traded against some other fitness traits, such as competitiveness, to defend colonies against intruders or to enhance their chances for successful dispersal (Dammann and Burda, 2006). A constitutively more activated HPA stress axis is expected to offer advantages for both family defense and dispersal but it carries the risk of status-specific aging symptoms, such as muscle weakness, lower GH/IGF1 axis activity, and lower bone density in the long run (Ferraù and Korbonits, 2015). This effect may become even more pronounced under laboratory conditions where grown non-breeders cannot decide to disperse even if they would like to.
However, even the down-regulation of the GH/IGF-axis may be adaptive in itself for non-breeders if it has the potential to protect them from further damage. Note that for today's conventional view that stronger activation of the GH/IGF1 axis accelerates aging (Table 2) it is generally challenging that decreasing activity is well documented to correlate with chronological age in a wide range of mammals, including mice, rats, dogs, and humans (Bartke, 2019). Also, decreasing activity correlates, under pathological conditions such as Cushing’s syndrome, with many symptoms akin to aging (Figure 5). It has been suggested that one solution to this apparent contradiction may be that the GH/IGF1 axis is adaptively down-regulated in aging organisms as a reaction to already accumulated aging-related symptoms so that additional damage can be avoided (Berryman et al., 2008; Milman et al., 2016).
Finally, recent findings of positively selected genes in African mole-rats (family Bathyergidae, containing also Fukomys and Heterocephalus) could partly explain some of the surprising results. It is striking that translation and oxidative phosphorylation were among the strongest differentially expressed molecular processes concerning the breeding status. Earlier, these processes were also reported to be the most affected by positive selection in the phylogeny of African mole-rats (Sahm et al., 2018b). Furthermore, IGF1 was one of 13 genes that were found to be under positive selection in the last common ancestor of the mole-rats. This may indicate that the corresponding mechanisms were evolutionarily adapted to be less detrimental and make their up-regulation more compatible with a long lifespan. Since the mere fact of positive selection does not permit to draw conclusions about the direction of the mechanistic effect, this hypothesis, however, needs to remain speculative.
Conclusions
We performed a comprehensive transcriptome analysis that, for the first time within mammalian species, compared naturally occurring cohorts of species with massively diverging longevities. The comparison of faster-aging Fukomys non-breeders with similar animals that were experimentally elevated to the slower-aging breeder status revealed by far the most robust transcriptome differences within endocrine tissue: adrenal gland, ovary, thyroid, and pituitary gland. Genes and pathways involved in anabolism, such as GH, IGF1, translation, and oxidative phosphorylation, were differentially expressed. Their inhibition is among the best-documented life-prolonging interventions in a wide range of short-lived model organisms (Table 2). Surprisingly, however, we found that the expression of these mechanisms was consistently higher in slower-aging breeders than in faster-aging non-breeders. This indicates that even basic molecular mechanisms of the aging process known from short-lived species cannot easily be transferred to long-lived species. In particular, this applies to the role of the GH/IGF1 axis, which has in recent years been predominantly described as having a negative impact on lifespan (López-Otín et al., 2013; Pitt and Kaeberlein, 2015; Bartke, 2017). In addition, special features of the mole-rats could also contribute to the explanation of the unexpected result that genes and processes differentially expressed between reproductive statuses were also strongly altered during the evolution of the mole-rats (Sahm et al., 2018b). Another intriguing possibility is that, in line with the hormesis hypothesis (Ristow, 2014), moderate harmful effects of anabolic processes can be hyper-compensated for by up-regulation of pathways such as proteasomes, P53 signaling, and antioxidant defense against ROS that we observed in slower-aging breeders as well.
Furthermore, our work provides evidence that the HPA stress axis is a key regulator for the observed downstream effects, including the lifespan difference. The effects are likely to be triggered by differential expression of the gene MC2R coding for the ACTH receptor, resulting in an altered stress response in breeders vs. non-breeders. This is supported by the facts that (1) cortisol levels in the shorter-lived non-breeders are elevated and (2) life expectancy of non-breeders rises when cortisol levels are reduced. Furthermore, the set of direct and indirect target genes of the glucocorticoid receptors is strongly affected by differential expression, and numerous known downstream effects of glucocorticoid excess have been demonstrated for non-breeders, such as muscle weakness, weight gain, and GH/IGF1 axis impairment. Overall, this evidence suggests that MC2R and other genes along the described signaling pathway are promising targets for possible interventions in aging research.
Materials and methods
Animal care and sampling
Request a detailed protocolAll animals were housed in glass terraria with dimensions adjusted to the size of the family (min. 40 cm × 60 cm) in the Department of General Zoology, Faculty of Biosciences, University of Duisburg-Essen. The terraria are filled with a 5 cm layer of horticultural peat or sawdust. Tissue paper strips, tubes, and solid shelters were provided as bedding/nesting materials and environmental enrichment. Potatoes and carrots are supplied ad libitum as stable food, supplemented with apples, lettuce, and cereals. Fukomys mole-rats do not drink free water. Temperature was kept fairly constant at 26 ± 1°C and humidity at approximately 40%. The daily rhythm was set to 12 hr darkness and 12 hr light.
New breeder pairs (new families) were established between March and May 2014. Each new family was founded by two unfamiliar, randomly chosen adult non-breeders of similar age (min/max/mean: 1.56/6.5/3.58 years in F. mechowii; 1.8/5.4/3.1 years in F. micklemi) and opposite sex and were taken from already existing separate colonies. These founder animals were moved to a new terrarium in which they were permanently mated. In both species, more than 80% of these new pairs reproduced within the first 12 months (total number of offspring by the end of the year 2016, 82 F. micklemi and 81 F. mechowii). Only founders with offspring were subsequently assigned to the breeder cohort; founders without offspring were excluded from the study. Non-breeders remained in their natal family together with both parents and other siblings.
F. mechowii were sampled in five distinct sampling sessions between March 2015 and winter 2016/2017. F. micklemi were sampled in three distinct sampling sessions between November 2016 and July 2017. In both species, females were killed 4–6 months later than their male mates to ensure that these breeder females were neither pregnant nor lactating at the time of sampling in order to exclude additional uncontrolled variables. In total, females/males selected as breeders spent on average 513/353 (F. mechowii) and 1095/997 (F. micklemi) days in this state before sampling (Supplementary file 1d, e).
Before sampling, animals were anesthetized with 6 mg/kg ketamine combined with 2.5 mg/kg xylazine (Garcia Montero et al., 2015). Once the animals lost their pedal withdrawal reflex, 1–2 ml of blood was collected by cardiac puncture, and the animals were killed by surgical decapitation. Blood samples (100 µl) were collected in RNAprotect Animal Blood reagent (Qiagen, Venlo, Netherlands). Tissue samples – hippocampus, hypothalamus, pituitary gland, thyroid, heart, skeletal muscle (M. quadriceps femoris), lateral skin, small intestine (ileum), upper colon, spleen, liver, kidney, adrenal gland, testis, and ovary – were transferred to RNAlater (Qiagen, Venlo, Netherlands) immediately after dissection and, following incubation, were stored at −80°C until analysis.
Animal housing and tissue collection were compliant with national and state legislation (breeding allowances 32-2-1180-71/328 and 32-2-11-80-71/345; ethics/animal experimentation approval 84-02.04.2013/A164, Landesamt für Natur-, Umwelt- und Verbraucherschutz Nordrhein-Westfalen).
RNA preparation and sequencing
Request a detailed protocolFor all tissues except blood, RNA was purified with the RNeasy Mini Kit (Qiagen) according to the manufacturer’s protocol. Blood RNA was purified with the RNeasy Protect Animal Blood Kit (Qiagen). Kidney and heart samples were treated with proteinase K before extraction as recommended by the manufacturer. Library preparation was performed using the TruSeq RNA v2 kit (Illumina, San Diego, USA), which includes selection of poly-adenylated RNA molecules. RNA-seq was performed by single-end sequencing with 51 cycles in high-output mode on a HiSeq 2500 sequencing system (Illumina) and with at least 20 million reads per sample, as described in Supplementary file 1i. We have taken care to concentrate the samples of a tissue on a few sequencing runs and to create a balance of breeders and non-breeders within each run (Supplementary file 1j). Read data for F. mechowii and F. micklemi were deposited as European Nucleotide Archive study with the ID PRJEB29798 (Supplementary file 1i, j).
Read mapping and quantification
Request a detailed protocolIt was ensured for all samples that the results of the respective sequencing passed ‘per base’ and ‘per sequence’ quality checks of FASTQC (Andrews, 2020). The reads were then mapped against previously published and with human gene symbols annotated F. mechowii and F. micklemi transcriptome data (Sahm et al., 2018a; Sahm et al., 2018b). For both species, only the longest transcript isoform per gene was used; this is the method of choice for selecting a representative variant in large-scale experiments (Ezkurdia et al., 2015; Source data 1b, c). This selection resulted in 15,864 reference transcripts (genes) for F. mechowii and in 16,400 for F. micklemi. After mapping and quantification, we further analyzed only those reference transcripts whose gene symbols were present in the transcript catalogs of both species – this was the case for 15,199 transcripts (the size of the union was 17,065). As mapping, algorithm ‘bwa aln’ of the Burrows-Wheeler Aligner (Li and Durbin, 2009) was used, allowing no gaps and a maximum of two mismatches in the alignment. Only those reads that could be uniquely mapped to the respective gene were used for quantification. Read counts per gene and sample can be found in Source data 1d. As another check, we ensured that all samples exhibited a Pearson correlation coefficient of at least 90% in a pairwise comparison based on log2-transformed read counts against all other samples of the same experimental group as defined by samples that were equal in the tissue as well as the species, sex, and reproductive status of the source animal.
DEGs analysis
Request a detailed protocolp-Values for differential gene expression and fold-changes were determined with DESeq2 (Love et al., 2014) and a multifactorial design. The DESeq2 algorithm also includes strict filtering based on a normalized mean gene count that makes further pre-filtering unnecessary (Love et al., 2014). Therefore, those genes whose read count was zero for all examined samples were removed before further analysis, thereby reducing the number of analyzed genes to 15,181. The multifactorial design means that, separately for each tissue, we input the read count data of samples across species, sex, and reproductive status into DESeq2 for each sample. This allowed DESeq2 to perform DEG analysis between the two possible states of each of the variables by controlling for additional variance in the other two variables. This approach resulted in a four-times higher sample size than with an approach that would have been based on comparisons of two experimental groups, each of which would be equal in tissue, species, sex, and reproductive status. It is known that the statistical power in RNA-seq experiments can increase considerably with sample size (Ching et al., 2014). p-Values were corrected for multiple testing with the Benjamini–Hochberg correction (Benjamini and Hochberg, 1995) (FDR).
The results of the DEG analysis can be found in Source data 1e–g.
Enrichment analysis on pathway and cross-tissue level
Request a detailed protocolLet ) represent the p-values obtained from differential gene expression analysis in the tissue corresponding with index and the indices corresponding to the examined genes. Furthermore, let , with and , represent the p-values of genes with the indices belonging to a corresponding pathway that is tested for enrichment of differential expression signals. To determine the enrichment p-values at the pathway level, we calculated the test statistic for the gene indices in tissue index according to Fisher’s method for combining p-values:
If the underlying test statistics are independent, a combined p-value may be determined by Fishers’s method using a χ2 distribution. However, since this assumption might be violated in gene expression studies (Väremo et al., 2013), for example, due to gene regulation cascades, we have empirically estimated the null distributions, as recommended by the literature (Fridley et al., 2010; Evangelou et al., 2012; Poole et al., 2016), using a resampling approach. For this purpose, for each KEGG/MSigDB pathway with m elements, we repeatedly drew p-values from the total set of p-values ) and calculated the test statistics. Specifically, this was done by calculating for 1000 random draws, each without replacement, , with . When required for numerical precision, that is, , sampling was performed again with 10,000 and 100,000 random draws. The combined p-values for our hypothesis tests can now be obtained by determining the probability of a test statistic in the empirically estimated null distribution. In addition, the indices were divided into and , depending on whether their fold-change was or in breeder vs. non-breeder comparison, and and calculated. The ratio was used as an indicator for functional up- or down-regulation of the corresponding pathway (Figure 4, Figure 4—figure supplements 1 and 2). Using this approach, enrichment p-values were estimated for all KEGG pathways (Kanehisa et al., 2017) and MSigDB hallmarks (Liberzon et al., 2015), as well as across all examined tissues (Source data 2a, b). In addition, the procedure was applied to test whether the known 300 direct and 1300 indirect glucocorticoid receptor target genes (Phuc Le et al., 2005) were enriched for status-dependent differential expression signals (Supplementary file 1g,h).
Similarly, cross-tissue DEG-p-values were weighted with a modified test statistic (c.f. Heard and Rubin-Delanchy, 2018). The test statistic weights the p-values of the various tissues by the respective expression levels in those tissues. This ensures that, for example, for a ubiquitously expressed gene such as TP53 all tissues contribute relatively equally to the cross-tissue p-value, whereas for typical steroid hormone biosynthesis genes such as CYP11A1 the endocrine tissue results almost exclusively determine the weighted cross-tissue p-value. Given the definitions from above, we calculated the weighted cross-tissue test statistic for the gene as follows:
with
where is logarithmic fold-change between reproductive states and is the normalized mean expression (across sexes, species, and reproductive status) for the gene with index and tissue with index – both calculated by DESeq2 (Love et al., 2014). Furthermore, is the signum function, and is the number of examined tissues. The test statistic rewards, based on the mentioned assumption, consistency in the direction of gene regulation throughout tissues. All calculated values for , , , as well as the resulting and p-values can be found in Source data 2c.
Finally, weighted cross-tissue enrichment p-values at the pathway level were estimated by applying the above-described method at the pathway level (based on test statistic ) to the gene-level weighted cross-tissue p-values. p-Values were corrected for multiple testing with the Benjamini–Hochberg correction (Benjamini and Hochberg, 1995) (FDR).
To not solely rely on p-value based statistics, we corroborated the pathways identified as differentially expressed by this procedure on the cross-tissue level by applying an alternative fold-change-based test statistics. The test statistics , for single tissues, and , for the weighed cross-tissue level, were evaluated by resampling via 1000–100,000 random draws the same way as and (see above), to obtain p-values for the examined gene sets:
Weighted gene co-expression network analysis
Request a detailed protocolWe used the WGCNA R package to perform weighted correlation network analysis (Langfelder and Horvath, 2008) of all 636 samples at once. We followed the authors’ usage recommendation by choosing a soft power threshold based on scale-free topology and mean connectivity development (we chose power=26 with a soft R2 of 0.92 and a mean connectivity of 38.6), using biweight midcorrelation, setting maxPOutliers to 0.1, and using 'signed' both as network and topological overlap matrix type. The maximum block size was chosen such that the analysis was performed with a single block and the minimum module size was set to 30. The analysis divided the genes into 26 modules, of which 5 were enriched for reproductive status DEGs based on Fisher’s exact test and an FDR threshold of 0.05. Those five modules were tested for enrichment among KEGG pathways (Kanehisa et al., 2017) with the same test and significance threshold (Supplementary file 1k). In addition, module eigengenes were determined and clustered (Supplementary file 1k). Then the topological overlap matrix that resulted from the WGCNA analysis (, where the row indices correspond to genes and the column indices correspond to samples) was used to determine pairwise connectivity between all KEGG pathways that showed differential expression at the weighted cross-tissue level (Figure 4—figure supplement 3). Based on the definition of connectivity of genes in a WGCNA analysis (Langfelder and Horvath, 2008), we defined the connectivity between two sets of indices and each corresponding to genes as . p-Values for the connectivities were determined against null distributions that were empirically estimated by determining for each pair and the connectivities of 10,000 pairs of each and || randomly drawn indices (without replacement), respectively. Since 'signed' was used as the network and topological overlap matrix type, the tests were one-sided.
Other analysis steps
Request a detailed protocolHierarchical clustering (Figure 2—figure supplement 1) was performed based on Pearson correlation coefficients of log2-transformed read counts between all sample pairs using the complete-linkage method (Defays, 1977). The principal variant component analysis (Figure 2A) was performed with the pvca package from Bioconductor (Bushel, 2013) and a minimum demanded percentile value of the amount of the variabilities, which the selected principal components needs to explain, of 0.5. Enrichments of DEGs among genes enlisted in the Digital Aging Atlas database (Craig et al., 2015; Source data 2d) were determined with Fisher’s exact test, the Benjamini–Hochberg method (FDR) (Benjamini and Hochberg, 1995) for multiple test correction, and a significance threshold of 0.05. Pathway visualization (Figure 4—figure supplements 4–7) was performed with Pathview (Luo and Brouwer, 2013). For the direction analysis of Fukomys reproductive status DEGs in previous experiments in naked mole-rats and guinea pigs, we examined those 10 tissues that were examined in all species (Supplementary file 1f, Source data 2e). Separately for each tissue and combination of species – naked mole-rat or guinea pig – and sex, we determined how many Fukomys reproductive status DEGs were up-regulated or down-regulated. We also performed two-sided binomial tests on each of these number pairs with a hypothesized success probability of 0.5. Furthermore, for each combination of species and sex, two-sided exact binomial tests using 0.5 as parameter were performed based on the sums of up-regulated and down-regulated genes across tissues (Supplementary file 1f). For enrichment analysis of direct and indirect glucocorticoid receptor target genes, mouse mRNA RefSeq IDs from Phuc Le et al., 2005 were translated to human Entrez IDs and gene symbols via Ensembl Biomart (Source data 2f). To statistically analyze the weight gain of the animals during the experiment, we used a type II ANOVA with status, species, sex, and age as independent variables (Supplementary file 1b); the weight gain, defined as the difference in weights at beginning and the end of the experiment, as dependent variable (Figure 6B); and no interaction terms. If interaction terms were also used for the model, the p-value for the difference in means between breeders and non-breeders changed from 7.49*10−3, as reported above, to 7.46*10−6. To compare our gene expression with Cushing data, we correlated the respective log2fold-changes from our skin (as computed by DESeq2) and cross-tissue results with those of human subcutaneous adipose tissue (Hochberg et al., 2015). The comparison directions were breeder/non-breeder and control/Cushing patient, respectively. We performed weighted Pearson correlation using the human expression, specifically the baseMean values given in the human study, as weights. Ensembl IDs used in the human study were translated to entrez IDs via the org.Hs.eg.db Bioconductor package (Source data 2g).
Bone density measurements
Request a detailed protocolFrozen carcasses of all F. micklemi that had been part of the transcriptome study were scanned with a self-shielded desktop small-animal computed tomography scanner (X-CUBE, Molecubes, Belgium). The x-ray source was a tungsten anode (peak voltage, 50 kVp; tube current, 350 µA; 0.8 mm aluminum filter). The detector was a cesium iodide (CsI) flat panel, building up a screen with 1536 × 864 pixels. Measurements were carried for individual 120 ms exposures, with angular sampling intervals of 940 exposures per rotation, for a total of seven rotations and a total exposure time of 789.6 s.
First, we first performed a calibration of the reconstructed CT data in terms of equivalent mineral density. For this purpose, we used a bone density calibration phantom (BDC; QRM GmbH, Moehrendorf, Germany) composed of five cylindrical inserts with a diameter of 5 mm containing various densities of calcium hydroxyapatite (CaHA) surrounded by epoxy resin on a cylindrical shape. The nominal values of CaHA were 0, 100, 200, 400, and 800 mg HA/cm3, corresponding to a density of 1.13, 1.16, 1.25, 1.64, and 1.90 g/cm3 (certified with an accuracy of ±0.5%). The BDC was imaged and reconstructed with the same specifications as each probe. From the reconstructed Hounsfield units, a linear relationship was determined against the known mineral concentrations.
Reconstruction of the acquired computed tomography data was carried out with an Image Space Reconstruction Algorithm, and spatial resolution was limited to the 100 µm voxel matrix reconstruction. Spherical regions of interest (radius, 0.7 mm) were drawn on the sagittal plane of vertebrae T12, L1, and L2. Care was taken to include all cancellous bone, excluding the cortical edges. Average Hounsfield unit values were computed on the calibration curve to finally retrieve equivalent densities of the regions of interest.
Statistical analysis was performed using general linear models with bone density (Hounsfield units) as dependent variable, age (in days) as continuous covariate, and reproductive status and sex as nominal cofactors. Models were calculated for each vertebra individually (individual models) and across all three vertebrae (full model); in this latter case, vertebral number was added as additional categorial cofactor. In all models, only main effects were calculated, no interactions. Analyses were performed with IBM SPSS version 25 (Figure 6C, Supplementary file 1l).
Analysis of mortality rates
Request a detailed protocolIncreasing mortality with age is defined as actuarial senescence (e.g., Møller, 2006). To calculate status-specific mortality rates across age classes, we subdivided the first 10 years of observation time (corresponding to 1.5–11.5 years of calendar age) into intervals of 6 months. We then recapitulated the times each individual spent as conventional non-breeder (both parents present; cortisol levels high, see Begall et al., 2021), orphaned non-breeder (at least one parent absent since more than 6 months; cortisol levels low, see Begall et al., 2021), or breeder (cortisol levels low, see Begall et al., 2021), and calculated status-specific annual mortality rates for each interval by dividing the number of death events of the respective status during a given interval by the sum of observation years for all individuals holding this status in the same interval. In each interval, individuals that were alive across the entire interval contributed 0.5 years observation time, whereas individuals who died or were censored during that interval contributed the time (expressed in years) up to the death or censorship event. Only animals whose exact age was known were included (i.e., wild-caught animals were excluded). To account for decreasing sample sizes with advancing age, mortality calculations per status were restricted to those intervals containing at least 10 individuals of that status (i.e., around the first 6 years of life for conventional non-breeders and about 11.5 years for orphaned non-breeders and breeders). Status-specific mortality rates were compared by repeated measures ANOVA followed by Tukey post-hoc testing (for the first 6 years) or by paired t-test (breeders vs. all non-breeders, and breeders vs. orphaned non-breeders across all 10 years of observation time) after confirming normal distribution using Shapiro–Wilk normality test using GraphPad Prism.
Control analyses
Request a detailed protocolFor practical reasons, samples collected for this project over several years were sequenced progressively in different batches. Therefore, we analyzed our data and sample scheme across all tissues and both species using BatchQC (Manimaran et al., 2016). In summary, the extensive search did not reveal any possible batch effects in the data set from F. micklemi, although they may not be ruled out completely for F. mechowii (Source data 2f). A source for the minor batch effects in F. mechowii may be attributable to the repeated sequencing of the same sample(s) in different sequencing experiments – an approach avoided during generation of the F. micklemi data set. However, because we followed a robust multifactorial analysis, in which differences must occur consistently across both species to be considered (see DEGs analysis in Materials and methods), false-positive results due to batch effects are substantially reduced (Manimaran et al., 2016).
In addition to detecting DEGs with DESeq2 (Love et al., 2014), we also examined the central contrast, that is, the juxtaposition of workers and breeders, using the combination of limma (Smyth, 2004) and voom (Law et al., 2014). The absolute numbers of DEGs found in each tissue, and thus the highly uneven distribution, were nearly identical for all tissues except thyroid (Figure 2—figure supplement 2). For the latter, however, a closer look shows that 85% of the DESeq2 DEGs (threshold FDR < 0.05) are confirmed if the FDR threshold for limma/voom is raised from 0.05 to 0.1. If an FDR threshold of 0.05 is applied for both approaches across all tissues, a total of 58% of the DESeq2-DEGs are confirmed by limma/voom. Notably, all genes mentioned in the paper, including the candidates from Table 1, are among this set. Using a less strict threshold of 0.1 for limma/voom, 91% of DESeq2-DEGs are confirmed.
We checked the DEG enrichment for potential detection biases with regard to expression level of genes/pathways or size of pathways at the cross-tissue level. For this, we looked at the correlations of the respective p-values with gene expression, pathway expression, and pathway size (Figure 4—figure supplements 8–10). In summary, we find a slightly enhanced sensitivity for highly expressed pathways at the cross-tissue level.
In order to avoid potential biases by using solely p-values as indicator of differential gene expression, we additionally resorted to fold-changes as a test statistic for our pathway enrichment analyses. This exercise roughly confirms half of the pathways/hallmarks identified as differentially expressed at the cross-tissue level. Notably the set includes, consistent with the HPA hypothesis, steroid hormone biosynthesis and ribosomes as a central part of up-regulated anabolism in breeders (Figure 4—figure supplement 11). Furthermore, myogenesis, the P53 pathway, and coagulation were confirmed as cross-tissue differentially expressed by this approach. On the other hand, proteasome and oxidative phosphorylation gene sets were not identified by the fold-change method.
Data availability
Read datasets generated during the current study are available in the European Nucleotide Archive, study ID: PRJEB29798.
-
European Nucleotide ArchiveID PRJEB29798. Transcriptome signatures of fast vs. slow aging in Fukomys mole-rat breeders vs. non-breeders.
References
-
DHEA treatment: myth or reality?Trends in Endocrinology & Metabolism 13:288–294.https://doi.org/10.1016/S1043-2760(02)00617-3
-
SoftwareFastQC a Quality Control Tool for High Throughput Sequence DataFastQC a Quality Control Tool for High Throughput Sequence Data.
-
Checkpoints in TNF-Induced cell death: implications in inflammation and CancerTrends in Molecular Medicine 24:49–65.https://doi.org/10.1016/j.molmed.2017.11.002
-
Comparative biology of agingThe Journals of Gerontology Series A: Biological Sciences and Medical Sciences 64:199–201.https://doi.org/10.1093/gerona/gln060
-
Candidate bird species for use in aging researchILAR Journal 52:89–96.https://doi.org/10.1093/ilar.52.1.89
-
Somatic growth, aging, and longevityNpj Aging and Mechanisms of Disease 3:14.https://doi.org/10.1038/s41514-017-0014-y
-
Growth hormone and aging: updated reviewThe World Journal of Men's Health 37:19–30.https://doi.org/10.5534/wjmh.180018
-
Dehydroepiandrosterone (DHEA): a fountain of youth?The Journal of Clinical Endocrinology & Metabolism 81:3147–3151.https://doi.org/10.1210/jcem.81.9.8784058
-
Basic and clinical pharmacology of glucocorticosteroidsAnesthesia Progress 60:25–32.https://doi.org/10.2344/0003-3006-60.1.25
-
Life expectancy, family constellation, and stress in giant mole-rats (Fukomys mechowii)Philosophical Transactions of the Royal Society B 376:20200207.https://doi.org/10.1098/rstb.2020.0207
-
Controlling the false discovery rate: a practical and powerful approach to multiple testingJournal of the Royal Statistical Society: Series B 57:289–300.https://doi.org/10.1111/j.2517-6161.1995.tb02031.x
-
Role of the GH/IGF-1 Axis in lifespan and healthspan: lessons from animal modelsGrowth Hormone & IGF Research 18:455–471.https://doi.org/10.1016/j.ghir.2008.05.005
-
ACTH-receptor expression, regulation and role in adrenocortial tumor formationEuropean Journal of Endocrinology 144:199–206.https://doi.org/10.1530/eje.0.1440199
-
Unravelling mechanisms of p53-mediated tumour suppressionNature Reviews Cancer 14:359–370.https://doi.org/10.1038/nrc3711
-
Somatotrope GHRH/GH/IGF-1 Axis at the crossroads between immunosenescence and frailtyAnnals of the New York Academy of Sciences 1351:61–67.https://doi.org/10.1111/nyas.12857
-
Negligible senescence in the longest living rodent, the naked mole-rat: insights from a successfully aging speciesJournal of Comparative Physiology B 178:439–445.https://doi.org/10.1007/s00360-007-0237-5
-
DHEA(S)--a novel marker in Cushing's diseaseActa Neurochirurgica 155:479–484.https://doi.org/10.1007/s00701-012-1596-6
-
SoftwarePvca: Principal Variance Component Analysis (PVCA)Pvca: Principal Variance Component Analysis (PVCA).
-
The GH/IGF-1 Axis in growth and development: new insights derived from animal modelsAdvances in Pediatrics 57:331–351.https://doi.org/10.1016/j.yapd.2010.09.003
-
Biology of healthy aging and longevityRevista De Investigacion Clinica; Organo Del Hospital De Enfermedades De La Nutricion 68:7–16.
-
The crucial impact of lysosomes in aging and longevityAgeing Research Reviews 32:2–12.https://doi.org/10.1016/j.arr.2016.04.009
-
Impairment of GH/IGF-1 Axis in the liver of patients with HCV-Related chronic hepatitisHormone and Metabolic Research 50:145–151.https://doi.org/10.1055/s-0043-118911
-
Dehydroepiandrosterone - is the fountain of youth drying out?Physiological Research 52:397–407.
-
The endocrine system and ageingThe Journal of Pathology 211:173–180.https://doi.org/10.1002/path.2110
-
The digital ageing atlas: integrating the diversity of age-related changes into a unified resourceNucleic Acids Research 43:D873–D878.https://doi.org/10.1093/nar/gku843
-
Advanced glycation end-products as markers of aging and longevity in the long-lived Ansell's mole-rat (Fukomys anselli)The Journals of Gerontology. Series A, Biological Sciences and Medical Sciences 67:573–583.https://doi.org/10.1093/gerona/glr208
-
Slow aging in mammals-Lessons from african mole-rats and batsSeminars in Cell & Developmental Biology 70:154–163.https://doi.org/10.1016/j.semcdb.2017.07.006
-
Sexual activity and reproduction delay ageing in a mammalCurrent Biology 16:R117–R118.https://doi.org/10.1016/j.cub.2006.02.012
-
An efficient algorithm for a complete link methodThe Computer Journal 20:364–366.https://doi.org/10.1093/comjnl/20.4.364
-
Revisiting an age-old question regarding oxidative stressFree Radical Biology and Medicine 71:368–378.https://doi.org/10.1016/j.freeradbiomed.2014.03.038
-
Morbidity and mortality in Cushing's disease: an epidemiological approachClinical Endocrinology 40:479–484.https://doi.org/10.1111/j.1365-2265.1994.tb02486.x
-
Most highly expressed protein-coding genes have a single dominant isoformJournal of Proteome Research 14:1880–1887.https://doi.org/10.1021/pr501286b
-
Metabolic comorbidities in Cushing's syndromeEuropean Journal of Endocrinology 173:M133–M157.https://doi.org/10.1530/EJE-15-0354
-
Chemical restraint of african mole-rats (Fukomys sp.) with a combination of ketamine and xylazineVeterinary Anaesthesia and Analgesia 42:187–191.https://doi.org/10.1111/vaa.12180
-
No evidence for hepatic conversion of dehydroepiandrosterone (DHEA) sulfate to DHEA: in vivo and in vitro studiesThe Journal of Clinical Endocrinology & Metabolism 90:3600–3605.https://doi.org/10.1210/jc.2004-2386
-
Aging: a theory based on free radical and radiation chemistryJournal of Gerontology 11:298–300.https://doi.org/10.1093/geronj/11.3.298
-
Aging: overviewAnnals of the New York Academy of Sciences 928:1–21.https://doi.org/10.1111/j.1749-6632.2001.tb05631.x
-
Choosing between methods of combining $p$-valuesBiometrika 105:239–246.https://doi.org/10.1093/biomet/asx076
-
On why decreasing protein synthesis can increase lifespanMechanisms of Ageing and Development 128:412–414.https://doi.org/10.1016/j.mad.2007.03.002
-
Gene expression changes in subcutaneous adipose tissue due to Cushing's diseaseJournal of Molecular Endocrinology 55:81–94.https://doi.org/10.1530/JME-15-0119
-
Extending life-span in C. elegansScience 305:1238–1239.https://doi.org/10.1126/science.305.5688.1238c
-
Eusociality has evolved independently in two genera of bathyergid mole-rats — but occurs in no other subterranean mammalBehavioral Ecology and Sociobiology 33:253–260.https://doi.org/10.1007/BF02027122
-
The endocrinology of ageing: a Mini-ReviewGerontology 61:291–300.https://doi.org/10.1159/000367692
-
Regulation of IGF -1 signaling by microRNAsFrontiers in Genetics 5:472.https://doi.org/10.3389/fgene.2014.00472
-
The GH/IGF-1 Axis in ageing and longevityNature Reviews Endocrinology 9:366–376.https://doi.org/10.1038/nrendo.2013.67
-
KEGG: new perspectives on genomes, pathways, diseases and drugsNucleic Acids Research 45:D353–D361.https://doi.org/10.1093/nar/gkw1092
-
Social insects as a model to study the molecular basis of ageingExperimental Gerontology 41:553–556.https://doi.org/10.1016/j.exger.2006.04.002
-
Recent advances in the diagnosis of Cushing's syndrome in dogsTopics in Companion Animal Medicine 27:21–24.https://doi.org/10.1053/j.tcam.2012.06.001
-
Expression of SIRT1 and SIRT3 varies according to age in miceAnatomy & Cell Biology 48:54–61.https://doi.org/10.5115/acb.2015.48.1.54
-
The SCARB1 gene is associated with lipid response to dietary and pharmacological interventionsJournal of Human Genetics 53:709–717.https://doi.org/10.1007/s10038-008-0302-2
-
A practical guide to the monitoring and management of the complications of systemic corticosteroid therapyAllergy, Asthma & Clinical Immunology 9:30.https://doi.org/10.1186/1710-1492-9-30
-
Effects of cushing disease on bone mineral density in a pediatric populationThe Journal of Pediatrics 156:1001–1005.https://doi.org/10.1016/j.jpeds.2009.12.027
-
The relevance of coagulation in cardiovascular disease: what do the biomarkers tell Us?Thrombosis and Haemostasis 112:860–867.https://doi.org/10.1160/TH14-03-0199
-
Equine Cushing's diseaseVeterinary Clinics of North America: Equine Practice 18:533–543.https://doi.org/10.1016/S0749-0739(02)00038-X
-
Cushing's syndrome due to adrenocortical adenoma in a catTijdschrift Voor Diergeneeskunde 103:1048–1051.
-
Sociality, age at first reproduction and senescence: comparative analyses of birdsJournal of Evolutionary Biology 19:682–689.https://doi.org/10.1111/j.1420-9101.2005.01065.x
-
Protein homeostasis and aging: taking care of proteins from the cradle to the graveThe Journals of Gerontology Series A: Biological Sciences and Medical Sciences 64:167–170.https://doi.org/10.1093/gerona/gln071
-
Why is aging conserved and what can we do about it?PLOS Biology 13:e1002131.https://doi.org/10.1371/journal.pbio.1002131
-
Long-lived species have improved proteostasis compared to phylogenetically-related shorter-lived speciesBiochemical and Biophysical Research Communications 457:669–675.https://doi.org/10.1016/j.bbrc.2015.01.046
-
Skeletal growth and bone mineral acquisition in type 1 diabetic children; abnormalities of the GH/IGF-1 AxisGrowth Hormone & IGF Research 34:13–21.https://doi.org/10.1016/j.ghir.2017.04.003
-
Effects of human growth hormone in men over 60 years oldNew England Journal of Medicine 323:1–6.https://doi.org/10.1056/NEJM199007053230101
-
Age related dehydroepiandrosterone decrease: clinical significance and therapeutic interest]Revue Medicale Suisse 11:321–324.
-
Extraordinary life spans of naked mole-rats (Heterocephalus glaber)Journal of Zoology 258:307–311.https://doi.org/10.1017/S0952836902001437
-
Linear models and empirical bayes methods for assessing differential expression in microarray experimentsStatistical Applications in Genetics and Molecular Biology 3:1–25.https://doi.org/10.2202/1544-6115.1027
-
The role of mitochondria in reactive oxygen species metabolism and signalingAnnals of the New York Academy of Sciences 1147:37–52.https://doi.org/10.1196/annals.1427.015
-
Human ageing genomic resources: integrated databases and tools for the biology and genetics of ageingNucleic Acids Research 41:D1027–D1033.https://doi.org/10.1093/nar/gks1155
-
Dehydroepiandrosterone (DHEA)--a precursor steroid or an active hormone in human physiologyThe Journal of Sexual Medicine 8:2960–2982.https://doi.org/10.1111/j.1743-6109.2011.02523.x
-
Rapid intra-adrenal feedback regulation of glucocorticoid synthesisJournal of the Royal Society Interface 12:20140875.https://doi.org/10.1098/rsif.2014.0875
-
Relative roles of ABCG5/ABCG8 in liver and intestineJournal of Lipid Research 56:319–330.https://doi.org/10.1194/jlr.M054544
-
The biological actions of dehydroepiandrosterone involves multiple receptorsDrug Metabolism Reviews 38:89–116.https://doi.org/10.1080/03602530600569877
-
Cushing's syndrome in a guinea pigVeterinary Record 160:878–880.https://doi.org/10.1136/vr.160.25.878
Decision letter
-
Jing-Dong Jackie HanReviewing Editor; Chinese Academy of Sciences, China
-
Jessica K TylerSenior Editor; Weill Cornell Medicine, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
This study analyzed 636 RNA-seq samples across 15 tissues, which is a massive effort and a rich resource. Changes in the hypothalamic-pituitary-adrenal stress axis is found to be associated with the extended life expectancy of reproductive vs. non-reproductive mole-rats, which is consistent with existing knowledge, whereas other associations are contrary to existing knowledge, such as the up-regulation of the IGF1/GH axis and several other anabolic processes in the long lived species. The authors should tune down their claims of causality such as "play a key role", which was not demonstrated without genetic perturbations, by only concluding associations, which are essentially what were described.
Decision letter after peer review:
Thank you for submitting your article "The HPA stress axis shapes aging rates in long-lived, social mole-rats" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Jessica Tyler as the Senior Editor. The reviewers have opted to remain anonymous.
While we cannot invite revision of the work in its current form, we acknowledge the potential and uniqueness of the data resource you have generated. If in the future you could improve the computational analyses to arrive at more definitive conclusions and include at least limited experimental validation of hypotheses you are advancing, we would be interested to receive such a thoroughly revised submission. We would treat it as a new submission, but we would try to recruit the same reviewers.
Summary:
Differential gene expression between breeding and non-breeding individuals of the two long-lived Fukomys species investigated here is a very interesting topic, and the scale of collected data is impressive. Unfortunately, however, both the data analyses and experimental support for the conclusions are rather superficial and descriptive. Given the size and rarity of the dataset, this could be an important resource for the aging research community, if data quality can be clearly demonstrated by robust and extensive computational analyses and experimental validation.
Without demonstrating causality, the authors should refrain from making assertive mechanistic claims, such as "the HPA stress axis shapes aging rates".
Please see the detailed reviews for specific concerns.
Reviewer #1:
This is a primarily descriptive study reporting gene expression differences between breeding and non-breeding individuals of two long-lived Fukomys species. These mole rat species are interesting for this type of analysis because it has been previously shown that breeders live ~2-fold longer than non-breeders. Thus, it may be possible to learn about mechanisms that determine longevity by studying differences between breeders and non-breeders. Although the manuscript is primarily descriptive, there is value in observational data sets such as this, which can be hypothesis generating and can spur future, more mechanistic studies. I found the presentation to be a bit of a "data dump" with many different comparisons of differentially expressed genes, GO terms, etc., but without a lot in terms true biological insight. As it was largely limited to gene expression, it is also a bit one-dimensional. This work will likely be of interest and value to scientists who study the comparative biology of aging and systems biology of aging, but in its current form seems somewhat incremental and probably will not capture the interest of a broad scientific audience or even the broader gerontology community.
My largest concerns with the content of this manuscript are related to interpretation of the data. While this type of comparative gene expression analysis can generate hypotheses, it cannot strongly support or refute causal relationships without additional experimentation. The authors repeatedly appear to interpret correlative observations with causation and make claims that overreach the data. The title itself is a good example of this where the authors claim "the HPA stress axis shapes aging rates". There is no direct evidence to support this claim or others similar in nature throughout the manuscript.
Another area of overinterpretation is with respect to the relevance of these findings as a test of the validity of prior aging theories or mechanisms. It seems very likely that this is a somewhat unique evolutionary case where upon the switch from non-breeder to breeder there is a dramatic rewiring of physiology at many levels (transcriptional, translational, post-translational, metabolome, epigenome, etc.). Simply because this switch does not match patterns seen in longevity interventions such as caloric restriction, reduced GH/IGF-1 signaling, or mTOR inhibition does not refute or call into question literature describing potential mechanisms for how those interventions act. It likely simply reflects that this is a different path to achieve longevity.
Specific comments:
– The title is problematic as described above.
– The first sentence in the Abstract is problematic as it implies causality between sexual activity/reproduction and longevity. Sexual activity/reproduction are associated with life expectancy. It is interesting to consider whether these two things could perhaps be uncoupled in this animal, as has been shown for reduced fecundity and longevity in invertebrate models.
– The phrase "oppose crucial findings" in the Abstract is also problematic. First, the word "crucial" does not seem to makes sense in this context. What makes them crucial? Second, it is intuitive and expected, indeed perhaps required, that reproduction would be associated with anabolic processes, and this study does not show causality between the observed changes in these processes and longevity, so it does not actually "oppose" the prior findings in genetic models with reduced IGF-1/GH signaling where lifespan is extended.
– I have a problem with the premise that it is possible to "confirm or falsify" results from cross-species or intra-species studies through the type of approach taken here as implied in the text. Confirming or falsifying results implies something about the quality of the prior data itself. I assume what the authors mean is confirming or falsifying the underlying hypotheses or assumptions. Even that is questionable, however, since the mechanisms by which longevity are determined could simply be different within species versus across species versus cases like this where you have dramatically different life expectancies in breeders versus non-breeders.
– I thought the discussion of the GH/IGF-1 results was fairly balanced, but I would encourage the authors to consider more deeply the within species versus across species observations in the literature. The evidence for reduced GH/IGF-1 increasing lifespan comes from within species studies and appears to hold true from worms to dogs and likely in humans as well, although the correlation between body size and lifespan is a bit more complicated in humans for obvious reasons. Within species, smaller individuals who have reduced GH/IGF-1 tend to live longer – that's a correlation. In worms, flies, and mice a reduction in growth signaling has been shown to be sufficient – and very likely causal – for enhanced longevity through genetic and pharmacological studies. The comment that these interventions are all performed during development is not exactly true – in mice rapamycin at least works in adulthood and even when only given transiently or intermittently in adulthood. Across species, larger species tend to live longer. Perhaps the mechanisms going from non-breeders to breeders more resemble the evolutionary longevity strategies that have been taken at the species level rather than the mechanism that appear to determine longevity within species. Personally, I don't see this as a contradiction or a controversy within the field.
– I find the "short-lived" versus "long-lived" species argument to be overly speculative and arbitrary. Interventions such as CR, mTOR inhibition, reduced IGF-1, etc. extend lifespan at least from worms to mice, which is a >50-fold difference in lifespan. Mice to people is ~30-fold. Compared to mice, worms are shorter-lived (by a fold difference metric) than mice are compared to humans.
– The authors state that breeders and non-breeders have "massively diverging aging rates", but I think it is important to keep in mind that differences in lifespan do not necessarily imply differences in aging rate. Especially going from the long-lived state to the short-lived state. Perhaps there is good evidence that functional and molecular declines and diseases/pathologies of aging are accelerated in the non-breeders and delayed in the breeders, which would support this assertion, but this is unclear from the manuscript.
– The phrase "unilaterally described as harmful" in the conclusion is simply not true. GH/IGF-1 signaling limits lifespan in worms, flies, and mice but none of the papers cited unilaterally claim that it is harmful. In fact, some of those same papers note that high GH/IGF-1 signaling often confers a selective advantage in terms of faster maturation and reproduction. So, at the species level, this is beneficial. Even at the individual level, high GH/IGF-1 may be associated with better outcomes in the wild where predation is a factor. High GF/IGF-1 in these species is detrimental (only?) in the context of aging/longevity in a relatively safe laboratory environment at the individual level.
Reviewer #2:
Sahm et al. have analyzed gene expression of 22 samples across 15 tissues of Fukomys mole rat, and with these resources, they try to explain why breeders and non-breeders have different lifespan in those species. While I was intrigued by the idea and design of the research, I found the bioinformatic approaches, as well as the overall result of the study, fell far short of being acceptable for publication. Below I highlight key aspects of the approach (and related conclusions, or lack thereof) that I believe represent serious issues. One positive note, again, the questions and data of this paper that I believe are highly interesting and important if the author can pursue it in a more focused and solid way.
1) First, in my opinion, the discussion of the paper is not well synthesized, and the content does not help the reader for the aim of study and of what is discussed in the Title (e.g., HPA) and Abstract. While the entire discussion fail to build logics between HPA stress and aging in mole rats, it says story of intervention and hypothesis testing instead, with the last of results immediately jumps into comparisons positive selection genes and DEG (It is unclear why it is relevant even one gene is found under positive selection).
2) Given hundreds of samples have been sequenced in the paper, there is no extensive examination of batch effect, which could ultimately, in its present form, put all the results (e.g., DEG analysis, pathways enrichment, multifactors analysis) and conclusions at high risk. Particularly, the PCA analysis has indicated the species, tissue, and the combination of both variables accounted for 98.4 % of the total variance in the data set, it becomes more important to know how the author have organized the sequencing strategy.
3).Unsophisticated use of enrichment analysis on pathway likely lead to misleading conclusions – a large portion of the analyses presented use an unreasonable approach to identify genes with expression shift across tissues and species to make conclusions that I believe are largely misleading about genes important for longevity. The authors identify pathways that have an overall (gene-wide) changes in gene expression as p value looks like an potential indication of expression shift. While this approach MAY be lucky enough to catch a few of these “true positives", the VAST majority of what it will identify will be pathways with global accumulated changes (which is what the point would be), but rather small pathways or pathways fully of with wired p values or undetectable expression (and thus perhaps some of the least important genes). This overly coarse approach seems to me unreasonable and publishable for this purpose. Then, the extension of these approaches to pathways further muddies the waters of any discussions to the extent that I believe it is largely nonsense that overwhelms the results. The approaches for detecting affected global pathways have been the subject of a very large body of literature, and there are well developed hierarchical/empirical models for testing these hypotheses that the authors have completely and inexcusably ignored. And, of course, it is strong encouraged the author could formulate an new model/index to re-evaluate pathways analysis as the data and experimental design is unique to other study.
4) Last but not the least, small KEGG set, i.e., KEGG that with very few gene could be also confusing such whole categories analysis. The author should check this potential bias using random sampling “pseudo-KEGG set” with same gene number and/or identify a threshold to filer the each true KEGG categories that considered.
Reviewer #3:
The manuscript by Sahm et al. describes the transcriptomic comparison of breeders and non-breeders in two species of mole rats from genus Fukomys. The remarkable aspect of Fukomys mole rats is that the breeders live significantly longer than workers. The authors produced new breeder couples by pairing animals from different family groups and then compared their transcriptomes to non-breeders of the same age from the original colonies.
There were very few differences identified between breeders and non-breeders. Most transcriptomic changes were confined to gonads and endocrine glands. This is somewhat unsurprising and these organs become active for breeding. The pathways that became activated were related to ribosome biogenesis, protein translation, and MYC signaling, which all reflect physiological activation of these tissues in breeders.
Interestingly, some activation of these signaling pathways was also observed in non-gonadal tissues, which contradicts the common dogma that downregulation of these pathways is associated with lifespan extension. This is a remarkable result that suggests that short-lived model organisms may not correctly reflect signaling effects required for longevity in long-lived species. The authors also point out to higher glucocorticoid levels in workers and speculate that may be showing Cushing's like syndrome.
Overall, this is an important study of high interest.
Items to address:
1) It is not clear how long were the animals maintained in a breeder status prior to analysis.
2) Figure 3A, the rational for the comparison of changes in young breeder versus non breeder animals to changes that occur with age over time is unclear.
3) Figure 4, were the changes driven by gonads mainly? Hoe many were non-gonadal? For example, in Figure 2, muscle showed 0 DEGs by status, but in pathway analysis we see upregulation of pathways. Please explain.
4) Table 2: TOR is downregulated in breeders while ribosome processed are upregulated, please explain.
5) The finding of Cushing's syndrome needs more support. Please consider comparing Cushing's related transcriptome changes as a whole to the changes observed in mole rats, otherwise this conclusion may need to be toned down.
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Thank you for submitting your article "Slower aging due to sexual activity in mole-rats is associated with transcriptional changes in the HPA stress axis" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Jessica Tyler as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Vera Gorbunova (Reviewer #3).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity, data analysis and presentation.
We are glad to see that you smoothed your HPA stress hypothesis according to our suggestions, but further revisions are still necessary to address our previous concerns about data quality and robustness of results have and remove the overinterpretation from the title.
Specifically you need to:
1) Control the batch effect using established packages, such as "BatchQC" or "BEclear";
2) Detect DEG using limma with voom to account for large variances and compare the results with DESeq2's output (see PMID: 23497356);
3) Redo the pathway enrichment with changes of expression values or other static indicators (see PMID: 23625889 and PMID: 32515539 for examples) but not p-value; to remove the bias and validate the HPA hypothesis.
4) Reconsider the title, which is still an over-interpretation of the data and therefore needs to be fixed as follows. The words "slower aging" imply that the same biological aging process is occurring in the breeders and non-breeders, which is not actually clear from the data. Indeed, it seems quite plausible that the non-breeders die prematurely due to too much stress (as you speculate). Whether the rate of aging is different or not remains unclear. Replacing "slower aging" with "Increased longevity" or something like that would address this issue.
https://doi.org/10.7554/eLife.57843.sa1Author response
Reviewer #1:
This is a primarily descriptive study reporting gene expression differences between breeding and non-breeding individuals of two long-lived Fukomys species. These mole rat species are interesting for this type of analysis because it has been previously shown that breeders live ~2-fold longer than non-breeders. Thus, it may be possible to learn about mechanisms that determine longevity by studying differences between breeders and non-breeders. Although the manuscript is primarily descriptive, there is value in observational data sets such as this, which can be hypothesis generating and can spur future, more mechanistic studies. I found the presentation to be a bit of a "data dump" with many different comparisons of differentially expressed genes, GO terms, etc., but without a lot in terms true biological insight. As it was largely limited to gene expression, it is also a bit one-dimensional. This work will likely be of interest and value to scientists who study the comparative biology of aging and systems biology of aging, but in its current form seems somewhat incremental and probably will not capture the interest of a broad scientific audience or even the broader gerontology community.
My largest concerns with the content of this manuscript are related to interpretation of the data. While this type of comparative gene expression analysis can generate hypotheses, it cannot strongly support or refute causal relationships without additional experimentation. The authors repeatedly appear to interpret correlative observations with causation and make claims that overreach the data. The title itself is a good example of this where the authors claim "the HPA stress axis shapes aging rates". There is no direct evidence to support this claim or others similar in nature throughout the manuscript.
Another area of overinterpretation is with respect to the relevance of these findings as a test of the validity of prior aging theories or mechanisms. It seems very likely that this is a somewhat unique evolutionary case where upon the switch from non-breeder to breeder there is a dramatic rewiring of physiology at many levels (transcriptional, translational, post-translational, metabolome, epigenome, etc.). Simply because this switch does not match patterns seen in longevity interventions such as caloric restriction, reduced GH/IGF-1 signaling, or mTOR inhibition does not refute or call into question literature describing potential mechanisms for how those interventions act. It likely simply reflects that this is a different path to achieve longevity.
We thank you for your outright comment that we have made too strong statements in light of the evidence presented. Your review was the first we received for this manuscript, and considerably helped us to realize the pitfalls of over-interpretation. We implemented your constructive criticism in two ways. First, we have, in accordance with your suggestions, toned down or deleted statements and generally made the correlative character of the facts presented more transparent. In addition, we are now discussing the previously collected and new, complementary evidence for the HPA stress hypothesis in more detail and with more differentiation than before.
Specific comments:
– The title is problematic as described above.
We agree and use a new title that reflects the correlative character of the underlying evidence: Slower aging due to sexual activity in mole-rats is associated with transcriptional changes in the HPA stress axis
– The first sentence in the Abstract is problematic as it implies causality between sexual activity/reproduction and longevity. Sexual activity/reproduction are associated with life expectancy. It is interesting to consider whether these two things could perhaps be uncoupled in this animal, as has been shown for reduced fecundity and longevity in invertebrate models.
Yes, we modified the Abstract according to your suggestions. Specifically, we made the following replacements:
“Sexual activity and/or reproduction are associated with a doubling of life expectancy in the long-lived rodent genus Fukomys.”
“This analysis suggests that changes in the regulation of the hypothalamic-pituitary-adrenal stress axis play a key role regarding the extended life expectancy of reproductive vs. non-reproductive mole-rats. This is substantiated by a corpus of independent evidence.”
– The phrase "oppose crucial findings" in the Abstract is also problematic. First, the word "crucial" does not seem to makes sense in this context. What makes them crucial? Second, it is intuitive and expected, indeed perhaps required, that reproduction would be associated with anabolic processes, and this study does not show causality between the observed changes in these processes and longevity, so it does not actually "oppose" the prior findings in genetic models with reduced IGF-1/GH signaling where lifespan is extended.
We agree and modified the Abstract according to your suggestions:
“On the other hand, several of our results are not consistent with knowledge about aging of short-lived model organisms.”
Moreover, we added to the Introduction:
“Nevertheless, both Fukomys and naked mole-rats can be seen as counter-examples to the disposable soma theory of aging, which assumes that the aging process is linked to an evolutionary tradeoff in which the preservation of the organism is essentially sacrificed to reproductive success (Kirkwood, 1977).”
– I have a problem with the premise that it is possible to "confirm or falsify" results from cross-species or intra-species studies through the type of approach taken here as implied in the text. Confirming or falsifying results implies something about the quality of the prior data itself. I assume what the authors mean is confirming or falsifying the underlying hypotheses or assumptions. Even that is questionable, however, since the mechanisms by which longevity are determined could simply be different within species versus across species versus cases like this where you have dramatically different life expectancies in breeders versus non-breeders.
Yes, we agree and removed the respective lines from the manuscript.
– I thought the discussion of the GH/IGF-1 results was fairly balanced, but I would encourage the authors to consider more deeply the within species versus across species observations in the literature. The evidence for reduced GH/IGF-1 increasing lifespan comes from within species studies and appears to hold true from worms to dogs and likely in humans as well, although the correlation between body size and lifespan is a bit more complicated in humans for obvious reasons. Within species, smaller individuals who have reduced GH/IGF-1 tend to live longer – that's a correlation. In worms, flies, and mice a reduction in growth signaling has been shown to be sufficient – and very likely causal – for enhanced longevity through genetic and pharmacological studies.
That is, if we see it correctly, so far completely in line with the manuscript, which describes that today’s literature based on evidence in different species predominantly assumes that increased IGF1/GH signaling shortens life and decreased signaling prolongs it (Table2). If we have misunderstood you on this point and in your opinion important arguments or literature evidence is still missing, we would be happy to add it.
The comment that these interventions are all performed during development is not exactly true – in mice rapamycin at least works in adulthood and even when only given transiently or intermittently in adulthood.
Yes, we have supplemented and qualified the statement as follows:
“However, there are some exceptions: for example, late treatments with rapamycin prolong the life of mice almost as much as those that are started early (Johnson et al., 2013). Therefore, it is still under debate…”
Across species, larger species tend to live longer. Perhaps the mechanisms going from non-breeders to breeders more resemble the evolutionary longevity strategies that have been taken at the species level rather than the mechanism that appear to determine longevity within species. Personally, I don't see this as a contradiction or a controversy within the field.
We understand your concerns, however, here we are making the case that there is no need to postulate two general strategies/mechanisms for longevity determination at the inter- and intraspecies levels. For example, in both intra- and interspecies comparisons, the growth rate is clearly anticorrelated with lifespan (https://pubmed.ncbi.nlm.nih.gov/11868658/, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3574304/, https://pubmed.ncbi.nlm.nih.gov/12954479/). In respect to the well accepted positive correlation of body weight and longevity this means, that large long-lived species grow relative longer than small short-lived ones. Noteworthy, in mammals maximum life span is even stronger correlated to time to maturity than either to adult weight or growth rate (see, e.g., Figure 2 of https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4406664/). In respect to the manuscript, however, we believe that such a speculation and discussion would extend too far from the topics of our work.
– I find the "short-lived" versus "long-lived" species argument to be overly speculative and arbitrary. Interventions such as CR, mTOR inhibition, reduced IGF-1, etc. extend lifespan at least from worms to mice, which is a >50-fold difference in lifespan. Mice to people is ~30-fold. Compared to mice, worms are shorter-lived (by a fold difference metric) than mice are compared to humans.
We would like to point out that the concept of comparing short- versus long-lived has been applied successfully (see e.g. Austad SN, Comparative biology of aging, PMID: 19223603 and Gorbunova V, Comparative genetics of longevity and cancer, PMID: 24981598). Particular in species with a similar biology like mammals, where a strong correlation exists between body-weight and longevity, it is quite rational and powerful. Based on this correlation, short- and long-lived outliers can be objectively identified. Classification of Fukomys as a long-lived outlier is a basic motivation of our research efforts.
However, we understand your concern with respect to the fact that the mentioned treatments work in principle consistently in the same direction in a wide range of species. Therefore, we have shortened the paragraph considerably and concentrated on what we consider a more stringent argument:
“A third point is the question of the transferability of knowledge obtained in one species to other species. Most insights into current aging research originate from very short-lived model organisms (Table 2). […]Therefore, it may also be possible that those gene expression patterns caused by our lifespan-extending intervention in Fukomys mole-rats highlight particularly the differences in aging mechanisms between short-lived and long-lived species.”
– The authors state that breeders and non-breeders have "massively diverging aging rates", but I think it is important to keep in mind that differences in lifespan do not necessarily imply differences in aging rate. Especially going from the long-lived state to the short-lived state. Perhaps there is good evidence that functional and molecular declines and diseases/pathologies of aging are accelerated in the non-breeders and delayed in the breeders, which would support this assertion, but this is unclear from the manuscript.
Sorry for being insufficiently clear in this respect in the submitted manuscript. We agree of course that differences in lifespan do not necessarily imply differences in aging rates, as organisms can die for various reasons, many of which have little to do with physiological aging, and as a matter of fact, we have not actually measured aging rates in our study species. However, the survival data for Fukomys non-breeders and breeders (https://pubmed.ncbi.nlm.nih.gov/16488857/, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3076438/) have been obtained in protected captive environments where both groups lived essentially under the same conditions, and where starvation, predation or dispersal could be excluded as mortality causes. Accidental death cases were very rare. Kaplan-Meier-survival curves and analyses of mortality rates indicate that both groups show actuarial senescence, but with obviously earlier onset of aging and significantly higher annual mortality hazards in non-breeders across all age classes (https://pubmed.ncbi.nlm.nih.gov/16488857/, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3076438/, https://cloud.leibniz-fli.de/index.php/s/d2Wob9jFAeyRHQZ [accepted at Philos. Trans. R. Soc. B]). Hence, we are confident to say that the longer mean and absolute lifespans of breeders result from a combination of delayed onset of aging and henceforth lower age-specific mortality rates in Fukomys breeders as compared to non-breeders. As, in the absence of extrinsic mortality causes, both parameters are parameters of aging, we consider it justified to speak of “slower aging breeders” and/or “faster aging non-breeders” whenever appropriate.
However, we modified the text to avoid misunderstandings:
“… compared naturally occurring cohorts of species with massively diverging longevities.”
plus numerous other passages where the phrase “aging-rate” has been substituted likewise.
– The phrase "unilaterally described as harmful" in the conclusion is simply not true. GH/IGF-1 signaling limits lifespan in worms, flies, and mice but none of the papers cited unilaterally claim that it is harmful. In fact, some of those same papers note that high GH/IGF-1 signaling often confers a selective advantage in terms of faster maturation and reproduction. So, at the species level, this is beneficial. Even at the individual level, high GH/IGF-1 may be associated with better outcomes in the wild where predation is a factor. High GF/IGF-1 in these species is detrimental (only?) in the context of aging/longevity in a relatively safe laboratory environment at the individual level.
Yes, we phrased that wrongly. In our discussion, we even deal with the possible evolutionary advantages of increased IGF1/GH signalling in non-breeders in such a tradeoff view.
“In particular, this applies to the role of the GH/IGF1 axis, which has in recent years been predominantly described as having a negative impact on lifespan (Lopez-Otin et al., 2013; Pitt and Kaeberlein, 2015; Bartke, 2017).”
Reviewer #2:
Sahm et al. have analyzed gene expression of 22 samples across 15 tissues of Fukomys mole rat, and with these resources, they try to explain why breeders and non-breeders have different lifespan in those species. While I was intrigued by the idea and design of the research, I found the bioinformatic approaches, as well as the overall result of the study, fell far short of being acceptable for publication. Below I highlight key aspects of the approach (and related conclusions, or lack thereof) that I believe represent serious issues. One positive note, again, the questions and data of this paper that I believe are highly interesting and important if the author can pursue it in a more focused and solid way.
1) First, in my opinion, the discussion of the paper is not well synthesized, and the content does not help the reader for the aim of study and of what is discussed in the Title (e.g., HPA) and Abstract. While the entire discussion fail to build logics between HPA stress and aging in mole rats, it says story of intervention and hypothesis testing instead, with the last of results immediately jumps into comparisons positive selection genes and DEG (It is unclear why it is relevant even one gene is found under positive selection).
In response to the comment, we now discuss the hypotheses, results and verifications concerning the HPA axis in greater detail. Furthermore, we have shortened those parts of the manuscript describing the up-regulation of various anabolic components such as the IGF1/GH axis, which are predominantly associated with a shortening of lifespan (see answer 1.10 to Reviewer #1). You also assert that we jump “immediately” into a comparison of positively selected genes and DEGs. This is very surprising to us, since this issue is addressed only once in the entire manuscript: in the last of about fifteen paragraphs of our Discussion – by far not the main pillar of our manuscript.
2) Given hundreds of samples have been sequenced in the paper, there is no extensive examination of batch effect, which could ultimately, in its present form, put all the results (e.g., DEG analysis, pathways enrichment, multifactors analysis) and conclusions at high risk. Particularly, the PCA analysis has indicated the species, tissue, and the combination of both variables accounted for 98.4 % of the total variance in the data set, it becomes more important to know how the author have organized the sequencing strategy.
We have added the rational one for our sequencing strategy as well as a new supplement table that details the organization of sequencing experiments together. In brief, we have taken great care to distribute samples of breeders and non-breeders across the sequencing runs to avoid such systematic errors.
“ We have taken care to concentrate the samples of a tissue on a few sequencing runs and to create a balance of breeders and non-breeders within each run (Supplementary file 1J).”
3) Unsophisticated use of enrichment analysis on pathway likely lead to misleading conclusions – a large portion of the analyses presented use an unreasonable approach to identify genes with expression shift across tissues and species to make conclusions that I believe are largely misleading about genes important for longevity. The authors identify pathways that have an overall (gene-wide) changes in gene expression as p value looks like an potential indication of expression shift. While this approach MAY be lucky enough to catch a few of these “true positives", the VAST majority of what it will identify will be pathways with global accumulated changes (which is what the point would be), but rather small pathways or pathways fully of with wired p values or undetectable expression (and thus perhaps some of the least important genes). This overly coarse approach seems to me unreasonable and publishable for this purpose. Then, the extension of these approaches to pathways further muddies the waters of any discussions to the extent that I believe it is largely nonsense that overwhelms the results. The approaches for detecting affected global pathways have been the subject of a very large body of literature, and there are well developed hierarchical/empirical models for testing these hypotheses that the authors have completely and inexcusably ignored. And, of course, it is strong encouraged the author could formulate an new model/index to re-evaluate pathways analysis as the data and experimental design is unique to other study.
Since we have reason to believe that this could be a methodological misunderstanding (see below), we would like to take the opportunity to explain the basic aspects of our approach again with references to the relevant literature. In order to combine different p-values in a threshold-free way into a test statistic, we have used the Fisher’s method.In fact, Fisher's method is one of the most widely used methods in statistics (https://academic.oup.com/biomet/article/105/1/239/4788722). The provided reference also illustrates the use of weights. This is an appropriate approach to account for different levels of gene expression in the individual tissues in the context of cross-tissue analysis. In addition, the method is also specifically recommended for the identification of differentially expressed pathways in gene expression analyses (e.g. https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0012693,
https://academic.oup.com/biomet/article/105/1/239/4788722, https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0041018, https://www.tandfonline.com/doi/abs/10.1080/19466315.2014.888013, https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-14-368, https://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-8-96, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3632109/).
“Simulation results demonstrated that overall Fisher's method and the global model with random effects have the highest power for a wide range of scenarios” (first linked article)
Unfortunately, you did not specify any instance of the "large body of literature" that you consider better suited for threshold-free, weighted enrichment analyses of gene expression comparisons across several tissues.
Strategies used to determine the combined p-value based on Fisher's method, are described, for example, in the first, third and seventh linked article. In general Fisher’s method assumes independence of m individual test statistics and uniformity of the associated p-values under the null. Consequently, the test statistic obtained with Fisher’s method is approximately chi square distributed with 2m degrees of freedom. In gene expression analyses, however, the assumption of independence may be violated (https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0012693, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3632109/). Therefore, a Monte-Carlo approach is often used instead to determine the p-values with empirically determined null distributions.
We have taken your comments as an opportunity to describe the procedure in more detail:
“If the underlying test statistics are independent, a combined p-value may be determined by Fishers’s method using an χ2-distribution. […] The combined p-values for our hypothesis tests can now be obtained by determining the probability of a test statistic in the empirically estimated null-distribution.”
Regarding the assertions about validity and integrity of our results we would like to refer to our cross-tissue-level analysis. We are able to falsify the hypothesis that our method preferentially identifies poorly expressed genes/pathways or small pathways:
“We checked whether this method may have a detection bias with regard to expression level of genes/pathways or size of pathways at the cross-tissue level. For this, we looked at the correlations of the respective p-values with gene expression, pathway expression, and pathway size (Figure 4—figure supplement 8-10). In summary, we find a slightly enhanced sensitivity for highly expressed pathways at the cross-tissue level.”
Furthermore, a preference for pathways in terms of size can already be excluded on the basis of the empirical determination of the null distribution described above. There, the number of elements in the respective pathway is taken into account (see above).
4) Last but not the least, small KEGG set, i.e., KEGG that with very few gene could be also confusing such whole categories analysis. The author should check this potential bias using random sampling “pseudo-KEGG set” with same gene number and/or identify a threshold to filer the each true KEGG categories that considered.
Confusingly, the description of this “pseudo-KEGG set” approach, which you suggest as a possible solution, matches exactly the method used in our manuscript. Given this last statement, we believe that the fundamental criticism our method may have been the result of a misunderstanding. We sincerely hope to have resolved this with our explanations.
Reviewer #3:
The manuscript by Sahm et al. describes the transcriptomic comparison of breeders and non-breeders in two species of mole rats from genus Fukomys. The remarkable aspect of Fukomys mole rats is that the breeders live significantly longer than workers. The authors produced new breeder couples by pairing animals from different family groups and then compared their transcriptomes to non-breeders of the same age from the original colonies.
There were very few differences identified between breeders and non-breeders. Most transcriptomic changes were confined to gonads and endocrine glands. This is somewhat unsurprising and these organs become active for breeding. The pathways that became activated were related to ribosome biogenesis, protein translation, and MYC signaling, which all reflect physiological activation of these tissues in breeders.
Interestingly, some activation of these signaling pathways was also observed in non-gonadal tissues, which contradicts the common dogma that downregulation of these pathways is associated with lifespan extension. This is a remarkable result that suggests that short-lived model organisms may not correctly reflect signaling effects required for longevity in long-lived species. The authors also point out to higher glucocorticoid levels in workers and speculate that may be showing Cushing's like syndrome.
Overall, this is an important study of high interest.
Items to address:
1) It is not clear how long were the animals maintained in a breeder status prior to analysis.
Yes, we added the detailed information for every animal in the study to the supplement and summarized it in the manuscript as follows:
“In total, females/males selected as breeders spent on average 513/353 (F. mechowii) and 1095/997 (F. micklemi) days in this state before sampling (Supplementary file 1D,E).”
2) Figure 3A, the rational for the comparison of changes in young breeder versus non breeder animals to changes that occur with age over time is unclear.
We thank you for the question and introduced the two corresponding paragraphs with the rationale:
“To gain insights into the molecular basis of the bimodal aging pattern in Fukomys, it would be of course interesting to compare sufficiently large cohorts of old breeders and non-breeders both at similar chronological and biological ages. But as old animals, particularly of a long-lived species, are extremely valuable due to obvious logistical, timely and financial reasons, such an approach was impossible for us. As an alternative, we decided to study young breeders and non-breeders, that are much easier to obtain, and hypothesized that the status change marks the beginning of a slowdown in the aging process.”
“Since the change from non-breeder to breeder status apparently marks the beginning of a slowdown in the aging process, we first wanted to find out whether and where there are intersections of reproductive status DEGs with those whose expression level is known to change during aging. Therefore, we determined overlaps by using…”
3) Figure 4, were the changes driven by gonads mainly? Hoe many were non-gonadal? For example, in Figure 2, muscle showed 0 DEGs by status, but in pathway analysis we see upregulation of pathways. Please explain.
The enrichment analysis is not based on the differentially expressed genes only but was performed in a threshold-free manner. This was done in order not to lose the ability to detect differences at the pathway level for which the statistical power at the single gene level is insufficient due to the available sample size. An enrichment analysis that would focus only on DEGs weights all genes with a p-value below a certain threshold (e.g. 0.05) with 1 and all genes above this threshold with 0. Instead, we have adopted an approach that takes each gene of a pathway into account in the test statistics – weighted according to its p-value, providing test sensitivity in tissues with no or few detected DEGs. We have addressed this approach in the Results section of our manuscript as follows:
“For this, we used Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Kanehisa et al., 2017), and Molecular Signatures Database (MSigDB) hallmarks (Liberzon et al., 2015) as concise knowledge bases. We used an established method that combines all p-values of genes in a given pathway in a threshold-free manner. The advantage of this approach is that it bundles the p-values from test results of individual gene expression differences at the level of pathways (see Materials and methods section for details).”
To avoid confusion when viewing the figure in isolation or when the paragraph above has been overlooked, we have also added the following to the figure description:
“This enrichment analysis was carried out threshold-free, which is why tissues without differentially expressed individual genes (see Figure 2B) can also show differentially expressed pathways.”
4) Table 2: TOR is downregulated in breeders while ribosome processed are upregulated, please explain.
It should be noted that the mTOR downregulation, as shown in Table 2, is limited to the adrenal gland only. Nevertheless, some interesting questions do indeed arise with regard to mTOR regulation due to the scenario on the one hand (massive life span differences) and the concrete results on the other hand (mTOR largely unchanged, but protein biosynthesis up-regulated). Although we cannot conclusively answer these questions, we will discuss them in a new paragraph:
“An important positive regulator of ribosome biogenesis and protein synthesis is mTOR (Johnson et al., 2013). Based on this, one could expect an upregulation of the corresponding gene in our scenario. At the same time, however, the inhibition of mTOR is one of the best documented life-prolonging interventions from invertebrates to mammals (Table 2). Given these conflicting premises, we find mTOR in almost all Fukomys tissues as not significantly altered; the exception is the adrenal gland, where mTOR is significantly downregulated in breeders. It is obvious that, during evolution of Fukomys mole-rats, both the extension of breeders’ lifespan and an increase in their anabolic processes provided fitness benefits. Consequently, in the underlying organismic framework mTOR may have acquired expression patterns and functions, different from those in organisms studied before. Thus, in the future it may be worth to study Fukomys mTOR biology in more detail, particularly in respect to its potential role in naturally evolved ways towards lifespan extension and healthy aging.”
5) The finding of Cushing's syndrome needs more support. Please consider comparing Cushing's related transcriptome changes as a whole to the changes observed in mole rats, otherwise this conclusion may need to be toned down.
Many thanks for this really good suggestion! We were able to find one single data set that performs an RNA-seq comparison of Cushing vs. controls. We have correlated this data on a global level with our non-breeder/breeder comparison and have extended the Materials and methods and the last section of the Results accordingly. We discuss the new results and the limitations of this comparison together with our other hypothesis tests about the HPA axis and Cushing in a new section. The part corresponding to your suggestion is this:
“Furthermore, according to our hypothesis, the expression changes in the comparison of control individuals against patients with Cushing's syndrome should correspond to those of the comparison of breeders vs. non-breeders. A global comparison of fold-changes from human subcutaneous adipose tissue (Hochberg et al., 2015) against our cross-tissue data showed a significant correlation (R=0.11, p=6.6*10-36). A very similar result was obtained by comparing the human data with our data from skin as the tissue closest to subcutaneous adipose tissue that was examined in our study (R=0.11, p=1.1*10‑40). Despite the clear significance of these correlations, the reported values of the correlation coefficients could be considered to be only modest. It should, however, be taken into account that two evolutionarily widely separated species and not exactly matching tissues were compared.”
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
We are glad to see that you smoothed your HPA stress hypothesis according to our suggestions, but further revisions are still necessary to address our previous concerns about data quality and robustness of results have and remove the overinterpretation from the title.
Specifically you need to:
1) Control the batch effect using established packages, such as "BatchQC" or "BEclear";
As suggested, we analyzed potential batch effects systematically using BatchQC, added a respective source data set 2h with all reports for full transparency, summarized them in a paragraph in the Control analyses section and referenced the issue in the first paragraph of the Results.
“The number of samples collected and processed for this work required splitting them into different batches for sequencing. To analyze the degree this procedure potentially biased our results, we systematically investigated potential batch effects. Our results suggest that the sequencing strategy has little effect on the outcomes reported in the following (see Control analyses).”
“Control analyses: For practical reasons, samples collected for this project over several years were sequenced progressively in different batches. Therefore, we analyzed our data and sample scheme across all tissues and both species using BatchQC (Manimaran et al., 2016). In summary, the extensive search did not reveal any possible batch effects in the dataset from F. micklemi, although they may not be ruled out completely for F. mechowii (Source data 2h). A source for the minor batch effects in F. mechowii may be attributable to the repeated sequencing of the same sample(s) in different sequencing experiments – an approach avoided during generation of the F. micklemi dataset. However, because we followed a robust multifactorial analysis, in which differences must occur consistently across both species to be considered (see “Differentially expressed genes analysis” in Materials and methods), false-positive results due to batch effects are substantially reduced.”
2) Detect DEG using limma with voom to account for large variances and compare the results with DESeq2's output (see PMID: 23497356);
As suggested, we searched for DEGs using limma/voom. The respective p- and FDR-values obtained with this method are juxtaposed to the respective DESeq2 values in the tables in Source data 1e-g. We compare DESeq2 and limma/voom results across single tissues in the new Figure 2—figure supplement 2 and summarized this cross-validation in the following paragraph as part of the new Control analyses section.
“Control analyses: In addition to detecting DEGs with DESeq2 (Love et al., 2014), we also examined the central contrast, i.e., the juxtaposition of workers and breeders, using the combination of limma (Smyth, 2004) and voom (Law et al., 2014). The absolute numbers of DEGs found in each tissue, and thus the highly uneven distribution, were nearly identical for all tissues except thyroid (Figure 2—figure supplement 2). For the latter, however, a closer look shows that 85% of the DESeq2 DEGs (threshold FDR <0.05) are confirmed if the FDR threshold for limma/voom is raised from 0.05 to 0.1. If an FDR threshold of 0.05 is applied for both approaches across all tissues, a total of 58% of the DESeq2-DEGs are confirmed by limma/voom. Notably, all genes mentioned in the manuscript, including the candidates from Table 1, are among this set. Using a less strict threshold of 0.1 for limma/voom, 91% of DESeq2-DEGs are confirmed.”
Discussion: “For further validation, we applied the combination of limma (Smyth, 2004) and voom (Law et al., 2014) as an alternative DEG detection method to DESeq2 to account for potential dispersion problems in our data. All single DEGs mentioned in the manuscript were confirmed by this validation. This specifically includes all genes supporting our HPA hypothesis, e.g., SULT2A1 (liver), MC2R (adrenal gland), CYP11A1 (adrenal gland and ovary), GH (pituitary gland) and IGF1 (ovary). In a similar way, we checked the main findings obtained with our standard bioinformatic p-value based enrichment approach (Fridley et al., 2010; Evangelou et al., 2012; Poole et al., 2016), using a second method that aggregates fold-changes instead. In particular, this step confirmed our initial finding of an enrichment of genes involved in the steroid hormone biosynthesis and ribosome metabolism as a central part of up-regulated anabolism in breeders, consistent with our HPA hypothesis (Figure 4—figure supplement 11) (see Control analyses).”
3) Redo the pathway enrichment with changes of expression values or other static indicators (see PMID: 23625889 and PMID: 32515539 for examples) but not p-value; to remove the bias and validate the HPA hypothesis.
The requested analysis was added to the Control analyses section and illustrated in a new supplemental figure. We refer to this additional check in our Results section. However, we would like to emphasize at the aggregation of p-values for the purpose of pathway enrichment in the context of gene expression studies has been used or even recommended in numerous previous studies (e.g. https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0012693, https://www.tandfonline.com/doi/abs/10.1080/19466315.2014.888013, https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-14-368, https://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-8-96, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3632109/).
Results: “As standard approach, we used an established method (Fridley et al., 2010; Evangelou et al., 2012; Poole et al., 2016) combining all p-values of genes in a given pathway in a threshold-free manner (Figure 4). […] In addition, we applied a second enrichment method that aggregates fold-changes instead (Figure 4—figure supplement 11; for a comparison of the results of both approaches, see Control analyses section).”
Discussion (see issue 2, above).
Control analyses: “In order to avoid potential biases by using solely p-values as indicator of differential gene expression, we additionally resorted to fold-changes as a test statistic for our pathway enrichment analyses. This exercise roughly confirms half of the pathways/hallmarks identified as differentially expressed at the cross-tissue level. Notably the set includes, consistent with the HPA hypothesis, steroid hormone biosynthesis and ribosomes as a central part of up-regulated anabolism in breeders (Figure 4—figure supplement 11). Furthermore, myogenesis, the P53 pathway, and coagulation were confirmed as cross-tissue differentially expressed by this approach. On the other hand, proteasome and oxidative phosphorylation gene sets, were not identified by the fold-change method.”
We have extended the corresponding Materials and methods section with the respective formulas. Also, other analyses pertaining to pathway enrichments were moved to the “Control analyses” section.
4) Reconsider the title, which is still an over-interpretation of the data and therefore needs to be fixed as follows. The words "slower aging" imply that the same biological aging process is occurring in the breeders and non-breeders, which is not actually clear from the data. Indeed, it seems quite plausible that the non-breeders die prematurely due to too much stress (as you speculate). Whether the rate of aging is different or not remains unclear. Replacing "slower aging" with "Increased longevity" or something like that would address this issue.
As suggested, we changed the title to “Increased longevity due to sexual activity in mole-rats is associated with transcriptional changes in the HPA stress axis”.
https://doi.org/10.7554/eLife.57843.sa2Article and author information
Author details
Funding
Deutsche Forschungsgemeinschaft (PL 173/8-1)
- Matthias Platzer
Deutsche Forschungsgemeinschaft (DA 992/3-1)
- Philip Dammann
Deutsche Forschungsgemeinschaft (Research Training Group 1739)
- Magdalena Staniszewska
Wiedenfeld-Stiftung, Stiftung Krebsforschung Duisburg
- Magdalena Staniszewska
Joachim Herz Stiftung
- Arne Sahm
Leibniz Association (Open Access Fund)
- Arne Sahm
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Ivonne Görlich, Christiane Vole, and Klaus Huse for excellent assistance in the preparation of biological samples. We thank Konstantin Riege for fruitful discussions on the analysis of the data. We thank Debra Weih and Flo Witte for proofreading the manuscript.
Ethics
Animal experimentation: Animal housing and tissue collection were compliant with national and state legislation (breeding allowances 32-2-1180-71/328 and 32-2-11-80-71/345; ethics/animal experimentation approval 84-02.04.2013/A164, Landesamt für Natur-, Umwelt- und Verbraucherschutz Nordrhein-Westfalen). Before sampling, animals were anaesthetized with ketamine combined with xylazine (Garcia Montero et al. 2015). Every effort was made to minimize suffering.
Senior Editor
- Jessica K Tyler, Weill Cornell Medicine, United States
Reviewing Editor
- Jing-Dong Jackie Han, Chinese Academy of Sciences, China
Version history
- Received: April 15, 2020
- Accepted: March 15, 2021
- Accepted Manuscript published: March 16, 2021 (version 1)
- Version of Record published: March 31, 2021 (version 2)
- Version of Record updated: July 21, 2021 (version 3)
Copyright
© 2021, Sahm et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 2,481
- Page views
-
- 254
- Downloads
-
- 16
- Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Chromosomes and Gene Expression
- Neuroscience
Mathys et al. conducted the first single-nucleus RNA-seq (snRNA-seq) study of Alzheimer’s disease (AD) (Mathys et al., 2019). With bulk RNA-seq, changes in gene expression across cell types can be lost, potentially masking the differentially expressed genes (DEGs) across different cell types. Through the use of single-cell techniques, the authors benefitted from increased resolution with the potential to uncover cell type-specific DEGs in AD for the first time. However, there were limitations in both their data processing and quality control and their differential expression analysis. Here, we correct these issues and use best-practice approaches to snRNA-seq differential expression, resulting in 549 times fewer DEGs at a false discovery rate of 0.05. Thus, this study highlights the impact of quality control and differential analysis methods on the discovery of disease-associated genes and aims to refocus the AD research field away from spuriously identified genes.
-
- Chromosomes and Gene Expression
- Genetics and Genomics
Spermatogenesis in the Drosophila male germline proceeds through a unique transcriptional program controlled both by germline-specific transcription factors and by testis-specific versions of core transcriptional machinery. This program includes the activation of genes on the heterochromatic Y chromosome, and reduced transcription from the X chromosome, but how expression from these sex chromosomes is regulated has not been defined. To resolve this, we profiled active chromatin features in the testes from wildtype and meiotic arrest mutants and integrate this with single-cell gene expression data from the Fly Cell Atlas. These data assign the timing of promoter activation for genes with germline-enriched expression throughout spermatogenesis, and general alterations of promoter regulation in germline cells. By profiling both active RNA polymerase II and histone modifications in isolated spermatocytes, we detail widespread patterns associated with regulation of the sex chromosomes. Our results demonstrate that the X chromosome is not enriched for silencing histone modifications, implying that sex chromosome inactivation does not occur in the Drosophila male germline. Instead, a lack of dosage compensation in spermatocytes accounts for the reduced expression from this chromosome. Finally, profiling uncovers dramatic ubiquitinylation of histone H2A and lysine-16 acetylation of histone H4 across the Y chromosome in spermatocytes that may contribute to the activation of this heterochromatic chromosome.