The cell rejuvenation atlas: leveraging network biology to identify master regulators of rejuvenation strategies

A unified processing and analysis pipeline for single-cell based rejuvenation studies

To overcome the abiding issue of heterogeneous processing and analysis approaches between different studies, we propose a unified multiscale analysis pipeline that allows for the characterization and comparison of the effects of rejuvenation interventions. Starting from quantified expression profiles of single-cell RNA-seq experiments from treated and untreated donor samples, our pipeline first filters low quality cells based on dynamic thresholds for the percentage of mitochondrial and ribosomal reads as well as the relationship of read counts to detected genes. Next, the expression profiles of all cells in a dataset after regressing out the effect of cell cycle induced differences and normalization using scTransform [13]. Finally, the optimal clustering of cells is automatically identified by maximizing the Calinski-Harabasz Index [14].

After processing the data, our pipeline analyzes each dataset of treated and untreated samples at different levels of biological organization. Initially, it characterizes differentially expressed genes and the cellular processes they belong to. Subsequently, transcriptional regulatory networks (TRN) among differentially expressed genes by following a previously published method [15]. In brief, assuming an “inhibition dominant” regulatory logic in which one upregulated inhibitor is sufficient to cause the downregulation of a gene (no matter the number of activating relationships), a prior knowledge network (pkn) of TF - gene interactions is pruned to remove interactions in which the gene activity is behaving differently to what is seen in the differential expression profile. As a next step, our pipeline integrates signaling and transcriptional regulation to reconstruct sustained signaling cascades and identify their key molecules using SigHotSpotter [16]. Eventually, we employ InterCom to interrogate intercellular communication by reconstructing cell-cell interactions mediated by ligands and their cognate receptors [17]. In brief, InterCom infers ligand-receptor interactions by modeling intracellular signaling and downstream TF expression to ensure compatibility.

We collected 9 previously published single-cell RNA-seq datasets of heterochronic parabiosis (3 datasets), calorie restriction (1), exercise (1), metformin (2), rapamycin (1) and in vivo partial reprogramming (2) [7–10, 18–22] (Figure 1A, Supplementary Table 1, Supplementary Figure 1). As expected, we observed substantial technical variability in these datasets evidenced by large differences in sequencing depth, which further underscores the need for a homogeneous data processing pipeline (Figure 1B). Altogether, the employed datasets span a total of 74 cell types across 18 organs. Notably, tissues from the central nervous system, adipose tissue, liver and bone marrow could be found in multiple intervention datasets. To begin with, we set out to characterize the differential expression patterns of each cell type in each organ in response to the individual rejuvenation strategies and identified a considerable heterogeneity (Figure 1C). While systemic interventions such as calorie restriction and heterochronic parabiosis consistently exert large effects on the transcriptome of multiple organs, metformin has little to no effect on the organs it has been examined in. Interestingly, although exercising is directly affecting the muscles by diverting blood to them, the largest transcriptional effects were observed in the liver, artery and spinal cord.

Identification of transcriptional master regulators that mediate rejuvenation effects

In order to gain insights into the regulatory relationships explaining the observed differential expression profile, we reconstructed TRNs among differentially expressed genes for each cell type in different tissues. As a result, we obtained 317 TRNs of cell types that were affected by a rejuvenation intervention. On average, TRNs are composed of 72 genes (range: 2–867) although most networks contain less than 35 genes (Figure 2A). Interestingly, the size of the TRNs is only weakly related to the number of differentially expressed genes (Pearson correlation, r = 0.21, p < 0.001), which suggests that the transcriptional response to rejuvenation interventions is dependent on other regulatory mechanisms. In this regard, we hypothesized that signaling dependent TFs, whose activity is not only mediated by their expression level but also extracellular signals, may regulate the genes that cannot be explained in the TRNs. Indeed, using a previously curated collection of signal-dependent TFs [23], we found between 77.8% and 100% (median 87.97%) of genes that are differentially expressed, have a known potential regulator but are not part of the TRNs can be regulated by signaling dependent TFs. In addition, we assessed how hierarchical each network is using the Krackhardt Hierarchy Score and found the TRNs to be highly hierarchical (average: 0.994, range: 0.932–1) (Figure 2B). This indicates the presence of a few ‘master regulators’, i.e., TFs that explain a large fraction of gene expression changes (see Methods for details).

