The HPA stress axis shapes aging rates in long-lived, social mole-rats

Sexual activity and/or reproduction doubles life expectancy in the long-lived rodent genus Fukomys. To investigate the molecular mechanisms underlying this phenomenon, we analyzed a total of 636 RNA-seq samples across 15 tissues. This analysis suggests that the differences in life expectancy between reproductive and non-reproductive mole-rats are mainly caused by critical changes in the regulation of the hypothalamic-pituitary-adrenal stress axis, which we further substantiate with a series of independent evidence. In accordance with previous studies, the up-regulation of the proteasome and several so-called “anti-aging molecules”, such as DHEA, is also linked with enhanced life expectancy. On the other hand, several our results oppose crucial findings in short-lived model organisms. For example, we found the up-regulation of the IGF1/GH 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.

. Therefore, non-breeders of opposite 111 sexes were taken from different families -labeled as "Family A" (blue) and "Family B" (red) -and permanently 112 housed in a separate terrarium. The two unrelated animals mated with each other, thus producing offspring and 113 becoming breeders of the new "Family C" (green). In addition to the animals that were promoted to be slower-114 aging breeders, age-matched controls that remained in the faster-aging non-breeders of "Family A" and "Family B" 115 were included in our study -in most cases full siblings (ideally litter mates) of the respective new breeders. After 116 at least two years and two pregnancies in "Family C", breeders from "Family C" and their controls from Colonies A 117 and B were put to death, and tissues were sampled for later analysis, that included identification of differentially 118 expressed genes (DEGs). The shown experimental scheme was conducted with 5 to 7 (median 6) specimens per 119 cohort (defined by breeding status, sex, species) and 12 to 15 tissues (median 14) per specimen: in total, 46 120 animals and 636 samples.

122
We measured gene expression differences between breeders and non-breeders in two African mole-rat 123 species, F. mechowii and F. micklemi. Altogether, we performed RNA-seq for 636 tissue samples 124 covering 16 tissue types from both species, sexes, and reproductive states (breeders and non-breeders). 125 Each of the four groups (male/female breeders/non-breeders) of each species consisted of 5 to 7 126 animals (see Tables S1 -S5 for sample sizes, animal data, and pairing schemes). For each tissue, we 127 conducted a multi-factorial analysis of differentially expressed genes (DEGs): the analysis was based on 128 the variables reproductive state, sex, and species. During this exercise, we focused on the differences 129 between slower-aging breeders and faster-aging non-breeders. This approach increases our statistical 130 power by giving us a four-fold increase of sample size in comparison to species-and sex-specific breeder 131 vs. non-breeder analyses. At the same time, we can additionally reduce the number of false-positive 132 DEGs by restricting the analysis to those breeding status-related genes that show the same direction in 133 both sexes and both species. We deliberately focused on those genes to concentrate our study on 134 universal mechanisms that hold for both sexes and species. 135 To globally quantify the transcriptomic differences between the reproductive states, we performed 136 three analyses: clustering of the samples based on pairwise correlation, principal variant component 137 analysis, and an overview of the number of DEGs between reproductive states in comparison to DEGs 138 between species and sex. Clustering of the samples based on pairwise correlations showed a full 139 separation of the two species at the highest cluster level (Fig. S1). Below that level, an almost complete 140 separation according to tissues was observed. Within the tissue clusters, the samples did not show a 141 clear-cut separation between sex or breeder/non-breeder status. Accordingly, a principal variance 142 component analysis showed that species, tissue, and the combination of both variables accounted for 143 98.4 % of the total variance in the data set; individual differences explained 1.4 % of the variance, and 144 only 0.004 % was explained by breeder/non-breeder status ( Fig. 2A). Regarding the numbers of DEGs, 145 we found -unsurprisingly considering the aforementioned facts -by far the highest number of DEGs in 146 the species comparison (Fig. 2B). Although in almost every examined tissue the numbers of detected 147 DEGs were also high between sexes, most tissues exhibited very few DEGs due to breeder/non-breeder 148 status. Exceptions were liver, spleen, ovary and, especially, tissues of the endocrine system (adrenal 149 gland, pituitary gland, thyroid), in which the number of DEGs between breeders and non-breeders 150 ranged from more than sixty to several thousand.

