A universal transcriptomic signature of age reveals the temporal scaling of Caenorhabditis elegans aging trajectories

We collected 60 age-dependent transcriptomes for C. elegans strains including four exceptionally long-lived mutants (mean adult lifespan extended 2.2- to 9.4-fold) and three examples of lifespan-increasing RNAi treatments. Principal Component Analysis (PCA) reveals aging as a transcriptomic drift along a single direction, consistent across the vastly diverse biological conditions and coinciding with the first principal component, a hallmark of the criticality of the underlying gene regulatory network. We therefore expected that the organism’s aging state could be characterized by a single number closely related to vitality deficit or biological age. The “aging trajectory”, i.e. the dependence of the biological age on chronological age, is then a universal stochastic function modulated by the network stiffness; a macroscopic parameter reflecting the network topology and associated with the rate of aging. To corroborate this view, we used publicly available datasets to define a transcriptomic biomarker of age and observed that the rescaling of age by lifespan simultaneously brings together aging trajectories of transcription and survival curves. In accordance with the theoretical prediction, the limiting mortality value at the plateau agrees closely with the mortality rate doubling exponent estimated at the cross-over age near the average lifespan. Finally, we used the transcriptomic signature of age to identify possible life-extending drug compounds and successfully tested a handful of the top-ranking molecules in C. elegans survival assays and achieved up to a +30% extension of mean lifespan.

