Abstract

   Emerging evidence has unveiled heterogeneity in phenotypic and
   transcriptional alterations at the single-cell level during oxidative
   stress and senescence. Despite the pivotal roles of cellular
   metabolism, a comprehensive elucidation of metabolomic heterogeneity in
   cells and its connection with cellular oxidative and senescent status
   remains elusive. By integrating single-cell live imaging with mass
   spectrometry (SCLIMS), we establish a cross-modality technique
   capturing both metabolome and oxidative level in individual cells. The
   SCLIMS demonstrates substantial metabolomic heterogeneity among cells
   with diverse oxidative levels. Furthermore, the single-cell metabolome
   predicted heterogeneous states of cells. Remarkably, the pre-existing
   metabolomic heterogeneity determines the divergent cellular fate upon
   oxidative insult. Supplementation of key metabolites screened by SCLIMS
   resulted in a reduction in cellular oxidative levels and an extension
   of C. elegans lifespan. Altogether, SCLIMS represents a potent tool for
   integrative metabolomics and phenotypic profiling at the single-cell
   level, offering innovative approaches to investigate metabolic
   heterogeneity in cellular processes.

   Subject terms: Metabolomics, Senescence, Metabolomics
     __________________________________________________________________

   Integrated analysis of metabolome and oxidative stress at single-cell
   level is challenging. Here, the authors develop SCLIMS, enabling
   simultaneous profiling of metabolome and oxidative stress levels and
   discoveries of key metabolites regulating oxidative stress, senescence,
   and healthy aging.

Introduction

   With the advancements in single-cell analysis techniques, the
   heterogeneity of single-cell genetics, proteomics, and metabolomics has
   captivated the attention of scientists for decades. The phenotypic and
   transcriptional heterogeneity of individual cells has been extensively
   elucidated in diverse biological processes, encompassing development,
   cancer, and aging. For instance, through the lens of single-cell
   transcriptomic studies, developmental trajectories are meticulously
   unraveled, unveiling the exquisite tapestry of single-cell
   heterogeneity and sub-populations during embryonic
   development^[46]1,[47]2. The genetic and epigenetic heterogeneity
   exhibited by cancer cells^[48]3 poses significant challenges in the
   field of cancer therapy^[49]4. Delving into the single-cell
   heterogeneity of these cells illuminates promising avenues for
   therapeutic targets and strategies^[50]5,[51]6. It has been reported
   that geriatric organs are mosaics of cells with varying ages, thereby
   highlighting the remarkable heterogeneity exhibited by individual cells
   during the process of aging^[52]7. Another study uncovers the
   heterogeneity of transcriptional landscape in aged ovaries where genes
   undergo specific changes during aging that are exclusive to different
   oocyte subtypes^[53]8. In essence, comprehending cellular heterogeneity
   is paramount in unraveling the intricate cellular mechanisms that
   underlie biological processes, including development, aging, and
   diseases.

   A groundbreaking advancement in single-cell research lies in the
   seamless integration of multiple features or functions within
   individual cells, known as cross-modality analysis. This includes
   integrating single-cell omics with cellular function or phenotype,
   resulting in a truly remarkable and sophisticated
   approach^[54]9,[55]10. For instance, with the Patch-seq technique,
   single-cell transcriptome and electric activities are combined, linking
   transcriptome to cellular function^[56]11. In a study investigating the
   exocytosis and excretion function of pancreatic β-cells in diabetes,
   the Patch-seq technology was employed to establish a connection between
   β-cell functionality, as represented by electrical activities, and gene
   expression profiles. This innovative approach successfully identified
   genes that are closely associated with β-cell dysfunction^[57]11.
   Another study on extracellular vesicles (EV) links lipidomics and
   proteomic with EV membrane function and crosstalk of discrete tissues
   in different stages of COVID-19 infection^[58]12. The integrative
   analysis of metabolism and transcriptome reveals a rare subpopulation
   of hyperactive T cells and subtypes of monocyte which links to disease
   severity in COVID-19, suggesting an association between multi-omics and
   immune cell functionality^[59]13. The combination of single-cell RNAseq
   and fluorescence-activated cell sorting (FACS) illustrates the relation
   between gene expression and stem cell function represented by specific
   surface markers, discovering key molecules associated with self-renewal
   in stem cells^[60]14. Thus, the cross-modality analysis enables a more
   comprehensive view of single-cell omics and a deeper understanding of
   cellular function, shedding light on mechanisms underlying complicated
   biological processes.

   Metabolism serves as both a reflection and regulator of cellular
   function^[61]15. Currently, most metabolomics research remains confined
   to homogenate-level analysis involving the preparation of samples using
   bulk tissue or cell suspension. Consequently, the enigmatic metabolic
   heterogeneity exhibited by cells during biological processes and the
   underlying mechanisms at the single-cell level remain shrouded in
   ambiguity. Recently, our laboratory and others have pioneered the
   development of cutting-edge single-cell mass spectrometry techniques,
   enabling us to unveil the intricate metabolic heterogeneity at an
   unprecedented resolution^[62]16–[63]18. These new techniques mainly
   focused on the interpretation of the metabolome without the phenotype
   of single cells, omitting cross-modality features important for the
   better understanding of cellular metabolism and function. However, the
   integration of metabolic heterogeneity with phenotypic heterogeneity,
   such as distinct levels of oxidative or senescent states within
   individual cells, remains a formidable technical challenge.

   In the present study, we employ a combination of single-cell mass
   spectrometry (SCMS) and live-cell imaging techniques to simultaneously
   capture the metabolomic features and cellular identities of individual
   cells, thereby establishing a link between the metabolome and cellular
   status at a single-cell resolution. By utilizing this approach, we
   design our study using a cellular model of oxidative stress-induced
   senescence, establishing the correlation between single-cell metabolome
   and their oxidative or senescent status. This technique unlocks
   possibilities in the realm of cross-modality analysis, integrating
   single-cell metabolome with fluorescent labeling techniques, thereby
   paving the way for innovative discoveries at the single-cell level.

Results

Integration of live-cell imaging and single-cell mass spectrometry

   To simultaneously capture the metabolome and phenotypic features of a
   single cell, we combined single-cell live imaging with mass
   spectrometry technique^[64]17,[65]18 (SCLIMS) and setup a
   cross-modality analysis platform. We employed a cellular oxidative
   stress (OS) model^[66]19–[67]24 which can be readily induced and
   labeled with dichlorodihydrofluorescein diacetate (DCFDA), a widely
   utilized live-cell probe for the detection of cellular
   OS^[68]8,[69]25,[70]26. Briefly, HEK293T cells were incubated with
   H[2]O[2] in culture medium at a final concentration of 80 μM for 1 h,
   and the medium was then replaced by fresh medium. The cells were
   allowed to recover for 48 h prior to subsequent testing. Cellular
   viability remained uncompromised in the model and was deemed suitable
   for SCMS analysis (Supplementary Fig. [71]1a). The OS level was
   evaluated through DCFDA staining^[72]8,[73]25,[74]26 (Supplementary
   Fig. [75]1b, c). Initially, the cells were incubated with DCFDA for a
   duration of 25 min and then imaged using a microscope. Subsequently,
   cellular sampling was performed via patch clamp technique utilizing
   micro-pipettes, followed by SCMS analysis. Finally, the fluorescent
   intensity was calculated and paired with metabolomic features in single
   cells, yielding a pairwise dataset of metabolome and oxidative levels
   (Fig. [76]1a).

Fig. 1. The cross-modality analysis platform integrating single-cell
metabolome and cellular phenotype.

   [77]Fig. 1
   [78]Open in a new tab

   a A workflow and experimental setup of the cross-modality analysis.
   Cells were first labeled with DCFDA and photographed with a fluorescent
   microscope, followed by sampling and single-cell MS analysis. The
   oxidative levels were reflected by DCFDA fluorescent intensity and the
   metabolic information was acquired by SCMS. b Heatmap showing relative
   abundance of representative metabolites corresponding to single-cell
   DCFDA intensity. DCFDA intensity was indicated by color: dark green,
   relative high intensity; light green, relative low intensity.
   Metabolite abundance was represented by color: red, relative high
   abundance; blue, relative low abundance. c–e PCA score plot (c), UMAP
   analysis (d), and tSNE analysis (e) showing no significant difference
   in metabolome of DCFDA incubated (n = 257, pink) and non-incubated
   cells (n = 325, blue). f Correlation heatmap illustrating Pearson’s
   correlation coefficient (r) between metabolites in non-incubated (left,
   n = 325) and DCFDA incubated (right, n = 257) cells. Glutamate was
   correlated with GABA and glutamine (inset). Two-sided Pearson’s
   correlation analysis was performed. P values were not adjusted. For
   (c–e), Source data are provided as Source Data files.

   In m/z ranging from 67 to 1000, a total of more than 500 ion signals,
   with signal-to-noise ratio greater than 3 (refs. ^[79]27,[80]28) and
   detected at a frequency greater than 20% in all single cells were used
   in subsequent analysis. Among these signals, 162 matched annotations in
   the HMDB database through MS spectra. The annotated metabolites
   underwent further scrutiny via MS/MS, resulting in confirmation of 83
   metabolites that were subsequently employed for pathway enrichment
   analysis (Supplementary Data File [81]1). Fluorescent intensities of
   sampled cells were calculated and matched to every single cell. The
   heatmap illustrated the metabolite abundance and corresponding DCFDA
   intensity in single cells, showing the gradual alteration of metabolism
   as DCFDA intensity differs (Fig. [82]1b). A variety of metabolites,
   exhibiting diverse changes in response to DCFDA intensification,
   including glutathione (GSH), phosphocreatine, hypotaurine and adenosine
   triphosphate (ATP), were unearthed, suggesting a profound impact of OS
   on cellular metabolism.

   To rule out the possibility that DCFDA incubation may disturb cellular
   metabolism, we analyzed the metabolome of the cells treated with or
   without DCFDA by using analysis methods including Principal component
   analysis (PCA), uniform manifold approximation and projection (UMAP),
   and t-distributed stochastic neighbor embedding (t-SNE). As a result,
   no significant difference in cellular metabolome between DCFDA
   incubated and non-incubated cells was observed (Fig. [83]1c-e).
   Consistently, the abundance of 31 common metabolites including GSH,
   ATP, adenosine diphosphate (ADP), glutamate, glutamine, creatine,
   oxidized glutathione (GSSG), phosphorylcholine, and phosphoserine were
   almost identical in DCFDA incubated and non-incubated cells
   (Supplementary Fig. [84]1d). We further performed a correlation
   analysis of metabolites in both DCFDA incubated and non-incubated
   cells, revealing the similarity in the correlations among metabolites.
   As shown in the heatmap, the correlation coefficients between
   metabolites remain virtually indistinguishable in both DCFDA incubated
   and non-incubated cells (Fig. [85]1f), suggesting that the metabolic
   landscape was not significantly altered by DCFDA incubation. For
   instance, the correlation and conversion of metabolites in the
   classical glutamine-glutamate-GABA metabolic pathway remained
   unaffected by the incubation of DCFDA, as evidenced by a robust
   association between glutamate and GABA, as well as between glutamine
   and glutamate, observed in both DCFDA-incubated and non-incubated
   cells. The difference of correlation coefficients (r) between
   non-incubated and DCFDA incubated cells was calculated and visually
   represented using a heatmap (Supplementary Fig. [86]1e), illustrating
   minimal divergence in the correlations of metabolites within the
   non-incubated and DCFDA incubated cells. Altogether, the cross-modality
   analysis platform successfully integrated single-cell metabolome and
   cellular phenotype with high stability and reliability, while
   preserving cellular viability and metabolism.

The SCLIMS reveals correlation between cellular metabolism and oxidative
levels

   The Multi-modal properties of SCLIMS facilitate its capacity to explore
   the link between cellular metabolism and their OS levels, by means of
   conducting correlation analysis between single-cell metabolomics data
   and DCFDA intensity in each cell (Fig. [87]2a). The metabolomics data
   was confirmed with no obvious batch effect by PCA and HCA
   (Supplementary Fig. [88]2a, b). For each cell, the DCFDA intensity and
   the abundance of metabolites were concomitantly collected. Then the
   correlation between DCFDA intensity and metabolite abundance was
   analyzed, with the calculation of Pearson’s r coefficient to assess the
   strength of the correlation. For example, the intracellular GSH
   abundance exhibited an inverse correlation with the DCFDA intensity in
   the cells (Fig. [89]2a). The consistency with previous studies, which
   highlight the pivotal role of GSH as a crucial metabolite in
   maintaining redox balance^[90]29,[91]30, further confirms our method’s
   reliability.

