Universal transcriptomic signature of age reveals temporal scaling of Caenorhabditis elegans aging trajectories

To detect overlap or convergence among the diverse genetic pathways that can extend lifespan, we collected a dataset of 60 C.elegans age-dependent transcriptomes by RNA-seq technique for worm strains with vastly different lifespans. We selected four exceptionally long-lived mutants and three examples of the most successful life-extending RNAi treatments (which increased mean lifespan by 35% rather than 120% as reported). We used the dataset augmented with publicly available gene expression datasets to produce a transcriptomic signature of biological age. We introduced a transcriptomic measure of biological age and observed that its dependence on chronological age is modulated by a single parameter, the rate of aging. We hypothesized that the scaling revealed in the gene expression kinetics underlies the recently observed scaling of the survival curves in C.elegans, and the stochasticity in gene expressions leads to deceleration of mortality with age, reaching a plateau at advanced ages. Using experimental survival data, we confirm that the plateau mortality agrees closely with the estimate of Gompertz exponent at the cross-over age near the mean lifespan. The genes associated with aging in our data are enriched with the targets of transcription factors such as DAF-16, ELT-2, ELT-6, NHR-10, ZTF-9, NHR-86, and miRNAs including miR-57, -59, and -244, which is in agreement with previous studies. Overall, our meta-analysis results are consistent with a concept of aging based on critical dynamics of molecular level variables (e.g., gene expression), and support our view of aging as arising from dynamic instability of a single (critical) mode.


I. INTRODUCTION
To date, the strongest lifespan extension effect is achieved in C.elegans and corresponds to almost 10-fold increase of lifespan by the mg44 nonsense mutation of age-1 gene [1,2]. That extreme hyperlongevity, however, requires two generations of homozygosity for a mutation, and thus total pre-embryonal genetic disruption. In human subjects any reasonable therapy against aging would instead be applied at a considerably later stage (in adulthood, ideally at advanced ages). Sadly, the best "therapeutic" models yield significant, but considerably smaller reported increase of lifespan (up to +150% by let-363 RNAi) by interventions at embryonic and postembryonic stages [3]. Pharmacological inhibition of the gene homologues yielded even smaller effects on lifespan in flies, worms, yeasts [4][5][6] and mice [7][8][9]. It seems challenging that a single nonsense mutation can dramatically extend the lifespan of the animals, whereas an RNAi of the same gene is not able to produce the same level of effect, especially if administered later in life. Should we think that such a mutation changes the molecular machinery of the whole organism during development in such a dramatic way that the course of aging of the super-long-living strains is qualitatively different both in terms of rates and form, and hence cannot be reproduced therapeutically? Or is the Gene Regulatory Network (GRN) sufficiently robust, so that the pace of aging can be reduced without qualitative alterations of the relevant molecular mechanisms?
In order to address these questions, we prepared a dataset of C.elegans age-dependent transcriptomes by RNA-seq technique for worm strains with vastly different lifespans. We started by studying isogenic strains carrying the most long-lived C.elegans mutations: daf-2 (e1370 ) [strain SR806], age-1 (mg44 ) [SR808, at the first and second generations of homozygosity], and the longest-lived daf-2; daf-12 double mutant [strain DR1694] [1,2]. The average lifespans in the series range from twofold to tenfold longer than that of the wild type. To model the effects of therapeutic interventions in adult animals, we have chosen five RNAi treatments with the largest reported increase of lifespan (daf-4 [10], che-3 [11], cyc-1 [12,13], cco-1 [14], eat-4 [15]). We confirmed longevities of worm strains subjected to these five RNAi interventions. The lifespan modification effects in our hands, however, proved to be generally smaller than originally reported: the maximum increase of average lifespan in the series turned out to be +35%, instead of +120%. We further selected the three RNAi treatments representing, apparently, the most diverse modes of action and targeting che-3, daf-4 and cyc-1 for a mRNA collection follow-up. The overall variation in lifespan among the strains in both series together ranges from 17 up to 160 days. For each of the mutants or intervention examples we characterized the transcriptomes at several ages by RNA-seq technique in order to cover the whole lifespan range, 60 transcriptomes in total (9 different biological conditions). The experimental data were augmented by a meta-analysis of publicly available gene expression measurements in C.elegans (roughly 4000 samples).
We analyzed the gene expression signatures in the experiments within the context of the "aging at criticality" framework [16], borrowed from dynamical systems and bifurcation theory [17]. To summarize briefly, we expect that GRNs of most species operate near an orderdisorder bifurcation point [18] and are intrinsically unstable [19]. The order parameters associated with the instability undergo a stochastic evolution in an effective potential determined by the underlying genetic interactions and hence can be used to produce a molecular (transcriptomic in our case) signature of age. The curvature of the potential imposes a natural time scale, which turns out to be the Mortality Rate Doubling Time (MRDT), the most important characteristic of the aging trajectory. The pattern of molecular changes associated with age over time was robust across diverse biological conditions tested in our experiments. We also showed that life-extending mutations and knockdown therapies extend lifespan by stabilizing the regulatory network and thus reducing the rate of aging. The universality of the aging trajectories is therefore the molecular basis for the universal character of the survival curves in previous experiments [20]. We supported these theoretical propositions by confirming the predicted quantitative relation between the MRDT and the mortality plateau at late ages. Therefore the aging-at-criticality model provides the best mechanistic explanation of all observations hence is expected to be useful for understanding of aging dynamics and in future anti-aging target identification efforts in C.elegans and other species.