such that the course of aging of the super-long-living strains is qualitatively different both regarding rates and form, and hence could not be easily reproduced therapeutically. Alternatively, perhaps the gene regulatory network is sufficiently robust that a therapy can reduce the pace of aging without qualitative alterations of the relevant molecular mechanisms.
To address these alternatives, we compiled an RNA-seq dataset of age-dependent transcriptomes produced from C. elegans isogenic strains and populations that have vastly different lifespans. Among them are three long-lived isogenic strains carrying mutations: daf-2(e1370), age- 1(mg44) [at the first and second generations of homozygosity], and daf-2(e1391); daf-12(m20) double mutant 1,11 ; three RNAi treatments (daf-4, che-3 and cyc-1); and two independent control runs represented by wild-type (Bristol-N2, DRM stock). The overall range of adult lifespans across all the experiments extends from 17 to 160 days. For each of the mutants or interventions, we measured gene-expression levels over time, across their lifespans, collecting 60 transcriptomes in total (9 different biological time-series, each in duplicate).
Principal Component Analysis (PCA) of gene expression reveals aging in all strains and treated groups as a transcriptomic drift in a single direction, consistent across the vastly diverse biological conditions and coinciding with the first principal component of the combined dataset, which is a hallmark of the criticality of the underlying regulatory network 12 . We therefore expected that the organism's physiological aging state can be characterized by a single stochastic variable having the meaning of biological age and coinciding approximately with the first principal component score. The aging trajectory, i.e., the dependence of the biological age variable on chronological age, is then universally determined by the underlying regulatory interactions and the experimental conditions through a single phenomenological property describing the effective regulatory-network stiffness. The quantity imposes a natural time scale proportional to the mortality rate doubling time, the fundamental dynamic characteristic of the aging process 12 .
To identify a set of genes universally associated with aging across many different biological conditions, the aging signature, and to evaluate the theoretical model, we performed a meta-analysis of publicly available gene-expression measurements in C. elegans (more than 4000 samples in total). The identified aging signature comprises a set of genes, many of which have no known role in the regulation of aging or longevity. We used the same data to introduce a robust transcriptomic biomarker of aging, as a read-out or predictor of "biological age", and demonstrated its utility across the datasets. The biological age dynamics in our experiments reveal a universal "aging trajectory": the rescaling of age by lifespan simultaneously brings together the time-dependent trajectories of the transcriptomic biomarker on age and the survival curves. Throughout the paper, "age" means the chronological adult age (post-L4/adult molt for C. elegans). Therefore, the universality of aging trajectories may provide a natural molecular basis for the scaling universality of survival curves recently observed 13 and independently confirmed in the survival data of all the strains and treatments in our experiments. We investigated the relationship between the stochastic evolution of the biological age variable and mortality using the survival data from an independent experiment. We also experimentally confirmed the model prediction of the equivalence between the mortality rate doubling exponent (inferred at the cross-over age, corresponding to the average lifespan) and the limiting mortality value (corresponding to the mortality plateau). Finally, we used the transcriptomic signature of age to identify possible life-extending drug compounds and successfully tested a handful of them in C. elegans survival assays.
We confirmed the longevity of worm strains subjected to these five RNAi interventions at 20 °C (see Table 1 and Fig. 1). The extensions of mean lifespan ranged from 5-16%, often well below those previously reported, which was likely due to the different mode of action (mutation or RNAi treatment), or variation in the study protocols, including RNAi-efficacy, exposure times and/or culture temperatures. Neuronal resistance to RNAi may also contribute, and could also impinge on non-neuronal tissues through an inter-tissue feedback loop [23][24][25] . Current and previous lifespan data are summarized in Table 2. The most substantial relative effect obtained by RNAi corresponds to an increase in mean lifespan from 20.1 days to 23.3 days (+16.2%) by RNAi of che-3.
www.nature.com/scientificreports www.nature.com/scientificreports/ We started by performing an exploratory analysis of all the gene-expression data with the help of principal component analysis (PCA). The first principal component (PC1), along which the variance of the data is maximal, is the only component significantly correlated with age (r = 0.75, p < 10 −10 , accounting for r 2 = 56% of total variance). PC1 simultaneously arranges the mutants, Fig. 2(a), and the RNAi treated strains, Fig. 2(b), according to their respective values of chronological age. The total amplitude of change from youngest to oldest ages is approximately the same for all 9 groups despite their wide range of longevities.
Our gene expression data suggest that the aging signature, i.e., the set of genes transcriptionally associated with age, is robust and remains consistent under a reasonably broad range of experimental conditions and genetic or epigenetic interventions. No other principal component score has a statistically significant correlation with age after correction for multiple comparisons. It is therefore plausible to assume that the state of the organism concerning development and aging can be characterized by a single number, such as PC1 score, indicating normalized biological age. This, however, by no means implies that the system state dynamics along the first principal component alone can explain all transcript-level variation caused by altered biology, mutations or gene silencing. Therefore, the positions of the strain-representing markers in Fig. 2(a,b) Table 1 for the summary of survivals). ( ours, yield a massive number of transcript levels, measured in a relatively small number of samples. The accurate PCA inference of the aging signature and the biomarker of age is challenging due to the lack of consistency of PCA in high dimensions 26,27 . The procedure is only appropriate for exploratory purposes since it is prone to over-fitting and false-positive errors 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 further complicated by divergence in experimental procedures, leading to uncontrollable batch effects requiring extensive normalization [28][29][30] . To address this hurdle, we proposed that a robust transcriptomic biomarker of age could also be obtained from a sufficiently large collection of publicly deposited "shallow" datasets from small, independent experiments (a dozen samples each, on average) since all the transcriptomes would have to share the same essential biology of aging. In contrast, the specific experimental conditions and uncontrollable batch differences should be mostly uncorrelated, and thus would comprise "noise" in a combined dataset of sufficient overall size.
To investigate such a possibility, we compiled a comprehensive transcriptomic collection for C. elegans by combining almost all publicly available gene-expression data for aging cohorts into a single database. The resulting "MetaWorm" dataset contains, in total, more than 400 different transcriptomic experiments with N = 3724   Table 1 for lifespans). The markers in (a,b) are colored according to age rescaled to the average lifespan in the groups: 0 and 1 correspond to the L4/adult molt and the mean adult lifespan respectively (see the color bar).
www.nature.com/scientificreports www.nature.com/scientificreports/ nematode samples characterized by G = 4861 of the most commonly expressed/detected genes in the samples (see Methods for more detail on the composition of the dataset and its normalization).
The gene-expression variance in the combined MetaWorm dataset is dominated by batch effects, and hence we do not expect PCA to reveal aging in association with the first principal component in an entirely unsupervised way. Instead, we attempted to identify the aging signature by testing many individual genes for differential expression during aging. Our MetaWorm dataset is sufficiently large to generate the cross-validation ensemble of single-gene association tests using exhaustive random resampling. We further reduced the number of candidate genes by thresholding the transcripts based on the frequency of significant associations in the resampling; we estimate that the chosen cross-validation threshold corresponds to p < 10 −6 uncorrected, or p < 0.005 after Bonferroni correction. The final list of genes robustly associated with aging in the MetaWorm study consists of 327 genes (7% of all genes in MetaWorm): 260 up-and 67 down-regulated with age. We suggest using this gene set as the transcriptomic signature of aging. It is noteworthy that approximately 4000 out of 4861 genes never showed a significant association with aging during the resampling (see Electronic Supplementary Materials for the full list of genes implicated in aging in our calculations).
The transcriptomic signature of age may not be exhaustive, and yet by design, it was reproducible across independent experiments and hence should be useful for future C. elegans aging studies. In our experimental RNA-seq dataset, for example, 902 genes are significantly associated with age rescaled by lifespan (for the same threshold as for MetaWorm, p < 0.005 after Bonferroni correction) out of 4861 genes most commonly detected in the MetaWorm samples. Even though the selection using Bonferroni correction is conservative, the list of significant gene-associations in our dataset is larger than for MetaWorm. The "extra" genes in our list may reflect an association with any laboratory-specific external parameter changing monotonically in time; meta-analysis of data obtained under different external conditions in different laboratories would strongly suppress such age-associations, thus leaving a smaller number of significant hits. There may also be strong correlations among genes governed by a transcription factor whose abundance varies with time within any one laboratory, but which may be lost from the MetaWorm database due to inter-laboratory variation in its induction. The prominence of transcription factors among the genes that are age-dependent would inevitably lead to a gene-set with high internal cross-correlation, and a far higher-than-expected fraction of age-associated genes. The MetaWorm list of 327 candidate is congruent to the list of 902 genes from our experimental dataset, since the Mann-Whitney U test shows that 327 MetaWorm candidate genes are significantly enriched with the genes having the most significant correlation with age rescaled by lifespan among 902 most significant ones in our RNA-seq data. The corresponding area-under-curve (AUC) statistic for the receiver operating characteristic (ROC) curve is AUC = 0.610 ± 0.015, at the significance level p < 10 −30 . This implies that the MetaWorm set of "aging signature" genes very likely includes the same genes that determined PC1 in our RNAseq data, among many more co-varying (and hence partially redundant) genes with age-dependent expression.
A correlation with age does not necessarily imply a causal relation to aging, yet genes correlated with age are usually the primary target in aging studies. As a first approach to inference of the regulators of aging, we checked whether the transcriptomic signature of aging is enriched for the targets of known gene-expression regulators (see Table 3). We used four databases for the enrichment analyses: a curated database for transcription factors and RNA-binding proteins from published high-throughput expression studies in C. elegans, WormExp 31 ; a high-quality protein-DNA interaction network 32 ; and two databases of miRNA-target interactions: the in silico predicted TargetScan 33 and the experimentally validated MirTarBase 34 . Enrichment analysis of the list detected ten hits already experimentally characterised as regulators of aging: DAF-16 35 , ELT-2 36 , ELT-6 37 , PMK-1 38,39 , PQM-1 40 , NHR-1, NHR-10, NHR-86 41 , let-7 42,43 , and miR-60 44 .
Universal transcriptomic biomarker of age. The robust nature of the gene-expression signature of age across widely varying experimental conditions suggests that at any given time, the organism's aging state can be characterized by a single number representing biological age. A typical approach for biological age modeling relies on linear regression of measured parameters, gene expression in our case, on chronological age. Naturally, due to the high internal cross-correlation among the gene expression levels and a limited number of samples, the multivariate regression problem is ill-defined, and any number of convenient biomarkers of age could be obtained by applying additional constraints to the regression problem. In this work we impose a requirement of sparsity on the transcriptomic biomarker of age, i.e., the number of genes used in the biomarker model should be minimized while preserving its predictive power. This is possible to do by performing a cross-validated lasso regression of gene expression in the MetaWorm dataset on age rescaled by lifespan in the transcriptomic signature of aging comprising 327 genes. To ensure that the obtained transcriptomic biological age model is not over-fitted and hence retains its predictive power, we have not used our experimental RNA-seq data during training. The final version of the sparse transcriptomic biological age predictor comprises the contributions of only 71 genes (see Electronic Supplementary Materials for the list of regression weights). simultaneous temporal scaling of survival curves and aging trajectories. The biological age predictor can now be used to transform our multi-dimensional experimental RNA-seq data representing every sample as a single number, the biological age. In Fig. 3(a,b) we plotted the aging trajectories (the dependence of the biological age on chronological age) and the survival plots. We choose to plot the age-dependent quantities not as a function of age, but as a function of age rescaled by lifespan, in contrast to Fig. 4(a,b) where the same data are plotted as a function of age without rescaling. The scaling transformation works exceptionally well and simultaneously brings together the survival curves and the aging trajectories of animals with drastically different average adult lifespans: from 17 days for wild-type N2-DRM control worms to 160 days for age-1(mg44) mutants (the Pearson correlation of the biological age with age rescaled by lifespan is r = 0.86 (p = 2 · 10 −18 ), cf. r = 0.54 (p = 8 · 10 −6 ) for the correlation of the biological age with chronological age (see Fig. 4(a)).
It is worth noting that the biological age naively inferred from the small dataset as a first principal component score from Fig. 2(a,b) is not sufficiently accurate to reveal the temporal scaling of the aging trajectories in the same experiment.
Mortality deceleration. The scaling property of the aging trajectories and the survival functions can be naturally explained using the "aging-at-criticality" model, providing a coarse-grained description of the biological age variable and gene expression dynamics with the help of a simple stochastic Langevin equation and allowing an analytic solution for mortality and the survival fraction 12 .
Early in life, up to approximately the average lifespan, mortality increases exponentially with age, where M 0 is the initial mortality rate (IMR). The Gompertz exponent α is determined by the regulatory network stiffness and is inversely proportional to the mortality rate doubling time (see Methods for a summary of the model). We predicted, however, that the exponential rise in mortality rates would cease at late ages, approaching a plateau determined at the value of α 12 .
High-quality mortality records 13 were used to test the theoretical prediction. In that study, nematodes were subjected to various life-shortening stresses and had their lifespans reduced in such a way that any two survival curves could be transformed one into another by a simple rescaling of age. We computed the approximate values of the mortality rate doubling exponent using the data in mid-life and the mortality plateau estimates later in life for all the reported conditions (see Methods for details of the calculations). The results, summarized in Fig. 5, demonstrate a remarkably tight correlation between the parameters, in good semi-quantitative agreement with the theoretical calculation, across a life-span range of almost two orders of magnitude.
Using the signature of aging to identify novel life-extending pharmacological interventions.
The results presented so far confirm our conjecture of an association between aging and critical dynamics of the underlying regulatory network. Aging appears to be a consequence of intrinsic instability manifesting itself as lack of dynamic control over the expression of genes comprising the signature of aging. We therefore predicted that interventions exerting perturbations opposing the aging change in the animals would reduce the rate of aging and extend lifespan.  Table 3. Transcription factors, RNA-binding proteins and miRNAs whose targets are significantly enriched in the list of potential transcriptomic biomarkers of age (327 genes). a Longevity-related regulators are highlighted in bold. b The top six transcription-factor hits are reported. c,d The top five miRNA hits from each of the databases.
www.nature.com/scientificreports www.nature.com/scientificreports/ To illustrate the predictive value of our gene-expression signature of age, we attempted to identify novel life-extending pharmacological interventions by comparing the signature of aging from this study with gene expression profile changes in response to pharmacologic perturbations from the Connectivity Map (CMAP) database 45,46 . This is no trivial task since most of the available transcriptomic studies represent the results of experiments characterizing the effects of drug compounds in human cancer cell strains. We transformed the list of genes associated with aging in C. elegans into the form recognizable by CMAP and obtained a list of prospective medicines with gene expression signatures opposing the aging direction and thus presumably capable of reversing the progression of aging in these animals (see Methods for the necessary details). A similar approach for drug repurposing against aging has been already demonstrated, e.g., using human brain tissue transcriptomics as the input 47 and in 48 We observed, however, that the list of the predicted compounds turned out to be sensitive to minor variations in our worm-to-human gene conversion pipeline and the difference between the latest CMAP versions 49 .
From the top-10 list of predicted compounds (see Table 4) we selected 5 drugs, spanning the entire range of p-values: anisomycin, camptothecin, alsterpaullone, azacytidine and metamizole. Because camptothecin did not pass our initial screen for a sparing effect on aggregation (see Methods), it was not tested for lifespan effects. The  Table 1 for the lifespans): red for small and green for large lifespans. Overall, the scaling spans almost tenfold variation in median adult lifespan, ranging from 17 to 160 days.  Table 1 for the lifespans): red for small and green for large lifespans. Overall, the adult lifespan ranges from 17 to 160 days.
www.nature.com/scientificreports www.nature.com/scientificreports/ other four compounds were assessed in the C. elegans lifespan assay using two concentrations (1 μM and 10 μM) at 20 °C; see Table 5 for a summary of the experimental results. All four compounds turned out to be more effective at the lower concentration, which suggests toxicity at the higher dose, probably due to off-target effects. In Fig. 6, we show the survival curves in the respective cohorts. Remarkably, temporal rescaling of survival curves accounts for all variation in survival of the stocks treated with these drugs at both concentrations (compare Fig. 6  to 7 and 8).