Fig. 2. The correlation analysis of intracellular metabolome and cellular OS
level.

   [92]Fig. 2
   [93]Open in a new tab

   a The diagram illustrates the pairing of metabolic data and OS levels,
   as well as the subsequent correlation analysis. The DCFDA intensity
   indicated the specific level of OS in each cell, while single-cell MS
   was employed to determine the abundance of metabolites. Here, the
   metabolite GSH in Cells 1–4 was used as an illustrative example. Data
   from Cell 1–4 was normalized to Cell 1. The vectors for x and y were
   constructed based on single-cell OS levels and metabolite abundances,
   followed by calculation of Pearson’s coefficient (r) and significance
   (P). In correlation analysis, the data was z scored. Scale bar, 15 μm.
   b Metabolites that were inversely and positively correlated with the
   level of OS. Metabolites are shown as dots by color: red, positively
   correlated (P < 0.05); blue, inversely correlated (P < 0.05). The
   representative metabolites are annotated adjacent to the corresponding
   dots. P values were not adjusted. Two-sided Pearson’s correlation
   analysis was performed. c–h The correlation analysis between the scaled
   levels of OS and the scaled intensities of representative metabolites
   including ATP (c), PCr (d), UTP (e), GTP (f), Hypt (g), and energy
   charge (h). The data were standardized by calculating z score. P values
   were not adjusted. Two-sided Pearson’s correlation analysis was
   performed. i The metabolic pathways enriched by the MSEA analysis of
   metabolites exhibiting a reverse correlation with OS levels. Only
   metabolites with MS/MS confirmation were included in the analysis. Only
   pathways with P < 0.05 in the MSEA analysis were included. P values
   were not adjusted. Related metabolic processes are annotated on the
   left. Dot size represents enrichment ratio, and dot color represents
   significance of the enrichment (P value). Yellow, relatively high
   significance; Blue, relatively low significance. n = 190 cells in (a)
   and (c–h). For (a–h), Source data are provided as Source Data files.
   ATP adenosine triphosphate, PCr phosphocreatine, UTP uridine
   triphosphate, GTP guanosine triphosphate, Hypt hypotaurine, GSH
   glutathione, CTP cytosine triphosphate, O-PE O-phosphoethanolamine,
   G-3-P glycerol 3-phosphate, AMP adenosine monophosphate.

   Among the whole metabolome, we observed a total of 254 metabolites
   significantly correlated with OS level (P < 0.05). The majority (61.4%)
   of the metabolites correlated with OS level (P < 0.05) exhibited an
   inverse correlation, which is nearly double the number of metabolites
   that showed a positive correlation with cellular OS. This suggests that
   the downregulation of multiple metabolites may serve as a crucial
   hallmark of OS (Fig. [94]2b). For instance, the abundance of key
   metabolites associated with energy metabolism, such as ATP
   (Fig. [95]2c) and phosphocreatine (Fig. [96]2d), exhibited a
   progressive decline as OS level increased. The linear downregulation
   was also observed in other high energy compounds such as uridine
   triphosphate (UTP) and guanosine triphosphate (GTP) (Fig. [97]2e, f).
   The abundance of hypotaurine (Hypt), which generates NADH as a
   by-product upon conversion to taurine^[98]31 and plays crucial roles in
   redox homeostasis and DNA protection^[99]32, gradually declines during
   OS (Fig. [100]2g). The energy metabolism is composed of numerous
   metabolites and is detrimentally affected by OS^[101]33. Single-cell
   energy charge, a metabolic parameter used to assess cellular energy
   supply by measuring ATP, ADP and adenosine monophosphate (AMP) levels,
   reflects the energy homeostasis within a cell^[102]34. The inverse
   correlation observed between the single-cell energy charge and the
   cellular OS level (Fig. [103]2h) indeed suggests an association between
   OS and dysfunctionality in energy metabolism.

   By performing the metabolite set enrichment analysis (MSEA)^[104]35
   using the metabolites inversely correlated with OS level, we found a
   variety of metabolic pathways related to different biological processes
   were disturbed in OS (Fig. [105]2i), such as the mitochondrial and
   energy metabolism including “Mitochondrial Electron Transport Chain”
   and “Citric Acid Cycle”, and the redox metabolism including
   “Glutathione Metabolism Pathway”. These findings were in consistence
   with previous reports observing the disruption of glutathione
   metabolism^[106]36,[107]37 and energy metabolism^[108]33 during OS.
   Lipid metabolism has been reported to be altered in OS^[109]38,[110]39.
   This was also verified by our finding of the downregulated pathways
   such as “Phosphatidylethanolamine Biosynthesis”, “Phosphatidylcholine
   Biosynthesis”, and “Sphingolipid metabolism”. The downregulation of
   purine and pyrimidine metabolic pathways revealed by the SCLIMS was
   also supported by previous studies reporting a depletion of purine and
   pyrimidine nucleotide during OS^[111]40. The SCLIMS also revealed an
   interference of vitamin metabolic pathways including “Thiamine
   Metabolism Pathway” and “Riboflavin Metabolism Pathway” during OS,
   which is in line with previous studies illustrating the crucial role of
   thiamine and riboflavin in regulating OS^[112]41,[113]42. The
   consistency of these discoveries by SCLIMS and previous reports again
   suggests the reliability of the technique.

   The SCLIMS additionally detected numerous previously undisclosed
   metabolic alterations during OS. For instance, pathways related to
   amino acids metabolism such as “glutamate metabolism”, “alanine
   metabolism”, “arginine and proline metabolism” were all disturbed in
   the cells with OS. Moreover, pathways related to sugar and derivatives
   metabolism such as “Amino sugar metabolism”, “Fructose and mannose
   degradation”, “Lactose synthesis”, “Lactose degradation”, and
   “Nucleotide sugars metabolism” were also discovered to be downregulated
   in OS. Other pathways like “threonine and 2-oxobutanoate degradation”
   were also found to be disturbed in OS (Fig. [114]2i). The discoveries
   made by SCLIMS have shed light on the pivotal role of metabolism in
   cellular oxidative stress, unveiling a diverse array of metabolic
   pathways that undergo significant alterations (P < 0.05). These
   findings indicate that the mechanisms governing protein synthesis and
   degradation, glycosylation reactions, and energy metabolism may undergo
   drastic transformations under OS.

   The SCLIMS technique was also employed for multi-modal analysis in
   MEFs. Firstly, the metabolomics data were validated for the absence of
   any noticeable batch effect through PCA and HCA (Supplementary
   Fig. [115]2c, d). In MEFs, the key metabolites found to be inversely
   correlated with OS levels in HEK cells, including GSH, O-PE, CTP, ATP,
   UTP, Hypt, GTP, and PCr, were similarly observed to exhibit reverse
   correlation with single-cell OS levels (Supplementary Fig. [116]3a).
   Additionally, MSEA analysis was conducted using these metabolites
   exhibiting inverse correlation with OS levels to unravel the metabolic
   pathways involved in oxidative stress of MEFs. Similar to the results
   of HEK cells, a variety of metabolic pathways were enriched
   (Supplementary Fig. [117]3b). For example, “Pyruvate metabolism”,
   “Citric acid cycle”, “Beta oxidation of very long chain fatty acids”,
   and “Mitochondrial electron transport chain”, in mitochondria and
   energy metabolism were also involved in the OS of MEFs. Lipid
   metabolism, including “Phosphatidylcholine Biosynthesis”,
   “Phosphatidylethanolamine Biosynthesis”, “Sphingolipid metabolism”,
   were altered during OS in MEFs as well. The downregulation of purine
   and pyrimidine metabolic pathways were also discovered in MEFs under
   OS. Similar to the interruption of vitamin metabolism in HEK cells
   under OS, “Pantothenate and CoA Biosynthesis”, “Thiamine Metabolism
   Pathway” were also discovered to be downregulated in OS-stressed MEFs.
   The outcomes in other pathways also exhibited a resemblance to those
   observed in HEK cells. These results confirmed the robustness of the
   multi-modal analysis of SCLIMS across diverse cellular phenotypes.

   Taken together, revealing a profound interplay, the SCLIMS-based
   cross-modality analysis has unveiled strong connections between
   single-cell metabolism and OS, thereby highlighting the paramount
   importance of integrating metabolome and cellular phenotype for
   comprehensive insights across different cell types.

Cell types identified by the SCLIMS exhibit divergent OS levels

   Cluster analysis was conducted using k-medoids
   algorithm^[118]43,[119]44 based on the single-cell metabolome acquired
   from the SCLIMS. As a result, cells were clustered into six subtypes
   (C1–C6) with distinct metabolic features according to their metabolome
   profiles, highlighting the diversity of metabolic characteristics
   within these cells (Fig. [120]3a). Notably, the six subtypes of cells
   displayed distinct levels of OS as indicated by varying DCFDA
   intensities (Fig. [121]3b). We then performed a pseudotime analysis
   using single-cell metabolomics data and generated a trajectory of the
   six subtypes. The trajectory originated from subtype C1, which
   exhibited the lowest level of OS, and gradually progressed towards
   other subtypes. Notably, three distinct branches were observed that led
   to C2/3 (representing a low oxidative level), C4 (representing a medium
   oxidative level), and C5/6 (representing a high oxidative level),
   respectively (Fig. [122]3c). These findings suggest a metabolism-guided
   stepwise progression of cellular OS.