158
Next, we evaluated the relevance of reproductive status DEGs for aging and aging-related diseases. For 159 this analysis, we first determined overlaps by using the Digital Aging Atlas (DAA) -a database of genes 160 that show aging-related changes in humans (Craig et al. 2015). Across species and sexes, significant 161 overlaps (FDR < 0.05, Fisher's exact test) with the DAA were found in three tissues: adrenal gland, ovary 162 and pituitary gland (false discovery rate [FDR] = 0.005, each; Fig. 3A). Among these three endocrine 163 tissues, the DEGs of the ovaries overlapped significantly with those from adrenal (p=2.8*10 -27 ) and 164 pituitary glands (p=0.005), but there was no significant overlap between the two glands (Fig. 3A). Thus, 165 together, we found indications for aging-relevant expression changes after the transition from non-166 breeders to breeders in three tissues of the endocrine system, which presumably affect separate 167 aspects of aging in adrenal and pituitary glands. 168 Moreover, we compared the DEGs with respect to the reproductive status that we identified in Fukomys 169 with regard to their direction to transcript-level changes observed in similar experiments using naked 170 mole-rats and guinea pigs (Bens et al. 2018). The direction of the status-dependent DEGs regulation in 171 Fukomys, as found in this study, was significantly more often the same rather than opposite as in the 172 naked mole-rat (females, 60 %, p=5.2*10 -58 ; males, 62 %, p=10 -44 for females, Fig.3B, Table S6). In the 173 guinea pig, on the contrary, the Fukomys reproductive status DEGs were significantly more often 174 regulated in the opposite direction (females, 57 %, p=9*10 -25 ; males, 59 %, p=4.7*10 -23 , Fig. 3B, Table  175 S6). Thus, at the single-gene level, the expression changes linked to reproductive status may affect 176 lifespan differently in long-lived African mole-rats than in shorter-lived guinea pigs. 177

179
For each tissue, we separately tested whether the identified differentially expressed genes between status groups significantly 180 overlapped with the genes within the Digital Aging Atlas database (Fisher's exact test, FDR < 0.05). Significant overlaps were 181 found for three tissues: adrenal gland, ovary and pituitary gland. The Venn-diagram depicts the overlaps of these three tissues 182 with the Digital Aging Atlas and with each other (** : FDR < 0.01; *** : FDR < 0.001). B) A similar experiment comparing the 183 transcriptomes of breeders versus non-breeders was recently conducted in naked mole-rats (NMRs) and guinea pig (GPs) (Bens 184 et al. 2018). For NMR there is also evidence that breeders have a (slightly) longer lifespan than non-breeders, whereas for GP 185 the opposite is assumed (Bens et al. 2018;Ruby et al. 2018). Across ten tissues that were examined in both studies, the analysis 186 determined whether status-dependent differentially expressed genes identified in the current study were regulated in the 187 same or opposite direction in NMR and GP (Table S6). The listed p-values (two-sided binomial exact test; hypothesized 188 probability, 0.5) describe how extremely the ratio of genes expressed in the same and opposite directions deviates from a 50:50 189 ratio. Green and orange indicate the majority and minority of genes within a comparison, respectively. Figure Fig. S3 and S4). Because the individual interpretation of each of these pathways/hallmarks would 200 go beyond the scope of this study, we focus here on those 14 pathways and 13 hallmarks that were 201 significantly different between breeders and non-breeders (FDR < 0.1) in a global analysis across all 202 tissues (Fig. 4). Because many pathways are driven mainly by gene expression in subsets of tissues, we 203 weighted gene-wise the differential expression signals from the various tissues by the respective 204 expression levels in the tissues. For instance, the expression level of the growth hormone (GH) gene GH1 205 is known to be almost exclusively expressed in the pituitary glands. In our data set the GH1 level of the 206 pituitary gland accounted for 99.96 % of the total GH1 across all tissues. Accordingly, in pathways that 207 contain GH, our weighted cross-tissue differential expression signal for this gene is almost exclusively 208 determined by the pituitary gland. On the contrary, a differential expression signal of this gene in 209 another tissue with a very low fraction of the gene's total expression would have almost no impact on 210 the weighted cross-tissue level -even if that signal were very strong (see Methods section for details). 211 respectively. White indicates a pathway/hallmark that is significantly affected by differential expression and whose signals for 218 up-and down-regulation are approximately balanced. Dark colors up to black mean that there is little or no evidence that the 219 corresponding pathway/hallmark is affected by differential gene expression. Figures S3 and S4 provide detailed overviews of all 220 pathways/hallmarks that are enriched in at least one tissue for status-dependent differential expression signals.  Fig. 4A,B). We also found strong indications for increased 230 protein degradation (hsa03050 PROTEASOME, Fig. 4A). This weighted cross-tissue signal was, in contrast 231 to the situation regarding ribosomes, driven mainly by two tissue types: the adrenal gland and the 232 gonads. 233 To examine whether the simultaneous up-regulation of the ribosome, proteasome, and oxidative 234 phosphorylation is a coordinated regulation, we performed a weighted gene co-expression network 235 analysis (WGCNA) (Langfelder and Horvath 2008) from our gene count data and examined the 236 With xenobiotic metabolism and TNF-α-signaling, two defense hallmarks were also found to be up-265 regulated in breeders by the weighted cross-tissue analysis (HALLMARK_XENOBIOTIC_METABLISM and 266 HALLMARK_TNFA_SIGNALING_VIA_NFKB, Fig. 4B). The up-regulation of the reactive oxygen species 267 (ROS) hallmark comprising genes coding for proteins that detoxify ROS 268 (HALLMARK_REACTIVE_OXYGEN_SPECIES_PATHWAY, Fig. 4B) falls into a similar category. 269 Another interesting aspect that was found to be significantly altered in breeders is steroid hormone 270 biosynthesis (hsa00140 Steroid hormone biosynthesis, Fig. 4A). In this case, both up-and down-271 regulated genes were in the pathway, and their absolute fold-changes were roughly balanced. Steroid 272 hormones, on the one hand, comprise sex steroids -because these hormones are important players in 273 sexual reproduction, such differences should be expected given the experimental setup. On the other 274 hand, the class of steroid hormones -corticosteroids -has regulatory functions in metabolism, growth, 275 and the cardiovascular system, as well as in the calibration of the immune system and response to stress 276 (Liu et al. 2013). The most influential contributor to the differential pathway signal by far on the 277 weighted cross-tissue level was CYP11A1, which codes for the (single) enzyme that converts cholesterol 278 to pregnenolone. This is the first and rate-limiting step in steroid hormone synthesis (Miller and Auchus 279 2011). Because CYP11A1 was found to be up-regulated in breeders in its main places of synthesis -the 280 gonads and the adrenal gland (Table 1) -it can be assumed that the total output of steroid hormone 281 biosynthesis in breeders is increased. The pattern of up-and down-regulation on the KEGG pathway, 282 however, suggests that sex steroids especially are produced at a higher rate in breeders whereas 283 circulating levels of glucocorticoids -such as cortisol -should be lower than in non-breeders ( Fig. S6; 284 see also our discussion on ACTH-R below). 285 Finally, several pathways flagged by the weighted cross-tissue analysis seem to be derivatives of the 286 above-mentioned differentially expressed pathways instead of representing altered functions on their own. For example, Huntington's (hsa05016), Parkinson's (hsa05012), and Alzheimer's (hsa05010) 288 diseases could, in principle, be interpreted as highly relevant for aging and lifespan. A closer inspection 289 of these pathways reveals, however, that the genes of the mitochondrial respiratory chain -which is the 290 core of the oxidative phosphorylation pathway -are in all three cases the main contributors to the 291 respective differential expression signals ( Fig. S7-S9, Data S1). Similarly, we see in the case of fat 292 digestion (hsa04975 Fat digestion and absorption) that two of the three largest contributors to the 293 differential expression signal of that pathway -ABCG8 and SCARB1 -are directly involved in the 294 transport of cholesterol (Liu et al. 2008;Wang et al. 2015). Therefore, it seems likely that this signal is an 295 expression of the altered steroid hormone biosynthesis rather than indicating altered fat digestion. 296 Note -tissues are listed if they contribute at least 10% to cross-tissue expression 299 Interestingly, three genes, which we had already mentioned as potential regulators during the pathway 300 analysis, also appeared among the ten most clearly altered genes on the weighted cross-tissue level: GH1, IGF1, and CYP11A1 (Table 1). The top two among these ten are sulfotransferase family 2A member 302 1 (SULT2A1) and melanocortin 2 receptor (MC2R). SULT2A1 is the main catalyzer of the sulfonation of 303 the steroid hormone dehydroepiandrosterone (DHEA) to its non-active form DHEA-S (Hammer et al. 304 2005). DHEA has repeatedly been proposed to be an "anti-aging hormone" because its levels are 305 negatively associated with chronological aging (Baulieu 1996 transcriptional pattern eventually leads to the listed symptoms. Our hypothesis is that the permanent, high expression of the 328 ACTH-receptor in Fukomys non-breeders causes effects similar to those known from overproduction of the hormone. In line 329 with this hypothesis, i) cortisol levels are increased in non-breeders and ii) target genes of the glucocorticoid-receptor are 330 highly enriched for status-dependent differential gene expression. Furthermore, the animals were examined for common 331 symptoms of chronic glucocorticoid excess: iii) non-breeders gained more weight during the experiment than breeders, iv) 332 exhibited lower bone density at the end of the experiment, and v) lower gene expression in the GH-/IGF1 axis than breeders. 333 We tested this hypothesis ( about three hundred direct target genes of the receptor that were identified by 340 chromatinimmunoprecipitation (i), and about 1300 genes that were found to be differentially expressed 341 depending on the presence or absence of exogenous glucocorticoid (ii). Both gene lists were found to be 342 significantly affected by differential expression at the weighted cross-tissue level (I, p=0.001; ii, p<10 -9 ) 343 as well as in 5 (i) and 8 (ii) single tissues (Table S7, S8). In line with our hypothesis, we observed that the 344 weight gain in non-breeders was, on average, twice as strong compared to the weight gain in breeders 345 during the experiment (p=7.49*10 -3 ,type II ANOVA, Fig. 6B). In addition, we found a subtle but 346 significant influence of reproductive status on the density of the vertebrae': the vertebrae of breeders 347 were slightly denser than those of age-matched non-breeders (p=0.03 for vertebra T12 only, and p=0.01 348 across all examined vertebrae L1, L2, and T12; ANOVA, Fig. 6C, Table S11). 349

357
The vast lifespan differences between Fukomys breeders and non-breeders are, according to our RNA-358 seq data, associated with only subtle global pattern shifts in transcript levels. Concerning the tested 359 explanatory variables (Fukomys species, sex, breeding status), we found by far the highest number of 360 DEGs at the level of the species comparison. Although the number of DEGs between the sexes was 361 comparably high in almost each examined tissue, only very few DEGs were found comparing breeders 362 and non-breeders. One exception is the ovary, whose high number of DEGs corresponds well with the 363 disparity in reproductive activity. Other exceptions are liver, spleen, and, especially, the tissues of the 364 endocrine system (adrenal gland, pituitary gland, thyroid), in which the number of DEGs between 365 breeders and non-breeders ranged from more than 100 to more than 2,500. Regarding corticosteroid synthesis, we hypothesize that the regulation of adrenal gland MC2R, coding 389 for the ACTH receptor, is likely to cause a considerable proportion of the overall observed expression 390 patterns and of the lifespan extension in breeders (Fig. 5). As a critical component of the HPA axis, ACTH 391 is a stress hormone that is produced by the pituitary gland and transported by the the blood to the 392 adrenal cortex, where it binds to the ACTH receptor (Fridmanis et al. 2017). Subsequently, the ACTH 393 receptor induces the synthesis of glucocorticoids (Walker et al. 2015), e.g., cortisol, which in turn cause 394 immunosuppressive and various metabolic effects throughout the organism (Becker 2013). Although in 395 many mammals (such as humans, dogs, and guinea pigs) glucocorticoid excess disorder, Cushing's 396 syndrome, is caused by overproduction of the hormone ACTH. Our results for Fukomys suggest, that the 397 increased levels of the ACTH receptor may lead to symptoms and expression patterns that could 398 resemble some of molecular and phenotypic aspects of this pathological condition (Fig. 5). Our 399 hypothesis is supported by a number of confirmed downstream effects: the target genes of the 400 glucocorticoid receptor are highly significantly affected by differential gene expression; non-breeder 401 symptoms were similar to those of humans with abnormally elevated glucocorticoid levels: weight gain, 402 decreases in bone density, and impairment of the GH/IGF1 axis (Fig. 5, Fig.6). Interestingly, blocking of 403 the ACTH receptor has been suggested as a treatment for human Cushing's syndrome (Newfield 2010). 404 We looked for evidence of a regulation upstream of the ACTH receptor. Given the known positive 405 feedback loop between the ACTH receptor and its own ligand, decreased ACTH synthesis in breeders 406 would have been an obvious explanation (Imai et al. 2001). We found that the expression of the ACTH 407 polypeptide precursor gene (POMC) is unchanged. However, also other post-transcriptional or post-408 translational mechanisms such as cleaving may influence the ACTH levels in breeders and non-breeders. The fact that the expression of proteasomal genes is significantly linked to the genes of ribosome 437 biogenesis and oxidative phosphorylation indicates that those processes influence, or even trigger, each 438 other and hence are regulated in a coordinated manner in Fukomys mole-rats. 439 The results of differential expression of anabolic components such as the GH/IGF1 axis are surprising.  (Table 2). Therefore, it is astonishing that massive up-regulation of these anabolic key components accompanies a 457 lifespan extension of approximately 100 % in long-lived mammals and potentially even contributes to it. 458 Several points could help to resolve this apparent contradiction: First, the up-regulation of anabolic 459 pathways and key genes is at least partially accompanied by the regulation of other mechanisms that 460 could plausibly compensate for deleterious effects. For example, it is, widely assumed that the negative 461 impact of enhanced protein synthesis on lifespan is to a large extent caused by the accumulation of 462 damaged or misfolded proteins, that is also known to contribute to aging-associated neurodegenerative   The transition from non-breeder to breeder takes place in adulthood, when by far the largest portion of 499 the growth process has already been completed. In contrast, genetic interventions aimed at prolonging 500 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 507 clear, on the other hand, that the observed effects of lifespan-prolonging interventions listed in Table 2  508 are by far the smallest in the model organisms with the relatively longest lifespans: mice and rats. 509 Compared to most other mammals, however, even mice and rats are short-lived. Given their body 510 weight and an often-used correlation between body weight and lifespan across mammals, their 511 observed maximum lifespan of about four years corresponds to only 51 % (for mice) and 32 % (for rats) 512 of the expected maximum lifespan. In contrast, humans can live as much as 463% as long as they would 513 be expected based on their body weight. This disparity makes them, given this particular model, the 514 most extreme known non-flying mammal species regarding maximum longevity residual (Tacutu et al. 515 2013). According to current data, Fukomys mole-rats reach values of ca. 200 % and thus can be regarded 516 to be closer to humans than to mice or rats in this respect. It is currently unclear to what degree aging 517 mechanisms and lifespan-affecting interventions that were discovered in short-lived model species 518 apply to organisms that are far more long-lived, as in our experiment (e.g. only a marginal effect in laboratory mice, whose primary cause of natural death is cancer (Seluanov et 524 al. 2018). Therefore, it may also be possible that those gene expression patterns caused by our lifespan-525 extending intervention in Fukomys mole-rats but contradicting the current knowledge obtained from 526 short-lived organisms may highlight differences in aging mechanisms between short-lived and long-lived 527 species. 528 As a fourth perspective, one could interpret the down-regulation of the GH/IGF1 axis in non-breeders as 529 a byproduct of the apparent up-regulation of the HPA axis in non-breeders, that may well be adaptive in 530 itself. In the wild, non-breeding Fukomys mole-rats can maximize their inclusive fitness by either 531 supporting their kin in their natal family or by founding a new family elsewhere. It has therefore been 532 suggested that the shorter lives of non-breeders could be adaptive on the ultimate level if longevity 533 were traded against some other fitness traits, such as competitiveness, to defend colonies against 534 intruders or to enhance their chances for successful dispersal (Dammann and Burda 2006). A 535 constitutively more activated HPA stress axis is expected to offer advantages for both family defense 536 and dispersal but it carries the risk of status-specific aging symptoms, such as muscle weakness, lower 537 GH/IGF1 axis activity, lower bone density etc. in the long run (Ferrau and Korbonits 2015). This effect may become even more pronounced under laboratory conditions where grown non-breeders cannot 539 decide to disperse even if they would like to. 540 However, even the down-regulation of the GH/IGF-axis may be adaptive in itself for non-breeders if it 541 has the potential to protect them from further damage. Note that for today's conventional view that 542 stronger activation of the GH/IGF1 axis accelerates aging (Table 2), it is generally challenging that 543 decreasing activity is well documented to correlate with chronological age in a wide range of mammals, 544 including mice, rats, dogs, and humans (Bartke 2019). Also, decreasing activity correlates, under 545 pathologic conditions such as Cushing's syndrome, with many symptoms akin to aging (Fig. 5). It has 546 been suggested that one solution to this apparent contradiction may be that that the GH/IGF1 axis is 547 adaptively down-regulated in aging organisms as a reaction to already accumulated aging-related 548 symptoms so that additional damage can be avoided (Berryman et al. 2008;Milman et al. 2016). 549 Finally, recent findings of positively selected genes in African mole-rats (family Bathyergidae, containing 550 also Fukomys and Heterocephalus) could partly explain some of the surprising results. It is striking that 551 translation, and oxidative phosphorylation were among the strongest differentially expressed molecular 552 processes concering the breeding status. Earlier, these processes were also reported to be the most 553 affected by positive selection in the phylogeny of African mole-rats (Sahm et al. 2018b). Furthermore, 554 IGF1 was one of thirteen genes that were found to be under positive selection in the last common 555 ancestor of the mole-rats. This may indicate that the corresponding mechanisms were evolutionarily 556 adapted to be less detrimental and make their up-regulation more compatible with a long lifespan. Since 557 the mere fact of positive selection does not permit to draw conclusions about the direction of the 558 mechanistic effect, this hypothesis, however, needs to remain speculative.

560
We performed a comprehensive transcriptome analysis that, for the first time within mammalian 561 species, compared naturally occurring cohorts of species with massively diverging aging rates. The 562 comparison of faster-aging Fukomys non-breeders with similar animals that were experimentally 563 elevated to the slower-aging breeder status revealed by far the most robust transcriptome differences 564 within endocrine tissue: adrenal gland, ovary, thyroid, and pituitary gland. Genes and pathways involved 565 in anabolism, such as GH, IGF1, translation and oxidative phosphorylation, were differentially expressed. 566 Their inhibition is among the best-documented life-prolonging interventions in a wide range of short-567 lived model organisms (Table 2). Surprisingly, however, we found that the expression of these 568 mechanisms was consistently higher in slower-aging breeders than in faster-aging non-breeders. This 569 indicates that even basic molecular mechanisms of the aging process known from short-lived species 570 cannot easily be transferred to long-lived species. In particular, this applies to the role of the GH/IGF1 571 axis, which has in recent years been unilaterally described as harmful Furthermore, our work provides evidence that the HPA stress axis is a key regulator for the observed 580 downstream effects, including the lifespan difference. The effects are likely to be triggered by 581 differential expression of the gene MC2R coding for the ACTH receptor resulting in an altered stress 582 response in breeders vs. non-breeders. This is supported by the fact that cortisol levels in the non-583 breeders are elevated. Furthermore, the set of direct and indirect target genes of the glucocorticoid 584 receptors is strongly affected by differential expression, and numerous known downstream effects of 585 glucocorticoid excess have been demonstrated for non-breeders, such as muscle weakness, weight gain, 586 and GH/IGF1 axis impairment. Overall, this evidence suggests that MC2R and other genes along the 587 described signalling pathway are promising targets for possible interventions in aging research. 588 Kidney and heart samples were treated with proteinase K before extraction, as recommended by the 626 manufacturer. Library preparation was performed using the TruSeq RNA v2 kit (Illumina, San Diego, USA) 627 which includes selection of poly-adenylated RNA molecules. RNA-seq was performed by single-end 628 sequencing with 51 cycles in high-output mode on a HiSeq 2500 sequencing system (Illumina) and with 629 at least 20 million reads per sample, as described in Table S9. Read data for F. mechowii and F. micklemi 630 were deposited as European Nucleotide Archive study with the ID PRJEB29798 (Table S9). 631 was zero (meaning ∀ 1 ≤ ≤ 1000: ℱ > ), the procedure was repeated with 10,000 and 100,000 678 random drawings, respectively. In addition, the indices were divided into and , depending 679 on whether their fold-change was > 1 or < 1 in breeder vs. non-breeder comparison, and ℱ and 680 ℱ calculated. The ratio ℱ ℱ was used as an indicator for functional up-or down-regulation of 681 the corresponding pathway (Fig. 4, Fig. S3 hsa050160 Alzheimer's disease at the cross-tissue level. 835   Table S1. Overview of number of samples that were examined with regard to status, sex, species, and 836 tissue. 837 Table S2. Animal description. 838 Table S3. Overview of the tissues that were successfully sampled for each type of animal. 839  Table S10. WGCNA module clustering and functional enrichment analysis regarding these modules.

Supplement Data 856
Data S1. (zip) For each KEGG pathway and MSigDB hallmark that was detected to be significantly (FDR < 857 0.1) enriched for status-dependent differential expression at the weighted cross-tissue level, the data 858 set contains a *.tsv file with the genes that form the respective pathway/hallmark sorted by their 859 individual contribution to the enrichment. The files also provide an overview of the p-values and fold-860 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 876 weighted cross-tissue level results. 877 Data S10. (tsv.gz) Overview of the weighted cross-tissue differential gene expression analysis. Contains 878 all p-values, fold-changes, and mean expression values for all genes across all tissues as well as the 879 weighted, combined cross-tissue test statistics, p-values, and fold-changes for all genes. 880 Data S11. (tsv) Genes of the Digital Aging Atlas used in this study. 881 Data S12. (zip) Comparison of status-dependent differentially expressed genes with results of an earlier, 882 similar study involving naked mole-rats (NMRs) and guinea pigs (GPs). The data set contains one *.tsv 883 file for each of the ten tissues that were examined in both studies. Each *.tsv file lists the differentially 884 expressed genes for the respective tissue as identified in this study with the determined FDRs and fold-885 changes, as well as the fold-changes determined in the earlier study using naked mole-rats and guinea 886 pigs. 887 Data S13. (zip) Glucocorticoid receptor target genes that were determined by Phuc Le et al. and tested 888 in this study for enrichment of status-dependent differential gene expression. The data set contains one 889 *.tsv file each for target genes determined via chromatin immunoprecipitation (CHIP) and differential 890 gene expression analysis. Phuc Le et al. determined differentially expressed genes between mice that 891 were treated with exogenous glucocorticoids and untreated controls. The relevant table columns from 892 Phuc Le et al. were added by human gene symbol and Entrez IDs that were used for enrichment analysis. 893