Discussion
The key findings from this study are reduction of gene-expression dynamics to a one-dimension manifold revealed by principal component analysis (PCA), the stability of the aging signature across biological conditions, scaling self-similarity of both aging transcript trajectories and survival curves, plateauing of experimental mortality at the predicted level of the Gompertz exponent, and identification of new potential life-extending pharmaceutical treatments.
The observed features of aging dynamics can be explained with the help of an "aging-at-criticality" hypothesis 12 . This hypothesis proposes that the gene regulatory networks of most species operate near an order-disorder bifurcation point 50 and are intrinsically unstable. In close proximity to the bifurcation, the dynamics of an organism's physiological state are effectively one-dimensional. Such a reduction of physiological-state vector trajectory in the multidimensional gene-transcript space lets us quantify aging progress by a single stochastic variable representing biological age. This natural indicator of the organism's aging is directly associated with mortality (see Methods for further details). This property of the underlying gene regulatory network is a common feature of complex networks; no matter how complex and large the network is, one can characterize the system by its natural state and control variables, thus effectively describing the system by a one-dimensional nonlinear equation 51,52 . In doing so,   www.nature.com/scientificreports www.nature.com/scientificreports/ one introduces an effective organism-level parameter combining all microscopic features of a network topology into a single number, stiffness or resilience of the network (a counterpart of the rate of aging). The effective state variable is another macroscopic parameter (the order parameter), that plays a role of the aging process indicator and can be ascribed the meaning of biological age.
The apparent stability of the aging signature across vastly different biological conditions is not surprising from the theoretical perspective since the lifespan and aging signatures are related to the smallest eigenvalue and the corresponding eigenvector of the gene-gene interaction matrix, respectively. Small perturbations, such as effects of mutations or RNAi, thus produce small shifts in the already small eigenvalue and hence may yield very large variations in lifespan. At the same time, alterations of the aging direction by a weak perturbation are expected to be small. The conclusion is rather general and applies to aging in other species, see, e.g., PCA of gene expression levels in normally fed and calorically restricted flies 53,54 .
Given the robust and effectively one-dimensional character of changes during aging, a sufficiently large dataset could be used to produce a universal transcriptomic biological age model, such as, in its simplest form, a regression of the gene expression levels on the chronological age, suitable for future aging studies in C. elegans. The magnitude and sign of contributions of individual transcripts to the biological age are not unique due to high    Table 5) from red to green for the smaller and for the larger values of lifespans, respectively.
Scientific RepoRts | (2019) 9:7368 | https://doi.org/10.1038/s41598-019-43075-z www.nature.com/scientificreports www.nature.com/scientificreports/ covariance within each subset of coordinately expressed genes. The model should be fixed by any reasonable additional requirement, such as sparsity, and hence has no direct physical or biological meaning other than providing a convenient tool for experimental data analysis.
The biological age as a function of chronological age defines the aging trajectory and can thus be used to distinguish the progression of aging across strains. Early in life, up to approximately the average lifespan, the age-dependent rise 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 interpret if the influence of stochastic forces is small beyond a certain age, 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 theory, this happens whenever the average lifespan greatly exceeds the mortality rate doubling time. In practice, the assumption is only qualitatively correct, but still provides a reasonable explanation of the experimental situation. According to the model, the temporal scale is defined by the underlying gene regulatory network stiffness and thereby mechanistically relates the organism-level properties, such as lifespan, with potentially modifiable molecular-level interaction properties of the underlying regulatory network, such as the characteristic molecular and genetic damage and repair rates 54 .
We expect that the dynamics of gene expression and mortality should increasingly depend on non-linearity of the gradually disintegrating gene regulatory network, as the aging drift and stochastic forces perturb it. Large deviations of gene expression levels from the youthful state are incompatible with survival. Hence the stochastic dynamics of the biological age variable provide a mechanistic coupling to mortality. The scaling universality of the variation in gene transcript levels, along the aging trajectory exemplified by Fig. 3(a), should in turn be the molecular basis for the temporal scaling of survival curves 13 . In that study, nematodes were subjected to various life-shortening stresses, and had their lifespans reduced in such a way that any two survival curves could be superimposed by a simple rescaling of age. Our survival data with life-extending mutations, RNA interference, and pharmacological interventions, follow the same pattern. The scaling transformation works exceptionally well and brings together the aging trajectories of animals with drastically different average adult lifespans: from a mere 17 days for wild-type N2-DRM control worms to 160 days for age-1(mg44) mutants. Whether the temporal scaling of aging trajectories can be generalized to life-shortening interventions has not been investigated yet, and may be complicated by the multiplicity of pathways whose disruption curtails lifespan.
The temporal scaling of aging trajectories and survival curves is, of course, an approximate statement, since gene expression and lifespans also depend on random environmental and endogenous factors. This may be an explanation behind the deviations from the temporal scaling of survival curves in different replicates of the same  Table 5 for the lifespans): red for small and green for large lifespans. (2019) 9:7368 | https://doi.org/10.1038/s41598-019-43075-z www.nature.com/scientificreports www.nature.com/scientificreports/ strain in ref. 13 . Certain external conditions or therapeutic interventions, in principle, may produce a developmental delay or acceleration without a change in the rate of aging, and thus produce deviations from the universal scaling behavior. The DR1694 strain, demonstrating the most discordant survival curve in our analysis, may be an example of such behavior and deserves an independent study. DR1694 is the only strain harboring two mutations (to genes encoding the insulin/IGF1 receptor and a nuclear hormone receptor). Those two mutated genes interact in allele-specific ways, such that some combinations are short lived while others are quite long-lived. This suggests a delicate balance between the gene products as they impact lifespan.
The biological age should plateau at roughly the average lifespan, which is indeed observed in all our experimental cohorts across a 10-fold range of lifespan difference ( Fig. 3(a)). As the biological age approaches the limiting value, mortality also decelerates and reaches a plateau at the level of the Gompertz exponent obtained from an exponential fit in the age range close to the mean lifespan 12 . Using high-quality survival data 13 , we fully confirmed the theoretical prediction and showed that age-dependent mortality in C. elegans indeed decelerates and reaches a plateau at late ages near the expected level. This phenomenon is not limited to experiments with nematodes 55 and presumably underlies the plateau in mortality rates observed previously in very large populations of medflies and fruit flies [56][57][58][59][60][61] , which we have extended here to relatively small and homogeneous populations of C. elegans (Fig. 5). The results match expectations and, together with the scaling universality of the aging trajectory, both in transcriptomes, Fig. 3(a), and corresponding survival curves, Fig. 3(b), support our theoretical model.
For humans, the existence of mortality plateaux is a subject of debate [62][63][64][65][66][67][68] , mostly due to lack of survivors to ages largely exceeding the average human lifespan. The lack of such data results in high sensitivity of a human mortality plateau to demographic errors 66,68 . In contrast, for C. elegans high-quality mortality data are available for ages up to double the average nematode lifespan, thus, in this particular case, the artifactual nature of mortality plateaux is unlikely.
Multiple explanations have been proposed, for plateauing mortality at advanced ages [69][70][71] , all involving multi-parametric models. The main advantage of our approach is that, at least in C. elegans, the exponential increase of mortality earlier in life and the saturation of mortality late in life are explained within the framework of a simple reaction-kinetics theory dependent only on a single parameter. This parameter is identifiable with the mortality rate doubling exponent measured at midlife on the population level, and with the underlying regulatory network stiffness on the microscopic molecular levels.
Quantification of aging progress using a single number, such as a regression on age of physiological variables representing an organism state, is gaining traction in the aging research community. One of the most successful models of the kind is "DNA methylation age", which is a weighted sum of DNA methylation features, trained to  . (b,d) The survival curves from (a,c) respectively, as a function of age rescaled by lifespan. In both panels, all markers are colored according to the lifespans of the strains (see Table 5 for the lifespans): red for small and green for large lifespans. (2019) 9:7368 | https://doi.org/10.1038/s41598-019-43075-z www.nature.com/scientificreports www.nature.com/scientificreports/ "predict" chronological age in humans [72][73][74] and mice 75,76 . We note, however, that practical utility of the biological age concept implicitly depends on the assumed lack of variability of the rate of aging in a study. The scaling universality of the aging trajectory reported here suggests that this assumption is not necessarily true, at least in the experiments with C. elegans. Nevertheless, the rate of aging is apparently stable in mice 77 and in humans [78][79][80] , suggesting that lifespan is much more tightly regulated in mammals than in nematodes-consistent with the dearth of spontaneous, very long-lived mutants among "higher" organisms.
The increasing stochastic heterogeneity effects (including the leveling-off of mortality and bioage) help explain when an anti-aging 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 developmental program. At more advanced ages, the stochasticity of the gene regulatory network kinetics takes its toll and leads to increasing phenotypic heterogeneity at every level. Accordingly, we expect that anti-aging interventions at the early stages have a broader and more generic effect on aging across diverse species. In contrast, interventions applied at late ages should be precisely selected to help treat specific conditions of an individual at a given point along the aging trajectory; a consequence of life history in the form of stochastically accumulated errors.
The universal aging signature consists of relatively few genes (less than 7% of all the available transcripts), and these are enriched with the targets of gene-expression regulators that promote longevity via disparate pathways, such as DAF-16 (a well-studied FOXO transcription factor that mediates key longevity pathways) 81 , ELT-2 (a transcription factor regulating downstream components of the intestinal DAF-2/DAF-16 pathway; overexpression extends lifespan 15-25% 36 ), ELT-6 (RNAi extends lifespan 37 ), NHR-1, NHR-10, NHR-86, ZTF-9, let-7, and miR-60 (see Table 3 for the complete analysis of over-representation). It would thus be of interest to experimentally test some of the uncharacterized hits from our lists. We tested whether the pharmaceutical interventions (azacytidine, metamizole, alsterpaullone and anisomycin) predicted to exert perturbations opposing the aging change would reduce the rate of aging and extend lifespan, and showed that they indeed prolong lifespan. A version of the extensive LOPAC compound-database with 1280 entries was already screened for lifespan-extending effects in C.elegans 6 . Nevertheless, three of our four hits (metamizole, alsterpaullone and anisomycin) were not tested there, and the fourth (azacytidine) was tested with a negative outcome, which is probably due to toxicity at the higher dose of 33 μM used for the primary screening. In our experiment, all drugs were more effective at 1 μM than at 10 μM, suggesting some toxicity at the higher dose. This is most evident for anisomysin, which is neutral or deleterious at 10 μM.
Alsterpaullone is an ATP-competitive inhibitor of cyclin-dependent kinases (Cdk1/cyclin B, Cdk2/cyclin A, and Cdk5/p25), and with even greater potency, of glycogen synthase kinase GSK-3β. Through the latter activity, it inhibits pathogenic phosphorylation of tau in Alzheimer's disease, and may have other pathogenic targets. Metamizole, or dipyrone, is an inhibitor of cyclooxygenase-3 (Cox-3), observed to activate opioid and cannabinoid receptors; however, it is not considered to be either an opioid or an NSAID. Clinically, it is employed as an analgesic with antipyretic and spasmolytic properties, but only minimal anti-inflammatory effects. It reduces lipopolysaccharide-induced fever (via prostaglandin-dependent and -independent pathways), and disrupts biosynthesis of inositol phosphate. Anisomycin, also known as flagecidin, is a bicyclic derivative of tyrosine that is produced by Streptomyces griseolus and inhibits peptidyl transferase activity of eukaryotic ribosomes. It secondarily interferes with DNA synthesis, induces apoptosis in diverse cell types, and is also used as an anti-protozoan agent. It would be interesting to see if it preferentially induces apoptosis in senescent cells, like azithromycin. It activates stress-and mitogen-activated protein kinases (SAP and MAP kinases) including Jnk and p38/Mapk. Azacytidine is an analog of cytidine, which upon incorporation into DNA (and possibly RNA) irreversibly binds and inactivates DNA methyltransferases. We note that it may inhibit additional targets, e.g. enzymes or transcription factors that bind cytidine or deoxycytidine. Azacytidine and its deoxy derivative, 5-aza-2′-deoxycytidine, are used in the treatment of myelodysplastic syndrome, and of numerous cancers in which anti-oncogenes have been epigenetically silenced. Given the diverse mechanisms of these drugs, they are quite likely to complement one another in a multiple-drug "cocktail". Moreover, each drug has known, deleterious side effects, which might be avoided or minimized at the low doses evidently required for gero-protection, and especially in drug-combination formulations.
The observed temporal scaling of survival curves and aging trajectories, together with the robust pattern of gene-expression changes associated with aging, appeared to be universal across extremely diverse biological conditions tested in our experiments. From this, we deduce that life-extending effects are achieved by stabilizing the gene regulatory network and by slowing the rate of aging, rather than by qualitatively changing the molecular machinery of the whole organism. This means that the course of aging of the super-long-lived strains can be potentially mimicked therapeutically, and hence eventually would lead to dramatically increased lifespan without detrimental effects. The "aging at criticality" hypothesis emerges as a robust theoretical and practical framework for the understanding of a broad range of aging-dynamic and survival properties helpful for future efforts to identify anti-aging interventions in C. elegans and other species.
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 after that. 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.
Our survival study confirms longevity of the worm strains subjected to the treatments targeting genes known to affect aging in previous studies. The relative lifespan modification effects in some cases proved to be somewhat smaller, which can probably be attributed to the use of mutation instead of RNAi or a different developmental temperature in the original studies, or neuronal resistance to RNAi, which may be required for full life extension (see Table 2 for comparison).
We note that the SR806/SR808 survivals and their controls from Table 1 were run quite a few years ago 1,2 , and while all results with those strains were consistent over a 4-5 year period preceding publication, the adult lifespans have not been repeated recently under our current lab conditions. The control lifespans for DR1694 were anomalously long, but still indicate a roughly 2.5-fold life extension by the double mutant. We are currently planning additional experiments to see if it is the bacterial stock we use to feed "normal" controls, but it deserves an independent study and is left for future research.
Drugs were prepared in small volumes (60-100 μl per 10-cm plate), at levels calculated to achieve the indicated concentrations upon equilibration with the full agar-medium volume. Plates were overlaid with drug solutions and rocked with rotation as liquid was absorbed into agar, 24 h prior to use. Worms were transferred to fresh drug-equilibrated plates daily for 12 days and after that, on alternate days (M-W-F).