Fig. 3. Subtype-specific metabolic signatures and dynamic change of
metabolome revealed by cross-modality analysis.

   [123]Fig. 3
   [124]Open in a new tab

   a UMAP visualization of six cellular subtypes based on the single-cell
   metabolome. Each point represents a single cell, color represents
   different subtypes. n = 55, 34, 14, 41, 29, and 17 for C1, C2, C3, C4,
   C5, and C6 respectively. b Average OS levels (indicated by the
   normalized DCFDA fluorescence intensity) of the six metabolic
   subpopulations. n = 55, 34, 14, 41, 29, and 17 for C1, C2, C3, C4, C5,
   and C6, respectively. F(5, 184) = 5.880, P = 4.57e−5 in One-way ANOVA
   (labeled in the plot). Data was normalized to values of the respective
   control (Cluster C1). Data is presented as mean ± s.e.m. Color
   represents different subtypes and each dot represents a cell. c
   Single-cell trajectory of pseudotime analysis showing temporal
   progression of cell subtypes originating from C1 and gradually
   transitioning towards C2/3, C4, and C5/6. Each subtype is represented
   by a distinct color. For (a–c), Blue: C1; yellow: C2; purple: C3;
   green: C4; brown: C5; red: C6. d Heatmap of characteristic metabolites
   and their corresponding relative intensities in each subtype. Color
   indicates z scores of metabolite abundance. Red: relatively high
   abundance; blue: relatively low abundance. e Representative metabolic
   pathways significantly (P < 0.05) enriched in Cluster C1 and C6 through
   the MSEA analysis. Only metabolites with MS/MS confirmation were
   included in the analysis. Dot size represents enrichment ratio, while
   dot color indicates significance (−log[10] P value) of the enrichment.
   Red: relatively high significance; blue: relatively low significance. P
   values were not adjusted. For (a–c), Source data are provided as Source
   Data files. Data was collected from at least three biological
   replicates.

   We subsequently conducted a more detailed analysis of the six metabolic
   subtypes by utilizing a Wilcox rank sum test to compare metabolite
   abundance between each cluster and other clusters. Each subtype was
   characterized with specific metabolic markers as shown in the heatmap
   (Fig. [125]3d). For instance, cells of Cluster C1 exhibited enrichment
   of energy rich phosphate compounds such as ATP, GTP, UTP,
   phosphocreatine along with antioxidants including GSH; cells of Cluster
   C2 were enriched with dimethyldithiophosphate and pyrophosphate; cells
   of Cluster C3 showed an abundance of catecholamine derivatives such as
   homovanillic acid and sugar derivatives such as deoxyribose
   5-phosphate; cells of Cluster C4 were enriched with intermediates of
   nucleotide metabolism such as guanine monophosphate (GMP), along with
   acetylated compounds such as N-acetylneuraminic acid and
   UDP-N-acetylglucosamine; cells of Cluster C5 were enriched with amino
   acids such as proline, threonine, asparagine, alanine and taurine,
   along with glucose and ribitol; Nucleotide monophosphates, including
   AMP, cytidine monophosphate (CMP), and uridine monophosphate (UMP), as
   well as intermediates in the glycolysis pathway (i.e. glyceraldehyde
   3-phosphate and fructose 1,6-bisphosphate), were found to be enriched
   in cells belonging to Cluster C6.

   Next, we conducted a comparison of the metabolome at the single-cell
   level between Cluster C6 (with the most oxidative level) and Cluster C1
   (with the least oxidative level). A total of 148 metabolites were
   downregulated and only 64 metabolites were upregulated from C1 cells to
   C6 cells (Supplementary Fig. [126]4a). To investigate the alteration of
   the metabolic process in cellular OS, we performed the MSEA analysis of
   the characteristic metabolites enriched in C1 and C6 cells
   (Fig. [127]3e), and other clusters of cells (Supplementary
   Fig. [128]5). Specifically, pathways related to lipid metabolism such
   as “Phosphatidylethanolamine Biosynthesis”, “Phosphatidylcholine
   Biosynthesis”, “Sphingolipid Metabolism”, and “Phosphatidylinositol
   Phosphate Metabolism” were less enriched or completely depleted in C6
   cells, consistent with previous reports^[129]38,[130]39. The pathways
   associated with mitochondrial function and redox balance, such as
   “Citric Acid Cycle”, “Nicotinate and nicotinamide metabolism” and
   “Glutathione metabolism” were also found to be depleted or less
   enriched in cells of C6, thus confirming previous reports on the
   disruption of glutathione metabolism^[131]36,[132]37 and energy
   metabolism^[133]33 under OS. The pathways related to metabolism of
   sulfinic acids and organosulfonic acids, such as hypotaurine and
   taurine metabolism, were depleted in C6 cells. This was supported by
   studies illustrating the antioxidant effect of these
   metabolites^[134]45,[135]46. In addition, the glycine and serine
   metabolism exhibited greater enrichment in C6 cells, which may be a
   compensatory response of cells under OS as glycine and serine
   metabolism are reported to attenuate OS in C. elegans^[136]47. The
   nucleotide metabolism pathways, such as “Purine metabolism” and
   “Pyrimidine metabolism”, were depleted in C6 cells. The metabolism of
   vitamins, including “Pantothenate and CoA metabolism”, “Thiamine
   metabolism”, and “Folate metabolism”, exhibited reduced enrichment or
   depletion in C6 cells. The alteration of nucleotide and vitamin
   metabolism in OS were well described by previous
   researches^[137]40–[138]42. Other pathways such as “Plasmalogen
   Synthesis” was enriched exclusively in the cells of C6 (Fig. [139]3e).
   Plasmalogen synthesis is reported to enhance the resistance to OS in E.
   Coli.^[140]48, which could be served as a protective response in cells
   under OS. Together, the convergence of these consistent findings
   between the SCLIMS and the existing literature further substantiates
   the dependability and steadfastness of this technique.

   Moreover, the SCLIMS has also unearthed previously unreported perturbed
   pathways during OS. For instance, several metabolic pathways, such as
   “Fructose and mannose degradation”, “Amino sugar metabolism”,
   “Nucleotide sugars metabolism” along with amino acids metabolism
   including “alanine metabolism”, “glutamate metabolism”, and “arginine
   and proline metabolism”, were identified to be implicated in OS, which
   has not been elucidated in previous studies. Another finding gleaned
   from the SCLIMS is the transition from the “Malate-Aspartate Shuttle”
   pathway to the “Glycerol Phosphate Shuttle” pathway, observed between
   C1 and C6, implying a reduced efficiency in ATP generation during OS.

   The metabolic heterogeneity of MEFs under OS was meticulously
   investigated using SCLIMS, revealing a parallel to HEK293T cells where
   six distinct subtypes (M1–M6) with specific metabolic characteristics
   were identified (Supplementary Fig. [141]6a). These subtypes in MEFs
   also exhibited varying levels of OS (Supplementary Fig. [142]6b).
   Furthermore, pseudotime analysis utilizing single-cell metabolic
   features unveiled a trajectory of the six subtypes originating from M1
   (with the lowest level of OS) and progressing towards M3 (with moderate
   OS), followed by M5 and finally converging at M4/6 (with higher levels
   of OS) (Supplementary Fig. [143]6c). Notably, a metabolism-guided
   progression of cellular oxidative stress was similarly observed in
   MEFs.

   The metabolite abundance was also analyzed between each cluster and
   other clusters, revealing subtype-specific metabolic markers
   (Supplementary Fig. [144]6d). In comparison to HEK293T cells, MEFs
   exhibited strikingly similar metabolic characteristics. Similar to C1
   in HEK293T cells, M1 in MEFs shared crucial metabolites such as Hypt,
   PCr, O-PE, ATP, GTP, UTP, and GSH. Both C1 and M1 displayed the lowest
   level of OS. A comparison of the metabolites in M6 (with the highest OS
   level) and M1 (with the lowest OS level) yielded a result akin to that
   observed in HEK293T cells. A total of 195 metabolites were
   downregulated, while 134 metabolites were upregulated from M1 to M6
   (Supplementary Fig. [145]4b). The metabolic pathways enriched in M1 and
   M6 were also analyzed using MSEA and compared (Supplementary
   Fig. [146]6e). Consistent with the findings in C1 and C6 of HEK293T
   cells, common lipid metabolic pathways such as
   “Phosphatidylethanolamine Biosynthesis”, “Phosphatidylcholine
   Biosynthesis”, and “Sphingolipid Metabolism” exhibited depletion in M6
   MEFs characterized by the highest level of OS. Additionally, key
   pathways involved in mitochondrial function and redox balance, namely
   “Citric Acid Cycle” and “Glutathione metabolism”, were also found to be
   depleted in M6 MEFs. Furthermore, several other pathways including
   “taurine and hypotaurine metabolism”, “alanine metabolism”, “glutamate
   metabolism”, and “arginine and proline metabolism” were identified as
   being implicated specifically in the OS subtype of MEFs at stage M6.
   Similarly, nucleotide metabolism pathways including “Purine metabolism”
   and “Pyrimidine metabolism” showed depletion or less significant
   enrichment in the same group of cells. Moreover, there was a depletion
   observed in the vitamin-related metabolic processes such as
   “Pantothenate and CoA metabolism”, “Thiamine metabolism”, and “Folate
   metabolism” within the context of OS progression among these specific
   type of cells. The obtained results exhibited consistency with those
   observed in HEK293T cells, as well as previous studies that have
   dissected cellular OS metabolism mentioned above. Consequently, SCLIMS
   has been demonstrated to possess robustness in analyzing metabolic
   alterations at the single-cell level.

   Taken collectively, the SCLIMS not only validate metabolic changes
   observed in previous studies but also unveil metabolic modifications
   associated with OS and metabolic heterogeneity within individual cells.
   These findings emphasize the reliability and universality of this
   technique, establishing it as a powerful tool for comprehensive
   analysis at the single-cell level.

The SCLIMS reveals capability of the single-cell metabolome in predicting
cellular OS status

   Although the metabolic alterations and heterogeneity of cells under OS
   were explored in detail, the strength of the link between metabolism
   and cellular phenotype is not fully established. Machine learning is
   utilized to ascertain whether the metabolic profile of individual cells
   can accurately predict heterogeneous subtypes with distinct OS levels.
   We employed discriminant analysis algorithms to train classification
   models using single-cell metabolic features aiming to distinguish the 6
   subtypes with specific OS levels identified in the above-mentioned
   clustering analysis (Fig. [147]4a). The data was randomly divided into
   independent training and testing datasets, with 2/3 of the original
   data composing the former and 1/3 of it composing the latter. Both
   datasets included all m/z signals meeting our criteria (S/N > 3 and
   detected in greater than 20% single cells) and no variable selection is
   performed before the model was trained. The model was trained on the
   training dataset to acquire the ability to classify metabolic subtypes
   based on full features of the single-cell metabolome, and its
   performance was assessed using the testing dataset. The performance of
   the model was evaluated using receiver operating characteristic (ROC)
   curve and confusion matrix. In multiclassification, the ROC curve had
   an average area under curve (AUC) of 0.98 (Fig. [148]4b). In cluster
   prediction, the model achieved an accuracy ranging from 77.8% to 100%
   (Fig. [149]4c). These results demonstrate that the single-cell
   metabolic profiles can directly predict metabolic subtypes.

Fig. 4. Machine learning-guided prediction of OS levels based on single-cell
metabolome.

   [150]Fig. 4
   [151]Open in a new tab

   a Flowchart of classification analysis with machine learning. The
   training and testing dataset were randomly assigned according to a
   ratio of 2:1. The model was trained and built with the training dataset
   with 5-fold cross validation. The testing dataset was held out and used
   for the evaluation of model accuracy independently. The performance of
   the model was evaluated with ROC curve and confusion matrix. b ROC
   curve of model testing. AUC for each cluster was determined separately
   by the classification model. Higher AUC value indicates a better
   performance of the model in predicting the clusters. The average AUC
   represented an overall performance of the model. The color represents
   classification of a certain subtype. Blue: C1; yellow: C2; purple: C3;
   brown: C4; green: C5; red: C6. c Confusion matrix of model testing,
   illustrating the distribution of errors in multi-class prediction. The
   average accuracy was used to evaluate the overall performance of the
   model. The color depth indicates the proportion of correct (blue) and
   incorrect (red) predictions, as displayed in the bar chart. d Flowchart
   of building a regression model with machine learning. The training and
   testing datasets were randomly assigned according to a ratio of 2:1.
   The model was trained and built with the training dataset with 5-fold
   cross validation. The testing dataset was held out and used for the
   evaluation of the model independently. e The correlation of real values
   and values predicted by the regression model (n = 61). Dash line
   represents the perfect fit (predicted values = real values). The model
   performance was evaluated by the correlation coefficient (r) and P
   value (P). Two-sided Pearson’s correlation analysis was performed.
   P = 1.45e−20. P value was not adjusted. For b and e, Source data are
   provided as Source Data files. AUC area under curve.

   The classification model is capable of predicting metabolic subtypes
   rather than precise levels of OS. We further trained a regression
   model^[152]49,[153]50 based on neuronal network algorithms
   (Fig. [154]4d), which may enable direct prediction of single-cell OS
   levels utilizing metabolic features. Similarly, the data was
   independently and randomly divided into training and testing datasets,
   with 2/3 of the original data comprising the former and 1/3 comprising
   the latter. All m/z signals meeting our criteria (S/N > 3 and detected
   in greater than 20% single cells) were included in building the
   regression model without variable selection. The model was trained on
   the training dataset to predict single-cell OS levels, as manifested by
   DCFDA intensity, based on single-cell metabolic features. The
   performance of the model was evaluated with the testing dataset, and
   the predicted values were plotted against real values (Fig. [155]4e).
   The correlation coefficient (r) was 0.88, indicating a good predictive
   power of the model. Furthermore, the metabolic profile of single-cells
   was further validated in MEFs, demonstrating a high predictive power
   with an AUC of 0.99 and an average accuracy of 88.42% in the
   classification model (Supplementary Fig. [156]7a, b). Additionally,
   there was a significant correlation between predicted and true OS
   levels (r = 0.6, P < 0.0001) as shown in the correlation analysis
   (Supplementary Fig. [157]7c).

   Thus, through the implementation of SCLIMS technique and its
   multi-modal integration, we have successfully demonstrated a link
   between single-cell metabolome and cellular OS states in two different
   cell types for the very first time, suggesting the potential role of
   metabolome in determining cellular phenotype.

The SCLIMS unveils the causal relationship between metabolic heterogeneity
and OS status

   Although there appears to be a strong correlation between the
   heterogeneity of OS levels and metabolic heterogeneity, the causal
   relationship between them remains unclear. We wonder if the baseline
   metabolic profile determines the phenotypic or metabolic heterogeneity
   of cells after they are induced to OS. We define the cells prior to OS
   induction with hydrogen peroxide as “initial cells”, representing the
   baseline state of both metabolism and phenotype. With SCLIMS, we
   analyzed the heterogeneity of the OS levels and the single-cell
   metabolome in the initial cells untreated with H[2]O[2]. We
   surprisingly found that while there was minimal heterogeneity in
   cellular OS levels among the initial cells (Fig. [158]5a and
   Supplementary Table [159]1), there existed heterogeneity in their
   single-cell metabolome (Fig. [160]5b, c). Specifically, the k-medoids
   clustering analysis revealed two major metabolic subtypes (Cluster-I
   and Cluster-II) in the initial cells based on their metabolomics
   features (Fig. [161]5b). Through machine learning, the single-cell
   metabolome can directly predict the subtype of a cell, providing
   further evidence for the robust heterogeneity of metabolism in the
   initial cells (Supplementary Fig. [162]8a, b). To investigate metabolic
   properties of the two subtypes of initial cells, we performed a
   differential analysis by comparing metabolite abundance between
   Cluster-I and Cluster-II cells. A series of characteristic metabolites,
   such as hypotaurine, GSH, ATP, UTP, O-phosphoethanolamine, and GSSG,
   distinguishing the two subtypes were discovered (Fig. [163]5c).