Based on network statistics, we sought to identify these master regulators and simulated the downstream effects of activating a single TF in the network to assess the number of genes whose differential expression could be determined by this gene alone. Thus, a TF with a score of ‘1’ determines all genes in a network while a TF with a score of ‘0’ determines no other gene. Following this approach, we detected 493 TFs with a non-zero score across all cell types, organs and interventions. However, the majority of these TFs only act as master regulators in less than 5 conditions (Figure 2C). Moreover, the master regulator score of many TFs is low across the majority cell types and interventions. For example, the rejuvenation response in adipocytes after exercising is orchestrated by the co-expression of Clock and Arntl, which induce different downstream factors depending on the organ of origin. On the other hand, Nfkb and Esr1 regulates varying fractions of differentially expressed genes depending on the intervention (Supplementary Figure 2A). Indeed, it is not uncommon that in different conditions both shared and distinct mechanisms are found, suggesting therapeutic approaches to be equally promising through similar mechanisms. In basal cells of the Skin, for instance, Srf, Cebpb, Atf4, Jun and Myc shared the majority of their downstream regulatory target genes whereas other TFs mostly acted in a non-overlapping manner (Supplementary Figure 2B). Similar patterns could also be found in different cell types of the same intervention, with Ddit3, Spib and Cebpb mediating the effect in the granulocyte lineage while the remainder of the transcription factor response is determined by the maturity of the cell (Supplementary Figure 2C). Not surprisingly, we also observed distinct mediators of the intervention response. For instance, Ybx1, Klf4 and Ets1 were found to be master regulators of exercise and calorie restriction in hepatocytes, whereas only Foxo3 attained a high master regulator score in case of heterochronic parabiosis (Supplementary Figure 2D).

Next, we aimed at interrogating the most common intervention mediators and selected the 30 TFs that have the highest average master regulator score across all cell types (Figure 2D). Surprisingly, when contrasted against previously existing analysis that documented substantial declines or increases in TF activity with ageing [12], the overlap with those TFs is limited; with only 4 of our 30 mater regulators appearing in such an analysis (Nfkb1, Irf1, Arntl and Id3). Moreover, the sign of the change in TF activity varied depending on cell type, rather than being consistently positive or negative. This would suggest a marked distinction between the regulatory agents associated with age and those able to orchestrate the rejuvenation response.

Interestingly, our master regulator TFs have been previously associated with diverse cellular functions, including differentiation, proliferation, immune response and cell migration. Intriguingly, when grouping these TFs by their master regulator score in every cell type, we observed the presence of several clusters. As a general observation, we conclude that our master regulators rather group by intervention instead of cell type. In light of the diverse set of enriched cellular processes, this suggests the induction of distinct signaling pathways that differentially activate broad TFs. However, there exists one cluster that almost exclusively contains immune cells after treatment of heterochronic parabiosis or calorie restriction (Figure 2D; right part). Although all of these TFs contribute to the mediation of the intervention effects, Jun, Junb, Jund, Atf4 and Fos, all of which belong to the AP-1 transcription factor complex, display a consistent involvement (Figure 2E). In addition, we observed that many smaller clusters are formed that predominantly contain a single master regulator even though the AP-1 complex TFs mostly co-occur. The presence of the AP-1 complex across multiple interventions prompted us to dissect the involvement of all known subunits [24] separately. Consistent with the known dimerization patterns of the AP-1 complex [25], Fos and Jun clearly emerged as the most common master regulators across interventions and cell types. In contrast, other Jun-, Fos- and Atf-family TFs act more selectively as master regulatory, which leads us to hypothesize that the cell type and intervention specific effects are exerted by distinct AP-1 dimers. Interestingly, although previous studies have documented the influence of this complex in promoting age-related inflammation (“inflammaging”) [26], our analysis strongly suggests their action as anti-aging mediators depending on their dimerization.