RNA isolation.
Synchronized strains of C. elegans were grown on 100-mm 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 83 .
Experimental RNA-seq dataset. RNA-seq reads were mapped to the C. elegans genome (WBcel235, Ensembl annotation) using TopHat 2.1.1 (with-b2-very-sensitive and-GTF options) 83 and gene-level read counts were obtained using the htseq-count software 84 . Low-expressed genes with at least one zero read count per sample were removed from subsequent analysis. Raw read counts were normalized using the upper quartile method 85 and converted to RPKM values using the edgeR library 86 .
MetaWorm dataset. The "MetaWorm" dataset was compiled from almost all publicly available information on gene-expression profiles for aging cohorts of C. elegans from GEO database 87 and annotated with the corresponding worms' ages and strain lifespans. For individual genes represented by multiple probesets, the probeset with the largest signal was used. Gene expression in all datasets was normalized using the YuGene 28 algorithm, which was specifically developed for gene expression comparisons among different platforms. The final dataset represents a 3724 × 4861 matrix (samples-x-genes) and incorporates more than 400 transcriptomic experiments (see Electronic Supplementary Materials). To visualize the composition of the MetaWorm dataset, in Table 6, we specify the top-13 largest datasets comprising in total more than 1000 samples; in Fig. 9, we plot the distribution of the datasets according to the number of samples in them.
Critical dynamics of gene regulatory networks. We focus on transcriptomic data and describe time-evolution of gene expression by a matrix x i n , 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 . We characterize the gene-expression kinetics by a differential equation 12 : t), where all the kinetic properties of an organism representing the gene-gene interactions are encapsulated into the function f(x i , t). The coarse-grained description of aging dynamics can be obtained from the linearized version, = + ∑  x f Kx i j ij j 0 , assuming small deviations from the steady state. Here K ij = df i /dx j is the matrix of interactions among the genes. The stability properties of this matrix determine, whether the corresponding gene regulatory network would be stable (all eigenvalues of K ij are negative), or unstable (at least one eigenvalue is positive).
In systems without evident symmetries, the system's dynamics phenomenology near the critical point, separating stable and unstable regimes (critical behaviour), is that of co-dimension-one bifurcation. More specifically, there is only one of all negative eigenvalues of the interaction matrix K ij approaching zero and becoming small and positive, α > 0. The system's kinetics are mostly associated with fluctuations along the right eigenvector b i of the matrix K ij , corresponding to this eigenvalue. The gene expression variation is dominated by the critical mode associated with the singular eigenvalue of the interaction matrix K. Therefore, the transcriptome of aging animals can be accurately described by a one-factor model i n i i n i n www.nature.com/scientificreports www.nature.com/scientificreports/ where x i is the initial system state. To consider higher-order factor models, a trade-off between model complexity and a dearth of experimental data must be reached: a more complex model would require more experimental data to resolve true model vectors from stochastic noise. The data covariance matrix is highly singular, 〈δx i (t)δxj (t)〉 ~ b i b j exp (αt) and hence the mode vector b i coincides with the first principal component, see Fig. 2(a,b).
The mode variable z n is the order parameter with the meaning of biological age, which follows the stochastic Langevin equation: here the random variable η represents the stochastic effects of external and endogenous stress factors. Over time, on average, the solution of the stochastic equation describes exponential deviations from the initial point, subsequent disintegration of the gene regulatory network, and, eventually, death of the organism. The characteristic time scale in Eq. (2) is defined by the underlying gene regulatory network stiffness α and thereby mechanistically relates the organism level properties, such as the mortality rate doubling time and lifespan, with the regulatory network topology quantified by potentially modifiable molecular-level interactions variables encoded in K ij and characterizing molecular and genetic damage and repair rates 54 . Small perturbations modify the gene-gene interactions and produce a change of already small eigenvalue, α, and hence may result in huge variations in the lifespan. At the same time, alterations of the aging direction, b i by the very same weak perturbation would remain small. Therefore, we expect that aging trajectories corresponding to different lifespans are self-similar and different by a single time scale factor α.
Practically, one calculates biological age by projecting the gene expression data into it with the help of a transcriptomic biomarker of age, a i GEO accession number Number of samples  Table 6. The list of 13 GEO datasets comprising the largest number of samples from MetaWorm dataset, in total more than 1000 samples. The definition of the transcriptomic biomarker a i is not unique, since 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 the transcriptomic biomarker of age a i is the left eigenvector of the interaction matrix K ij corresponding to the eigenvalue α.
Mortality analysis. 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. In the mortality data of appropriate quality 13 , a temporal scaling law of aging in C. elegans was observed, similar to that inferred for D. melanogaster 88 . 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 13 , 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 89 . This is fortunate, since a gene regulatory network'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 90 with the help of the Lifelines package 91 .
Since the mortality rate reaches a plateau at late ages for C. elegans 55 , 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 Fig. 10, 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).
Preparation of the signature of aging for the Connectivity Map screening. To transform the list of genes associated with aging in C. elegans into the form appropriate for the Connectivity Map (CMAP) database 45,46 , we first identified human orthologs for the genes from this list using OrthoList database 92 comprising information from four other databases: Ensembl Compara 93 , InParanoid 94 , NCBI HomoloGene Database, OrthoMCL 95 . Since, CMAP requires human genes to be presented by HG-U133A tags (Affymetrix Human Genome U133A Array), the g:Profiler database 96 was used to map human Ensembl gene IDs to HG-U133A tags. Finally, the lists of up-and down-differentially expressed with age genes were formed and used to predict the list of prospective drugs-perturbagens using CMAP. These drugs are expected to reverse the gene-expression profiles to a younger state. www.nature.com/scientificreports www.nature.com/scientificreports/ Aggregation test. The aggregation tests were done in the AM141 strain, a model of huntingin-like aggregation in which a tract of 40 glutamines is fused in frame to YFP. These worms have only diffuse YFP fluorescence as larvae, but as adults progressively accrue punctate aggregates over about 6 days. Drugs are usually introduced just before the start of adulthood (late L4 stage). We take pictures of the fluorescent signal and plot either aggregate count (using imageJ) or total YFP intensity within foci.