Fig. 5. Metabolic heterogeneity in initial cells.

   [164]Fig. 5
   [165]Open in a new tab

   a Distribution and box plot (inset) of DCFDA intensity in initial cells
   and cells under OS. Data in distribution plot was z scored. Variance
   was indicated by IQR and MAD values in Supplementary Table [166]1.
   W = 58627, P < 2.2e-16 in unpaired two-tailed Wilcox rank sum test. The
   data presented in the inset was normalized to the values of initial
   cells. Box plots extend from 25th to 75th percentiles; central lines
   represent medians; whiskers extend over 1.5 times the interquartile
   range (IQR, the distance from 25th to 75th percentile); dots represent
   outliers. For single-cell DCFDA intensity, a total of 1740 initial
   cells (blue) and 960 oxidative stressed cells (gray) from 3 independent
   experiments were analyzed. b UMAP visualization of metabolic subtypes
   in initial cells. Green: cells of Cluster-I (n = 43). Purple: cells of
   Cluster-II (n = 183). c Heatmap of potential metabolic markers in
   Cluster-I and Cluster-II. The data was z score scaled. The color
   represents relative abundance of metabolites. Yellow: relatively high
   abundance; green: relatively low abundance. Each row represents a
   metabolite and each column represents a cell. The representative
   metabolites are labeled on the left and the clusters are labeled on the
   top. ATP: adenosine triphosphate; GSH: glutathione; GSSG: oxidized
   glutathione; Hypt: hypotaurine; O-PE: O-phosphoethanolamine; UTP:
   uridine triphosphate. d A heatmap illustrating the metabolic similarity
   between subtypes of initial cells (Cluster-I/II) and subtypes of
   oxidative stressed cells (C1 to C6). Each row represents an initial
   cell (subtypes are labeled on the left) and each column represents an
   oxidative stressed cell (subtypes are labeled on the top). The heatmap
   was plotted with similarity (1/distance) and the data was z score
   scaled. The distance between cells was calculated based on single-cell
   metabolome. The color represented the relative similarity: red,
   relative high similarity; blue, relative low similarity. Cells of
   Cluster-I is more similar to the cells with lower OS levels (cells of
   C1 and part of C2). e Enrichment of GSH in cells of Cluster-I
   visualized on the UMAP plot. Each dot represents a cell, the color of
   the dots represents the relative GSH abundance. Red: relatively high
   abundance; blue: relatively low abundance. Data was z score scaled. f
   Quantification of GSH abundance in Cluster-I (n = 43, green) and
   Cluster-II (n = 183, purple). Data was normalized to values in
   Cluster-I. Data is represented as mean ± s.e.m. Data were collected
   from at least three biological replicates. W = 7843, P = 5.85e−43 in
   unpaired two-tailed Wilcox rank sum test. g The correlation of every
   two metabolites were calculated (reflected as Pearson’s r) and the
   network was constructed based on the correlation data for initial cells
   belonging to Cluster-I (left) and Cluster-II (right) subtypes. In the
   network, each node represents a metabolite while an edge connecting two
   nodes indicates their correlation. The big green dot denotes GSH and
   the small blue dots denote GSH correlated metabolites. Blank-colored
   dots indicate metabolites without any correlation to GSH. h Rewiring
   score calculated with DyNet algorithm. Metabolites with higher scores
   were more rewired in topology in the correlation network of initial
   cells. For (a, b, e, f), Source data are provided as Source Data files.
   GSH: glutathione.

   Interestingly, despite no initial cells being subjected to OS
   stimulation, the metabolic disparity between the two cell subtypes in
   these initial cells already demonstrated a consistent pattern with the
   metabolic variation observed in cells with varying levels of OS.
   Specifically, the metabolomic similarity between single cells was
   determined by calculating the statistical distance of their metabolomic
   data. The initial cells (from Cluster-I and Cluster-II) were paired
   with each cell under OS (from C1-C6), and subsequently, the distance
   was computed based on the metabolite abundance in the two cells. A
   larger value of statistical distance indicates a lower degree of
   similarity^[167]6. The heatmap visualization (Fig. [168]5d) was
   generated based on the reciprocal value of statistical distance, which
   represents the level of similarity among single cells. From a
   metabolomic perspective, cells in Cluster-II exhibit greater similarity
   to cells with higher levels of OS (i.e., C2-4), while cells in
   Cluster-I are more similar to cells with lower levels of OS (i.e., C1)
   (Fig. [169]5d). The abundance of GSH, ATP and hypotaurine were 82%,
   42%, and 87% lower in cells of Cluster-II (Fig. [170]5e, f and
   Supplementary Fig. [171]8c), which was in line with the above-mentioned
   downregulation of GSH, ATP, and hypotaurine discovered in cells with
   high OS levels (Supplementary Fig. [172]4). Metabolic pathways
   including “Glutathione metabolism”, “Phosphatidylcholine Biosynthesis”,
   “Phosphatidylethanolamine Biosynthesis”, “Sphingolipid Metabolism”,
   “Pyrimidine Metabolism”, “Purine Metabolism”, “Thiamine Metabolism”,
   “Pantothenate and CoA Biosynthesis”, and “Taurine and Hypotaurine
   Metabolism” were enriched with marker metabolites in cells of Cluster-I
   (Supplementary Fig. [173]8d). Notably, these pathways were all
   downregulated in cells with high OS levels (i.e., cells in C6) when
   compared to cells with low OS levels (i.e., cells in C1)
   (Fig. [174]3e).

   A remarkable demonstration of the metabolic heterogeneity exhibited by
   initial cells is exemplified by GSH. The UMAP map vividly depicted the
   remarkable disparity in GSH abundance among individual cells of
   Cluster-I and Cluster-II (Fig. [175]5e). We further constructed
   metabolic networks based by leveraging the correlation among
   metabolites. The correlations between metabolites exhibited alterations
   in cells belonging to Cluster-I and Cluster-II, indicating a
   subtype-specific metabolic profile. For instance, the intracellular
   level of GSH was strongly correlated with 15 other metabolites in
   Cluster-I, whereas only 5 metabolites displayed a correlation (r > 0.8,
   P < 0.05) with GSH in Cluster-II (Fig. [176]5g), indicating a notable
   divergence in the GSH-centered network between the two cell subtypes.
   The DyNet algorithm was employed to calculate the rewiring
   score^[177]51 for each metabolite in Cluster-I and Cluster-II networks,
   which quantifies modifications in a specific metabolite’s correlation
   with other metabolites between the two clusters. Notably, the
   identification of GSH as the most rewired metabolite (Fig. [178]5h)
   highlights its pivotal role in connecting with other metabolites within
   the metabolome.

   The metabolic heterogeneity in initial MEFs was also investigated.
   Notably, MEFs were classified into two distinct subtypes (Cluster-1 and
   Cluster-2) exhibiting metabolic characteristics akin to the two
   clusters observed in initial HEK293T cells (Cluster-I and Cluster-II)
   (Supplementary Fig. [179]9a, b). The detailed metabolic properties of
   Cluster-1 and Cluster-2 initial MEFs were further studied. Several key
   metabolites, including Hypt, GSH, O-PE, and ATP exhibited enrichment in
   Cluster-1 within the initial MEFs, aligning with the characteristic
   metabolites in Cluster-I within the initial HEK293T cells. The Pathways
   including “Glutathione metabolism”, “Phosphatidylcholine Biosynthesis”,
   “Phosphatidylethanolamine Biosynthesis”, “Sphingolipid Metabolism”,
   “Pyrimidine Metabolism”, “Purine Metabolism”, “Thiamine Metabolism”,
   “Pantothenate and CoA Biosynthesis”, and “Taurine and Hypotaurine
   Metabolism” were enriched with marker metabolites in Cluster-1
   cells (Supplementary Fig. [180]10c), which was consistent with the
   pathways enriched in Cluster-I in HEK293T cells (Supplementary
   Fig. [181]8d). These pathways were also depleted or less enriched in
   MEFs with the highest OS levels (Cluster-M6). Machine learning analysis
   of the single-cell metabolome of initial MEFs confirmed the robust
   heterogeneity in cellular metabolism (Supplementary Fig. [182]10a, b).
   We next calculated the metabolic similarity between Cluster-1 and
   Cluster-2 initial MEFs in the absence of OS with the M1-M6 MEFs under
   OS conditions. Similarly, Cluster-1 initial MEFs exhibited a higher
   degree of similarity to MEFs with lower levels of OS, such as M1, M2,
   and part of M3 (Supplementary Fig. [183]9c). On the other hand,
   Cluster-2 initial MEFs showed a greater resemblance to MEFs with higher
   levels of OS including M4, M5, and M6 (Supplementary Fig. [184]9c).
   These findings were consistent with the results observed in HEK293T
   cells.

   Next, the role of GSH in distinguishing the metabolic landscape of
   Cluster-1 and Cluster-2 was confirmed in the initial MEFs. As shown in
   Supplementary Fig. [185]9d, e, Cluster-1 exhibited higher levels of GSH
   abundance, which could potentially serve as a discriminative factor
   between the two clusters. Metabolic networks were constructed for both
   Cluster-1 and Cluster-2 based on metabolite correlations (Supplementary
   Fig. [186]9f). Similar to the networks observed in HEK293T cells for
   Cluster-I and Cluster-II, it was found that compared to Cluster-2,
   Cluster-1 displayed a more intricate network connectivity among
   metabolites, indicating a subtype-specific metabolic profile within the
   initial MEFs. Furthermore, the key role of GSH in governing metabolism
   and correlation between metabolites, as well as its ability to
   distinguish between the two subtypes of initial MEFs, were also
   elucidated. These findings align with the observed role of GSH in
   initial HEK293T cells. Moreover, the divergence in the GSH-centered
   network between these two cell subtypes was confirmed in initial MEFs.
   By employing the DyNet algorithm to compare the two metabolic networks,
   GSH emerged as the fourth-highest rewired metabolite (Supplementary
   Fig. [187]9g). These results underscored both the generality and
   reliability of the SCLIMS technique while demonstrating its
   applicability across diverse cell types.

   In summary, utilizing the SCLIMS technique has revealed that while
   there was no obvious heterogeneity in initial cells at the OS level,
   the metabolome exhibited a high degree of heterogeneity. Interestingly,
   this metabolic diversity in initial cells mirrored that observed in
   cells with varying OS levels, suggesting that it may be the root cause
   of heterogeneous OS levels within cells.

The metabolic heterogeneity of cells determines their senescence fate under
OS

   A common fate of cellular OS is the oxidative stress-induced premature
   senescence, which is a crucial type of cellular
   senescence^[188]21–[189]24,[190]52. Subsequently, our objective was to
   identify a live-cell probe capable of distinguishing between the cells
   belonging to the two distinct metabolic subtypes, enabling us to
   further investigate their respective cellular fates in OS-induced
   senescence. The disparity in GSH levels between the two metabolic
   subtypes (Fig. [191]5e) prompted us to explore the potential of
   utilizing intracellular GSH levels as a means to distinguish between
   these two cell subtypes.

   The initial cells exhibiting high/low (ranging from 5% to 50%) GSH
   levels were selected and mapped onto the UMAP projection alongside the
   remaining initial cells, followed by calculating the proportion of
   accurate matches between cells with high/low GSH levels and
   Cluster-I/II (Fig. [192]6a). The fraction of correct matches exhibited
   a drastic decline when the selection included more than 15% of cells
   with top/bottom GSH levels (Fig. [193]6b–g), suggesting that only cells
   within the range of 5%-15% accurately represented the metabolome of
   Cluster-I and Cluster-II, respectively, and potentially reflected
   resistance to OS and sensitivity to OS. To ensure stability, we
   carefully selected cells with the highest and lowest 5% GSH levels for
   further analysis, designating them as “initial cells exhibiting a
   metabolome resembling that of OS-resistant cells (C1)” (MROR cells) and
   “initial cells exhibiting a metabolome resembling that of OS-sensitive
   cells (C6)” (MROS cells).

Fig. 6. Role of intracellular metabolic features in determining the cellular
fate in OS and OS induced senescence.

   [194]Fig. 6
   [195]Open in a new tab

   a A flowchart showing the process of the analysis. The initial cells
   were clustered based on their metabolomic features with unsupervised
   clustering, revealing Cluster-I (blue palette) and Cluster-II (red
   palette). Then cells with top (blue strip) and bottom (red strip)
   5%-50% GSH levels were projected to the UMAP scatter plot. The correct
   matches were defined as cells with top 5%-50% GSH levels to Cluster-I
   (blue dots), and cells with bottom 5–50% GSH levels to Cluster-II (red
   dots). Incorrect matches were labeled as green dots. Finally, the
   fraction of correct matches was calculated. b–f, The visualization of
   projection of cells with top/bottom 5% (b), 15% (c), 25% (d), 35% (e),
   and 50% (f) GSH levels into Cluster-I and Cluster-II. Correct matches
   were labeled as blue (cells with top 5–50% GSH levels vs Cluster-I) and
   red (cells with bottom 5%-50% GSH levels vs Cluster-II) dots. Incorrect
   matches were labeled as green dots. Cells with intermediate GSH levels
   were labeled as gray dots. g Quantification of fraction of accurate
   matches between cells with top/bottom 5%-50% GSH levels and Cluster-I
   and Cluster-II. Blue: correct matches; green: incorrect matches. h,
   Experimental setup of FACS separation of MROR and MROS cells and gating
   of the FACS: cells with top 5% (MROR, blue) and bottom 5% (MROS, red)
   GSH intensity were collected according to their fluorescent intensity.
   i, j Representative images (i) and quantification (j) of DCFDA staining
   of MROR and MROS cells before and after OS induction. F(1, 8) = 7.405,
   P = 0.0262 in two-way ANOVA. k, l Representative images (k) and
   quantification (l) of SA-β-Gal staining of MROR and MROS cells before
   and after induction of OS-induced senescence. F(1, 8) = 9.177,
   P = 0.0163 in two-way ANOVA. Scale bar, 50 μm. P values in two-way
   ANOVA with Turkey’s multiple comparisons were labeled in the plot.
   n = 3 for each group. All P values were reported as multiplicity
   adjusted P values for multiple comparisons. Data was normalized to the
   values of MROR group in control cells. Data is presented as mean ±
   s.e.m. For (j, l), Source data are provided as Source Data files. Blue:
   MROR cells; red: MROS cells. OS: oxidative stress. MROR cells: initial
   cells exhibiting a metabolome resembling that of OS-resistant cells.
   MROS cells: initial cells exhibiting a metabolome resembling that of
   OS-sensitive cells.

   To isolate the MROR and MROS cells, we employed a live-cell GSH
   fluorescent dye (monochlorobimane, mClB)^[196]53 and fluorescence
   activated cell sorting (FACS). The MROR cells (Top 5% fluorescent
   intensity) and the MROS cells (Bottom 5% fluorescent intensity) were
   collected by the FACS according to their mClB fluorescent intensity
   (Fig. [197]6h and Supplementary Figs. [198]11, [199]12). The DCFDA
   intensity (Fig. [200]6i, j) had only 3% difference between MROR and
   MROS cells. However, following OS induction, the MROS cells exhibited
   2.7-fold higher of DCFDA intensity compared to their MROR counterparts
   (Fig. [201]6i, j), indicating an elevated state of OS. Furthermore,
   SA-β-Gal staining confirmed the sensitivity of MROS cells in OS-induced
   senescence (Fig. [202]6k, l). Thus, the data from the SCLIMS suggest
   that the metabolic features of an initial cell may dictate its destiny
   under OS and OS-induced senescence.

