Metabolic Reprogramming Induced by Aging Modifies the Tumor Microenvironment

1. Introduction

Aging is the strong risk factor for tumorigenesis, with more than 60% of cancers occurring attributed to patients aged 60 and over [1]. One of interesting example is that the Werner syndrome, a premature-aging syndrome occurring in adults, had an increasing incidence of cancer [2]. Accumulating evidence indicates that metabolic reprogramming is a hallmark of both aging and tumor development [3,4,5]. Some evidence also indicates that the ability of neurons in the aging brain to utilize glucose was decreased [6]. Warburg effect, a well-known metabolic reprogramming event, takes up glucose by glycolysis rather than oxidized phosphorylation (OXPHOS) even in oxygen-rich environment to maximize biomass and fast provide energy for tumor cell growth and proliferation [7,8]. Furthermore, recent reports suggest that age-driven metabolic reprogramming in different cell types and significantly affect the cellular functions, thus, providing perfect environments for premalignant cells [9]. In particular, dysfunction of one-carbon metabolism led to adequate priming and effector differentiation of pre-effector-like T cells, thereby exhibiting resistance to anti-PD-1 therapy [10,11]. However, the comprehensive investigation into the association between the age-related metabolic alterations and tumorigenesis remains poorly explored.

Herein, we firstly comprehensively depicted the aging-related metabolic plasticity across 17 different cancer types using an in-silico framework developed by our previous study, termed as MMP³C (Modeling Metabolic Plasticity by Pathway Pairwise Comparison) [12,13,14]. MMP³C analysis facilitates us to characterize aged-related metabolic plasticity events associated with tumor initiation and progression, thereby laying the foundation for understanding the significance of metabolic reprogramming in the context of aging and cancer. MMP³C framework initially calculated the MP activity score based on gene expression profiles and built-in pathway PPI network weight value. Then, intra-sample comparison of two MPs was conducted to eliminate the batch effect of samples. Subsequently, significant metabolic switches were identified based on chi-square test among different aging stages. The tendency of metabolic switches between two MPs was evaluated using odd ratio value among different aging stages.

Thus, we investigated the association between metabo-plastic pairs derived from 84 KEGG MPs and aging based on MMP³C using bulk RNA-seq. The results exhibited increased fluxes of glucose-related pathways compared with ATP generation pathway. Meanwhile, single cell transcriptome profiling with approximately 180,000 individual cells were utilized to characterize the age-induced metabolic switches in TME. This revealed that the metabolic dysfunction of immune cells (T cells and macrophages) might lead to TME reprogramming and thus provide prefect environment for malignant cells in elder individuals. In summary, with the increase in life expectancy, understanding the interaction between age-related metabolic reprogramming and immune cells will aid in the development of precise treatment strategies for elderly patients.

2. Materials and Methods

2.1. Datasets

The pan-cancer gene expression, methylation, and clinic information datasets produced using the HiSeq Illumina platform, including 26 solid cancer types were downloaded from the UCSC Xena Database (https://xenabrowser.net/datapages/, accessed on 26 April 2023). Tumor type-specific age quartiles were considered as age group boundaries for TCGA cohorts, defining the younger group and the older group by the lower and upper quartiles, respectively.

The single-cell RNA-seq data (GSE163120) with seven samples of primary glioma was downloaded from gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi, accessed on 26 April 2023) from a previous study [15]. Additionally, cptac-3 is available at DATE from https://registry.opendata.aws/cptac-3. The age group limits of single cell samples were set by fixed thresholds and consistent with samples from bulk sequencing.

2.2. Immune Infiltration Analysis by ssGSEA

The 22 immune cell infiltration distribution was estimated by the immune cell signatures [16] by ssGSEA method. The Gaussian kernel was employed with TPM data.

2.3. DNAm Age Calculation