A. Selection of long-lived strains and life-extending interventions
Several mutations leading to exceptional longevity of C.elegans have been identified [21][22][23][24] and studied extensively for their remarkable elevations of both lifespan and stress-resistance traits [1,2] (see also Table Ia). To model the effects of anti-aging interventions we have tested five RNAi treatments targeting genes with the largest reported life-span extending effects (see Figure 1a and Table Ib for the summary of our confirmatory survival experiment and the lifespans reported in the original pub-lications). The strongest effect we obtained by RNAi corresponds to an increase in lifespan by +35%. We chose the three genes producing significant results in our survival experiment and representing apparently diverse modes of action, che-3, daf-4 and cyc-1, for the mRNA collection follow-up (see Table Ib, highlighted in bold).
B. Age-dependent transcriptomes of long-lived C.elegans strains To understand the nature of gene expression changes, associated with aging and lifespan, we produced an agedependent RNA-seq experimental dataset, consisting of the four mutant strains from Table Ia, three examples of life-extending RNAi from Table Ib and two control survivals representing C.elegans wild-type (Bristol-N2, strain DRM). We utilized a standard pipeline for RNAseq data processing (see Materials and Methods V D for details), yielding a dataset of 60 transcriptomes in total.
We perform the Principal Components Analysis (PCA) of all the acquired data and present its results for the mutants (Fig. 1b) and the RNAi treated strains (Fig. 1c) separately. Even though the number of samples in the data is dramatically smaller than the number of identified genes, the data covariance matrix is profoundly singular (the first four components explain nearly 75% of the overall variance in the data). The first principal component is universally associated with aging in all our experiments. The correlation with age rescaled by the average lifespan has Pearson correlation coefficient r = 0.75 (p < 10 10 ). This indicates that the dynamics of gene expression in relation to aging occurs along the single direction in the multidimensional space of all possible gene transcripts levels, coinciding with the first principal component loading vector. We also note that the variance of gene expression along PC1 remained quite consistent across all experiments, for samples ranging from youngest to oldest of each cohort. Because the first principle component comprised transcripts of the precisely the same genes in genetic mutants with large effects on lifespan (b), as in RNAi interventions with much more modest effects (c), these results strongly support their convergence on a common proximal pathway of life extension.
C. Meta-analysis of aging and signature of biological age in C.elegans transcriptomes We first asked whether the aging direction associated with the first principal component is indeed a property of aging dynamics, rather than an artifact due to inconsistency of PCA for relatively small sample size. In theory, this could be addressed by a meta-analysis of age-dependent transcriptomes from unrelated experiments using the publicly deposited gene-expression data. Typically, every single experiment ends up with a very   RNAi treatments correspond to different targeted genes/pathways affecting longevity (see Table Ib for the lifespans). Principal components analysis (PCA) of the experimental RNA-seq transcriptomic datasets for (b) four long-lived mutants and C.elegans wild-type (Bristol-N2, strain DRM), and (c) three life-prolonging RNAi-treated groups and C.elegans FV controls (fed bacteria harboring empty feeding-vector plasmid). The marker type denotes strains and RNAi groups (see Tables Ia and Ib for lifespans). The markers in (b) and (c) are colored according to age rescaled to lifespan (see the colorbar).
large number of measured transcripts' levels but very few samples. Analysis of such small datasets is very challenging [25,26] because it is prone to over-fitting even if all of the samples are collected in the same laboratory under the same conditions. A direct comparison of gene expression data obtained in different laboratories is notoriously hard due to the nature of experimental procedures, leading to batch effects requiring, but not necessarily cured by, extensive normalization [27][28][29].
On the horns of this dilemma, we theorized that a sufficiently large collection of "shallow" datasets (up to a dozen samples in a single experiment) would share the essential biology and include effects of experimental and batch differences as uncorrelated noise. To see if the "depth" can exceed the noise, we decided to create a comprehensive transcriptomic dataset for C.elegans by combining all the publicly-available gene expression experi-ments into a single database. This "MetaWorm" dataset contains, in toto, more than 400 different transcriptomic experiments with N = 3724 nematode samples characterized by the expression of G = 4861 genes most commonly expressed/detected in the samples. The composition of the dataset and the details of the normalization procedures are described in Section V G.
The gene expression variance in the combined dataset is dominated by batch effects and hence we do not expect PCA to reveal aging in association with the first principal component. Given the theoretical expectations, however, we can hope to identify the direction of aging by searching for genes associated with age. The MetaWorm dataset is sufficiently large to generate the cross-validation ensemble of single-gene association tests using random resampling. We further reduced the number of candidate genes by thresholding the transcripts by the frequency Geneset FDR DAF-16 targets [30] 4.6e-7 CEC-3 targets (spr-5 mutant, generation 10) [31] 2.0e-4 small RNAs decreased by starvation at P0 [32] 1.3e-3 PQM-1 L3 targets [33] 1.4e-3 Rb/E2F pathway (DPL-1,EFL-1,LIN-35), intestine [34] 2.7e-3 PMK-1 targets down in Day 15 vs. Day 6 [35] 3.8e-3 age-regulated ELT-2 targets [36] 3.8e-3 up by CSR-1 and w/out CSR-1-bound 22G RNA [37] 7.1e-3 proteins interacting with CEY-1 [38] 7.1e-3 MUT-14, SMUT-1 (DEAD box) [39] 1.6e-2 Hox gene targets (LIN-39, MAB-5, EGL-5) [40] 3.3e-2 Table II: Transcription factors or RNA-binding proteins for which targets are significantly enriched in the identified direction of aging (327 genes). Over-representation analysis was performed using WormExp database [41]. of significant associations in the resampling; we estimate that the chosen cross-validation threshold corresponds to p < 10 6 uncorrected, or p < 0.05 after Bonferroni correction. Finally, the set of genes associated with aging in our MetaWorm study consists of 327 genes: 260 upand 67 down-regulated with age, respectively. The aging genes are very few in number, whereas approximately 4000 out of 4861 genes were never observed in the top-200 gene-list after resampling, and thus do not seem to be significantly associated with age in the majority of experiments (see Electronic Supplementary Materials for the full list of genes implicated in aging in our calculations).