The key metabolites identified by SCLIMS mitigate OS and cellular senescence

   The aforementioned results suggest that the downregulation of critical
   metabolites, identified by the SCLIMS, may be a contributor to OS,
   indicating the potential for reversing OS and OS-induced senescence
   through targeted metabolite supplementation. To explore this
   possibility, we opted to conduct an experiment using three key
   metabolites: hypotaurine, phosphocreatine, and O-phosphoethanolamine,
   which all showed similar characteristics such as declining in OS
   (Fig. [203]2b), distributing heterogeneously across metabolic subtypes
   (C1-6) (Fig. [204]3d), and serving as metabolic markers distinguishing
   the initial cells with different fates (Fig. [205]5c).

   Treatment of cells with the three metabolites resulted in an average of
   67% reduction in OS levels, as evidenced by a decrease in DCFDA
   intensity within the cells (Fig. [206]7a, b). OS is a common inducer of
   senescence^[207]22,[208]54,[209]55 and also plays a vital role in
   natural aging^[210]56. We therefore explored the effects of these key
   metabolites on OS-induced senescence. The SA-β-Gal staining showed that
   hypotaurine, phosphocreatine and O-phosphoethanolamine considerably
   decreased the cellular senescence (Fig. [211]7c, d and Supplementary
   Fig. [212]13a–c). The staining intensity exhibited an average of 48%
   reduction in key metabolite-treated cells compared with vehicle-treated
   OS cells. The growth arrest of cells, a marker for cellular
   senescence^[213]55,[214]57, was also rescued by these metabolites
   (Fig. [215]7e, f and Supplementary Fig. [216]13d–f). The mitochondrial
   membrane potential (MMP), which is compromised in cellular
   OS^[217]58,[218]59, was recovered by targeted supplement of key
   metabolites (Fig. [219]7g, h). Furthermore, the examination of the
   three key metabolites on MEFs under OS yielded analogous outcomes,
   thereby indicating the generalized applicability of these key
   metabolites across diverse cell types (Supplementary Fig. [220]14).
   Collectively, these findings imply that the SCLIMS-identified
   metabolites exert a protective effect by reducing the OS level and
   alleviating senescence.

Fig. 7. Effects of key metabolites on OS and induced senescence.

   [221]Fig. 7
   [222]Open in a new tab

   a, b Representative images (a) and quantification (b) of DCFDA staining
   of control cells and oxidative stressed cells with indicated
   treatments. n = 5 for each group. Scale bar, 50 μm. F(4, 20) = 5.991,
   P = 0.0024 in One-way ANOVA. P values in one-way ANOVA with multiple
   comparison were labeled in the plot. Data was normalized to values of
   control group. c, d Representative images (c) and quantification (d) of
   SA-β-Gal staining of control and oxidative stressed cells with
   indicated treatments. n = 9, 9, 3, 3, and 3 for Control, Vehicle
   treated, Hypt treated, PCr treated, and O-PE treated group
   respectively. Scale bar, 50 μm. F(4, 22) = 17.29, P = 5.50e-9 in
   One-way ANOVA. P values in one-way ANOVA with multiple comparison were
   labeled in the plot. Data were normalized to values of control group.
   e, f Representative images (e) and growth curve (f) of control and
   oxidative stressed cells with indicated treatments at indicated time
   points. Growth curve was plotted with at least 15 random fields from 3
   independent biological replicates for each group at each indicated time
   point. Scale bar, 50 μm. F(4, 520) = 41.26, P < 2.2e−16 in two-way
   ANOVA. P values in two-way ANOVA with Turkey’s HSD comparison were
   labeled in the plot. Each group was compared with vehicle-treated
   oxidative stressed cells. Data were normalized to values of control
   group at 0 h. g, h Representative images (g) and quantification (h) of
   TMRE staining of control and oxidative stressed cells with indicated
   treatments. n = 9, 9, 5, 5, and 4 for Control, Vehicle treated, Hypt
   treated, PCr treated, and O-PE treated group respectively. Scale bar,
   50 μm. F(4, 27) = 26.12, P = 6.14e−9 in one-way ANOVA. P values in
   one-way ANOVA with multiple comparison were labeled in the plot. Data
   was normalized to values of control group. All data are presented as
   mean ± s.e.m. OS: oxidative stressed cells; Hypt: hypotaurine (1 mM);
   PCr: phosphocreatine (0.5 mM); O-PE: O-phosphoethanolamine (40 μM). All
   P values were reported as multiplicity adjusted P values for multiple
   comparisons. For (b, d, f, h), Source data are provided as Source Data
   files. Blue: control; red: OS+vehicle; green: OS+Hypt; purple: OS+PCr;
   yellow: OS + O-PE.

Treatment of key metabolites regulated the metabolome of single cells

   To investigate the precise alterations in the metabolome of individual
   cells upon treatment with key metabolites, we conducted metabolic
   analysis on cells from five groups: Non-OS, vehicle-treated OS,
   Hypt-treated OS, PCr-treated OS, and O-PE-treated OS. Subsequently, we
   examined the distinct metabolic characteristics exhibited by these cell
   populations (Supplementary Fig. [223]15a). Notably, cells treated with
   key metabolites were positioned between non-OS and vehicle-treated OS
   cells, indicating an intermediary metabolic state between the two
   groups. In consistency, cells treated with key metabolites exhibited
   intermediate intensity in DCFDA and SA-β-Gal staining along with
   intermediate growth rate when compared to non-OS and OS cells
   (Fig. [224]7a–f). This suggests that the treatment of metabolites may
   have partially restored the disrupted metabolism in oxidative stressed
   cells and further alleviated cellular oxidative stress and senescence.
   To quantify the similarity between the metabolome of metabolite-treated
   cells and non-OS/OS cells, we calculated the paired distance between
   these cell types. Interestingly, compared to non-OS and OS cells, the
   cells treated with key metabolites exhibited a smaller distance to
   non-OS cells (Supplementary Fig. [225]15b). These findings indicate
   that treatment with key metabolites modulated the metabolome of OS
   cells towards a state similar to that of non-OS cells.

   Next, the detailed alteration of metabolome in OS and
   metabolite-treated cells was further studied. By conducting a
   comparative analysis of metabolites in non-OS cells and those treated
   with OS, we identified a series of downregulated metabolites under OS
   conditions, which were subsequently restored upon treatment with key
   metabolites (Supplementary Fig. [226]15c). For instance, some key
   metabolites downregulated in OS, such as hypotaurine, phosphocreatine,
   O-phosphoethanolamine, GSH, ATP (Figs. [227]2b, [228]3d), were
   recovered under the treatment of hypotaurine, phosphocreatine, and
   O-phosphoethanolamine (Supplementary Fig. [229]15c). The metabolic
   pathways involved in the recovery of the metabolome were dissected
   through MSEA analysis using the recovered metabolites under treatment
   (Supplementary Fig. [230]15d). A subset of these metabolic pathways
   overlapped with those downregulated in OS. Specifically, the “Citric
   acid cycle” in mitochondrial metabolism; “Phosphatidylethanolamine
   Biosynthesis,” “Phosphatidylcholine Biosynthesis,” and “Sphingolipid
   Metabolism” in lipid metabolism; “Aspartate metabolism” and “Glutamate
   metabolism” in amino acid metabolism; “Pyrimidine metabolism” and
   “Purine metabolism” in nucleotide metabolism; and “Amino sugar
   metabolism,” “Galactose metabolism,” and “Nucleotide sugars metabolism”
   in sugar and derivatives metabolism were all enriched with recovered
   metabolites and overlapped with pathways downregulated in OS.
   Consistently, OS cells treated with key metabolites exhibited lower
   levels of OS, reduced SA-β-Gal intensity, and improved mitochondrial
   function (Fig. [231]7 and Supplementary Fig. [232]14).

   The collective findings suggest that the modulation of crucial
   metabolites governs the cellular metabolome under OS, leading to the
   restoration of a multitude of metabolites. Treatment with these pivotal
   metabolites induces a metabolic state resembling that of non-OS cells,
   thereby highlighting their significant role in regulating OS and
   senescence pathways. Consequently, it can be inferred that specific key
   metabolites exert control over the cellular metabolome. The observed
   extensive heterogeneity and metabolic disparities in OS cells may stem
   from variations in the abundance levels of these essential metabolites.

The SCLIMS-identified protective metabolites promote healthy aging and extend
lifespan

   The effect of key metabolites treatment across two different cell types
   confirmed the regulatory role of metabolism in OS and senescence,
   leading to the question whether key metabolites further regulate animal
   aging. To explore the regulatory effects of metabolism on animals, we
   introduced the C. elegans aging model^[233]60,[234]61 to investigate
   the impact of the key metabolites (hypotaurine, phosphocreatine, and
   O-phosphoethanolamine) on the process of natural aging (Fig. [235]8a).
   A serial of concentrations of the three metabolites were added into the
   Nematode Growth Medium (NGM) at the age of L4 and throughout the
   lifespan of the worms. The OS level, lifespan and healthspan of the
   worms were then evaluated. Strikingly, the supplement of hypotaurine,
   phosphocreatine and O-phosphoethanolamine at various doses caused a
   remarkable extension in the lifespan of C. elegans by approximately
   33%-50% (Fig. [236]8b and Supplementary Fig. [237]16a).

Fig. 8. Metabolic intervention extends lifespan and promotes healthy aging in
C. elegans.

   [238]Fig. 8
   [239]Open in a new tab

   a A flowchart of experimental setup of lifespan and healthspan assay of
   C. elegans. b Lifespan of C. elegans treated with vehicle (n = 221),
   0.4 mM hypotaurine (Hypt) (n = 320), 0.2 mM phosphocreatine (PCr)
   (n = 225) and 0.1 mM O-Phosphoethanolamine (O-PE) (n = 387). P values
   in two-tailed log rank test compared with vehicle-treated worms were
   labeled in the plot, P < 2.2e−16 for all comparisons indicated in (b).
   c, d Representative images (c) and quantification (d) of DHE staining
   of L4 (young adult) (n = 14) and Day-9 (aged) worms with indicated
   treatments (n = 21, 23, 19, and 24 for vehicle, Hypt, PCr, and O-PE
   treated worms respectively). Scale bar, 200 μm. F(4, 96) = 65.74,
   P < 2.2e-16 in One-way ANOVA. P values in One-way ANOVA with multiple
   comparison were labeled in the plot. P values were reported as
   multiplicity adjusted P values for multiple comparisons. Data was
   normalized to values of vehicle treated group at Day 9. Data is
   presented as mean ± s.e.m. e, f Representative images (e) and
   quantification (f) of L4 (n = 12) and aged C. elegans thrashing under
   treatment of vehicle (n = 19), 0.4 mM hypotaurine (Hypt) (n = 28),
   0.2 mM phosphocreatine (PCr) (n = 23) and 0.1 mM O-Phosphoethanolamine
   (O-PE) (n = 13). Arrows indicate immobilized worms. Scale bar, 1 mm.
   F(4,90) = 27.58, P = 6.00e−15 in One-way ANOVA. P values in One-way
   ANOVA with multiple comparison were labeled in the plot. P values were
   reported as multiplicity adjusted P values for multiple comparisons. g,
   h Representative traces (g) and quantification (h) of free moving C.
   elegans. The traces showed the track of free moving worms in 1 min. At
   day 1, data of tracks was derived from 294, 368, 317, and 391 worms for
   vehicle, Hypt treated, PCr treated, and O-PE treated group,
   respectively. At day 5, data of tracks was derived from 184, 284, 175,
   and 305 worms for vehicle, Hypt treated, PCr treated, and O-PE treated
   group, respectively. At day 9, data of tracks was derived from 176,
   309, 316, and 153 worms for vehicle, Hypt treated, PCr treated, and
   O-PE treated group, respectively. Scale bar, 1 mm. F(3,26296) = 73.15,
   P = 4.17e−47 in Two-way ANOVA. P values in Two-way ANOVA with Turkey’s
   HSD comparison (vehicle vs O-PE, Hypt, and PCr, respectively) were
   labeled in the plot. P values were reported as multiplicity adjusted P
   values for multiple comparisons. Data were collected from 3 independent
   biological replicates. For (b, d, f, h), Source data are provided as
   Source Data files. For (b and h), blue: Vehicle treated worms; green:
   Hypt treated worms; purple: PCr treated worms; yellow: O-PE treated
   worms. For d and f, gray: L4 worms; red: vehicle treated aged worms;
   green: Hypt treated aged worms; purple: PCr treated worms; yellow: O-PE
   treated aged worms. Hypt: 0.4 mM hypotaurine; PCr: 0.2 mM
   phosphocreatine; O-PE: 0.1 mM O-Phosphoethanolamine.

   The OS level of the aged C. elegans were assessed by dihydroethidium
   (DHE) staining^[240]62. The findings unveiled a remarkable elevation in
   OS levels among aged worms, as indicated by the intensity of DHE.
   However, this surge in OS was mitigated upon supplementation with the
   three metabolites in C. elegans (Fig. [241]8c, d). The intensity of DHE
   was reduced by an average of 16% in metabolite treated worms.

   It has been reported that the mobility of worms is drastically
   compromised in aging as a sign of deterioration of health^[242]63. We
   next examined the locomotion of worms treated with and without these
   metabolites from the L4 stage onwards. Aged nematodes (Day 9) were
   subjected to the “thrashing” assay, wherein their swimming movements
   were observed to assess their physical mobility^[243]63. The elderly
   worms exhibited a 43% decline in their thrashing rate, which was
   effectively ameliorated by the administration of hypotaurine,
   phosphocreatine and O-phosphoethanolamine (Fig. [244]8e, f and
   Supplementary Fig. [245]16b).

   Another indication of the deterioration in health of aged C. elegans is
   the reduction in their speed of movement^[246]64. The worms at three
   different ages were analyzed and a progressively decline in free moving
   speed was indeed observed (Fig. [247]8g, h). The supplementation of the
   three metabolites from L4 improved the free moving speed of the aged
   worms (Fig. [248]8h and Supplementary Fig. [249]16c, d). Free moving
   speed of metabolite-treated worms was 1.4–1.8 fold higher compared with
   vehicle-treated worms at Day 9, and 1.1–1.2 fold higher at Day 5.

   Thus, the key metabolites identified by SCLIMS, including hypotaurine,
   phosphocreatine, and O-phosphoethanolamine, possess remarkable
   potential to prolong lifespan and foster graceful aging in C. elegans
   by effectively preventing cellular OS and senescence.

