Currently, no generalizable tools exist for predicting the impact of interventions on aging using gene expression changes. Here, we build a model for predicting potential geroprotectors by deconstructing the aging transcriptome into co-expression modules and identifying those modulated by lifespan-extending interventions in Caenorhabditis elegans (C. elegans). Many of these modules exhibited adaptive changes during aging that are linked to longevity. Using these longevity-associated modules, we predicted and validated 10 compounds (7 previously unreported) that significantly extended lifespan in C. elegans. RNA-seq analysis of two top candidates also showed targeted modulation of our longevity-associated modules. Importantly, we identified conservation of several modules in mice that were also correlated with frailty. Together, our findings demonstrate that beneficial transcriptomic adaptations occur during aging and that targeting longevity-associated modules offers a novel framework for geroprotector discovery.

RNA-seq and preprocessing:

Freshly prepared bacteria (OP50-1 E. coli) were spotted on 94 cm NGM agar plates on the previous evening and allowed to dry. Synchronized, young adult worms were transferred to fresh plates. For maintaining synchronized populations, FUDR was added to NGM media. Worms are transferred every two days into fresh plates to avoid food exhaustion. Worms are harvested at day 6, 8, 11, 14, and 18 and were washed off the plates into 15 mL tubes, and then washed several times with fresh buffer until a clear solution was obtained. Clean worm pellets were then frozen and later used for RNA extraction. Total RNA was isolated using Qiagen RNeasy micro kit (Qiagen, Hilden, Germany) following the standard protocol. Afterward, RNA was quantified photometrically with a NanoDrop 2000 and stored at -80^oC until use.

For the drug treatment, synchronized young adult worms are transferred into plates supplemented with drugs dissolved in DMSO for 2 days and collected for RNA extraction.

For mouse soleus muscles, we collected tissues at different ages and treatments and extracted RNA using the Qiagen RNAeasy Mini kit (Qiagen, Hilden, Germany) following the standard protocol.

The integrity of total RNA was measured by Agilent Bioanalyzer 2100. For library preparation, an amount of 2mg of total RNA per sample was processed using Illumina’s RNA Sample Prep Kit following the manufacturer’s instructions (Illumina, San Diego, CA, USA). Libraries were sequenced using the Illumina HiSeq4000 sequencing platform (Illumina, San Diego, CA, USA) in a paired-end read approach at a read length of 150 nucleotides. Sequence data was extracted in FastQ format. Sequence data was subjected to quality control FastQC (78), then trimmed using Trimmomatic v0.32 (79) to remove adaptors. Processed/cleaned reads were then mapped with STAR (2.7.11b) (80) and counted via RSEM 1.3.3 (81). The raw count was then imported into R using the R package tximport 1.26.1 (82). To filter out low-count genes, we retained only genes with at least 0.5 count per million reads in at least 6 samples in C. elegans to include genes that are likely to be expressed in at least 1 time point or condition. The remained genes are then normalized using the variance stabilization method implemented in DESeq2 (83).

Statistical Analysis:

Preprocessing the aging and daf-16 RNA-seq data

We used the network construction methodology described in Esmaili et al. (22). To filter out outlier samples, we firstly transformed the raw counts with variance stabilization in DESeq2. Next, we calculated the inter-sample connectivity using Euclidean inter-sample distance using the scaled expression profiles of the top 9000 genes with the highest expression and transformed the inter-sample connectivity into Z-scores.

We removed the outlier in each sample dataset, generally using thresholds Z = -6. For our aging dataset, no sample is categorized into outliers. For the drug combination dataset, two samples are removed as outliers using the same method. To remove potential batch effect within the dataset, we used two different methods: the first one is to use the Surrogated Variable Analysis (SVA) from package sva (84) to detect the within-dataset technical differences. The second method is to use the RUVSeq (85) and detect unwanted variables using the replicate samples. The sva was performed on variance-stabilized data, while RUVSeq was performed directly on raw counts. We next compared the leading surrogated variables and the unwanted factors detected from RUVSeq by calculating the Pearson correlations between variables. The leading surrogated variables that are strongly correlated with at least one unwanted variable from RUVSeq are used for correcting the data. For the aging dataset, the data was corrected with the first two surrogated variables using the Empirical Bayes-moderated linear regression implemented in the EmpiricalBayesLM() function in the WGCNA package. For the drug combination dataset, the data was corrected with the first surrogated variables.

Construction of PCA-based transcriptomic clock