D. Biological interpretation of the direction of aging
We note that the genes associated with age may not include regulators of aging. The dataset we created by RNAseq, comprising 60 samples (section II A), is not large enough to produce a reliable set of relevant targets by a direct reconstruction of the gene-expression interaction subspace responsible for the progression of aging. The MetaWorm dataset is considerably larger, but is apparently noise-dominated and is not suitable, at least in our hands, for the identification of proximal mediators of aging. The solution of this problem would require a different methodology and we are leaving such an investigation for future research. In the meantime, we checked whether the 327 genes implicated by aging dynamics are enriched with the targets of natural gene expression regulators previously implicated in aging or longevity, such as Transcription Factors (TF) and miRNAs.

Targets of transcription factors or RNA binding proteins
We utilized WormExp [41], a curated database of gene sets built from published high-throughput expression studies in C.elegans, to predict the proteins whose targets are enriched within our "aging genes" set. The significant results of over-representation analysis against the known targets of transcription factors or RNA-binding proteins are presented in Table II. We observed considerable enrichment for the targets of DAF-16, a well-documented longevity-promoting transcription factor [42], that is also involved in antibacterial defense and is in part regulated by ELT-2 [43]. ELT-2 regulates downstream intestinal components of the DAF-2/DAF-16 pathway and controls p38-dependent gene induction [44], and is considered to be a key C.elegans transcription factor for intestinal development and function. Recently, normal aging in C.elegans was reported to be modulated by ELT-2, whose over-expression extends median lifespan by 15-25% [36]. The targets of ELT-2 were also significantly enriched.
We found enrichment for genes encoding protein targets of PMK-1 (p38 mitogen-activated protein kinase, a regulator of the innate immune response), reported to contribute to longevity and immunosenescence [35,45], and PQM-1, a stress-response transcription factor that interacts antagonistically with DAF-16, recently characterized as an important lifespan regulator [46].
We then turned to another listing of transcription factors and their target genes -a high-quality network obtained in [47]. We found that targets of the following transcription factors were the most over-represented among the genes comprising the identified direction of aging: NHR-1, ELT-6, NHR-10, MAB-3, ZTF-9, NHR-86. Three of these proteins are nuclear hormone receptors (NHRs) involved in lipid storage and/or catabolism, which may modulate aging via impaired autophagy and lysosomal lipolysis [48]. ELT-6 downregulates the expression of ELT-3 which in turn affects expression of multiple aging-related genes; ELT-6 knockdown increases C.elegans lifespan [49]. Finally, as annotated in Worm-Base [50], ZTF-9 is a zinc finger transcription factor regulated by several of the most potent longevity-limiting proteins including AGE-1, ISP-1, DAF-12 and DAF-2.

