Status-dependent aging rates in long-1 lived , social mole-rats are shaped by HPA 2 stress axis 3

20

glucocorticoids -such as cortisol -should be lower than in non-breeders ( Fig. S6; see also our discussion 240 on ACTH-R below). 241 Finally, several pathways flagged by the weighted cross-tissue analysis seem to be derivatives of the 242 above-mentioned differentially expressed pathways instead of representing altered functions on their 243 own. For example, Huntington's (hsa05016), Parkinson's (hsa05012), and Alzheimer's (hsa05010) 244 diseases could, in principle, be interpreted as highly relevant for aging and lifespan. A closer inspection 245 of these pathways reveals, however, that the genes of the mitochondrial respiratory chain -which is the 246 core of the oxidative phosphorylation pathway -are in all three cases the main contributors to the 247 respective differential expression signals ( Fig. S7-S9, Data S1). Similarly, we see in the case of fat 248 digestion (hsa04975 Fat digestion and absorption) that two of the three largest contributors to the 249 differential expression signal of that pathway -ABCG8 and SCARB1 -are directly involved in the 250 transport of cholesterol [44,45]. Therefore, it seems likely that this signal is an expression of the altered 251 steroid hormone biosynthesis rather than indicating altered fat digestion. 252 Interestingly, three genes, which we had already mentioned as potential regulators during the pathway 253 analysis, also appeared among the ten most clearly altered genes on the weighted cross-tissue level: 254 GH1, IGF1, and CYP11A1 (Table 1). The top two among these ten are sulfotransferase family 2A member 255 1 (SULT2A1) and melanocortin 2 receptor (MC2R). SULT2A1 is the main catalyzer of the sulfonation of 256 the steroid hormone dehydroepiandrosterone (DHEA) to its non-active form DHEA-S [46]. DHEA has 257 repeatedly been proposed to be an "anti-aging hormone" because its levels are negatively associated 258 with chronological aging [47][48][49]. We found that SULT2A1 is strongly down-regulated in breeders ', liver 259 which is also the main location of its enzymatic action. The second candidate, MC2R, encodes the 260 adrenocorticotropin (ACTH)-receptor, which is the main inducer of glucocorticoid synthesis and a crucial 261 component of the hypothalamic-pituitary-adrenal (HPA) axis [50]. In humans and many other 262 mammals, prolonged glucocorticoid excess leads to Cushing's syndrome. Affected individuals exhibit 263 muscle weakness, immune suppression, impairment of the GH/IGF1 axis, higher risk of diabetes, 264 cardiovascular disease (hypertension), osteoporosis, decreased fertility, depression, and weight gain [51, 265 52]. The large overlap of these symptoms with those of aging could explain to some extent that 266 Cushing's syndrome patients exhibit considerably higher mortality rates [53]. We hypothesized that the 267 increased expression of the ACTH receptor in Fukomys non-breeders can cause similar expression 268 patterns and consequences (Fig. 5). Note -tissues are listed if they contribute at least 10% to cross-tissue expression 272 We tested this hypothesis (Fig. 5) by checking five of its key predictions. Altered MC2R expression (Fig. 273 6A) coincides with higher cortisol levels in hair samples from non-breeding F. mechowii than in those 274 from breeders of the same species (Begall et al., in prep.). Furthermore, glucocorticoids such as cortisol 275 exert their effect by binding to the glucocorticoid receptor that, in turn, acts as a transcription factor for 276 many genes [54]. We tested whether the expression of targets of the glucocorticoid-receptor (NR3C1) 277 was significantly altered throughout our data using two gene lists [55]: about three hundred direct 278 target genes of the receptor that were identified by chromatinimmunoprecipitation (i), and about 1300 279 genes that were found to be differentially expressed depending on the presence or absence of 280 exogenous glucocorticoid (ii). Both gene lists were found to be significantly affected by differential 281 expression at the weighted cross-tissue level (I, p=0.001; ii, p<10 -9 ) as well as in 5 (i) and 8 (ii) single 282 tissues (Table S7, S8). In line with our hypothesis, we observed that the weight gain in non-breeders 283 was, on average, twice as strong compared to the weight gain in breeders during the experiment 284 (p=7.49*10 -3 ,type II ANOVA, Fig. 6B). In addition, we found a subtle but significant influence of 285 reproductive status on the density of the vertebrae': the vertebrae of breeders were slightly denser than 286 those of age-matched non-breeders (p=0.03 for vertebra T12 only, and p=0.01 across all examined 287 vertebrae L1, L2, and T12; ANOVA, Fig. 6C, Table S11). 288