To support the involvement of the identified master regulator TFs, we cross-referenced the 30 TFs having the highest master regulator score with aging-associated genes contained in GenAge [27]. As a result, we found 53% (16/30) master regulators to be linked to aging with varying degrees of evidence. In particular, Arntl, Cebpb, Foxo1 and Jun possess strong evidence and have been directly linked to aging in mammalian and non-mammalian model organisms. In addition, the gene products of Myc and Nfe2l2 have been directly linked to aging in a cellular model system and Foxo3 has been shown to be involved in human longevity. Interestingly, several genes, i.e., Ar, Egr1, Jun and Sp1, have been shown to regulate genes previously linked to aging. Less evidence is provided for Hif1a and Nfkb1, which are known to be involved in pathways or mechanisms linked to aging. Finally, Ddit3, Fos and Stat3 are known effectors of aging-related genes.

Despite the TFs that have been found in GenAge, we collected publicly available transcriptomic perturbation data of the top 30 master regulators and applied MultiTIMER [28], a multi-tissue transcriptional aging clock, to quantify potentially rejuvenating effects. Due to the nature of MultiTIMER as a predictor of transcriptional age in bulk data, we chose to validate the master regulators found in SINGULAR with experiments available in the Gene Expression Omnibus (GEO). In particular, we found knockdown/knockout experiments for Klf4, Irf1, Atf4, Myc, Hif1a and Esr1 in cell types where they have been identified as master regulators (Supplementary Table 3). Since master regulators are up-regulated upon rejuvenating interventions, we would consequently expect the age of normal cells to increase after their knockdown or knockout. Indeed, we observed increases in the predicted cellular age after perturbing Klf4 (9.1 years), Irf1 (3.9 years), Hif1a (2.6 years). In addition, we also collected transcriptional profiles after overexpression of Klf4 and Myc. Interestingly, the predicted cellular age after Myc overexpression is considerably younger (−9.7 years) whereas it slightly increased in case of Klf4 (1.6 years). These results suggest that master regulators act synergistically with other TFs to exert a rejuvenating effect in a cell type dependent manner. This idea is further supported by current partial reprogramming strategies that upregulate Klf4 and Myc in combination with Pou5f1 and Sox2 to achieve a significantly higher reduction in cellular age compared to what we observed for Klf4 or Myc alone [29].

Crosstalk between transcriptional master regulators and intracellular signaling response

Most of the TFs with the highest master regulator potential across interventions are well known to be activated or inhibited by multiple signaling pathways. Thus, we set out to identify the active signaling molecules that are likely to mediate the activation of master regulators in each cell type, tissue and intervention, as described before. Moreover, we compared the accordingly detected signaling molecules in treated and untreated samples to select those that are differentially active between both conditions. As a result, we identified 452 molecules in 33 cell types across conditions and organs. Of these molecules, 74 directly activated the downstream TFs of the corresponding TRNs after treatments. The full set of results can be viewed in the SINGULAR database.

To interrogate the function of the integrated signaling cascades and their induced TRNs, we performed Gene Set Enrichment Analysis (GSEA) of their constituent genes [30]. Strikingly, 23 integrated networks are negatively enriched in aging gene sets derived from Tabula Muris Senis, as found in The Molecular Signatures Database (MSigDB) [31]. This finding supports the rejuvenating effects of different interventions on cell types on the basis of an independent dataset (Figure 3A–3C). After heterochronic parabiosis however, neutrophils predominantly displayed a pro-aging signature (Figure 3D). This is consistent with a previous study reporting an increase of a gene signature of neutrophil activation in hematopoietic stem and progenitor cells after parabiosis treatment [8]. Of note, especially in the case of parabiosis, we observed a reduction of aging signatures derived from various tissues. Nevertheless, in other cases, the enriched signatures precisely matched the tested cell type, such as in case of lung B cells where Grb2 mediates the activation of the transcriptional master regulator Fos (Figure 3B, Supplementary Figure 3A).