Discussion

   Emerging evidence indicates a diversity of cellular subtypes in
   senescence^[250]7,[251]8,[252]65,[253]66, cancer^[254]6,[255]67,
   diabetes^[256]68,[257]69, and inflammatory diseases^[258]70, suggesting
   that there is vast heterogeneity among cells in various biological
   processes. However, previous studies on cellular metabolism have
   predominantly been performed at the homogenate level, potentially
   disregarding the metabolic heterogeneity and intricate metabolic
   changes of individual cells. Due to the challenges in obtaining
   metabolomic information from individual live cells, determining
   metabolic heterogeneity in single cells has proven to be a formidable
   task. The SCMS technique we have previously established excels at
   unraveling the metabolome at a single-cell resolution^[259]17,[260]18.
   However, it encounters challenges when it comes to correlating the
   single-cell metabolome with cellular function and phenotype. Recently,
   the integrative analysis across multiple modalities brought new
   insights into cellular heterogeneity and elucidated the underlying
   mechanisms governing biological processes^[261]11,[262]71–[263]74,
   thereby illustrating the crucial role of cross-modality analysis. In
   this study, we established an approach called SCLIMS by combining SCMS
   and live-cell imaging to the metabolome of individual cells with their
   cellular OS status, thereby enabling a cross-modality analysis of both
   metabolomics profiles and cellular phenotypes at single-cell
   resolution. This technology provides some notable contributions.
   Firstly, with SCLIMS we unveiled distinct metabolic signatures among
   the six subtypes of cells under OS, each with specific oxidative
   levels. The detailed single-cell metabolic profile of OS has been
   dissected, revealing the deterioration of metabolic processes
   associated with redox balance, energy metabolism, lipid metabolism and
   mitochondrial function. The SCLIMS has not only confirmed various
   alterations in metabolism under OS as reported in previous studies, but
   also unveiled discoveries of metabolic changes, including modifications
   in amino acid metabolism and the transition from the “Malate-Aspartate
   Shuttle” to the more sophisticated “Glycerol Phosphate Shuttle”.
   Secondly, the utilization of machine learning analysis on the
   single-cell metabolome substantiates the predictive capacity of
   individual metabolic characteristics regarding cellular heterogeneity
   and phenotype. Thirdly, the SCLIMS has led to the discovery that the
   ultimate destiny of cells following OS induction is determined by their
   initial heterogeneity in metabolomics. Lastly, the key metabolites
   identified by SCLIMS exhibit protective effects against OS, cellular
   senescence, and natural aging. Overall, the SCLIMS technique sheds
   lights into the study of metabolic changes in OS and therapeutic
   interventions in aging.

   In the present study, the SCLIMS shows that the heterogenous states of
   cells can be predicted directly with their intracellular metabolome,
   showing a tight link between metabolic features and cellular phenotype.
   In addition, a more remarkable discovery revealed by the SCLIMS is that
   the destiny of cells under OS can be determined by their initial
   metabolomics status. This technology has verified that GSH is strongly
   correlated with many metabolites in the metabolome and may play a key
   role in cell fate determination, as evidenced by the fact that GSH-rich
   cells exhibit greater resistance to OS. However, it should be noted
   that such resistance is not solely attributed to GSH; rather, our
   SCLIMS analysis reveals elevated levels of other key metabolites such
   as hypotaurine, and O-phosphoethanolamine in cells exhibiting enhanced
   oxidative resistance. Indeed, the anti-oxidant effect of these
   metabolites has been further confirmed in the aged C. elegans. In
   essence, the GSH-cored metabolome serves as the determinant of cell
   phenotype and fate in OS. Among the top 15 rewired metabolites, GSH is
   well correlated with the other 13 metabolites including glutamate,
   creatine, glutamine, taurine, threonine, N-acetyl-aspartic acid,
   aspartic acid, asparagine, UDP-N-acetylglucosamine, proline, GABA,
   cystathionine, and glucose. These metabolites participate amino acids
   metabolism as well as carbohydrate metabolism, which were compromised
   in cells with higher OS levels (Fig. [264]3e). Other studies reported
   the effects of these top-rewired metabolites in aging and senescence as
   well. For example, taurine was reported to decline in aging and
   supplement of taurine increases healthspan and lifespan in various
   species^[265]75. Taurine protected telomerase and mitochondrial
   function, along with decreasing inflammation and DNA damage. Creatine
   was reported to promote healthy aging by attenuating inflammation and
   preventing bone mineral loss^[266]76. Creatine improved neuronal
   function by antioxidant effect and exhibited therapeutic effect against
   age-related diseases including Alzheimer’s disease, Parkinson’s
   disease, and heart failure^[267]77. The level of threonine, aspartic
   acid, and proline were reported to be positively correlated with
   lifespan of yeasts^[268]78. Glutamine promotes autophagy via AMPKα
   lactylation and suppresses senescence^[269]79. N-acetyl-aspartic acid
   was reported to decline in aging and was related to brain
   atrophy^[270]80. Asparagine prevented stem cell aging by regulating the
   autophagy-lysosome pathway^[271]81. The underlying mechanism behind the
   notable disparity in metabolic characteristics among initial cells
   remains elusive and necessitates further investigations. However, a
   possible explanation is that it could be a result of the asymmetric
   division of cytoplasm during cell division^[272]82–[273]85. Another
   possibility may be attributed to differences in cellular contact,
   micro-environmental development, or potential transport of certain
   metabolites between neighboring cells that has gone undetected^[274]86.
   However, from a metabolic view, abundance of key metabolites determined
   the metabolomic profile of cells and recovered the disrupted metabolic
   pathways and the metabolome under OS along with reducing the OS and
   senescent level (Fig. [275]7 and Supplementary Figs. [276]14, [277]15),
   suggesting the role of key metabolites in determining the metabolic and
   phenotypic heterogeneity.

   The SCLIMS-identified key metabolites provides opportunities in
   senescence and aging intervention. The OS-induced senescence is one of
   the important models of cellular senescence^[278]22,[279]54,[280]55 and
   play crucial role in various diseases^[281]87,[282]88. The
   supplementation of the key metabolites, including hypotaurine,
   phosphocreatine, and O-phosphoethanolamine, can effectively mitigate
   the cellular OS and the OS-induced senescence. Moreover, these
   metabolites also prevent OS, prolong lifespan, promote healthy aging
   and delay the decline in mobility during aging in C. elegans. Combining
   the clues from the literature, we posit that these metabolites may
   regulate cellular OS and senescence through multiple mechanisms. For
   instance, hypotaurine serves as a hydrogen donor for NAD^+ during its
   conversion into taurine, thereby generating NADH as a
   by-product^[283]31. This process effectively restores redox balance in
   the presence of OS. Phosphocreatine serves as a direct catalyst for the
   conversion of ADP into ATP^[284]89, functioning as a quintessential
   cellular energy reservoir that possesses the inherent capability to
   reinstate equilibrium in energy levels. O-phosphoethanolamine has been
   reported to mitigate mitochondrial dysfunction^[285]90 and effectively
   restore membrane, as it serves as the fundamental precursor of membrane
   lipids^[286]91. Therefore, these metabolites may potentially impede
   senescence by modulating various cellular processes, including energy
   metabolism, mitochondrial function, and lipid metabolism. The
   deficiency of such metabolites may render cells more susceptible to
   senescence-inducing factors, such as OS. Furthermore, a more paramount
   consequence of the metabolic intervention lies in its ability to
   prolong healthspan, which assumes a relatively pivotal role in the
   realm of aging research when compared to lifespan^[287]63,[288]92. The
   decline of physical function is a common occurrence in both the early
   and late stages of aging, and bestowing longevity upon frailty offers
   minimal advantage to the individual^[289]93. Therefore, screening for
   potential metabolites from our SCLIMS database presents an opportunity
   for promoting health benefits during the aging process.

   The current study has certain limitations that we would like to
   address, along with potential solutions for future studies. Firstly,
   the identification of metabolites at the single-cell level poses a
   significant challenge due to the complexity of MS/MS analysis.
   Acquiring MS/MS spectra for hundreds of m/z in the metabolome is indeed
   arduous. However, there are promising avenues to enhance metabolite
   identification efficacy. For instance, ion mobility mass spectrometry
   can potentially differentiate metabolites sharing identical m/z values
   by considering collision cross section^[290]94–[291]96. Additionally,
   optimizing scan speed and extending sampling duration can facilitate
   acquiring comprehensive MS/MS spectra from single cells^[292]97.
   Secondly, the SCLIMS utilized in this current study was meticulously
   designed to incorporate cultured cells and cellular models of OS. The
   versatility of SCLIMS extends to tissue-embedded cells, as they can
   also be effectively labeled with fluorescence markers. However, the
   potential application of SCLIMS in tissue-embedded cells remains an
   intriguing area for future exploration and investigation. Thirdly, one
   must acknowledge the challenges associated with analyzing fixed cells
   when employing the SCLIMS technique. Nevertheless, by enabling analysis
   of fixed cells, a myriad of phenotypic features such as SA-β-Gal
   staining, immunofluorescence analysis, and immunohistochemistry can be
   seamlessly integrated into the cell metabolome profiling process.

   In summary, this study presents a cross-modality analysis integrating
   single-cell metabolomic profile and cellular phenotype enhancing our
   understanding of cell heterogeneity and subtype-specific metabolic
   signatures in a cellular model of OS. Most importantly, the
   cross-modality platform and analysis described in this study provide a
   way in single-cell research. The single-cell metabolome and the
   cellular phenotype such as OS status are directly linked and
   integrated. The heterogeneous states are explained with single-cell
   metabolome and metabolic pathways. Significantly, the insights into the
   metabolic regulation governing OS, cellular senescence, and natural
   aging serves as a valuable resource for future investigations into
   interventions targeting oxidative damage, aging and senescence.
   Furthermore, this cutting-edge platform possesses the remarkable
   capability to integrate single-cell metabolomics profiling with a
   diverse array of cellular phenotypes assessed by live-cell labeling.
   For instance, (1) The SCLIMS can be combined with live-cell
   mitochondrial probes such as probes for mitochondrial membrane
   potential (i.e. TMRE probe) and mitochondrial morphology (i.e.
   Mito-Tracker). This enables the study of the interaction between
   mitochondrial function and the metabolome at single-cell level; (2) The
   SCLIMS can be integrated with calcium imaging^[293]98 which labels
   neuronal activities and enables the performance of multi-modal analysis
   of heterogeneity in metabolome and neuronal functions; (3) The SCLIMS
   can be utilized to investigate the correlation between cellular
   metabolome and cell cycle by incorporating dynamic live-cell
   fluorescent probes for real-time monitoring of cell division and
   proliferation^[294]99. Thus, with any technique labeling live cells
   with fluorescent, this cross-modality platform will become a feasible
   way for integrative analysis and a powerful tool for the discovery of
   secrets in single cells.

Methods

Chemicals

   NaCl, KCl, CaCl[2], MgCl[2], HEPES, NaOH, sucrose, NH[4]HCO[3],
   Na[2]HPO[4], KH[2]PO[4], K[2]HPO[4], MgSO[4], cholesterol, hypotaurine,
   and O-phosphoethanolamine were purchased from Sigma-Aldrich. Phosphate
   buffer saline was purchased from Sangon Biotech. Dulbecco’s Modified
   Eagle’s medium (DMEM) was purchased from HyClone. Fetal bovine serum
   and trypsin-EDTA (0.25%) were purchased from Gibco. Trypan blue,
   penicillin and streptomycin were purchased from Biosharp. Hydrogen
   peroxide was purchased from Sinopharm. Phosphocreatine was purchased
   from Aladdin. The Senescence-Associated β-Galactosidase kit and
   Mitochondrial membrane potential assay kit were purchased from
   Beyotime. A cellular ROS assay kit and dihydroethidium (DHE) were
   purchased from Abcam. Live-cell GSH probe (mClB) was purchased from
   MedChemExpress.