289
The vast lifespan differences between Fukomys breeders and non-breeders are, according to our RNA-290 seq data, associated with only subtle global pattern shifts in transcript levels. Concerning the tested 291 explanatory variables (Fukomys species, sex, breeding status), we found by far the highest number of 292 DEGs at the level of the species comparison. Although the number of DEGs between the sexes was 293 comparably high in almost each examined tissue, only very few DEGs were found comparing breeders 294 and non-breeders. One exception is the ovary, whose high number of DEGs corresponds well with the 295 disparity in reproductive activity. Other exceptions are liver, spleen, and, especially, the tissues of the 296 endocrine system (adrenal gland, pituitary gland, thyroid), in which the number of DEGs between 297 breeders and non-breeders ranged from more than 100 to more than 2,500. 298 Changes in the gene expression of the endocrine system are well known to play an important role both 299 in sexual maturation an in aging and the development of aging-associated diseases, e.g., diabetes and 300 cardiovascular diseases [56][57][58]. This finding fits well with the observation of substantial changes in the 301 endocrine system after the transition from non-breeders to breeders in the related naked mole-rat -302 one of the key results of a recent study in that species [7]. The dominance of differential expression in 303 endocrine tissue is also plausible insofar as these tissues exert a strong control function for other tissues 304 via hormone release. 305 Steroid hormone biosynthesis exhibits a bipartite pattern in breeders, with up-regulated sex steroid 306 genes and a simultaneous down-regulation of corticosteroid synthesis genes. The former could be 307 expected as a consequence of sexual activity in breeders. For aging, it is, however, interesting that 308 SULT2A1, a gene that codes for the specialized sulfotransferase converting the sex steroid DHEA to its 309 non-active form DHEA-S was found to be heavily down-regulated in breeders. DHEA is the most 310 abundant steroid hormone; it serves as a precursor for sex steroid biosynthesis but also has various 311 metabolic functions on its own [59,60]. DHEA levels decrease continuously during the human aging 312 process to an extent that favors it as an aging biomarker [61,62]. As a result, an aggressive advertising 313 of DHEA as "anti-aging hormone" in the form of dietary supplements could be observed in recent years 314 [47][48][49]60]. Despite conflicting experimental data from various animal studies and clinical trials, a 315 positive effect of DHEA on human health is frequently considered to be likely in the literature. Large-316 scale and long-term studies, however, are still pending [47,62,63]. Interestingly, DHEA is not only 317 described as an aging marker but also as a marker for clinically relevant glucocorticoid excess [64]. 318 Regarding corticosteroid synthesis, we hypothesize that the regulation of adrenal gland MC2R, coding 319 for the ACTH receptor, is likely to cause a considerable proportion of the overall observed expression 320 patterns and of the lifespan extension in breeders (Fig. 5). As a critical component of the HPA axis, ACTH 321 is a stress hormone that is produced by the pituitary gland and transported by the the blood to the 322 adrenal cortex, where it binds to the ACTH receptor [65]. Subsequently, the ACTH receptor induces the 323 synthesis of glucocorticoids [50], e.g., cortisol, which in turn cause immunosuppressive and various 324 metabolic effects throughout the organism [66]. Although in many mammals (such as humans, dogs, 325 and guinea pigs) glucocorticoid excess disorder, Cushing's syndrome, is caused by overproduction of the 326 hormone ACTH. Our results for Fukomys suggest, that the increased levels of the ACTH receptor may 327 lead to symptoms and expression patterns that could resemble some of molecular and phenotypic aspects of this pathological condition (Fig. 5). Our hypothesis is supported by a number of confirmed 329 downstream effects: the target genes of the glucocorticoid receptor are highly significantly affected by 330 differential gene expression; non-breeder symptoms were similar to those of humans with abnormally 331 elevated glucocorticoid levels: weight gain, decreases in bone density, and impairment of the GH/IGF1 332 axis (Fig. 5, Fig.6). Interestingly, blocking of the ACTH receptor has been suggested as a treatment for 333 human Cushing's syndrome [67]. 334 We looked for evidence of a regulation upstream of the ACTH receptor. Given the known positive 335 feedback loop between the ACTH receptor and its own ligand, decreased ACTH synthesis in breeders 336 would have been an obvious explanation [68]. We found that the expression of the ACTH polypeptide 337 precursor gene (POMC) is unchanged. However, also other post-transcriptional or post-translational 338 mechanisms such as cleaving may influence the ACTH levels in breeders and non-breeders. Of those 339 genes known to be involved in the regulation of MC2R [69, 70], we found that only PRKAR1B was 340 differentially regulated in the adrenal gland. However, this gene codes for only one of 7 subunits of the 341 involved protein kinase A. Alternative explanations could be that epigenetic modifications, other still 342 unidentified regulators of the transcript level, or both, are responsible for the differential expression. 343 In the circulatory system, represented by blood and spleen, down-regulation of coagulation factors was 344 observed in breeders. Because coagulation factors are known to be up-regulated during aging in humans 345 mice, rats, and even fish, this can be interpreted as a sign of a more juvenile breeders' transcriptome 346 [71,72]. Furthermore, the activity of coagulation factors is associated with a higher risk of coronary 347 heart disease [73]. However, we found no obvious histopathological lesions in the hearts or other 348 organs (spleen, kidney, liver, lung) of non-breeders.. Therefore, if down-regulation of coagulation 349 attentuates the aging process in Fukomys, it seems to exert its effect only, if at all, at a latter age not 350 investigated here. 351 Two defense mechanisms were also found to be up-regulated in breeders: xenobiotic metabolism and 352 TNF-α-signaling. Increased TNF-α-signaling often leads to the induction of apoptosis [74]. In line with 353 this, we found that apoptosis and P53-signaling were also up-regulated in breeders. Apoptosis is 354 considered an important anti-cancer mechanism [75,76]. We hypothesize this to compensate 355 cancerogenic effects of the anabolic alterations described above (especially the up-regulation of the 356 GH/IGF1 axis). In line with this, our more than 30 years' breeding history with several Fukomys species in 357 Germany and the Czech Republic suggests that breeders are as "cancer-proof" as non-breeders despite 358 their much longer lifespan (own unpublished data and R. Šumbera, personal communication). 359 Many anabolic pathways are up-regulated in breeders across tissues: protein biosynthesis, myogenesis, 360 and the GH/IGF1 axis. In line with this finding and with the fact that protein synthesis consumes 30 to 361 40% of a cell's ATP budget [77], we observed increased expression of mitochondrial respiratory chain 362 components, implying an increase in the capacity for cellular ATP production. On the other hand, 363 protein degradation and clearance by the proteasome are also up-regulated in breeders. The fact that 364 the expression of proteasomal genes is significantly linked to the genes of ribosome biogenesis and 365 oxidative phosphorylation indicates that those processes influence, or even trigger, each other and 366 hence are regulated in a coordinated manner in Fukomys mole-rats. 367 The results of differential expression of anabolic components such as the GH/IGF1 axis are surprising. 368 They fall within a debate in aging research that has been highly controversial over time: based on the 369 well-known fact that the expression and secretion of GH and IGF1 decline with age in humans and other 370 mammals [78], Rudman et al. administered synthetic GH to elderly subjects in 1990, thereby reversing a 371 number of aging-associated effects such as expansion of adipose mass [79]. This led to GH being 372 celebrated as an anti-aging drug [36], including dubious commercial offers. Today's aging research, on 373 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 [80][81][82]. In addition to several studies of synthetic 375 GH in humans yielding less convincing results than those of Rudman et al., the main reasons for this turn 376 are the results of studies on short-lived model organisms. From worms to mice, the impairment of the 377 GH/IGF1 axis by genetic intervention consistently led to longer lifespans (Table 2) Table 2). 383 Therefore, it is astonishing that massive up-regulation of these anabolic key components accompanies a 384 lifespan extension of approximately 100 % in long-lived mammals and potentially even contributes to it. 385 Several points could help to resolve this apparent contradiction: First, the up-regulation of anabolic 386 pathways and key genes is at least partially accompanied by the regulation of other mechanisms that 387 could plausibly compensate for deleterious effects. For example, it is, widely assumed that the negative 388 impact of enhanced protein synthesis on lifespan is to a large extent caused by the accumulation of 389 damaged or misfolded proteins, that is also known to contribute to aging-associated neurodegenerative 390 diseases [15,85,86]. Up-regulation of the proteasome, as we observed in breeders in a weighted cross-391 tissue approach and especially in endocrine tissues, is known to counteract these effects by clearing 392 damaged proteins leading to lifespan extension in worm and fly [86]. Enhanced proteasome activity has 393 also been linked with higher longevity of the naked mole-rat compared to the laboratory mouse [87, 88] 394 and in comparative approaches across several mammalian lineages [89]. We hypothesize that the 395 simultaneously high anabolic synthesis and catabolic degradation of proteins will lead to a higher 396 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 398 chain (oxidative phosphorylation) in breeders is compensated for by simultaneous up-regulation of the 399 reactive oxygen hallmark: the mitochondrial respiratory chain is the main source of cellular ROS which 400 can damage DNA, proteins and other cellular components [90,91]; the reactive oxygen hallmark 401 consists by definition of genes that are known to be up-regulated in response to ROS treatments. 402 Unsurprisingly, at least 25 % of these genes code for typical antioxidant enzymes such as thioredoxin, 403 superoxide dismutase, peroxiredoxin, or catalase that can detoxify ROS. Furthermore, the known cancer 404 promoting effects of an enhanced GH/IGF1 axis [36] could, to some extent, be compensated for by up-405 regulation of apoptosis and p53 signaling, because these are major mechanisms of cancer suppression 406 [92]. More generally, potential lifespan-extending effects of moderate up-regulation of both the 407 GH/IGF1 axis and ROS production can also be viewed in the light of the hormesis hypothesis [93], which 408 postulates that mild stressors can induce overall beneficial adaptive stress responses. In line with these 409 arguments, we found higher resting metabolic rates in breeders compared to non-breeders in Fukomys 410 anselli [94], a species closely related to F. micklemi. 411 A second point that could help to resolve this apparent contradiction concerns the time of intervention. 420 The transition from non-breeder to breeder takes place in adulthood, when by far the largest portion of 421 the growth process has already been completed. In contrast, genetic interventions aimed at prolonging 422 the lifespan by inhibiting the GH/IGFH1 axis (Table 2)  Most insights into current aging research originate from very short-lived model organisms (Table 2). It is 429 clear, on the other hand, that the observed effects of lifespan-prolonging interventions listed in Table 2  430 are by far the smallest in the model organisms with the relatively longest lifespans: mice and rats. 431 Compared to most other mammals, however, even mice and rats are short-lived. Given their body 432 weight and an often-used correlation between body weight and lifespan across mammals, their 433 observed maximum lifespan of about four years corresponds to only 51 % (for mice) and 32 % (for rats) 434 of the expected maximum lifespan. In contrast, humans can live as much as 463% as long as they would 435 be expected based on their body weight. This disparity makes them, given this particular model, the 436 most extreme known non-flying mammal species regarding maximum longevity residual [107]. 437 According to current data, Fukomys mole-rats reach values of ca. 200 % and thus can be regarded to be 438 closer to humans than to mice or rats in this respect. It is currently unclear to what degree aging mechanisms and lifespan-affecting interventions that were discovered in short-lived model species 440 apply to organisms that are far more long-lived, as in our experiment (e.g. [1,2]). It seems reasonable to 441 hypothesize that specific interventions that prolong the lifespan of long-lived organisms may have no 442 major effect in short-lived species, and vice versa. For example, it could be that changes suitable to 443 prolong the life of species with a low cancer-risk (e.g., African mole-rats, blind mole rats) [26, 108, 109] 444 would have no or only a marginal effect in laboratory mice, whose primary cause of natural death is 445 cancer [110]. Therefore, it may also be possible that those gene expression patterns caused by our 446 lifespan-extending intervention in Fukomys mole-rats but contradicting the current knowledge obtained 447 from short-lived organisms may highlight differences in aging mechanisms between short-lived and 448 long-lived species. 449 As a fourth perspective, one could interpret the down-regulation of the GH/IGF1 axis in non-breeders as 450 a byproduct of the apparent up-regulation of the HPA axis in non-breeders, that may well be adaptive in 451 itself. In the wild, non-breeding Fukomys mole-rats can maximize their inclusive fitness by either 452 supporting their kin in their natal family or by founding a new family elsewhere. It has therefore been 453 suggested that the shorter lives of non-breeders could be adaptive on the ultimate level if longevity 454 were traded against some other fitness traits, such as competitiveness, to defend colonies against 455 intruders or to enhance their chances for successful dispersal [13]. A constitutively more activated HPA 456 stress axis is expected to offer advantages for both family defense and dispersal but it carries the risk of 457 status-specific aging symptoms, such as muscle weakness, lower GH/IGF1 axis activity, lower bone 458 density etc. in the long run [52]. This effect may become even more pronounced under laboratory 459 conditions where grown non-breeders cannot decide to disperse even if they would like to. 460 However, even the down-regulation of the GH/IGF-axis may be adaptive in itself for non-breeders if it 461 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 463 decreasing activity is well documented to correlate with chronological age in a wide range of mammals, 464 including mice, rats, dogs, and humans [78]. Also, decreasing activity correlates, under pathologic 465 conditions such as Cushing's syndrome, with many symptoms akin to aging (Fig. 5). It has been 466 suggested that one solution to this apparent contradiction may be that that the GH/IGF1 axis is 467 adaptively down-regulated in aging organisms as a reaction to already accumulated aging-related 468 symptoms so that additional damage can be avoided [111,112]. 469 Finally, recent findings of positively selected genes in African mole-rats (family Bathyergidae, containing 470 also Fukomys and Heterocephalus) could partly explain some of the surprising results. It is striking that 471 translation, and oxidative phosphorylation were among the strongest differentially expressed molecular 472 processes concering the breeding status. Earlier, these processes were also reported to be the most 473 affected by positive selection in the phylogeny of African mole-rats [113]. Furthermore, IGF1 was one of 474 thirteen genes that were found to be under positive selection in the last common ancestor of the mole-475 rats. This may indicate that the corresponding mechanisms were evolutionarily adapted to be less 476 detrimental and make their up-regulation more compatible with a long lifespan. Since the mere fact of 477 positive selection does not permit to draw conclusions about the direction of the mechanistic effect, this 478 hypothesis, however, needs to remain speculative. 479