Despite the reversal of pro-aging signatures in multiple cell types, we observed other enriched functions that are consistent with our current understanding of the interventions. For instance, E2f target genes, cell cycle control and the G2m checkpoint were negatively enriched after partial reprogramming of the muscle, which displays an expected decrease in cellular proliferation that is consistent with previous studies indicating an up-regulation of the cell cycle inhibitor Cdkn1a in the critical treatment window [32] (Supplementary Figure 3B). Moreover, in the same experiment, cardiac muscle organogenesis was positively enriched suggesting an ongoing re-commitment to a fully differentiated state after de-differentiation (Figure 3C, Supplementary Figure 3B). Interestingly, while Neutrophils in the bone marrow show signs of overall rejuvenation in response to heterochronic parabiosis, their inflammatory potential increases (Supplementary Figure 3C). This implies that parabiosis could have harmful effect on aging alongside its rejuvenating benefits. Finally, the response of T cells in the peripheral blood to Parabiosis2 was mediated by the kinase Pak2 and the master regulators Nfkb1 and Stat3. As a result, we observed a decline in the Pi3k/Akt/mTOR signaling pathway, which is a well-known aging determinant. In fact, hypomorphic Pi3k mice show an increased longevity (Supplementary Figure 3D).

Integration of gene network inference, signaling and intercellular communication analysis

Based on our finding that the transcriptional response to several interventions could be linked to sustained signaling cascades that get activated, we finally aimed at interrogating whether these effects are induced by ligand-receptor mediated cell-cell interactions. Thus, we employed InterCom [17] to reconstruct the cell-cell communication networks of each treated and untreated organ. Similar to our assessment of the intracellular signaling cascades, we focus in the remainder on the most significant interactions that are unique to the treated condition and that involve a receptor as a key signaling molecule. The complete information can be accessed in the SINGULAR database.

Deriving mechanistic insights and testable hypotheses from the response to different rejuvenation strategies could significantly accelerate the development of new anti-aging treatments. Therefore, we focus on two illustrative examples to demonstrate the potentiality of our unified analysis approach. First, we found our analysis to recapitulate known cell communication effects in macrophages after heterochronic parabiosis (Figure 4A). More specifically, our analysis revealed that Gnai2 activates the AP-1 complex genes Jun and Fos upon dissociation of Ccr2 upon recognition of its cognate ligand Ccl2 [33]. Activation of the AP-1 complex in turn leads to the up-regulation of known chemotaxis related genes, such as Cd14, Cxcl2 and Vegfa [34–36]. Indeed, positive enrichment of chemotaxis-related gene sets in a GSEA of the signaling cascades and downstream TFs underscored the chemotactic expression program induced by Ccr2 (Figure 4B). As a second example we chose to illustrate a novel, non-canonical signaling cascade that has not been reported before. In response to exercising, Purkinje cells in the cerebellum form an autocrine loop and interact with oligodendrocyte precursor cells via the Fgf10-Fgfr2 axis (Figure 4C). While our analysis recapitulates the downstream activation of Runx2, which in turn up-regulates and guarantees the expression of Fgfr2, we found that Pax6 is activated by Tcf12 in response to Fgfr2 activation. Although the function of Pax6 has not been reported in Purkinje cells, it is a known neuroprotective transcription factor [37].

Identification of potential drugs targeting key TFs and signaling molecules

In order to demonstrate the utility of SINGULAR, we asked whether we can determine drugs that can target the identified TF master regulators and key signaling molecules. For that, we collected all available drug-target relationships in DrugBank and searched for drugs that could activate our master regulators or mimic the effect of rejuvenation interventions on key signaling molecules. For this purpose, we classified TFs as master regulators if they determine at least 30% of the network TFs when activated according to our simulation studies (see online Methods). Unsurprisingly, of the 239 transcriptional master regulators across all cell types, organs and interventions, only 17 could be activated by drugs (Figure 4D). Moreover, these TFs predominantly belong to the class of nuclear receptors, including Nr3c1, Vdr, Nr1i2, Rxra and Ar. However, notable exceptions are the AP-1 complex proteins Jun and Fos as well Trp53. These further underscores the suitability to interfere with the AP-1 complex to mimic the effect of complex interventions. In order to determine whether any of these drugs possess known rejuvenating effects, we cross-referenced them with DrugAge, a database of aging related drugs [38]. As a result, we found several compounds with demonstrated effects on lifespan in model organisms. For instance, Curcumin, a Vdr agonist, extends the maximum lifespan of D. melanogaster on average by 19.5% at high concentrations and Vitamin D3 extends the average lifespan of C. elegans by 26.8% in a dose-dependent manner. Moreover, Bezafibrate, a partial agonist of Nr1i2, has been shown to increase the average lifespan of C. elegans by 13%. In contrast to TFs, the differentially active key signaling molecules between treated and untreated conditions are generally better druggable (Figure 4E). In particular, the microglia specific key signaling molecules after parabiosis App and Mapk14 are targeted by 24 and 56 drugs, respectively. However, none of the identified molecules target both genes.