miRNAs
To determine if there are any miRNAs with overrepresented targets in the direction of aging, we performed enrichment analysis. Two sources of data describing miRNA targets were used -TargetScan [51] and MirTarBase [52]. The former resource contains results of in silico prediction of targets of miRNAs, while the latter database contains experimentally validated interactions. The five top miRNAs enriched from each data source are presented in Table III Table III: The miRNAs whose targets are significantly enriched in the identified direction of aging. Significant findings (Bonferroni-corrected p < 0.05) are highlighted in bold. Data on miRNA targets were obtained from TargetScan [51] and MirTarBase [52].
ing, we found evidence that several of the other miRNAs are indeed involved in the regulation of aging. The miRNA let-7 is differentially expressed in several aging studies and is associated with DNA damage checkpoint genes and mitochondrial respiration genes [53]. A known target of let-7 is the nuclear hormone receptor daf-12, which serves as a key regulator of germline elimination-induced longevity [54]. In a recent paper, it was reported that miR-60 is expressed almost exclusively in the intestine and directly modulates genes functioning in intestinal endocytosis. Furthermore, C.elegans animals lacking miR-60 have a dramatically extended lifespan under chronic oxidative stress condition [55]. Despite 3-fold upregulation of miR-253 during aging, no significant change in lifespan was observed for miR-253 mutants [56], suggesting a redundant, parallel route to modulation of the same targets. Other miRNAs implicated here were not previously reported to affect (or be affected by) aging.