Typical transcriptomics contains a larger number of features than the number of samples, which makes the typical machine learning method unfeasible. We previously described the construction of aging clocks with lipidomic and clinical parameters using principal component analysis (PCA) to scale down the dimension (15, 16). PCA is a dimension down-scaling method based on the singular vector decomposition (SVD), which factorizes the transformed the data matrix into a new coordinate system of feature matrix with features replaced by the principal component (PCs). PCs are linear combination of the initial variables and are ranked based on the original information of variance explained. Therefore, the leading PCs can explain most of the variance of the original data matrix.

We firstly randomly selected 80% wild-type day 6 samples and 5 of 6 sample from each older age as training data, and the unselected samples are considered as test data. The training of the PCA-based transcriptomic clock was carried out on the variance-stabilized training data from the previous section before the correction for internal batch effect, since the latent batch effect or variance can be detected by PCA itself. The wild-type data across different ages were selected and projected into PC coordinates. The elastic net regression model was trained on training data using the top 20 PC coordinates as input and the chronological ages as dependent variables. Model parameters are optimized using 10-fold validation by grid search over an l1 ratio from 0 to 1 in 0.01 intervals. It’s worth noting that biological age is identical to the chronological age in this cohort because each transcriptome sample represents an average value of the cohort of animals. We used the test data, the daf-16(mu86) mutant in the dataset, the time-coursed sterile control fer-15(bs26ts) and fer-15(bs26ts);daf-16(mu86) dataset (GSE122892), and the drug treatment dataset (GSE108263) to test whether the clock can predict the effect of aging, longevity-associated genetic interventions, and drug treatment. For the time-coursed control (fer-15(bs26ts)) and daf-16(mu86) dataset, the PCA plot on the original dataset showed a strong separation between adult day 1 and the rest of the data on PC1 (data not shown), which is possibly due to contamination of worms in developmental stages and is excluded from the prediction. For the testing dataset, we first keep the genes that are used in the PCA training sets and normalized data using variance stabilization method. We next used the singular vectors derived from the training dataset to project each sample’s normalized feature-space coordinates into the same PC coordinate system of the training set.

BitAge application