Unified processing of the sc-RNAseq rejuvenation datasets

Data analysis was done in Seurat version 4.3.0 (R versions 4.2.3 and 4.3.0), following a unified, biologically informed approach in preprocessing similar to Subramanian et al. [45].

All datasets were processed from cellranger matrix, feature and barcodes when possible. The main exceptions were the second Parabiosis dataset [8], which was processed from FASTQ files with cellranger 6.1 with all argument set to default except an explicit call to –include-introns to get processed cellranger files, and the second Reprogramming dataset [18], in which the raw count matrix after quality control offered by the authors was used for all downstream analyses, due to the unavailability of the VectorBuilder sequences necessary to build the raw count matrices from the raw reads.

Every dataset was further processed with a unified preprocessing pipeline in which the median absolute deviation argument of the our filterCells function was set to 3 to filter outlying cells based on mitochondrial counts, ribosomal counts, number of features and number of counts, with the latter two further filtered based in a linear model with the formula = log10(nFeature_RNA) ~ log10(nCount_RNA)) for unified quality control.

Resulting Seurat objects were further normalized with the SCT transform function with the vars.to.regress argument with “C.C. difference” for cell cycle adjustment and vst.flavor set to “v2”. Doublets were removed using DoubletFinder version 2.0.3. Cell cycle scoring for adjustment relied on the cell marker database annotation retrieved from Ensembl and provided in the repository as a supplementary file.

Integration of within-experiment sc-RNA seq datasets and removal of unwanted sources of variation

Integration of the datasets was performed with SCTtransform [13] splitting every Seurat object by biological replicates of the same condition and sex (if any) sex (if only this was available) and condition in every case. rPCA dimensionality reduction was used together with the SCT normalization method.

Clustering of sc-RNAseq data and cell type identification

LSK and Skin in the second parabiosis dataset [8] were downsampled (keeping original files 1–4 for the former and 5–8 for the later) due to the full data causing extreme numbers of clusters with the Monocle3 clustering function.

Every other dataset was clustered in full after preprocessing with the cluster_cells argument of the Monocle3 library version 3.1, with 10 neighbors, leiden clustering, UMAP used for dimensionality reduction and iterating 5 over several orders of magnitude of the resolution parameter. Every iteration was scored with the Calinski-Harabasz Index and the clustering with the highest score kept, except if more than 30 clusters were detected after this optimizing process, in which case the partitions function of Monocle3 was used instead.

Preliminary cell type assignments were done using SCINA 1.2.0, but this was substantially supplemented by manual annotation with the support of the literature, the cell marker database and the panglao database of single-cell gene markers, as well as support from the markers provided by the original authors of each study, if it was provided.

Differential gene expression test for sc-RNAseq data