Cell culture

   HEK293T cell line were originally obtained from ATCC (CRL-3216). The
   cell line was authenticated by ATCC with STR profiling. Primary MEFs
   were a kind gift from Professor Chunlei Cang in University of Science
   and Technology of China. All cells were cultured in Dulbecco’s Modified
   Eagle’s medium (DMEM) (HyClone), supplemented with 10% fetal bovine
   serum (FBS, Gibco) and 100 U/ml penicillin 100 µg/ml streptomycin
   (Biosharp) at a temperature of 37 °C, with 5% CO[2] in a humidified
   atmosphere. The culture medium was refreshed every 2-3 days, and the
   cells were subcultured every 3-5 days when they reached approximately
   80% confluency.

C. elegans strain and maintenance

   The Caenorhabditis elegans (C. elegans) were cultured and maintained on
   Nematode Growth Medium (NGM) seeded with E. Coli OP50 at 20 °C. N2
   (wild isolate) strain was used in all C. elegans experiments. Plates
   were maintained by transferring the worms every 3 days. For DHE
   staining and thrashing analysis, worms at L4 and at Day 9 after L4 were
   used. For lifespan analysis, the survival of worms was observed
   throughout the whole lifespan. For activity analysis, worms at Day 1,
   Day 5 and Day 9 after L4 were used.

Oxidative stress model

   Cells were seeded and allowed to grow overnight. Subsequently, they
   were treated with hydrogen peroxide at a final concentration of 80 μM
   for HEK293T cells and 240 μM for MEFs in the culture medium for 1 h,
   followed by replacement with fresh medium. Finally, the cells were
   cultured for an additional 48 h to establish an oxidative stress model.

Treatment of cells with metabolites

   Cells were seeded and allowed to grow overnight. Subsequently, the
   cells were treated with hydrogen peroxide at a final concentration of
   80 μM for HEK293T cells and 240 μM for MEFs in culture medium for 1 h.
   After that, the medium was replaced with fresh medium supplemented with
   specific metabolites at indicated concentrations. The cells were then
   cultured for an additional 48 h before further assays.

Metabolite treatment and life span assay of C. elegans

   Life span assays of C. elegans were conducted at a temperature of
   20 °C. Metabolite treatment was administered by adding specific
   metabolites at the indicated concentrations to NGM plates, which were
   then incubated overnight prior to use. Following bleaching,
   age-synchronized eggs were washed with M9 buffer and subsequently
   placed on NGM plates. Late L4 larvae or young adult worms were
   subsequently transferred to NGM plates that had been seeded with
   heat-inactivated OP50 E. coli and supplemented with 0.1 mg/ml of
   5’-FUDR, as well as the indicated treatment of metabolites.
   Approximately 100 worms were placed on each plate, which was then
   inspected and scored every one to two days. The worms were moved to
   fresh plates every one to two days in order to ensure the efficacy of
   the drugs and metabolites. Worms that exhibited no response to
   mechanical stimulation were considered deceased. Worms displaying a
   “protruding vulva”, those that were lost, or had burrowed into the
   medium were censored. Statistical analysis was conducted using the
   MATLAB function “logrank”
   ([295]www.mathworks.com/matlabcentral/fileexchange/22317), and P values
   were calculated.