To apply BitAge to our drug treatment dataset, we firstly converted the gene-level counts to counts-per-million(CPM). We next binarized the data using the make_binary() function provided by BitAge (https://github.com/Meyer-DH/AgingClock/). The biological ages of samples were predicted using the test_prediction() function with the BitAge V2 coefficients.

Construction of the network

For C. elegans, WGCNA was carried out on the variance-stabilized data corrected for the first two surrogated variables (SVs) in the aging/daf-16 dataset and the variance-stabilized data corrected for the first surrogated variables in the drug combination dataset. To improve the robustness of the WGCNA analysis, we used the weighting method slightly modified from biweight midcorrelation, which is firstly described in Lee et al. (86). The weighting method is like the biweight midcorrelation, where the Tukey bi-square-like weight l for each observed gene is calculated as:

Where MAD is the median absolute deviation of the gene observation. The MAD is further adjusted to meet the conditions where (1) the 10^th percentile of the weight is at least 0.1, and the 40^th percentile of the weight is at least 0.9, and (2) the 40% percentile of the weight l for each individual group is at least 0.9.

            We used weighted correlation with individual sample weights determined above and the ‘signed-hybrid’ network which negatively correlated genes are considered uncorrelated. Scale-free topology on the aging/daf-16 dataset showed that the resulting networks are approximately scale-free with the power greater than 5. At power of 8, the mean connectivity is approximately 55, and the median connectivity is approximately 30, which are reasonable thresholds. For the drug intervention dataset, at power of 9, the mean connectivity is approximately 50, and the median connectivity is approximately 20, which are reasonable thresholds. Hence, we used the soft power of 8 for constructing the aging/daf-16 network and 9 for constructing the drug intervention network.

The topological overlap matrix (TOM) was calculated using the given soft threshold and the weights calculated above. The modules were constructed from the TOM with the following settings: deepSplit = 3, minimum module size of 30, merge cut height of 0.12. To improve the robustness of our modules, we also removed genes whose eigengene connectivity to their module eigengene lower than 0.2 from the modules. To improve the stability of our modules, we randomly subsampled approximately 2/3rds of samples in each analysis and created 50 sets of WGCNA module labels, which are used as the stability labels in the Dynamic Tree Cut to improve the stability of resulting modules. We used two different minimum stability criteria (0.55 and 0.35) to construct two different sets of modules. The higher criterion results in larger and more robust modules, which are noted as main modules, while the lower criterion produces finer split or those main modules, which are named submodules. Submodules are mainly used for subsequent analysis.

We next calculated the module eigengenes, which is the first singular vector of the standardized module expression matrix. The fuzzy module membership, which describes the similarity of gene expression to the eigengene of the given module, is calculated using the hierarchicalconsensusKME() function implemented in the WGCNA package. The top 30 genes with the highest module membership within each module are considered as hub genes.

We carried out the same approach on six mouse aging organs (liver, heart, lung, muscle, brain, and small intestine) from the Tabula Muris Senis cohort (https://tabula-muris-senis.sf.czbiohub.org) and the mouse aging soleus muscle data we obtained. For the he raw counts of mouse aging organs from bulk RNA-seq are obtained from the R package TabulaMurisSenisData (87). The 1-month samples are excluded from the analysis since mice are potentially in a developmental stage at the point. Each tissue was then pre-processed separately. We filtered out low-count genes and pre-processed the data using the same outlier-filtering and batch detection methodologies described above. The WGCNA network construction are performed using the same methodology above with an appropriate soft threshold power and a merge cut height of 0.15 for each dataset. For all aging organs, the resulting networks are approximately scale-free, with the mean connectivity lower than 60 and the median lower than 30. Both the main modules and submodules are obtained using the same minimum stability criteria. For the mouse aging soleus muscle, we chose power of 9, deepSplit = 3, minimum module size of 30, and a merge cut height of 0.1 to construct the mouse soleus muscle co-expression network.

Association of Modules with Ages, treatments, and mutations

            For the aging/daf-16 dataset, we used the linear models to relate module eigengenes to individual age contrasts, or mutant contrasts, with the first two singular vectors (SVs) as covariates. For the mouse soleus muscle dataset, we related module eigengenes to individual contrast of age and sex. We additionally treated ages as a continuous variable and related module eigengenes with it, with the first two SVs as covariates. For the drug treatment dataset, we related module eigengenes to individual contrast of drug treatment or mutation versus control, or median lifespan extracted from the original paper, with the first SV as a covariate.

DE analysis

DE analysis was performed using the R package DESeq2 (83) version 1.38.3 using default settings. For the design, we used the unwanted variables derived from the RUVseq to correct for the latent batch effect.

Association of Z statistics for each dataset:

To assess similarity of gene response across different conditions with aging-associated changes, we applied the following strategy: we firstly extracted genes differentially expressed during aging (FDR <=0.05) in our dataset by treating aging as a linear term. We next extracted the Z statistics of all aging-associated genes in each drug treatment versus control and computed the correlation of Z statistics of these genes between conditions using Spearman correlations. The Spearman correlation matrix were plotted as a heatmap with hierarchical clustering using average-linkage hierarchical clustering. For mouse data, we used GSE131901 to validate our hypothesis. We treated male and female mice separately and obtained differentially expressed genes between 12 months and 6 months for each sex following the DE analysis described above. Next, we extracted and pooled the Z-statistics of these genes for each condition versus control, and calculated Spearman correlation in males and females separately. 

Enrichment and annotation of gene sets of DE genes and modules:

We used the R package anRichment (https://github.com/plangfelder/anRichment) to calculate the enrichment of WGCNA modules in a large collection of the following gene sets:

1. Gene Ontology (GO);
2. BioSystems gene sets, including KEGG and Reactome;
3. WormEnrichr  Anatomic Association from Wormbase 2018;
4. WormEnrichr  TF2DNA 2018;
5. WormEnrichr  Human Disease from Wormbase 2018;
6. WormEnrichr  Gene Interaction Hubs BioGRID 2018;
7. WormEnrichr  PPI Hubs BioGRID 2018

Applying the same method, we build a collection of gene sets on mice with the following gene sets: 1. Gene Ontology; 2. BioSystems gene sets; 3. ChEA 2022 library; 4. ENCODE histone modification 2015 library; 5. ENCODE TF CHIP-seq 2015; 6. miRTarBase 2017 for annotating mouse modules. Fisher’s exact test was used to evaluate the overlap significance. The genes used in the network analysis are used as background genes for the enrichment analysis. To annotate the WGCNA modules, we combined the top enrichment from each gene sets that have a minimum fraction of the class size of 0.05 to form a merged annotation on each module. Modules with no enrichment with fraction of class size greater than 0.05 were left empty.

For enrichment of differentially expressed genes (DEG) across conditions, we performed a Fisher exact test in the collection of gene sets we described. We then kept gene sets enriched in DEGs with sizes between 20 and 500 for visualization.

PPI enrichment

The enrichment of the protein-protein network was performed using the method previously described (19). We firstly downloaded protein-protein interaction pairs from the STRING (version 12.0) database and kept only protein-protein pairs with evidence of binding or forming a physical complex. The network was then filtered to keep proteins in WGCNA analysis. We next performed hypergeometric testing with the null hypothesis that the frequency of both members of a protein interaction pair presenting inside a module is not greater that the frequency that at least one of the members of the PPI pair does not present in the module.

GSEA

Gene set enrichment analysis are carried out using the R package ClusterProfiler (88). DE Z statistics of genes were used as input and ranked with a decreasing order. For gene loading on PCs, we calculated the loadings of each genes on PC1 then rank with a decreasing order. GSEA was performed using the modules, or gene sets from GO or KEGG as input.

Module preservation analysis

Module preservation analysis was performed between the aging/daf-16 dataset and the drug combination dataset and between the aging/daf-16 worm datasets and six aging mouse organ datasets, as well as between aging/daf-16 worm datasets with mouse aging soleus muscles. For module preservation between aging/daf-16 worm modules and aging mouse organ modules, human ortholog of worm genes were downloaded from Ortholist 2 (OL2) (89). The mouse-to-human sequence homology are obtained from the Mouse Genome Informatics (https://www.informatics.jax.org/downloads/reports/HOM_MouseHumanSequence.rpt) and are merged with the human ortholog to worm genes to obtain the gene ortholog map between C. elegans and mice. We then filtered the genes to keep only genes with orthologs in both C. elegans and mice by mapping the Entrez identifiers between the two species using the function mapEntrez() provided by package anRichmentMethods (https://github.com/plangfelder/anRichment) with the ortholog map obtained as the reference. 6680 genes in C. elegans are eventually used for measuring module preservation.

We used the modulePreservation function in WGCNA with 500 permutations. The preservation statistics of aging/daf-16 and drug combination submodules in respect to each other are calculated. The Zsummary score is used to describe the overall preservation statistic of modules. Two thresholds of Zsummary are applied: a Zsummary score greater than 8 is considered as strongly preserved, while a Zsummary score greater than 4 is considered as moderately preserved.

Module GSVA association with median lifespan:

To calculate the association of module expression and longevity in the drug combination dataset and in the public datasets we analyzed, we applied Gene Set Variation Analysis (GSVA) (25) to compute relative module expression changes within datasets. GSVA is an algorithm that calculates an estimated variation of gene set expressions over a sample population. We firstly applied Gene Set Variation Analysis (GSVA) to compute GSVA scores based on the variance-stabilized data matrix for the co-expression modules we obtained from the aging/daf-16(mu86) dataset, with the kernel function set as ‘Gaussian’. Next, we obtained the expected median lifespan for each sample in the original published paper and associate the module GSVA score with expected median lifespan using the R package limma (90) version 3.54.2. Samples with no median lifespan reported in the original paper or not designed for lifespan analysis (such as heat shock) are omitted. The standard error of log2FC obtained from limma was calculated from t-statistics using the formula: The modules with Benjamini-Hochberg (BH) adjusted p-values <= 0.05 are considered as significantly associated with expected median lifespan.

For analyzing individual intervention versus control, we associate the module GSVA score with binarized comparison between individual treatment versus control using Limma, and obtained log₂FC for each module as well as the standard error.

Meta-analysis for conserved longevity-associated modules

To obtain conserved modules associated with lifespan across datasets, we expanded our data with publicly available RNA-seq datasets from Gene Expression Omnibus (GEO) (91), ArrayExpress (92), or NCBI Sequence Read Archive (SRA) (93) databases. All the public datasets are pre-processed using the same method described in the previous section. The public datasets are selected based on the following criteria: 1. The dataset contains at least 8 samples; 2. All samples are performed within 3 days of adulthood are kept for meta-analysis. We also extracted the expected median lifespan from the original paper. For the datasets with epistasis study, we primarily used the median lifespan from the epistasis lifespan analysis. The following datasets meet the criteria after filtering and are used for the meta-analysis: E-MTAB-7272, GSE54853, GSE55988, GSE63075, GSE86493, GSE93724, GSE94077, GSE95240, GSE100482, GSE107134, GSE110984, GSE112981, GSE115531, GSE119485, GSE123531, PRJNA523315, GSE146412, GSE150108, GSE72426, GSE173799, GSE239572, PRJNA508378, PRJNA479869.

Most of the RNA-seq experiments in the dataset focused on the mechanism of a certain lifespan-associated pathways, mutations, or interventions. We annotated the dataset with a single term describing the focus of the datasets, such as the main pathways, drugs, or mutant epistasis in the studies. For the dataset PRJNA479869, we split the dataset into three independent epistatic studies based on the epistasis of the mutant tested and used the daf-2 epistasis and glp-1 epistasis separately for the downstream analysis.

Each dataset was pre-processed using the same procedure described in the previous sections. For identification of longevity-associated modules, we firstly moved unwanted variables and potential outliers for each dataset as described in the previous section, and calculated the GSVA scores of each module. We next associated modules with expected median lifespan within each dataset using limma. We obtained the log2FC and standard error of each module with median lifespan for each dataset and aggregated the results using the restricted maximum likelihood (REML) random effects models using the R package metafor 4.6-0 (94). By applying this method, we obtained the aggregated log2FC of each module expression with corresponding p-values, which is further adjusted using the Benjamini-Hochberg (BH) adjustment. To improve the robustness of the meta-analysis, we further calculated the leave-one-out p-value (LOO p-value), which is defined as the highest p-value after removal of every dataset one by one, and the robust p-value, which is the lowest p-value after removal of each dataset one by one. The p-values of modules significantly associated with lifespan across datasets are then converted into Z scores, and a heatmap was generated with hierarchical clustering on modules.

To associating module expression with lifespan extension in individual conditions, we used a method similar to what reported in Tyshkovskiy et al. (14). We firstly collected log2 fold change and standard error of module GSVA scores for each intervention versus control in each dataset and aggregated 83 interventions with median lifespan reported. We next applied the mixed-effects model using the R package metafor to calculate the p-value of the association of module expression change with lifespan effect. We used both the mean and standard error of log2 fold change as input. To remove the batch effect, we considered interventions belonging to the same study and the same control group as a random effect. In addition, we categorized intervention types into three groups: Genetic intervention includes mutation, knock-down, knock-out, and overexpression; dietary restriction (DR) includes eat-2 mutants or dietary restriction; epistasis includes combination of genetic or drug interventions where the lifespan extension is compromised comparing with single genetic or drug intervention. The type of intervention was then used as a random effect for the model. We additionally calculated the LOO-pvalue and robust p-value from the meta-analysis. The p-values calculated are further adjusted using Benjamini-Hochberg (BH) adjustment, and modules with strong association with median lifespan change are obtained based on an appropriate threshold on adjusted p-values.

For visualizing the mutual organizations between longevity-associated modules we determined using the random-effect model as a network, we utilized Cytoscape. We calculated the Pearson correlations between module eigengenes for only day 6 wild-type samples and for all wild-type samples from day 6 to day 18. Only edges with an absolute Pearson correlation higher than 0.5 are shown. The width of the edges is determined by the absolute value of the Pearson correlation, with wider widths representing higher absolute values of correlation.

Querying Connectivity Map:

We selected genes in worm longevity-associated modules M1.1, M2.4, M27, M11, and M20 as input for querying CMap. The human orthologs of C. elegans genes in longevity-associated modules are predicted using the Ortholist 2 (89) and were applied to query CMap. For modules with size over 150, we used the top 100 genes with the highest module membership within the modules as representative genes (hub genes) and converted them into their human orthologs as input for CMap.

The CMap database used pattern-matching algorithms to compare input gene signatures with gene signatures in diverse perturbagens. CMap generated enrichment scores from +100 to -100 for these perturbagens and ranked depending on their similarities or dissimilarities comparing with the input gene signatures. We expect that a potential anti-aging drug would produce negative scores when using aging modules as input and positive scores when using longevity modules as input.

CMap provides perturbation scores on five cell lines. In this study, we used the summary score provided the CMap, which is obtained using the maximum quantile statistics from the vector of normalized connectivity scores of a perturbagen^^ (76), as the final score of a drug’s perturbation on a module. 

DrugAge

The DrugAge database (7) was accessed in November 2023 and was filtered to include all drugs that significantly increase lifespan in any model species.

GenAge gene selection:

The GenAge database was access in November 2023. Since many loss-of-function perturbations of genes results in lifespan shortening due to collapse of the biological system, we filtered longevity-associated genes so that genes are lifespan-associated only in the following condition: 1. genes whose loss-of-function mutation, RNA interference, or knock-out resulting in lifespan extension are considered as anti-longevity genes; 2. genes whose overexpression or gain-of-function resulting in lifespan extension are considered as pro-longevity genes; and 3. Genes whose gain-of-function result in lifespan shortening are considered as anti-longevity genes.