480
We performed a comprehensive transcriptome analysis that, for the first time within mammalian 481 species, compared naturally occurring cohorts of species with massively diverging aging rates. The 482 comparison of faster-aging Fukomys non-breeders with similar animals that were experimentally 483 elevated to the slower-aging breeder status revealed by far the most robust transcriptome differences 484 within endocrine tissue: adrenal gland, ovary, thyroid, and pituitary gland. Genes and pathways involved 485 in anabolism, such as GH, IGF1, translation and oxidative phosphorylation, were differentially expressed. 486 Their inhibition is among the best-documented life-prolonging interventions in a wide range of short-487 lived model organisms (Table 2). Surprisingly, however, we found that the expression of these 488 mechanisms was consistently higher in slower-aging breeders than in faster-aging non-breeders. This 489 indicates that even basic molecular mechanisms of the aging process known from short-lived species 490 cannot easily be transferred to long-lived species. In particular, this applies to the role of the GH/IGF1 491 axis, which has in recent years been unilaterally described as harmful [81,82,106]. In addition, special 492 features of the mole-rats could also contribute to the explanation of the unexpected result that genes 493 and processes differentially expressed between reproductive statuses were also strongly altered during 494 the evolution of the mole-rats [113]. Another intriguing possibility is that, in line with the hormesis 495 hypothesis [93], moderate harmful effects of anabolic processes can be hyper-compensated for by up-496 regulation of pathways such as proteasomes, P53-signaling and antioxidant defense against ROS that we 497 observed in slower-aging breeders as well. 498 Furthermore, our work provides evidence that the HPA stress axis is a key regulator for the observed 499 downstream effects, including the lifespan difference. The effects are likely to be triggered by 500 differential expression of the gene MC2R coding for the ACTH receptor resulting in an altered stress 501 response in breeders vs. non-breeders. This is supported by the fact that cortisol levels in the non-502 breeders are elevated. Furthermore, the set of direct and indirect target genes of the glucocorticoid 503 receptors is strongly affected by differential expression, and numerous known downstream effects of 504 glucocorticoid excess have been demonstrated for non-breeders, such as muscle weakness, weight gain, 505 and GH/IGF1 axis impairment. Overall, this evidence suggests that MC2R and other genes along the 506 described signalling pathway are promising targets for possible interventions in aging research. 507