E. Universality of aging trajectory
Thus far we have focused on the list of genes associated with aging in C.elegans. Another useful tool is the transcriptomic signature of biological age (cf. DNA methylation age, which is a weighted sum of DNA methylation features, trained to approximate chronological age in humans [57][58][59] and mice [60,61]). Since the progression of aging is inherently restricted to essentially a one-dimension manifold, the corresponding projecting operator is not unique. Any number of other convenient signatures can be obtained by applying additional requirements, such as for example, sparsity. We performed a cross-validated lasso regression of gene expressions in the MetaWorm dataset on age re-scaled by the average lifespan, in this 327-gene subset comprising the direction of aging (see Electronic Supplementary Materials for the summary of the respective regression weights). The pipeline is a variant of the distance correlation sure independent screening (DC-SIS), which is a fairly intuitive proposition given the extremely high dimensionality of the data [62].
We next applied the established regression model, the signature of biological age, to our own RNA-seq transcriptomic data. We observe that the biological age model yields a universal dependence of the transcriptomic bioage on the dimensionless measure of chronological age, rescaled to the average lifespan of a strain, for the mutants (Fig. 2a) and the RNAi treated strains (Fig. 2b). The universal scaling dependence appears to hold even for exceptionally long-lived strains, and altogether comprised animals with average adult lifespans ranging from 17 days for wild-type N2-DRM control worms to 160 days for age-1(mg44) mutants, and the trajectories appear identical for very long-lived mutants and wild-type worms with lifespans extended moderately by RNAi exposure. It thus appears that there is a single parameter describing the rate of aging.
We note that the direction of aging, as a list of transcripts physically associated with age, is a well defined object with a clear biological meaning. It is the right eigenvector of the gene-gene interaction matrix, characterizing the deterministic part of gene expression kinetics (see Section V E). On the other hand, the biological age signature is not unique and has no direct physical or biological meaning other than being a convenient tool for experimental data analysis.

F. Universal scaling in mortality law, mortality deceleration
Considering mortality as a consequence of molecular changes, the universality of the transcriptome aging trajectory ( Fig. 2a and Fig. 2b) should also manifest itself as a universality in survival curves. This is exactly what has been observed [20], when worms of different strains were subjected to different stresses and thus had their life-spans generally reduced. The probability of remaining alive to a certain age differed among the experimental conditions. At the same time, the shape of the survival curve was the same after a proper time-rescaling transformation.
The aging-at-criticality hypothesis has another testable prediction for mortality. Early in life, up to approximately the average lifespan, mortality can be approximated by the Gompertz equation where M 0 is the initial mortality rate (IMR) and ↵ is the Gompertz exponent, related to the MRDT. We predicted, however, that the exponential rise in mortality rates would cease at late ages, approaching a plateau at the mortality level determined as the value of the Gompertz exponent [16]. We obtained the mortality records from [20] and computed the Gompertz exponent and the mortality plateau estimates for all the reported conditions (Fig. 2c, see V F for the details of the calculations). We observe a fairly tight correlation between the parameters across a range of almost two orders of magnitude.