The epigenetic age of age-related cancer was computed via ‘methylclock’ R package [17]. This tool infers epigenetic age based on 353 CpG sites by Horvath method [18]. In addition, predicted ages with a correlation lower than 0.8 with the internal gold standard were discarded from downstream analysis.

2.4. RNA Age Calculation

The FPKM obtained from UCSC Xena database to determine transcriptional age using the ‘RNAAgeCalc’ R package [19]. The regression model for transcriptional age utilized the Dev signature, which encompassed coefficients for genes that had the largest variance across samples.

2.5. Trajectory Analysis Based on Bulk Sequencing

The pan-cancer cohort was leveraged to infer the tumor development with the age covariate via the ‘PhenoPath’ package [20]. Then, Principal Component Analysis (PCA) was processed via the ‘stats’ [21] R package to visualize the trajectory of tumor development with aging.

2.6. Single-Cell RNA-Seq Analysis

2.6.1. Data Preprocess, Cell Clustering, and Annotation

The ‘Seurat’ [22] R package (V 5.1.0) was used to preprocess single-cell gene expression profile. The VlnPlot() function was used to visualize three quality control indices and filter low-quality cells by: 1) gene number large than 4,500; 2) unique molecular identifier (UMI) count less than 200 or large than 20,000; 3) mitochondrial gene percentage more than 10% 4) the ratio between gene number and UMI count less than 0.8. The doublets also were identified and excluded via the ‘scDblFinder’ R package [23]. After quality control process, the filtered gene expression matrix is further conducted library size normalization by NormalizeData() function in Seurat. The “LogNormalize” method and 1,000 size.factor are selected. Then, top 2,000 variable genes were found using the “vst” method by the FindVariableFeatures() function to conduct PCA. The optimal number of Principal Components (PCs) for further analysis was selected by the PCElbowPlot() function. Subsequently, the FindNeighbors() and FindClusters() functions were utilized for cell clustering with optimal number of PCs. The RunUMAP() and RunTSNE function were performed for cell clustering visualization. The FindAllMarkers() or FindMarkers() function based on Wilcoxon rank-sum test detected differential expression genes (DEGs) in each cell cluster. The p value of 0.05 and log2(fold change) of 1 were selected as the cutoff for significantly DEGs. Then, the DEGs are used to annotate cell subtype based on the CellMarker V2.0 database and the previous study [15].

2.6.2. Integration of Various Single-Cell RNA-Seq Datasets Across Multiple Samples

The IntegrateLayers() function in ‘Seurat’ R package was utilized to integrate single-cell RNA-seq datasets from GEO and cptac3 database. Different parameters consisting of ‘harmony’, ‘CCA’, ‘RPCA’, and ‘scvi’ were performed and thus compared various integrated methods to find optimal method for the batch effect of datasets.

2.6.3. Inference of T Cell Fate by Trajectory Analysis

The trajectory inference analysis, also referred to pseudotime analysis, was implemented using the ‘Monocle2’ [24] R package (V2.3.6) to reveal the T cell-state fates. Before trajectory analysis, the T cell subgroups were divided and further performed T cell subgroup clustering and annotation. DEGs among T cell subtypes identified by Seurat were regarded as the ordering genes in the trajectory analysis. The naïve T cells were set as the root state and performed the orderCell() function in Monocle2.

2.6.4. Cell Type Enrichment of Various Age-Based Subgroups