Behavioral analysis of C. elegans

   The ‘thrashing’ assay was employed to assess the locomotion of C.
   elegans. On Day 9 post-adulthood, worms subjected to specific
   treatments were transferred to M9 buffer and allowed to acclimate for
   1 min before body bends were quantified using a dissecting microscope.

   The physical function of C. elegans was assessed by monitoring the
   locomotion of the worms on NGM plates. On Day 1, Day 5, and Day 9
   post-adulthood, freely moving worms subjected to specific treatments
   were recorded using a digital camera and analyzed in ImageJ (version
   1.54 g, [296]https://imagej.nih.gov/ij/) with the ‘wrMTrck’ plugin as
   per manual instructions. The traveling speed was then calculated using
   the plugin.

Cell viability assay

   Cell viability assay was performed according to the manufacturer’s
   manual. Cells were dissociated using trypsin-EDTA (0.25%) (Gibco) for
   1 min at 37 °C, followed by termination of the dissociation process
   through the addition of an equal volume of DMEM supplemented with 10%
   fetal bovine serum. Cells were suspended, centrifuged, and then
   resuspended with PBS. Subsequently, they were stained with trypan blue
   (Biosharp) at a final concentration of 0.04%. The cells were then
   enumerated under a microscope; the stained cells were designated as
   non-viable. Cell viability was determined by calculating the ratio of
   unstained cells to the total number of both stained and unstained
   cells.

SA-β-Gal assay

   The experimental procedure was performed in accordance with the
   protocols provided by the SA-β-Gal assay kit (Beyotime, China).
   Briefly, cells were washed with phosphate-buffered saline (PBS) and
   fixed with fixing reagents at room temperature for 15 min. After three
   washes with PBS for 3 min each, staining solution was prepared
   according to the manufacturer’s instructions prior to use. Cells were
   stained overnight at 37 °C, and images of five randomly selected fields
   were captured using a bright field setting for subsequent analysis.

Oxidative stress assay

   To evaluate the extent of cellular oxidative stress, a live-cell probe
   DCFDA was prepared according to the manufacturer’s protocol (Abcam).
   Subsequently, live cells were incubated with a 10 μM DCFDA solution for
   25 min at 37 °C and 5% CO[2], followed by PBS washing. Images were
   promptly captured using a fluorescent microscope (Leica), and five
   random fields were selected for analysis in each dish.

   To assess the level of oxidative stress in C. elegans, worms were
   exposed to a final concentration of 3 μM DHE (Dihydroethidium) in M9
   buffer at 20 °C for 30 min, followed by washing with M9 buffer.
   Subsequently, the worms were transferred onto glass slides and imaged
   using a fluorescent microscope (Leica).

Mitochondrial membrane potential assay

   TMRE was utilized to assess the mitochondrial membrane potential (MMP)
   of viable cells in accordance with the manufacturer’s instructions.
   Briefly, cells were incubated with 1X TMRE (Beyotime) in serum-free
   DMEM for 15 min at 37 °C. Subsequently, the cells were rinsed with warm
   serum-free DMEM and immediately imaged using a fluorescent microscope
   (Leica, version 4.6.2 build: 410). Images of five randomly selected
   fields were captured for analysis.

Analysis of SA-β-Gal intensity

   The SA-β-Gal-stained area in images of random fields was extracted
   using the ‘IHC toolbox’ plugin in ImageJ, resulting in a new image of
   the stained area. The color images of the extracted SA-β-Gal stained
   area were then converted to 8-bit grayscale and calibrated within
   ImageJ (version 1.54g). Finally, the gray value (integrated density) of
   the stained area was calculated using ImageJ. The intensity of SA-β-Gal
   was quantified by calculating the gray value (integrated density) of
   the stained area normalized to the total cell area in each image.

Analysis of fluorescence intensity

   For random field images, color images of cells stained with fluorescent
   live-cell probes (DCFDA and TMRE) were imported into ImageJ software.
   The images were then converted to 8-bit grayscale and calibrated before
   being thresholded. The fluorescence intensity was calculated as the
   mean gray value of cells in the field by measuring the integrated
   density of fluorescent positive area divided by the fluorescent
   positive area. For single cells, the image was converted to an 8-bit
   grayscale image and calibrated. Fluorescent intensity was analyzed
   using ROI manager in ImageJ, which allowed for analysis of individual
   cells through selection. The mean gray value for each cell was
   calculated by dividing the integrated density of a cell by its area.

   For C. elegans, the level of oxidative stress (indicated by DHE
   fluorescence) was quantified using ImageJ software. The color images
   were converted to 8-bit grayscale and calibrated prior to thresholding.
   The gray value (integrated density) of fluorescent positive areas in
   individual worms was then calculated and normalized by the
   corresponding area.

The workflow and experimental setting of SCLIMS

   The SCLIMS platform comprises two primary components: live-cell imaging
   and single-cell MS. Cells were initially stained with a 10 μM DCFDA
   solution for 25 min at 37 °C and 5% CO[2], followed by PBS washing.
   Subsequently, the cells were captured using a fluorescent microscope to
   record their spatial distribution. The acquired images were
   subsequently subjected to analysis, enabling the calculation of
   oxidative stress levels in individual cells.

   The cells were then transferred to the single-cell patch clamp platform
   and incubated in a bath solution containing 140 mM NaCl, 5 mM KCl, 2 mM
   CaCl[2], 1 mM MgCl[2] and 10 mM HEPES (pH adjusted to 7.4 with NaOH;
   ∼320 mosmol with sucrose) and approached by a borosilicate glass
   pipette filled with pipette solution (88 mM NH[4]HCO[3]) using a
   micromanipulator. The cells were selected based on the fluorescent
   images acquired in the preceding step. The cells were patched with a
   high-quality seal (>1 GΩ) and the cell membrane was disrupted by rapid
   application of negative pressure. Mild negative pressure was then
   applied to the pipette to obtain cytoplasmic chemical constituents,
   which were subsequently analyzed using mass spectrometry after removal
   of the pipette from the bath solution.

   Following the extraction of cellular cytoplasmic constituents, the
   capillary was connected to a MS system as described below. An AC
   voltage of 4 kV amplitude and approximately 500 Hz frequency was
   applied externally to the spray capillary micropipette, while
   maintaining a distance of approximately 5 mm between the tip of the
   spray micropipette and the orifice of the MS instrument.
   High-resolution mass measurements were analyzed using a Q Exactive Plus
   MS instrument (Thermo Fisher Scientific, San Jose, CA, USA). nanoESI
   source and Orbitrap mass analyzer were used. The main experimental
   parameters were established as follows: capillary temperature at
   275 °C, S-lens radio frequency (RF) level set to 50%, mass resolution
   of 70,000, maximum injection time of 10 ms, AGC target of 1e6, and
   microscan rate of 1. Negative ion mode was employed throughout all
   experiments. Data was acquired under full scan mode. For MS/MS,
   collision energy was set to hcd = 30. Data was collected with Thermo
   Scientific Exactive Tune software (version 2.9.0.2926). The MS data
   were then processed. Data of each single cell was paired with oxidative
   stress levels based on the fluorescent images.

Single-cell mass spectrometry data processing

   The spectral data of individual cells were initially stored in separate
   files by the instrument. Subsequently, all files were converted into
   mzML format using ProteoWizard (version 3.0.9870). The XCMS package
   (version version 3.9.1) in R statistical environment was utilized to
   process all single-cell MS data, encompassing peak calling, peak
   alignment, and quality control. A signal-to-noise (S/N) filter of 3 was
   applied to the m/z signal, and metabolic signals were identified as
   those with a frequency exceeding 20% across all tested cells. The data
   was organized into a matrix, where metabolites were represented by rows
   and cells were represented by columns. All intensities were normalized
   to the total ion current (TIC) ratio. Metabolites were annotated by
   comparing observed m/z values with theoretical values in the Human
   Metabolome Database v5.0 ([297]www.hmdb.ca), and m/z annotations were
   assigned if errors fell within 10 ppm. To ensure the precision of the
   annotations, metabolites underwent further identification employing
   MS/MS. For metabolites involved in treating OS cells and worms
   including hypotaurine, phosphocreatine, and O-phosphoethanolamine,
   their confirmation relied on matching MS/MS spectra between cells and
   standards (Supplementary Fig. [298]17). Other metabolites were
   identified by comparing MS/MS spectra between bulk cellular samples and
   the HMDB database. Batch effects of experiments were assessed through
   PCA and HCA analysis in the R statistical environment.

Flow cytometry

   The cells were seeded and cultured until they achieved a confluency of
   70-80%. Following a single wash with 1 mL of phosphate-buffered saline
   (PBS), the cells were dissociated using  0.25% trypsin solution
   containing EDTA, after which the reaction was halted by supplementing
   DMEM with FBS at a concentration of 10%. Subsequently, the cell
   suspension was collected in tubes measuring approximately 1.5 mL
   capacity and subjected to centrifugation at a speed of 600 × g for 5
   min before being resuspended in PBS containing mClB at a final
   concentration of 40 μM. The resuspended cells were then incubated at a
   temperature of precisely maintained at or around 37 °C for 20 min.
   Subsequently, the cells were washed and resuspended in PBS supplemented
   with 1% FBS. After filtration through a 40 μm cell strainer and
   transfer to clean tubes, all samples were vortexed for 5 s to ensure
   complete dissociation into single-cell suspension prior to flow
   cytometry analysis. Flow cytometric analysis was performed using a BD
   FACSAria III instrument (BD Biosciences), with initial gating based on
   FSC-A versus SSC-A parameters followed by measurement of whole-cell
   fluorescence. Unstained cells were utilized as a blank control to
   establish baseline correction. Finally, cells stained with mClB
   exhibiting the top and bottom 5% fluorescent intensity were isolated
   and collected independently for subsequent analysis. The FACS data was
   then analyzed and visualized using the ‘fca_readfcs’ function and the
   ‘Flow cytometry GUI for MATLAB’ plugin in MATLAB (version 2022b)
   ([299]www.mathworks.com/matlabcentral/fileexchange/9608-fca_readfcs and
   [300]www.mathworks.com/matlabcentral/fileexchange/38080-flow-cytometry-
   gui-for-matlab).

Identification of marker metabolites

   Wilcoxon rank-sum tests were used to compare metabolites in cells
   belonging to one cluster with those of all other cells, based on the
   results obtained from unsupervised clustering analysis. Ions exhibiting
   P < 0.05 and fold change (FC) > 1.5 were identified as marker
   metabolites for the specific cluster, which were subsequently z score
   scaled and visualized using MATLAB’s ‘heatmap’ function. The z score
   scales the data for each ion and is calculated as follows:
   [MATH: <mi>z</mi><mo>=</mo><mfrac><mrow><mi>x</mi><mo>−</mo><mover
   accent="true"><mrow><mi>X</mi></mrow><mo>¯</mo></mover></mrow><mrow><mi
   >S</mi></mrow></mfrac> :MATH]
   1

   where z represents z score, x represents a sample raw data,
   [MATH: <mover accent="true"><mrow><mi>X</mi></mrow><mo>¯</mo></mover>
   :MATH]
   represents the population mean, and S represents the population
   standard deviation.

Identification of metabolic subtypes

   The metabolic data matrix was imported into the MATLAB workspace, where
   it underwent z score scaling and unsupervised clustering using the
   k-medoids algorithm. The resulting clusters were visualized with
   Uniform Manifold Approximation and Projection (UMAP) via a MATLAB
   plugin
   ([301]https://www.mathworks.com/matlabcentral/fileexchange/71902). The
   data was dimensionally reduced to two dimensions and subsequently
   visualized using the ‘gscatter’ function in MATLAB (version 2022b),
   revealing distinct metabolic subtypes of cells within a bi-dimensional
   space.

Pseudotime analysis and single-cell trajectory construction

   The metabolite data matrix was initially filtered based on the
   correlation with oxidative stress level (indicated by DCFDA intensity).
   Metabolites exhibiting a correlation coefficient (r) greater than 0.2
   or less than -0.2 and P < 0.05 in correlation analysis were selected
   for pseudotime analysis. The refined data matrix was then imported into
   Monocle (version 2.16.0) in R statistical environment, followed by
   dimensional reduction using the ‘DDRTree’ algorithm according to the
   documentation of Monocle
   ([302]https://cole-trapnell-lab.github.io/monocle-release/). Finally,
   the ‘plot_cell_trajectory’ function was utilized to visualize the
   trajectory.

Supervised machine learning and model evaluation

   Supervised machine learning was conducted using MATLAB (version 2022b).
   The data matrix was randomly divided into a training dataset and a
   testing dataset at a ratio of 2:1. The testing dataset was exclusively
   used for evaluating the trained model and not exposed to the training
   phase. The training of both the classification and regression models
   was conducted without any feature selection, and all m/z signals that
   met our criteria (S/N > 3 and detected in greater than 20% single
   cells) were included in both the training and testing datasets. To
   train classification models, ensemble algorithm (function
   ‘fitcensemble’), discriminant analysis algorithm (function ‘fitcdiscr’)
   and neural network algorithm (function ‘fitcnet’) were employed as
   specified. The model underwent 5-fold cross-validation during training,
   with Bayesian optimization employed to optimize hyperparameters and
   minimize cross-validation loss (error). Parameters ‘Method’,
   ‘NumLearningCycles’, ‘MinLeafSize’, and ‘LearnRate’ were automatically
   optimized in the ensemble algorithm, while parameters ‘Delta’ and
   ‘Gamma’ were automatically optimized in the discriminant analysis
   algorithm. Parameters such as ‘Activations’, ‘Lambda’, ‘LayerSizes’,
   and ‘Standardize’ were automatically optimized for neural network
   algorithm. For the training of regression models, we utilized the
   neural network algorithm (function ‘fitrnet’). The data underwent a
   “log2” or “ln” transformation and was trained using 5-fold
   cross-validation and Bayesian optimization. Parameters such as
   ‘Activations’, ‘Lambda’, ‘LayerSizes’, and ‘Standardize’ were
   automatically optimized.

   After training the model with the training dataset, the testing dataset
   was utilized to evaluate the performance of the classification models
   through receiver operating characteristic (ROC) curve and confusion
   matrix analysis. For multiclass problem, ‘one versus rest’ strategy was
   used which transformed the problem into a two-classification task. The
   predicted vs real scatter plot and the Pearson’s r and P value were
   used to evaluate the performance of the regression model.

Evaluation of metabolic similarity

   The ‘pdist’ function in MATLAB (version 2022b) was utilized to compute
   the statistical distance based on the cellular metabolome, where a
   larger distance indicated a lower degree of similarity between two
   cells. Similarity was calculated as the reciprocal of distance. The
   heatmap was generated by plotting the reciprocal of the distance as a
   measure of similarity, with higher values indicating greater likelihood
   of shared metabolomic features between cells. We quantified the
   pairwise distances between cells in the initial and OS groups based on
   their metabolite abundance as variables. Spearman distance functions
   were employed.

   Spearman distance function:
   [MATH:
   <msub><mrow><mi>d</mi></mrow><mrow><mi>s</mi><mo>,</mo><mi>t</mi></mrow
   ></msub><mspace width="0.25em"></mspace><mo>=</mo><mspace
   width="0.25em"></mspace><mn>1</mn><mo>−</mo><mfrac><mrow><mfenced
   close=")"
   open="("><mrow><msub><mrow><mi>r</mi></mrow><mrow><mi>s</mi></mrow></ms
   ub><mo>−</mo><msub><mrow><mover
   accent="true"><mrow><mi>r</mi></mrow><mo>¯</mo></mover></mrow><mrow><mi
   >s</mi></mrow></msub></mrow></mfenced><msup><mrow><mfenced close=")"
   open="("><mrow><msub><mrow><mi>r</mi></mrow><mrow><mi>t</mi></mrow></ms
   ub><mo>−</mo><msub><mrow><mover
   accent="true"><mrow><mi>r</mi></mrow><mo>¯</mo></mover></mrow><mrow><mi
   >t</mi></mrow></msub></mrow></mfenced></mrow><mrow><mo>′</mo></mrow></m
   sup></mrow><mrow><msqrt><mrow><mfenced close=")"
   open="("><mrow><msub><mrow><mi>r</mi></mrow><mrow><mi>s</mi></mrow></ms
   ub><mo>−</mo><msub><mrow><mover
   accent="true"><mrow><mi>r</mi></mrow><mo>¯</mo></mover></mrow><mrow><mi
   >s</mi></mrow></msub></mrow></mfenced><msup><mrow><mfenced close=")"
   open="("><mrow><msub><mrow><mi>r</mi></mrow><mrow><mi>s</mi></mrow></ms
   ub><mo>−</mo><msub><mrow><mover
   accent="true"><mrow><mi>r</mi></mrow><mo>¯</mo></mover></mrow><mrow><mi
   >s</mi></mrow></msub></mrow></mfenced></mrow><mrow><mo>′</mo></mrow></m
   sup></mrow></msqrt><msqrt><mrow><mfenced close=")"
   open="("><mrow><msub><mrow><mi>r</mi></mrow><mrow><mi>t</mi></mrow></ms
   ub><mo>−</mo><msub><mrow><mover
   accent="true"><mrow><mi>r</mi></mrow><mo>¯</mo></mover></mrow><mrow><mi
   >t</mi></mrow></msub></mrow></mfenced><msup><mrow><mfenced close=")"
   open="("><mrow><msub><mrow><mi>r</mi></mrow><mrow><mi>t</mi></mrow></ms
   ub><mo>−</mo><msub><mrow><mover
   accent="true"><mrow><mi>r</mi></mrow><mo>¯</mo></mover></mrow><mrow><mi
   >t</mi></mrow></msub></mrow></mfenced></mrow><mrow><mo>′</mo></mrow></m
   sup></mrow></msqrt></mrow></mfrac> :MATH]
   2

   where d[s,t] represents the distance between two cells (s and t); r[sj]
   is the rank of x[sj] taken over x[1j], x[2j],…x[mj]; r[s] and r[t] are
   the coordinate-wise rank vectors of x[s] and x[t], as an example,
   r[s] = (r[s1], r[s2],… r[sn]);
   [MATH: <msub><mrow><mover
   accent="true"><mrow><mi>r</mi></mrow><mo>¯</mo></mover></mrow><mrow><mi
   >s</mi></mrow></msub><mo>=</mo><mfrac><mrow><mn>1</mn></mrow><mrow><mi>
   n</mi></mrow></mfrac><msubsup><mrow><mo>∑</mo></mrow><mrow><mi>i</mi><m
   o>−</mo><mn>1</mn></mrow><mrow><mi>n</mi></mrow></msubsup><msub><mrow><
   mi>r</mi></mrow><mrow><mi>s</mi></mrow></msub><mfenced close="]"
   open="["><mrow><mi>i</mi></mrow></mfenced> :MATH]
   and
   [MATH: <msub><mrow><mover
   accent="true"><mrow><mi>r</mi></mrow><mo>¯</mo></mover></mrow><mrow><mi
   >t</mi></mrow></msub><mo>=</mo><mfrac><mrow><mn>1</mn></mrow><mrow><mi>
   n</mi></mrow></mfrac><msubsup><mrow><mo>∑</mo></mrow><mrow><mi>i</mi><m
   o>−</mo><mn>1</mn></mrow><mrow><mi>n</mi></mrow></msubsup><msub><mrow><
   mi>r</mi></mrow><mrow><mi>t</mi></mrow></msub><mfenced close="]"
   open="["><mrow><mi>i</mi></mrow></mfenced> :MATH]
   . n represents the number of variables.

   Hamming distance was used in the comparison of metabolome of Hypt, PCr,
   O-PE treated OS cells with Non-OS and vehicle treated OS cells.

   Hamming distance function:
   [MATH:
   <msub><mrow><mi>d</mi></mrow><mrow><mi>s</mi><mo>,</mo><mi>t</mi></mrow
   ></msub><mo>=</mo><mfenced close=")" open="("><mrow><mfrac><mrow><mi
   mathvariant="normal">#</mi><msub><mrow><mi>x</mi></mrow><mrow><msub><mr
   ow><mi>s</mi></mrow><mrow><mi>j</mi></mrow></msub></mrow></msub><mo>≠</
   mo><msub><mrow><mi>x</mi></mrow><mrow><msub><mrow><mi>t</mi></mrow><mro
   w><mi>j</mi></mrow></msub></mrow></msub><mstyle
   mathsize="1.19em"><mfenced
   open=")"><mrow></mrow></mfenced></mstyle></mrow><mrow><mi>n</mi></mrow>
   </mfrac></mrow></mfenced> :MATH]
   3

   where d[s,t] represents the distance between two cells (s and t); x[sj]
   represents the variable j of the first cell s, x[tj] represents the
   variable j of the second cell t; n represents the number of variables.
   The function calculates the fraction of different variables in all
   variables (n) of the two cells.

Metabolite set enrichment analysis

   Metabolite set enrichment analysis (MSEA) was performed using
   MetaboAnalyst (v6.0), an online metabolomics analysis tool
   ([303]www.metaboanalyst.ca). The annotation of metabolites, including
   metabolic markers in clusters and those correlated with the oxidative
   stress level of single cells, was compiled into a list and uploaded to
   the website. Only metabolites identified with MS/MS confirmation was
   used in the MSEA. The algorithm processed the data and obtained
   enriched pathways, while also calculating two important parameters:
   Enrichment Ratio and P value. Enrichment Ratio indicates the degree of
   enrichment, while P value represents its significance. Pathways with a
   P value less than 0.05 were deemed significant. The results were
   downloaded from the website and visualized using MATLAB’s ‘bubblechart’
   function.

Construction of metabolic networks and rewiring analysis

   The annotated metabolite data was analyzed using the ‘corr’ function in
   MATLAB (version 2022b), resulting in a matrix of correlation
   coefficients (r). This matrix was then reorganized into three columns:
   metabolite-1, metabolite-2, and their respective correlation
   coefficient values. The reorganized matrix was imported into Cytoscape
   (version 3.10.1) to construct metabolic networks based on the
   correlations between metabolites. In analysis of Cluster-I and
   Cluster-II in initial cells, metabolite pairs with P > 0.05 or
   |r | <0.8 were excluded. In the network, nodes represent metabolites,
   and edges represent correlations. The networks of correlation were
   visualized using Cytoscape. By utilizing the DyNet algorithm to compare
   two different networks, changes in relationships with other metabolites
   for each individual metabolite were analyzed. Metabolite rewiring
   scores were calculated using the DyNet algorithm and then visualized
   through stem plot with MATLAB’s ‘stem’ function.

Statistical analysis

   The statistical analysis was conducted using Microsoft Excel
   (Microsoft), R statistical environment, MATLAB 2022b (Mathworks), and
   GraphPad Prism (version 8). The Wilcox rank sum test was employed to
   determine the significance of two-grouped data, while one-way ANOVA or
   two-way ANOVA were utilized for multi-group or multi-factor data,
   respectively. Distribution analysis (using the “Rayleigh” probability
   density function) and calculation of interquartile range (IQR) and
   median absolute deviation (MAD) were performed using MATLAB R2022b
   (MathWorks). Sample size was not predetermined by statistical methods
   in this study, but rather based on previous experience. The number of
   samples (n) is indicated in the figures or figure legends.

Reporting summary

   Further information on research design is available in the [304]Nature
   Portfolio Reporting Summary linked to this article.

Supplementary information

   [305]Supplementary Information^ (31.7MB, pdf)
   [306]Reporting Summary^ (450.3KB, pdf)
   [307]Supplementary Data 1^ (13.1KB, xlsx)
   [308]Supplementary Data 2^ (99.4KB, zip)
   [309]Transparent Peer Review file^ (13.3MB, pdf)

Source data

   [310]Source Data^ (4.9MB, zip)

Acknowledgements