III. DISCUSSION
The results of the present work support the concept of aging being an organism-level manifestation of critical dynamics of a single unstable mode of the underlying regulatory network [16]. A key observation, relevant to both theoretical and practical applications, is that the aging direction is extremely robust and is shared by the long-lived mutants characterized by as much as 10-fold increase in lifespan, even relative to long-lived controls. We inferred a transcriptomic signature of biological age for C.elegans and investigated aging trajectories, the dependence of the biological age on chronological age, to distinguish the progression of aging in one strain from another. Early in life, up to approximately average lifespan, the age-dependent rise in transcriptomic indices of biological age is a universal function of a dimensionless age variable, obtained by rescaling the chronological age to the strain life expectancy. This is easy to understand if the influence of stochastic forces is small in young animals and therefore the progression is almost deterministic with the same time scale defining the shape of the gene expression variations and the value of the average lifespan. In fact, the temporal scale is defined by the underlying gene regulatory network stiffness and hence can be modified by external conditions or therapies on the level of specific gene-gene interactions in a very well-controlled manner.
According to the theory, later in life the dynamics of aging are defined by the interplay of increasingly important non-linearities of the gradually disintegrating gene regulatory network, as it is perturbed by stochastic forces. First, large deviations from the youthful state are incompatible with survival and hence the biological age should also plateau at roughly the average lifespan of each cohort, across a 10-fold range of longevities ( Fig. 2a  and Fig. 2b). Second, the stochasticity of aging dynamics leads to mortality deceleration such that mortality reaches a plateau at advanced ages, at the level of the Gompertz exponent obtained from fitting the mortality to the Gompertz equation close to the mean lifespan. Using the best available experimental data from [20], experimental mortality in C.elegans decelerates and reaches a plateau at late ages near the level of the Gompertz exponent, related to the mortality rate doubling time in the vicinity of the average lifespan, thus fully confirming the theoretical prediction. This phenomenon evidently underlies the plateau in mortality rates observed previously in very large populations of medflies and fruitflies [63,64], which we have here extended to relatively small populations of C.elegans (Fig. 2c). The results match expectations and, together with the scaling universality of the aging trajectory ( Fig. 2a and Fig. 2b), are the key results of our study and serve as a good confirmation of the theoretical model.
We note that, in analyzing transcriptomic data, one confronts a ubiquitous dilemma. One option would be to work with a few samples and thousands of genes assessed in a single laboratory, which is in many respects an ill-defined mathematical problem and prone to overfitting. The other possibility is to combine samples from many different laboratories (thus tackling the mathematical problem), although this produces side effects that are notoriously difficult to address, due to the nature of experimental procedures, i.e. batch effects requiring exten-sive normalization. In all likelihood, these difficulties lead to additional problems such as lack of reproducibility and inconsistency of results obtained in different laboratories. In the present study, we assumed that batch effects in gene expression are uncorrelated from one laboratory to another, and thus can be considered as noise. The critical dynamics of gene expressions associated with aging are responsible for the most biologically relevant variation in any set of transcriptomes of the same species. Therefore, we expected that the aging signature has the best signal-to-noise ratio and therefore a simple "piling" of the gene expression datasets from multiple unrelated experiments into a single MetaWorm dataset, after thorough normalization, yields valuable and transferable insights into aging dynamics.
Using the combined MetaWorm dataset, we have identified a list of genes associated with aging. The genes are relatively few in number (less than 7% of the available transcripts) and are enriched with targets of transcription factors known for promoting longevity via disparate interacting pathways, namely DAF-16, ELT-2, ELT-6, NHR-10, ZTF-9, and NHR-86. They are also enriched for miRNAs including miR-57, -59, -244 and -256 (see Results and Appendix for the complete results of over-representation analysis). We theorize that targeting other regulators in our list, such as miR-59 and miR-256, in addition to the known regulators, may further influence aging in C.elegans. In the near future, we plan to test our predictions experimentally by inhibition of identified transcription factors and/or miRNAs associated with normal progression along the direction of aging.
We expect that the "aging at criticality" framework should be applicable beyond C.elegans studies, extending to many other species including humans. Recently, we observed evidence of critical dynamics in association with aging in a subset of the large-scale human 2003-2006 National Health and Nutrition Examination Survey (NHANES) dataset representing physical activity metrics [65]. Although aging dynamics in humans appears on a relatively slow (sub-exponential) time course, agerelated changes are responsible for most of the variance in the signal, and progress along the aging trajectory explains most of the variation in mortality. In the same work, we demonstrated that the collective mode variable associated with age has a simple meaning of biological age and is shown to drift (increase) with age and perform a random walk leading to irreversibly increasing heterogeneity in the population.
The evidence for stochasticity in aging dynamics at late ages helps answer a critical question of when an antiaging treatment should be applied to obtain the largest possible effect. We speculate that at pre-embryonal and embryonal stages in the simplest animals, or early in life in humans, the growth of an organism is to a large degree determined by a development program. At more advanced ages, the stochasticity of the GRN kinetics takes its toll over time and leads to increasing pheno-typic heterogeneity at every level. This may help to explain why anti-aging interventions at the early stages have a broader and more generic effect on aging in diverse species. In contrast, interventions applied at late ages have the potential to treat specific conditions of an individual trajectory of aging progression with all stochastically accumulated errors -a possibility that could not have been foreseen or implemented without the identification of a universal direction of aging, i.e. a specific set of genes departing from youthful expression levels with the advance of age.
We hypothesize that the universality of C.elegans aging trajectories, a hallmark of this study, might be extrapolated to other species. For humans this universality would mean that prospective anti-aging treatments would eventually lead to a dramatic increase of lifespan with almost no detrimental effect on well-being.
The criticality hypothesis provides a fruitful insight into the evolutionary process of adjusting GRN properties to meet temporal requirements on the level of individual organism development, and of ever-present pressure of natural selection on species, see, e.g. [66]. On the one hand, it is favorable for an individual GRN to be as stress-resistant and long-lived as possible, but, on the other hand, the long-term effects of such a network's rigidity would be disastrous for the whole species. Evolution would fine-tune and optimize the parameters of the GRN and finally stabilize the GRN by driving the network toward the critical point. We believe the concept of aging as criticality in a dynamic system will provide a systematic approach to map longevity-controlling mechanisms, and eventually bring about powerful lifespanmodifying interventions.