To quantify the enrichment of immune cell types among different age-based subgroups, we compared the observed and expected cell numbers in each age-based subgroup by calculating the R_o/e (the ratio between observed and expected cell numbers) value via the ‘epitools’ R package (https://cran.r-project.org/web/packages/epitools/index.html).

2.7. Statistical Analysis

All statistical analyses were conducted using R version 4.3.1. Chi-square or Fisher’s exact tests were used to test the difference between two binary variables. All P-values were two-sided and P < 0.05 was regarded as a statistically significant difference.

3. Results

3.1. Aging-Related Metabolic Alterations: Insights from Pan-Cancer Bulk Sequencing Transcriptome Analysis

Initially, the univariate Cox regression analysis was utilized to examine the correlation between aging and survival outcomes (time and status) across 26 solid cancer types from the cancer genome atlas (TCGA) database [25]. To delineate the effect of aging, the covariates such as sex, race, tumor stage, etc. were included the cox survival model. The results showed that 17 cancer types were correlated to age (adjusted p value < 0.05) after eliminating the influence of covariates. Among them, glioma (GBM) with significant p value (2.970067e-15) and low-grade glioma (LGG, 3.329830e-12; Figure 1A) were identified, which was consistent with many previous studies (refs). Then, these age-related cancers were utilized to investigate age-related metabolic plasticity events by comparison of adult and aged cancer patients using MMP³C framework (Figure 1B, C). Briefly, the 3,486 candidate MP pairs (yielded by pairwise pairing between 84 KEGG MPs; Table S1) were subjected to MMP³C toward global screening of metabolic plasticity events. The results displayed that distinct cancer type exhibited a diverse number of significant metabolic plasticity events (Figure 1B), implying an organized heterogeneity of metabolic plasticity across cancers. Nonetheless, many significant metabolic plasticity events were identified in several specific cancers (Figure 1B, C). For instance, a metabolic plasticity event composed of propanoate metabolism and fatty acid degradation was consistently observed in head and neck squamous cell carcinoma (HNSC), LGG, lung squamous cell carcinoma (LUSC), stomach adenocarcinoma (STAD), Uterine corpus endometrial carcinoma (UCEC), and thyroid cancer (THCA), whereas, this metabolic plasticity event reversed in kidney renal clear cell carcinoma (KIRC) and thyroid cancer (THCA; Figure 1B, C). In addition, steroid hormone biosynthesis prevalently increased flux in multiple cancer types. In GBM and LGG, the flux of fatty acid elongation and biosynthesis decreased in the older patients, but the flux of it increased in the younger patients. It might suggest that older patients had low immune in the glioma since lipid biosynthesis was required the activation of T cells [26]. Thereafter, to further identify consensus dysregulated MPs in the pan-cancer cohort, the top 200 significant metabo-plastic events for each cancer type ordered by the adjusted P-value yielded by MMP³C were selected to construct a metabo-plastic network (Figure 1D; Table S2). Subsequent network analysis via Gephi [27] disclosed many critical altered MPs (hub nodes in network; Figure 1D) in the process of senescence, involved the steroid hormone biosynthesis, oxidative phosphorylation, galactose metabolism, pantothenate and CoA biosynthesis, drug metabolism, and glycosylphosphatidylinositol (GPI)-anchor biosynthesis, implying their importance to senescence (Figure 1D). In addition, some of altered MP were validated by the comparison of the single MP activity score computed by MMP³C was compared between old and young groups via ‘limma’ [28] R package (Figure 1E and Table S3), i.e., glucose metabolism, energy metabolic pathways, etc., disclosing a high confidence of MMP³C analysis.

3.2. The metabolic switch correlates to molecular features of senescence in age-related cancers

To examinate the correlation between the metabolic switch and the molecular aging, as well as senescence, initially, the methylation age [17] and RNA age [19] were assessed using the ‘methylclock’ and ‘RNAAgeCalc’ R packages, respectively. Then the Pearson correlation analysis among them was performed. The result showed that many metabolic switches exhibited significant correlation with methylation age or RNA age (Table S4-S5), i.e., Nitrogen metabolism increased in contrast to Glutathione metabolism decreased with the methylation age (R = 0.45, p value < 0.01; Figure S1), fatty acid biosynthesis declined in contrast to fatty acid degradation increased with the RNA age (R = -0.47, p value < 0.01; Figure S2). Furthermore, on the basis of MP pair profile, pseudo-temporal analysis was conducted (chronological age was set as covariate) via ‘Phenopath’ R package [20]. As shown in Figure 2A, we identified distinct pattern of the evolving tumor metabolism among different age-based groups (the patients were categorized into three age-based groups: young, middle, old, based on the lower and upper quartiles of their age distribution [29]. Subsequently, correlation analysis between MP pair activities and pseudo-temporal scores disclosed several metabolic plasticity events associated to the changes in root time (Figure 2B). For example, Taurine and hypotaurine metabolism greatly decreased in the older group, while various types of N-glycan biosynthesis increased compared to the younger group, exhibiting opposing directions of metabolic transformation. Notably, distinct immunological activity among three age-based groups was detected through immune infiltration analysis based on single sample gene set enrichment analysis (ssGSEA) method. The infiltration scores increased with aging progression in most immune cell categories, except for Eosinophil and Monocyte (Figure 2C), in particular, the increased infiltration of effector memory CD4+ T cell was found in young group compared with the remaining two groups (Figure 2C). This result implies that close correlation between age and immunity in patients with tumor.

3.3. The Landscape of Tumor Microenvironment Showed Distinct Discrepancy among Aged-Based Subgroups

Above analysis indicated that aging exhibited a high degree of association with metabolism and immunity in cancer, thereby we seek to explore how interplay of immune and metabolism reprogramming during aging impacts on tumorigenesis using an integrated scRNA-seq data of GBM, which contains 24 primary tumor samples with different chronological ages. The reason that we chose GBM is that GBM has been demonstrated to exhibit a strong correlation with senescence [30]. Also, our previous Cox regression analysis showed that GBM has significant correlation to senescence. Conventionally, all GBM patients were divided into three age-related groups, young (age< 50, n=6), middle (50<=age<65, n=12), and aged (age>=65, n=6; Table S6). After filtering low-quality cells, we established the scRNA-seq atlas consisting of 181,286 qualified cells. The uniform manifold approximation and projection (UMAP) clustering mainly generated nine cell types including T cells (TRBC2, CD3E, and CD3G), oligodendrocytes (PLP1, MBP, and UGT8), neuroendocrine cells (CHGA, NRSN1, NDRG4), macrophages (C1QA, C1QB, and CSF1R), Fibroblasts (COL1A2 and COL3A1), Astrocytes (AQP4, CLU, CHI3L1), OPCs (ALDH1L1, BCAN, MEGF11), Neurons (TOP2A, RRM2, STMN2, DLX6-AS1) and B cells (MS4A1, MZB1, IGKC; Figure 3A, B and Figure S3-4). Younger groups showed an increased presence of the stromal cells (neuroendocrine cells and neurons) and OPCs (Figure 3C). Elder groups had a higher abundance of oligodendrocytes and exhibited increasing infiltration of T cells and macrophages in the tumor microenvironment (TME; Figure 3C), implying a close correlation of T cells or macrophages with aging This result is in line with the previous ssGSEA analysis using bulk RNA-seq data (Figure 2D). Therefore, different kinds of immune cells (T cells and Macrophages) were further curated for cell re-clustering and annotation based on their universal marker gene [15,31,32] and CellMarker2.0 database [33].

UMAP clustering of all 12,144 T cells generated 13 T cell subclusters were identified and nine main cell types were annotated, consisting of central memory T cells (Tcm), effector T cells, proliferative T cells (prol.T), NKT cells, CD8+ T cells, As T cells, helper T cells, naïve T cells and regulatory T cells (Treg) in Figure 3D and Figure S5-S6. CD8+ T cells expressed several pre-exhausted T cell markers (ITGAE, LAG3, GZMH, etc.), which indicated CD8+ T cells were pre-exhausted T cells (Figure S7). Additionally, one cluster of T cells expressed several Astrocytes’ markers (CLU, SLC1A3, CHI3L1, etc.) were termed As T cells (Figure S6). Afterwards, cellular composition analysis showed middle and old groups had a gradually higher abundance of CD8+ pre-exhausted T cells (Figure 3E). However, naïve T cells, activated T cells (NKT cells, effector T cells, Tcm, etc.) are enriched in younger populations (Figure 3E). Tumor-associated macrophages (TAMs) further divided into monocyte-derived brain macrophages (Mo-TAMs) expressed monocyte markers (TGFBI, VCAN, LYZ and CLEC12A) and cells termed as Mg-TAMs expressed microglial signature genes (TMEM119, P2RY12, and P2RY13; Figure S8). Cellular composition analysis revealed that brain TAMs derived from monocytes and microglia didn’t have significant difference (Figure S9). As a result, the younger group had more various immune cell types including naïve cells such as naïve T cells and monocyte (Figure 3E). By contrast, the older had the fewest of activated T cells (Figure 3E). This result is consistent with the previous ssGSEA analysis depicted in Figure 2C. These results suggest that elder individuals recruit more exhausted T cells and fewer activated T cells, leading to a decline in the ability of the immune response to tumor cells.

3.4. The Metabolic Reprogramming of Various Cell Types within TME Exhibits Increased Heterogeneity as the Aging Procession

Given GBM malignant cells originate from glial cells, then, CNV analysis was performed to distinguish malignant cells from glial cells (Astrocytes, oligodendrocytes, and OPCs) using ‘copykat’ R package (details see Methods). The CNV level of all glial cells were evaluated and significantly clustered into malignant and non-malignant cells (Figure S10). By choosing top age-related metabolic switches (adjusted p value < 0.05 based on MMP3C analysis), we found that the metabolic switches of malignant cells were mainly enriched in glucose metabolism and has certain commonalities with each other, such as Glycerophospholipid metabolism and Oxidative phosphorylation (Figure 3G). Some unique metabolic characteristics are also obvious, for example, OPCs were enriched in Ascorbate and aldarate metabolism (Figure 3G). MMP³C analysis showed that malignant OPCs mainly enriched in elder groups exhibited the increased flux of biosynthesis of unsaturated fatty acids and N-glycan in multiple metabolic switches (Figure 4H). By contrast, non-malignant OPCs exhibited Warburg effect and the decline in biosynthesis of unsaturated fatty acids and N-glycan (Figure 4H). The results disclosed the metabolic heterogeneity of GBM and indicate that non-malignant cells in GBM in elder populations might intake energy and thus promote tumor progression via reverse Warburg effect [34].

Next, to investigate which metabolic plasticity implicated in mainly immune cells in GBM, thereafter, T cells and macrophages were selected for further study. Initially, to dissect the evolutionary dynamics of T cell lineages in the aging progression, the pseudotime cell trajectory analysis of nine cell subclusters was constructed, which generated a three-branch trajectory representing differentiation process of naïve T cells into either CD8+ pre-exhausted T cells, Tcm cells or proliferative T cells (Figure 4A). Especially, two-branch trajectory mapped different age-related groups (Figure 4A) and evolved same proliferative T cells with evolution trajectory. Subsequently, we performed MMP³C analysis of each T cell subtypes among aged groups. We screened out numerous significant MP pairs through chi-square testing (adjusted P value < 0.05; Figure 4B). The activity of glucose metabolic pathways increased, and oxidative phosphorylation reversely decreased in CD8+ pre-exhausted T cells (Figure 4B, C). Furthermore, the correlation analysis between trajectory and the activity score of metabolic switches (Figure 4D). The results showed that a metabolic plasticity event composed of Inositol phosphate and oxidative phosphorylation between younger and older groups. Glycolysis and oxidative phosphorylation as the golden standard for tumor metabolism exhibited significance in the T cell’s pseudo-time series, which displayed indeed similar metabolic alterations (Figure 4B, D). The flux of oxidative phosphorylation increased in the young groups but decreased in older groups, which might result in the dysfunction of T cells in elder groups due to the impairment of ATP generation.

Next, we attempted to determine the association between metabolic reprogramming and the function of macrophages within TME. To this end, we initially conducted another round of cell annotation specifically on macrophages via the mRNA expression of inflammatory factors such as CXCR4, SEPP1, CXCL12, CCL2, CXCL8, CXCL10, and IL1β and anti-inflammatory factors such as CD163, FOLR2, and IL10 (Figure 4E, F and Table S7). Following this, we plotted a ternary diagram annotated with the percentage of mainly various cell types within each age group, divided by the total number of cells within that group (Figure 4C). Our analysis revealed that the proportion of CCL3+ TAM cells was significantly higher in the young group compared to both the middle-aged and elderly groups (Figure 4C, F). CCL3 is known for its proinflammatory and anticancer properties in GBM [15,31,35,36]. By contrast, FOLR+ TAM which expressed immune regulatory factors, for example, FOLR, CD163, etc. reversely enriched in aged groups (Figure 4 C, F). Similar to the findings observed in T cells, young patients exhibited superior immune response capabilities, enabling them to react to tumors by supporting proinflammatory processes and maintaining relatively normal energy metabolism. Energy generation pathway (oxidative phosphorylation) significantly increased than glycolysis, arginine and proline metabolism, etc. in the immune regulatory (FOLR+) TAM in adulted group (Figure 4G). TCA cycle and Terpenoid backbone biosynthesis decreased than primary bile acid biosynthesis and Glutathione metabolism in CCL3+ TAM in adulted group (Figure 4G). These metabolic switches, which demonstrate the distinct metabolic adaptations of macrophages, are presented in a way that highlights their potential functional significance.

4. Discussion

An increasing body of evidence indicates the significance of Warburg-like metabolic switches in the growth, survival, and treatment of cancer [9,37]. Aging drives novel metabolic reprogramming and gives rise to Warburg-like metabolic switches, which may help create a favorable environment for tumorigenesis and growth. At the same time, the differential changes relative to aging are to some extent unclear. Our research, through the MMP³C computational framework and rigorous statistical screening, has identified common Warburg-like metabolic switches and key metabolic pathway alterations in age-related cancers.

Aging can impair oxidative phosphorylation through one of several mechanisms. For instance, with advancing age, there is a gradual decline in phosphoenolpyruvate carboxykinase, accompanied by a corresponding increase in pyruvate kinase, leading to a shift in energy metabolism from oxidative metabolism towards anaerobic glycolysis [38]. Our pan-cancer MP profile of tumors in aged individuals reflects this trend compared to younger individuals, with an increase in glycolysis and a decrease in oxidative phosphorylation (Figure 1E). Also, MMP³C analysis of T cells at single-cell resolution supports that the impair of ATP generation pathway in CD8+ pre-exhausted T cells.

At the single-cell resolution, we observed distinct differentiation processes of T cells under the influence of aging and cancer (Figure 4B, D). It has been known that cancer suppresses glycolysis in activated T cells. Through single cell sequencing and pseudotemporal ordering, we found that glycolysis in proliferating T cells was significantly reduced in the middle and old groups compared to the young group. Naturally, the TCA cycle also exhibited a substantial decline. This led to a general decrease in biosynthesis and metabolism of proliferating T cells in both the middle and old groups (Figure 4B). In addition, the immune activated and regulatory TAMs also exhibited similar distribution among age-related groups.

Indeed, there are certain limitations to our research. First and foremost, while our methodology effectively captures metabolic plasticity events with statistical significance, it falls short in precisely dissecting the underlying causes of these events, necessitating additional molecular experimental validation. Secondly, the age-specific conditions under which our single-cell sequencing was conducted resulted in a relatively small sample size of T cells and their diverse subpopulations, potentially impacting the precision of our findings. But this discovery still could enhance our understanding of the complex interplay between aging and cancer metabolism and open new avenues for the development of targeted therapeutic strategies.