Differential gene expression analysis was in every case performed using the DElegate R package version 1.10 (https://github.com/cancerbits/DElegate) a wrapper to use the DESeq-2 differential gene expression analysis assigning cells to pseudo-replicates. This choice stems from its positive benchmarking against other methods in the metrics of precision, sensitivity and false discovery rate [46]. Moreover, random assignment of pseudo-replicates for pseudo-bulk analyses is the best option to mitigate the concern of drop-outs that is always present in scRNA-seq data. Final results were filtered to only keep those cases when the average log fold change in gene expression had an absolute value greater than 0.25 or the gene was seen in more than 10% of the cells in that cell type for both conditions. When more than one control condition was provided (such as isochronic parabiosis and old untreated) both were taken as the comparison group against the rejuvenation intervention. Comparisons were always made against the old condition, dropping the young data, if any, from further analysis, even if they were used during the integration.

Gene set enrichment analysis

Analysis of enriched gene ontology functions and KEGG pathways was realized with package WebGestaltR version 0.4.4 and the ORA (overrepresentation analysis) method. Differentially expressed genes at the p-value adjust < 0.05 significance level were used as a query, while all the genes in the original count matrix (i.e., the RNA assay Seurat object rownames) were used as background. As illustrated in the attached repository, the remaining arguments were overrepresentation analysis, false discovery rate corrected with the Benjamini-Hochberg with a threshold of 0.05.

Strategy to identify cell-type specific transcriptional master regulators

Inference of the master regulators uncovered in this study was done using R library GRNOpt, developed for a previous manuscript [15]. Differential gene expression results were booleanized by setting the log fold changes in the rejuvenation condition filtered as described above to 1 if they were positive and to 0 if they were negative. The inhibition dominance logic rule was used for building the transcriptional regulator networks with the prune gurobi function. This method has been utilized in prior studies to investigate cell state transitions, especially in the context of stem cell differentiation. Indeed, this approach is specifically designed to identify master regulators that orchestrate transcriptional shifts responsible for initiating and maintaining changes over time [47].

Every TF in a gene regulatory network was scored by exploring the consequences of initiating the state of every element in the network to undetermined and examining the consequences of activating each TF with a depth-first search algorithm. The fraction of the final state of the network (as described by the booleanized gene expression) that matched all the downstream changes due to this activation was then the TF score, implementing the same inhibition dominance logic rule that is used during the creation of gene regulatory network, in which any number of inhibitory relationships takes precedence (i.e., sets the state to ‘inhibited’) over any number of activating relationships. In this way, a TF in a completely hierarchical connected relationship to the rest of the network that could explain its entire state by being activated would get a score of 1.0, while sparse networks in which not all the elements are connected will always get a fractional score, an end node that was only scored based on predicting its own state would get a score of 1/(N), where N is the size of the network, and an inhibited gene in the booleanized gene expression would get a score of 0, being unable to explain its own state relative to the differential gene expression if activated.

The TF-regulon interactions database from CollecTRI [48] was used as a previous knowledge network. This previous knowledge network was retrieved using the omnipath python library with the omnipath.constants.InteractionDataset. COLLECTRI function with the 10090 argument (to get the interactions collected for mouse). Conversion of Protein IDs to gene symbols was performed with the getBM argument of biomart, with rows with non-matching NA values and those where the consensus stimulation and consensus inhibition entries were not opposites filtered out. Finally, for the Calorie Restriction study (performed on Rattus norvegicus) a further conversion of the mice previous knowledge network developed as described above was performed, from the mouse to the rat orthologs, using the “https://dec2021.archive.ensembl.org/” Esembl mirror during the call to the getLDS function.

Endpoint previous knowledge networks for mice and rat are provided in the repository.

Finally, during the development of the gene regulatory network, duplicate interactions and values in the TF “from” column that were not transcription factors were filtered out, both for the previous knowledge network and the booleanized differential gene expression input.

Discovery of stable signaling networks and signaling hotspots from sc-RNAseq data

For the cell communication results, we used Sighostpotter (https://gitlab.com/srikanth.ravichandran/sigHotSpotter) with a custom modification in the for_plotting_networks_functions.R function .trimResults. Lines 56 and 61 were replaced with res_trimmed <- res_trimmed(res_trimmed(,2)>0.7,) and res_trimmed <- res_trimmed(res_trimmed(,2)<0.3,) respectively, to facilitate downstream analysis. Everything else was as seen in the repository and the modified version is provided with the main repository for a custom install.

sigHotSpotter was selected for its ability to model both canonical and non-canonical signaling networks, which generate locally stable configurations. This modeling approach is ideal for studying the sustained, long-term changes involved in both aging and rejuvenation interventions. Moreover, sigHotSpotter identifies ‘hotspots’ in the cascade that sustain the new cell state.

Results were analyzed with the sigHotSpotter pipeline function with cutoff value set to 30, percentile set to 70, and the RNA assay counts as input matrix.

Signaling networks were kept after following two conditions. First, both the rejuvenation intervention and the control condition needed to have a non-NA results. Then, only signaling networks where the signaling intermediate was activated in one condition and inhibited in another (activation defined as a final score above 0.70, and inhibition as a value below 0.30, NA entries included) and present in the rejuvenation intervention were kept for downstream analysis.

Cross-talk between TRNs and signalling networks

To identify the cross talk between the master regulator transcriptional networks and the Sighostpotter results, filtered SigHotSpotter results were subset by signaling intermediate and the networks in which the same transcription factor shared a node with outgoing relationships both in the TRN and the SigHotSpotter edges object from column were concatenated to compute the connected components using igraph’s components function with the mode “weak” argument. If there was a single component, the elements of the entire shared network were taken, while if there was more than one, only the connected component with the signaling intermediate was used for downstream analysis.

The elements of this combined network were then queried for Gene Set Enrichment Analysis with the WebGestaltR library (version 0.4.4) using the online enrich method GSEA and three further queries from the Molecular Signatures Database retrieved at https://www.gsea-msigdb.org/gsea/msigdb. In particular, the M2 curated gene sets, M8 cell type signature gene sets and M5 ontology gene sets were used.

Cell to cell communication analysis in sc-RNA seq data

Intercom was used to model the communication between cells as it pertains to those events that are only seen in treated rejuvenation conditions. For a complete description of how the tool works, we refer to the original work that introduced Intercom [17] but we offer a summarized version as follows:

As a preliminary step, InterCom generates a scaffold of experimentally validated receptor-ligand cell-cell interactions, later integrated with intracellular signaling networks and gene regulatory interactions. For ligand-receptor interactions, a previously curated database [49] had the extent of its relationship narrowed down to only those ligands annotated as ‘Secreted’ in Uniprot.

For the intracellular signaling network scaffold, interactions from databases both publicly available (Omnipath, Reactome) and with limited access (MetaCore from Thomson Reuters) were collected, choosing those related to signal transmission (phosphorylation and ubiquitination events). For transcriptional regulator interactions, we again used Metacore, keeping only direct interactions that were known to entail activation or inhibition. These three elements compose the scaffold.

Regarding the analysis of the provided data, InterCom calculates an interaction score for each potential cell–cell interaction by multiplying the average receptor expression and average ligand expression in all cells of a population expressing the receptor or ligand, respectively.

The significance of these scores is then assessed by comparing the scores of all potential cell–cell interactions contained in the scaffold between the two interacting cell types. Interactions with scores in the top decile are considered significant and are the focus of our analysis.

As it pertains to specific parameters in our analysis, Intercom analysis was performed taking the SCT integration counts as an input matrix and with the sigcutoff and z.score.cutoff parameters set to 0. Every other argument was left as default.

Results were further filtered to keep only the interactions unique to the rejuvenation intervention and with a significance score above 0.90, to keep the top decile as usual, as described above.

Assessing the effect of MR perturbation on cellular age

To validate the discovered master regulators, we carried out a search of perturbation data the Gene Expression Omnibus (GEO), using ‘knockdown’, ‘knock-out', ‘shRNA’ ‘overexpression’ and ‘knock-in’ as keywords and selecting as ‘Study Type’ both ‘high throughput sequencing’ and ‘expression profiling by array. In addition, in order to keep comparisons informative, we selected datasets that meet the following criteria: (1) The relevant master regulator is perturbed. (2) This master regulator is the only gene perturbed in the experiment. (3) The study does not involve cancer, embryonic, or similarly altered cell lines that would confound the transcriptional age estimation. (4) The master regulator is found in our single-cell data at least once in a closely similar cell type to the one used in the bulk experiment.

All datasets that met these selection criteria and included both control and perturbed data were used as input for MultiTIMER [28], a transcriptional age clock able to generate predictions for any tissue. Identifiers as well as the difference in predicted age can be seen in Supplementary Table 3. Master regulators shown in Figure 2D that are not listed in Supplementary Table 3 either lacked perturbation data or did not meet the specified criteria described above.

Expression data were converted into log2 of the rank of the gene expression and further subset to the intersection of genes present in all datasets before comparing the difference in predicted biological age between experimental data and controls.

Data availability

The following publicly available datasets were used in the analysis of this study: GEO Accession numbers GSE137869 (Calorie Restriction), GSE176206 and GSE144600 (Reprogramming), GSE193093 and GSE222510 (Parabiosis) as well as GSA CRA004660 (Parabiosis) and CRA007207 (Exercise).

Code for the processing pipeline and auxiliary functions in the workflow is available at https://github.com/jarcoshodar/singularsource.

SINGULAR is available as a publicly available interactive database at https://singular.lcsb.uni.lu/. Source code for a local install and exploration of the data is available at https://git-r3lab.uni.lu/mohamed.soudy/singular.