IV. ACKNOWLEDGEMENTS
We thank Prof. Fontana and Nicholas Stroustrup for providing raw survival data and experimental lifetables from [20], G. Getmantsev and O. Burmistrova for proofreading of the work, I. Molodtsov and V. Kogan for invaluable help and advice on data analysis. The work was funded by Gero LLC, NIH Grant P01 AG012411-17A1 (WST Griffin, PI), and VA Merit Review Grant I01 BX001655 (RJ Shmookler Reis, PI).

B. Survivals
Lifespan assays were conducted at 20 C, as described previously [1]. Briefly, synchronous cultures were initiated by lysis of gravid hermaphrodites in alkaline hypochlorite. Worms were selected at the L4 larval stage, placed 50 worms per plate, and transferred at 1-to 2-day intervals onto fresh plates during days 1-7, and at 2-to 3-day intervals thereafter. A worm was scored as dead if it failed to move, either spontaneously or in response to a mechanical stimulus; lost worms were excluded (censored) from the survival analysis.

C. RNA isolation
Synchronized strains of C.elegans were grown on 100mm NGM plates, as above, and harvested for RNA extraction at the ages indicated. Worms were washed off plates and rinsed twice in survival buffer; after 30 min at 20 C (to allow digestion of enteral bacteria), they were flash frozen and stored at 80 C. Frozen worms were ground in a dry-ice-cooled mortar and pestle, and total RNA was extracted using RNeasy RNA extraction kits (Qiagen), followed by RNA purification for construction of transcript libraries using TruSeq RNA kits (Illumina, v.2). Sequences are generated as PE100 multiplexes, 100-bp paired-end reads from an Illumina HiSeq2500 or NextSeq instrument, producing 40 50 ⇥ 10 6 reads per sample. Paired samples are analyzed with DESeq2 (v1.4.5), and combined sequences are mapped to the C.elegans genome using TopHat 2.0.6.

D. Experimental RNA-seq dataset
RNA-seq reads were mapped to the human genome (WBcel235, Ensembl annotation) using TopHat 2.1.1 (with --b2-very-sensitive and --GTF options) [68] and gene-level read counts were obtained using the htseqcount software [69]. Low-expressed genes with at least one zero read count per sample were removed from subsequent analysis. Raw read counts were normalized using upper quartile method [70] and converted to RPKM values using edgeR library [71].

E. GRNs at criticality summary
We focus on transcriptomic data and describe timeevolution of gene expression by a matrix x n i , where indices n = 1 . . . N and i = 1 . . . G enumerate samples and gene transcripts, respectively, G is the total number of genes and N is the total number of samples. The measurements are taken at successive instances of time/age, t n . Following [16] we describe the dynamics of a transcriptome by a differential equation, representing the gene-expression kinetics: where all important properties of an organism are encapsulated into the function f (x i , t). The equations are still complex and for small deviations from the initial state we use a linearized version,ẋ i = f 0 + P i K ij x j . Here K ij = df i /dx j is the matrix, representing the simplest form of interaction in the model.
The matrix consists of kinetic rates corresponding to time scale of elementary biological processes, which are normally very high compared to characteristic rates associated with aging. The lifelong dynamics associated with aging is only possible whenever the interaction matrix K ij possesses a very small eigenvalue, ↵. Under the circumstances, the dynamics of the biological variables is said to be critical and can be accurately described with the help of the following one-factor model wherex i is the system state, corresponding to the animals early in life, b i is the right eigenvector of the matrix K ij corresponding to the singular eigenvalue ↵, and z n is the order parameter, associated with the GRN instability, a stochastic variable, associated with aging, and following the Langevin equation: z = ↵z + ⌘, where the stochastic variable ⌘ represents the effects of external and endogenous stress factors. In [16] we solved the equation and demonstrated, that ↵ is nothing else, but the Gompertz exponent, characterizing exponential increase of chances to die as a function of chronological age. The parameter defines a natural time scale, associated with the GRN instability, and depends on the balance between the damage accumulation and repair systems efficacy [19].
For the lack of a better word, the biological age z n , in contrast to the chronological age t n , is a universal parameter describing the progress of aging. Mutations or anti-aging therapies modify interactions in the GRN. Formally, this corresponds to alterations in the matrix elements K ij and change in its most relevant eigenvalue, ↵. The mortality in the model can be related to the effects of non-linearities: at sufficiently large values of z the higher order corrections can no longer be neglected and, starting at about average lifespan, the dynamics of the order parametre accelerates and the GRN disintegrates. Since normally the average lifespan exceeds the MRDT, the effects of the non-linearities and stochastic forces can be negletcted, at least early in life. Therefore, we expect that aging trajectories corresponding to different lifespans are self-similar and different by a single time scale factor ↵.
Practically, it would be interesting to derive an estimate of the biological age from the gene expression data using a linear model Since, however, the aging dynamics in biological signals is restricted to a one-dimensional manifold, the definition of biological age is not unique. This difficulty is not only a consequence of limited number of samples in a typical experiment. Any vector, orthogonal to b i can be added to the projector a i without changing the prediction results significantly, especially if the experimental noise (such as batch effects) is large. The best possible candidate for a i is the left eigenvector of the matrix K ij corresponding to the eigenvalue ↵.