508
Animal care and sampling 509 All animals were housed in glass terraria with dimensions adjusted to the size of the family (min. 40 cm × 510 60 cm) in the Department of General Zoology, Faculty of Biosciences, University of Duisburg-Essen. The 511 terraria are filled with a 5 cm layer of horticultural peat or sawdust. Tissue paper strips, tubes, and solid 512 shelters were provided as bedding/nesting materials and environmental enrichment. Potatoes and 513 carrots are supplied ad libitum as stable food, supplemented with apples, lettuce, and cereals. Fukomys 514 mole-rats do not drink free water. Temperature was kept fairly constant at 26± 1 °C, and humidity, at 515 approximately 40 %. The daily rhythm was set to 12 hours darkness and 12 hours light. 516

New breeder pairs (new families) were established between March and May 2014. Each new family was 517
founded by two unfamiliar, randomly chosen adult non-breeders of similar age (min/max/mean: 518 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 519 from already existing separate colonies. These founder animals were moved to a new terrarium in 520 which they were permanently mated. In both species, more than 80 % of these new pairs reproduced 521 within the first 12 months (total number of offspring by the end of the year 2016, 82 F. micklemi and 81 522 F. mechowii). Only founders with offspring were subsequently assigned to the breeder cohort; founders 523 without offspring were excluded from the study. Non-breeders remained in their natal family together 524 with both parents and other siblings. 525 Before sampling, animals were anaesthetized with 6 mg/kg ketamine combined with 2.5 mg/kg xylazine 531 [114]. Once the animals lost their pedal withdrawal reflex, 1 to 2 ml of blood was collected by cardiac 532 puncture, and the animals were killed by surgical decapitation. Blood samples (100 µl) were collected in 533 RNAprotect Animal Blood reagent (Qiagen, Venlo, Netherlands). Tissue samples -hippocampus, 534 hypothalamus, pituitary gland, thyroid, heart, skeletal muscle (M. quadriceps femoris), lateral skin, small 535 intestine (ileum), upper colon, spleen, liver, kidney, adrenal gland, testis, and ovary -were transferred 536 to RNAlater (Qiagen, Venlo, Netherlands) immediately after dissection and, following incubation, were 537 stored at -80°C until analysis. Kidney and heart samples were treated with proteinase K before extraction, as recommended by the 545 manufacturer. Library preparation was performed using the TruSeq RNA v2 kit (Illumina, San Diego, USA) 546 which includes selection of poly-adenylated RNA molecules. RNA-seq was performed by single-end 547 sequencing with 51 cycles in high-output mode on a HiSeq 2500 sequencing system (Illumina) and with 548 at least 20 million reads per sample, as described in Table S9. Read data for F. mechowii and F. micklemi 549 were deposited as European Nucleotide Archive study with the ID PRJEB29798 (Table S9). 550 It was ensured for all samples that the results of the respective sequencing passed "per base" and "per 552 sequence" quality checks of FASTQC [115]. The reads were then mapped against previously published 553 and with human gene symbols annotated F. mechowii and F. micklemi transcriptome data [24,113] quantification, we further analyzed only those reference transcripts whose gene symbols were present 558 in the transcript catalogs of both species -this was the case for 15,199 transcripts (the size of the union 559 was 17,065). As mapping algorithm "bwa aln" of the Burrows-Wheeler Aligner (BWA) [117] was used, 560 allowing no gaps and a maximum of two mismatches in the alignment. Only those reads that could be 561 uniquely mapped to the respective gene were used for quantification. Read counts per gene and sample 562 can be found in Data S4. As another check, we ensured that all samples exhibited a Pearson correlation 563 coefficient of at least 90% in a pairwise comparison based on log 2 -transformed read counts against all 564 other samples of the same experimental group as defined by samples that were equal in the tissue as 565 well as the species, sex, and reproductive status of the source animal. 566 Differentially expressed genes analysis 567 P-values for differential gene expression and fold-changes were determined with DESeq2 [118] and a 568 multi-factorial design. The DESeq2-algorithm also includes strict filtering based on a normalized mean 569 gene count that makes further pre-filtering unnecessary [118]. Therefore, those genes whose read count 570 was zero for all examined samples were removed before further analysis, thereby reducing the number 571 of analyzed genes to 15,181. The multi-factorial design means that, separately for each tissue, we input 572 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 574 variables by controlling for additional variance in the other two variables. This approach resulted in a 575 four-times higher sample size than with an approach that would have been based on comparisons of 576 two experimental groups, each of which would be equal in tissue, species, sex, and reproductive status. 577 It is known that the statistical power in RNA-seq experiments can increase considerably with sample size 578 [119]. P-values were corrected for multiple testing with the Benjamini-Hochberg correction [120] (false 579 discovery rate -FDR). 580 The results of the DEG-analysis can be found in Data S5-S7. 581

582
Let ( 1 , … ) represent the p-values obtained from differential gene expression analysis in the tissue 583 corresponding with index and the indices 1, … , corresponding to the examined genes. Forthermore, 584 let ( 1 , … . , ), with 1 ≤ ≤ and 1 ≤ ≤ , represent the p-values of genes with the indices 585 = ( 1 , … , ) belonging to a corresponding pathway that is tested for enrichment of differential 586 expression signals. To determine the enrichment p-values at the pathway-level, we calculated the test 587 statistic ℱ for the gene indices in tissue index according to Fisher's method for combining p-values 588 also known as Brown's method: 589 Because p-values at the gene level were, as is common in RNA-seq experiments [121], not equally 590 distributed, we empirically estimated the null distribution of the test statistic for each pathway instead 591 of using the χ 2 -distribution as frequently suggested in the literature [122][123][124]. This was done by 592 calculating ℱ for 1,000 random drawings, each without replacement, = ( 1 , … , ) , with 593 1 ≤ ≤ , 1 ≤ ≤ 1000. If the resulting p-value was zero (meaning ∀ 1 ≤ ≤ 1000: ℱ > ), the 594 procedure was repeated with 10,000 and 100,000 random drawings, respectively. In addition, the 595 indices were divided into and , depending on whether their fold-change was > 1 or < 1 in 596 breeder vs. non-breeder comparison, and ℱ and ℱ calculated. The ratio ℱ ℱ was used as an 597 indicator for functional up-or down-regulation of the corresponding pathway (Fig. 4, Fig. S3, S4). Using 598 this approach, enrichment p-values were estimated for all  and MSigDB-hallmarks 599 [31], as well as across all examined tissues (Data S8, S9). In addition, the procedure was applied to test 600 whether the known 300 direct and 1300 indirect glucocorticoid receptor target genes [55] were 601 enriched for status-dependent differential expression signals (Tables S7, S8). 602 Similarly, cross-tissue DEG-p-values were weighted with a modified test statistic (c.f. [125]). Given the 603 definitions from above, we calculated the weighted cross-tissue test statistic Ƒ for the gene as 604 follows: 605 where and ̅̅̅̅̅̅ are the logarithmized fold-change between reproductive states and 607 normalized mean expression (across sexes, species, and reproductive status) for the gene with index 608 and tissue with index -both calculated by DESeq2 [118] -, is the signum function, and is the 609 number of examined tissues. The reasoning for using this test statistic was as follows: For the weighted 610 cross-tissue analysis, we assumed that the gene serves the same function throughout the entire 611 organism. Therefore, the test statistic given above weights the p-values of the various tissues by the 612 respective expression levels in those tissues. This ensures that, weighting e. g., a ubiquitously expressed 613 genes such as TP53 is relatively equally across tissues, whereas for typical steroid hormone-biosynthesis 614 genes, such as CYP11A1, the endocrine tissue results determine almost exclusively the weighted cross-615 tissue p-value. Furthermore, the test statistic rewards, based on the mentioned assumption, consistency 616 in the direction of gene regulation throughout tissues. All calculated values for , ̅̅̅̅̅̅ , , as 617 well as the resulting Ƒ and p-values, can be found in Data S10. 618 Finally, weighted cross-tissue enrichment p-values at the pathway level were estimated by applying the 619 above-described method at the pathway level (based on test statistic ℱ) to the gene level weighted 620 cross-tissue p-values. 621

P-values were corrected for multiple testing with the Benjamini-Hochberg correction [120] (FDR). 622
Weighted gene co-expression network analysis 623 We used the WGCNA R package to perform weighted correlation network analysis [33] of all 636 624 samples at once. We followed the authors' usage recommendation by choosing a soft power threshold 625 based on scale-free topology and mean connectivity development (we chose power=26 with a soft R 2 of 626 0.92 and a mean connectivity of 38.6), using biweight midcorrelation, setting maxPOutliers to 0.1, and 627 using "signed" both as network and topological overlap matrix type. The maximum block size was 628 chosen such that the analysis was performed with a single block and the minimum module size was set 629 to 30. The analysis divided the genes into 26 modules, of which 5 were enriched for reproductive status 630 DEGs based on Fisher's exact test and an FDR threshold of 0.05. Those 5 modules were tested for 631 enrichment among  with the same test and significance threshold (Table S10). In 632 addition, module eigengenes were determined and clustered (Table S10). Then the topological overlap 633 matrix that resulted from the WGCNA analysis ( = [ , ], where the row indices 1 ≤ ≤ 634 | | correspond to genes and the column indices 1 ≤ ≤ | | 635 correspond to samples) was used to determine pairwise connectivity between all KEGG-pathways that 636 showed differential expression at the weighted cross-tissue level (Fig. S5). Based on the definition of 637 connectivity of genes in a WGCNA analysis [33], we defined the connectivity between two sets of indices 638 and each corresponding to genes as , = ∑ ∑ , ∈ Y\X ∈ X\Y . P-values for the connectivities 639 were determined against null distributions that were empirically estimated by determining for each pair 640 and the connectivities of 10,000 pairs of each | | and | | randomly drawn indices (without 641 replacement), respectively. Since "signed" was used as the network and topological overlap matrix type, 642 the tests were one-sided. 643 Other analysis steps 644 Hierarchical clustering (Fig. S2) was performed based on Pearson correlation coefficients of log 2 -645 transformed read counts between all sample pairs using the complete-linkage method [126]. The 646 principal variant component analysis ( Fig. 2A) was performed with the pvca package from Bioconductor 647 [127] and a minimum demanded percentile value of the amount of the variabilities, that the selected 648 principal components needs to explain, of 0.5. Enrichments of DEGs among genes enlisted in the Digital 649 Aging Atlas database ([29], Data S11) were determined with Fisher's exact test, the Benjamini-Hochberg 650 method (FDR) [120] for multiple test correction, and a significance threshold of 0.05. Pathway 651 visualization (Fig. S6-S9) was performed with Pathview [128]. For the direction analysis of Fukomys 652 reproductive status DEGs in previous experiments in naked-mole rats and guinea pigs, we examined 653 those 10 tissues that were examined in all species (Table S6, Data S12). Separately for each tissue and 654 combination of species -naked mole-rat or guinea pig -and sex, we determined how many Fukomys 655 reproductive status DEGs were up-regulated or down-regulated. We also performed two-sided binomial 656 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 Fukomys [12]. Therefore, non-breeders of opposite sexes were taken from different families -labeled as 738 "Family A" (blue) and "Family B" (red) -and permanently housed in a separate terrarium. The two 739 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-741 matched controls that remained in the faster-aging non-breeders of "Family A" and "Family B" were 742 included in our study -in most cases full siblings (ideally litter mates) of the respective new breeders. 743

Figure legends
After at least two years and two pregnancies in "Family C", breeders from "Family C" and their controls 744 from Colonies A and B were put to death, and tissues were sampled for later analysis, that included 745 identification of differentially expressed genes (DEGs). The shown experimental scheme was conducted 746 with 5 to 7 (median 6) specimens per cohort (defined by breeding status, sex, species) and 12 to 15 747 tissues (median 14) per specimen: in total, 46 animals and 636 samples. 748 guinea pig (GPs) [7]. For NMR there is also evidence that breeders have a (slightly) longer lifespan than 764 non-breeders, whereas for GP the opposite is assumed [7,27]. Across ten tissues that were examined in 765 both studies, the analysis determined whether status-dependent differentially expressed genes 766 identified in the current study were regulated in the same or opposite direction in NMR and GP (Table  767 S6). The listed p-values (two-sided binomial exact test; hypothesized probability, 0.5) describe how 768 extremely the ratio of genes expressed in the same and opposite directions deviates from a 50:50 ratio. 769 Green and orange indicate the majority and minority of genes within a comparison, respectively. Figure  770 created with BioRender.com. 771 within the cells give the FDR, i.e., the multiple testing corrected p-value, for the respective 776 pathway/hallmark and tissue. As indicated by the color key, red and green stand for up-or down-777 regulated in breeders, respectively. White indicates a pathway/hallmark that is significantly affected by 778 differential expression and whose signals for up-and down-regulation are approximately balanced. Dark 779 colors up to black mean that there is little or no evidence that the corresponding pathway/hallmark is 780 affected by differential gene expression. Figures S3 and S4 provide detailed overviews of all 781 pathways/hallmarks that are enriched in at least one tissue for status-dependent differential expression 782 signals. 783 Figure 5. Model of the stress axis as a key mechanism for status-dependent lifespan differences in 784 Fukomys. From a wide range of mammals, including humans [52], dogs [129], horses [130], cats [131], and guinea pigs [132], it is known that chronic glucocorticoid excess leads to a number of pathologic 786 symptoms that largely overlap with those of aging and result in considerably higher mortality rates for 787 affected individuals [53#808]. The most common cause of chronic glucocorticoid excess is excessive 788 secretion of the adrenocorticotropic hormone (ACTH) by the pituitary gland. ACTH is transported via the 789 blood to the adrenal cortex where it binds to the ACTH-receptor (ACTHR; encoded by the gene MC2R) 790 which induces the production of glucocorticoids, especially cortisol. Glucocorticoids are transported to 791 the various tissues, where they exert their effect by activating the glucocorticoid receptor (NR3C1) that 792 acts as a transcription factor and regulates hundreds of genes. The constant overuse of this 793 transcriptional pattern eventually leads to the listed symptoms. Our hypothesis is that the permanent, 794 high expression of the ACTH-receptor in Fukomys non-breeders causes effects similar to those known 795 from overproduction of the hormone. In line with this hypothesis, i) cortisol levels are increased in non-796 breeders and ii) target genes of the glucocorticoid-receptor are highly enriched for status-dependent 797 differential gene expression. Furthermore, the animals were examined for common symptoms of 798 chronic glucocorticoid excess: iii) non-breeders gained more weight during the experiment than 799 breeders, iv) exhibited lower bone density at the end of the experiment, and v) lower gene expression 800 in the GH-/IGF1 axis than breeders. 801   Figure S1. Overview of the tissues examined in this study with a schematic mole-rat representation. 811 hsa050160 Alzheimer's disease at the cross-tissue level. 828   Table S1. Overview of number of samples that were examined with regard to status, sex, species, and 829 tissue. 830 Table S2. Animal description. 831 Table S3. Overview of the tissues that were successfully sampled for each type of animal. 832  Table S7. Analysis of status-dependent differential gene expression enrichment on glucocorticoid 838 receptor target genes that were determined by Phuc Le at al. 2005 839 (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1186734/), using chromatin immunoprecipitation 840 (CHIP). 841 Table S8. Analysis of status-dependent differential gene expression enrichment on glucocorticoid 842 receptor target genes that were determined by Phuc Le at al. 2005 843 (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1186734/), using a differential gene expression 844 analysis. 845 Table S9. Sample description. 846 Table S10. WGCNA module clustering and functional enrichment analysis regarding these modules.

Supplement Data 849
Data S1. (zip) For each KEGG pathway and MSigDB hallmark that was detected to be significantly (FDR < 850 0.1) enriched for status-dependent differential expression at the weighted cross-tissue level, the data 851 set contains a *.tsv file with the genes that form the respective pathway/hallmark sorted by their 852 individual contribution to the enrichment. The files also provide an overview of the p-values and fold-853 changes of those genes in those tissues in which the genes are expressed most highly. The data set contains one *.tsv file per analyzed tissue, as well as an additional *.tsv file for the 869 weighted cross-tissue level results. 870 Data S10. (tsv.gz) Overview of the weighted cross-tissue differential gene expression analysis. Contains 871 all p-values, fold-changes, and mean expression values for all genes across all tissues as well as the 872 weighted, combined cross-tissue test statistics, p-values, and fold-changes for all genes. 873 Data S11. (tsv) Genes of the Digital Aging Atlas used in this study. 874 Data S12. (zip) Comparison of status-dependent differentially expressed genes with results of an earlier, 875 similar study involving naked mole-rats (NMRs) and guinea pigs (GPs). The data set contains one *.tsv 876 file for each of the ten tissues that were examined in both studies. Each *.tsv file lists the differentially 877 expressed genes for the respective tissue as identified in this study with the determined FDRs and fold-878 changes, as well as the fold-changes determined in the earlier study using naked mole-rats and guinea 879 pigs. 880 Data S13. (zip) Glucocorticoid receptor target genes that were determined by Phuc Le et al. and tested 881 in this study for enrichment of status-dependent differential gene expression. The data set contains one 882 *.tsv file each for target genes determined via chromatin immunoprecipitation (CHIP) and differential 883 gene expression analysis. Phuc Le et al. determined  1225 3305 1917 P-value 5.2*10 -58 1*10 -44 9*10 -25 4.7*10 -23