F. Mortality analysis
Since there exists and we have here defined a universal molecular signature of bioage, the universality at the molecular level should be apparent as a universal behavior of mortality. The discrepancy between the mortality behavior predicted by the Gompertz equation and experimental mortality for late ages in C. elegans is sufficiently large and thus can be used to test the aging theory predictions quantitatively with high-quality mortality data. Mortality data of appropriate quality were recently published in [20], where a temporal scaling law of aging in C.elegans was observed, similar to that inferred for D.melanogaster [72]. This scaling law states that under the influence of some intervention, survival curves are stretched along the age axis by a dimensionless factor.
To extract the Gompertz exponent ↵ from the mortality data [20], we used the corresponding survival curves and fitted them to the prediction of the Gompertz equation. The procedure is only sensitive to the behavior of the survival curves in the neighborhood of the average lifespan. This is fortunate, since a GRN's stiffness parameter ↵ coincides with the Gompertz exponent in this interval only. The value of the plateau mortality M (t t ) was then calculated from the tail of the cumulative hazard m(t), estimated from the raw mortality data by the well-defined Nelson-Aalen routine [73] with the help of the Lifelines package [74].
Since the mortality rate reaches a plateau at late ages, the behavior of the cumulative hazard for these ages is linear and the value of M (t t ) can be extracted by linear regression of the cumulative hazard on age. We calculated the cumulative hazard m(t) from the experimental data and as the prediction of the Gompertz equation and compared them in Figure 3, where the disagreement between the two is substantial and significant, both qualitatively (exponential growth for the Gompertz equation and linear growth for the plateau mortality) and quantitatively (the cumulative hazard for the Gompertz equation is several orders of magnitude larger for late ages).

G. MetaWorm transcriptome
The MetaWorm dataset was compiled from all publicly-available information on expression profiles C.elegans from GEO database [75] and annotated with the corresponding worms' ages. For individual genes represented by multiple probesets, the probeset with the largest signal was used. Gene expressions in all datasets were normalized using the YuGene [27] algorithm. The final MetaWorm dataset represents a 3724 ⇥ 4861 matrix (samples-x-genes) and incorporates more than 400 transcriptomic experiments (see Electronic Supplementary Materials).

H. Enrichment analysis
Functional annotation (GO and pathway) analysis of the gene components comprising the direction of aging b i against targets of transcription factors was performed with WormExp resource [41], a curated collection of genesets derived from analysis of C.elegans expression datasets. We also utilized transcription factor regulation data obtained in [47]. Over-representation analysis against targets of miRNAs was performed using the cluster-Profiler package [76]. The list of predicted regulatory targets of C.elegans miRNAs was downloaded from TargetScanWorm [51] and MirTarBase [52] databases.